Contrasting Phylogeographic Patterns Among Northern and Southern Hemisphere Fin Whale Populations With New Data From the Southern Pacific

Four fin whale sub-species are currently considered valid: Balaenoptera physalus physalus in the North Atlantic, B. p. velifera in the North Pacific, B. p. quoyi and B. p. patachonica in the Southern Hemisphere. The last, not genetically validated, was described as a pygmy-type sub-species, found in low to mid latitudes of the Southern Hemisphere. Genetic analyses across hemispheres show strong phylogeographic structure, yet low geographic coverage in middle latitudes of the Southern Hemisphere impeded an assessment within the area, as well as evaluating the validity of B. p. patachonica. New mtDNA sequences from the Southeastern Pacific allowed an improved coverage of the species’ distribution. Our phylogenetic analyses showed three main lineages and contrasting phylogeographic patterns between Northern and Southern Hemispheres. Absence of recurrent female mediated gene flow between hemispheres was found; however, rare dispersal events revealing old migrations were noted. The absence of genetic structure suggests the existence of one single taxa within the Southern Hemisphere. Thus, until further evidence supporting this subspecies can be produced, such as genetic, ecological, behavioral, or morphological data, we propose that all fin whales from the Southern Hemisphere, including those from middle latitudes of the Southeastern Pacific belong to B. p. quoyi subspecies. This information is important for the current assessment of fin whales, contributing to the evaluation of the taxonomic classification and the conservation of the species.

The taxonomic classification of the Northern (B. p. physalus) and Southern Hemisphere (B. p. quoyi) subspecies was proposed in the mid twentieth century (Tomilin, 1946), primarily based on size differences between fin whales from the Northern Hemisphere (NH) and Antarctic waters. This was later verified by comparisons of body measurements and organ weights of North Atlantic and Antarctic individuals (Lockyer and Waters, 1986). The fourth subspecies, B. p. patachonica, was defined as a pygmy-type (approximately 18-24 m total length) and dark in coloration, possibly with black baleen. Its distribution was described as being restricted to low to mid-latitudes in the SH, not extending further south than approximately 55 • S, hence feeding at lower latitudes than most fin whales (Clarke, 2004). Aside from the original description, there has been no further evidence supporting this subspecies.
Strong genetic differentiation was detected in the Northern Hemisphere, between whales of the North Pacific (NP) and North Atlantic (NA) (Bérubé et al., 1998;Archer et al., 2013;Cabrera et al., 2019). Based on new genetic evidence, Archer et al. (2019) recently formally described fin whales in the North Pacific as a new subspecies, B. p. velifera, which was accepted by the Committee on Taxonomy (Taxonomy, 2020). Additionally, Jiménez López et al. (2019) presented evidence supporting previous findings that fin whales from the Gulf of California (GoC) are a small year-round resident population isolated from the rest of the NP (Bérubé et al., 2002;Urbán et al., 2005;Nigenda-Morales et al., 2008;Rivera-León et al., 2019). Genetic analyses between the Northern and Southern Hemisphere show a strong phylogeographic structure among NP, NA, and SH fin whale populations (Archer et al., 2013(Archer et al., , 2019Cabrera et al., 2019). However, the reduced sample size, restricted geographic coverage in the SH, and inability to positively associate morphological characteristics of putative pygmy fin whales with those sampled impeded the assessment of potential phylogeographic structure within this area and evaluation of the validity of the putative B. p. patachonica sub-species (Archer et al., 2013(Archer et al., , 2019. The occurrence of fin whales has recently become regular in mid-latitudes of Southeastern Pacific (SEP), in particular along the Chilean coast between 23 • S and 29 • S in spring and summer (Pérez et al., 2006;Pacheco et al., 2015;Toro et al., 2016;Sepúlveda et al., 2018). Foraging areas for the species have been identified within this zone (Pérez et al., 2006;Toro et al., 2016;Sepúlveda et al., 2018). If pygmy fin whales effectively inhabit the Southern Ocean from around 35 • S (Burmeister, 1865) to ∼55 • S, and from the Equator southwards along the Pacific coast of South America (Clarke, 2004), fin whales found in these lower SH latitudes may potentially correspond to the B. p. patachonica subspecies.
Here, we present a global phylogeographic study of fin whales, increasing the geographical coverage of previous studies by adding 37 new sequences from the Southeastern Pacific as well as longer sequences from the Gulf of California. We evaluate the genetic structure among fin whale populations within the Southern Hemisphere, the existence of two SH subspecies, the degree of genetic structure between North Pacific and Southeastern Pacific populations, and differentiation the GoC and the wider North Pacific. Finally, we discuss the existence of the putative pygmy fin whale sub-species B.p. patachonica.

Sample Collection and Laboratory Procedures
Skin samples of fin whales were obtained from free-swimming whales by biopsy darting (Harlin et al., 1999). Biopsy samples were collected under permits 1502/2013, 1803/2015, and 2982/2015 from the Chilean Undersecretariat of Fisheries and Aquaculture and with approval and permits issued by the Mexican Wildlife Agency (Dirección General de Vida Silvestre, Subsecretaría de Gestión para la Protección Ambiental, Secretaría del Medio Ambiente y Recursos Naturales). All methods were carried out in accordance with relevant guidelines and regulations. Experimental protocols were approved by the Bioethical Committee at the Institute of Ecology and Biodiversity, Chile (IEB). Potentially re-sampled individuals among the GoC samples were ruled out by microsatellite analysis comparisons (13 loci, Cervus v3.0, Kalinowski et al., 2007), while genetic sampling in the SEP was guided by simultaneous photoidentification of individuals.
mini kit (Qiagen, Valencia, CA, United States) according to the manufacturer's protocol. The primers used were the same as previously mentioned, but amplifications conditions and cycling were as follows: denaturation at 95 • C for 3 min, followed by 34 cycles at 95 • C for 30 s, 50 • C for 45 s, 72 • C for 45 s and a final extension at 72 • C for 10 min. Reactions were performed in a total volume of 50 µl, with 2.5 mM MgCl 2 , 600 mM each deoxynucleotide triphosphate (dNTP), 4 pmol each primer, 1 unit of Taq DNA polymerase (Sigma-Aldrich, St Louis, MO, United States), 1 µl 10 mg/ml BSA and approximately 50-100 ng of template DNA. Amplified products were sequenced on an automated sequencer (ABI 3100, Applied Biosystems, Foster City, CA, United States). Sequence data were edited in SEQUENCHER v4.1 (Gene Codes Corporation, Ann Arbor, MI, United States).
The eight mitogenome sequences from Cabrera et al. (2019) were trimmed to our D-loop segment length of 624 bp, and included in the construction of the haplotype network and the phylogenetic tree. However, these samples were not used for the genetic diversity and structure analyses. As explained by the authors, these sequences were deliberately selected for mitogenome amplification, which is inconsistent with a random sampling scheme.

Genetic Analyses
For all strata, we estimated the following genetic diversity indices: number of segregating sites (k), number of haplotypes (h), haplotype diversity (Hd), mean number of differences between two random sequences ( ) and nucleotide diversity (π). Analyses of phylogeographic structure ( ST ) that evaluate genetic structure considering both haplotype frequencies and genetic distance between haplotypes were conducted in Arlequin v3.5.2 (Excoffier et al., 2005) with 1,000 permutations and a significance level of 0.05. The haplotype network was constructed in Hapview 1 , a software package that uses phylogenetic trees to construct haplotype genealogies. For this, a neighbor-joining tree was constructed in MEGA v7 (Kumar et al., 2016) to serve as input file for the reconstruction of the haplotype network in Hapview.
To estimate the pairwise migration rates between geographic strata [NP, GoC, SH (SEP + SO) and NA], the software IMa2 was used (Hey, 2010). This was conducted only between groups that either shared haplotypes or showed discrepancy between genetic and geographic clusters (i.e., SH-NA, SH-NP, and GoC-NP). We carried out 10 preliminary runs in the M-mode of the software (Markov Chain Monte Carlo; MCMC mode), in order to select the set of priors that best improved mixing and convergence of the MCMC. These priors were then used to determine migration patterns based on calculations of effective population size ( 1, 2, and ancestral a, = 1,000) and splitting time (t = 100). We performed 10 × 10 8 MCMC steps using the HKY mutation model and sampled every 100 iterations, with an initial burnin of 10%. Once MCMC convergence was achieved in the test runs (i.e., same results obtained with different runs through visual inspection of trace-plots and histograms of each parameter), simulated genealogies were used in the Load Tree mode (L-mode) to calculate the log of the maximum-likelihood and 95% HDP intervals for migration parameters. These were compared to a no migration scenario (null model), using a likelihood ratio-test.
For the phylogenetic reconstruction, we used multiple intraand extra-specific calibration points, as proposed by Ho et al. (2008). This was done under the premise that if an estimate of evolutionary rate is calibrated using only an external calibration point, the intraspecific rate will be underestimated and the resulting divergence time will be overestimated. Thus, we performed a phylogenetic reconstruction of the fin whale haplotypes combining a phylogenetic mutation rate with a population substitution rate (Ho et al., 2008). For the former, we estimated divergence times among the three main clades using the phylogenetic rate, as it relates to the rate at which mutations are fixed among clades. For the latter, we used a population substitution rate to estimate the TMRCA of each clade, as it considers the rate at which new haplotypes appear in the clade.
Divergence times among the main lineages were first estimated in BEAST v2.6.0 (Bouckaert et al., 2014), using a relaxed molecular clock with a lognormal distribution and a phylogenetic substitution rate of 1.5% per million years (Ho et al., 2005). Additionally, given their close phylogenetic relationship (Sasaki et al., 2005;Archer et al., 2013), we used the humpback whale, Megaptera novaeangliae, as outgroup (GenBank accession mitochondrial DNA, complete genome Sequence ID: AP006467.1 Length: 16,398 bp), with an estimated divergence time from fin whales of 11 million years (Marx and Fordyce, 2015) as the primary calibration point. Bayesian Inference posterior probabilities were estimated using the Metropolis coupled MCMC mode, running four chains for 3 × 10 8 iterations, with trees sampled every 1,000 steps. The first 20% of the trees were discarded as burn-in and posterior probabilities were estimated as the fraction of trees in the clades of interest. Finally, posterior probability density was summarized as a maximum clade credibility tree using TREEANNOTATOR v.2.6.1 2 and visualized using FIGTREE v.1.4 3 . The time to the Most Common Recent Ancestor (TMRCA) for each clade was estimated BEAST. For this, three independent Bayesian MCMC analyses with 100 × 10 6 iterations were carried out, considering a population substitution rate of 5.2% per million years (Alter and Palumbi, 2009) and a relaxed molecular clock analysis with an uncorrelated lognormal (ucln) distribution. Results of multiple runs were combined using LogCombiner v2.6.0 (Rambaut and Drummond, 2010), and convergence and TMRCA values were estimated with TRACER v1.6 (Rambaut et al., 2018). The final tree was constructed using the results of both approaches with the Chronos function in ape package (Paradis, 2013) in R v 4.0.5 (R Core Team, 2020), using a sequential secondary age calibration with the time of divergence of the main clades and TMRCA of each lineage.

Genetic Diversity
All analyses were performed considering a final fragment of 624 bp of the mitochondrial DNA control region. We obtained 37 sequences from the SEP samples, which were composed of 26 different haplotypes. Of these haplotypes, 17 (65%) occurred only once, 7 (27%) were shared by two SEP individuals, and one haplotype (8%) was shared by three SEP individuals. Additionally, these samples presented high haplotype (Hd = 0.980) and nucleotide diversities (π% = 0.81). Compared to the rest of the SH, in a similar number of samples, the SEP had a lower number of haplotypes (SEP: n = 37, h = 26; rest of SH: n = 42, h = 40). SEP shared six haplotypes with the rest of SH and had similar haplotype diversities, only slightly lower in the SEP (Hd = 0.980 vs. 0.997) as well as nucleotide diversity (π% = 0.81 vs. 0.88). Overall, sequences from the SEP were mixed with those from other SH localities, with no separation visually evident in the haplotype network (Figure 2).
Comparing these results with the global control region diversity patterns of the species, similar high diversity was found in four of the five ocean basins (NP, NA, SEP, and SO, Hd range = 0.957-0.997) ( Table 1). Haplotypic diversity in the GoC was approximately three times lower than in these regions (Hd = 0.314). This low diversity is the result of a single haplotype representing 93 of the 112 (83%) samples.

Phylogeographic Structure
Three main haplogroups can be distinguished in the resulting haplotype network (Figure 2). Each group is associated with a Hemisphere, haplotypes sampled in GoC were associated with the NP haplogroup. However, most of the GoC haplotypes were found at higher frequencies than in other areas, contrasting with the high diversity patterns observed in the NP, NA and SH. Fortyeight haplotypes were identified among the 93 NP individuals, while only five haplotypes with few mutational steps among them (from 1 to 4) were found among the 112 GoC samples (Figure 2). High and significant phylogeographic structure was detected among all geographic areas considered here (Figure 3), except within the Southern Hemisphere ( ST = 0.007, p = 0.157). The highest ST values were obtained comparing the GoC population to NA and SH ( ST = 0.800 and 0.656, respectively). Phylogeographic structure indices among NP, NA and SH were lower, ranging from 0.305 to 0.441, while ST between NP and GoC showed the lowest value (0.278). This value was highly significant, probably due to the higher frequency of three GoC haplotypes (Figure 3).
Regarding migration patterns, GoC-NP was the only pair that showed a statistically significant and asymmetrical migration rate. A statistically significant migration rate was obtained from the GoC to the NP, yet even though the migration rate in the opposite direction (from the NP to the GoC) showed non-zero values, it was not significantly different from a non-migration model (Figure 4). For the SH-NP and SH-NA pairs, near zero and non-statistically significant migration rate values were obtained.
The phylogenetic tree ( Figure 5) shows a clear assortment of haplotypes according to their geographic origin (NA, NP, SH), and as previously observed in the haplotype network (Figure 2), five haplotypes sampled in the NP and three in the NA were placed in the predominantly SH clade. Using sequential calibration points, the divergence times estimated among the main clades of fin whales suggest the NP linage as the first cladegroup to diverge, around 1.8 ± 0.4 million years ago (Mya), followed by the divergence between NA and SH linages around 1.2 ± 0.3 Mya. TMRCA estimates for the main clades were NA: 0.14 Mya (95% HDP = 0.06-0.24 Mya); NP: 0.13 Mya (95% HDP = 0.06-0.25 Mya) and SH: 0.17 Mya (95% HDP = 0.09-0.29 Mya, Figure 5).

DISCUSSION
The addition of new control region sequences from the waters off of Chile to the global data set of fin whale (Balaenoptera physlaus) sequences has allowed us to contribute new insights to our growing understanding of phylogeographic structure within this species. These sequences provide the first estimates of genetic diversity for a previously unsurveyed area in the Southern Hemisphere. Our analysis showed that SEP fin whales have high levels of genetic diversity similar to, albeit slightly lower than, that in the rest of the Southern Hemisphere. The haplotypes, of which six were shared, appeared homogeneously distributed in the network with no clear structure within the SH. This was further supported by the statistically non-significant phylogeographic structure observed between samples from the SEP and the rest of the Southern Hemisphere. Diversity within the SEP was similar to that in the neighboring North Pacific. Nevertheless, except for five haplotypes nested within the Southern Hemisphere, the NP haplotypes were more clearly separated in the network.
Our expanded dataset also confirms previously described patterns of global fin whale differentiation (Archer et al., 2013(Archer et al., , 2019Cabrera et al., 2019). The three previously described fin whale maternal linages, each corresponding to an ocean basin (NA, NP + GoC, and SH) were also recovered in our phylogenetic analyses. This is notable because in previous analyses, the SH was primarily represented by a set of 40 samples collected from the same area in the Southern Ocean over a period of 2 years (Archer et al., 2013). The addition of a comparable number of SEP samples, from the opposite side of the Southern Hemisphere does not appreciably change the structure of previously published phylogenetic trees (Archer et al., 2013;Cabrera et al., 2019). Estimated divergence times suggest that the diversification of the three maternal linages (1,8 ± 0,4 Mya) occurred during the Lower (Early) Pleistocene (dated from 2.58 to 0.77 Mya; Gibbard, 2015). This estimate is concordant with Archer et al. (2013), who estimate a divergence time of approximately 2 Mya. This overall pattern of divergence may be part of a complete oceanographic re-organization and climate variation that occurred at that time, including the gradual closure of the Panama Seaway, between 13 and 2.5 Mya (Lunt et al., 2008). The complete process promoted species diversification not only in cetaceans (Rosenbaum et al., 2000;Pastene et al., 2007), but also in many pelagic marine organisms (Jones and Hasson, 1985;Lima et al., 2020). The non-continuous distribution of the species across Southern and Northern Hemispheres may relate to the equatorial warm-water belt that operates as a variable but efficient barrier to dispersal, leading to antitropical distribution patterns (Davies, 1963).
Strong genetic differentiation was observed between the NA and NP, with each ocean basin being characterized by different haplogroups as previously described (Archer et al., 2013(Archer et al., , 2019. Given that continents separate ocean basins and prevent interoceanic dispersal in the NH, this pattern of structure has been reported for many taxa, including humpback whales (Jackson et al., 2014), sei whales (Balaenoptera borealis; Baker et al., 2004;Huijser et al., 2018), as well as odontocetes, including bottlenose dolphins (Tursiops truncatus; Dowling and Brown, 1993) and orcas (Orcinus orca; Morin et al., 2015). Additionally, absence of recurrent gene flow was observed between the two Northern Hemisphere fin whale subspecies (NP and NA) and the SH populations. However, few haplotypes sampled in the North Pacific nested genetically in the Southern Hemisphere clade, suggesting the existence of discrete dispersal events from SP to NP as previously proposed by Archer et al. (2013). Similarly, the finding of three haplotypes in NA that nested in the SH clade may indicate the existence of dispersal events from SH to NA, as suggested by Cabrera et al. (2019). However, asymmetrical migration was not supported by our IMa2 analyses and may indicate that migrations from SH to NA and NP are the result of rare or even exceptional dispersal events from the Southern to Northern Hemisphere. This is also supported by close genetic relationship between some of these geographically misplaced haplotypes, corresponding to small post-migration lineages. Although we note that incomplete lineage sorting cannot be discarded as a possible explanation of this pattern. Nonetheless, other occasional gene flow between baleen whale populations across hemispheres has been reported at least for southern minke whales (Balaenoptera bonaerensis; Glover et al., 2010) and humpback whales. These species perform synchronous seasonal migrations to feeding areas, which differ in timing between hemispheres and thus precluding them from mixing (Jackson et al., 2014).
Within the North Pacific, analyses of additional Gulf of California (GoC) fin whale sequences in this study were consistent with the known significant genetic differentiation and isolated nature of this population (Bérubé et al., 2002;Rivera-León et al., 2019). Haplotype diversity within the GoC was the lowest of all regions examined (Hd = 0.314). However, this value was twice as high as that previously reported by Rivera-León et al. (2019) (Hd = 0.15). Notwithstanding the higher number of samples used in this previous study (n = 253 vs. n = 112), this difference can be explained by the longer sequences used in our study (624 vs. 280 base pairs), which also generated a new GoC haplotype. We estimate the founding of the GoC from the North Pacific around 9,000 years ago 550 years;ESS 12,688). These results refine the TMRCA estimates and confidence interval between NP and GoC obtained by Rivera-León et al. (2019) of approximately 1,300 years ago (95% HPD: 0-2,700 kya; ESS: 559,148). This new estimate is now also more consistent with a sea-level rise associated with the deglaciation process, about 15,000 years ago (Keigwin and Jones, 1990). Paleohabitat modeling has shown strong variations in estuarine habitats since the Last Glacial Maximum (LGM), approximately 20 kya, when the sea level was approximately 130 m lower. The colonization of the Gulf of California by marine species may be related to this sea level rise, as mentioned by Dolby et al. (2016), who identified and assessed a sea-level driven recolonization process after the LGM along the coast of California (United States) and Baja California (Mexico). The colonization of the GoC by fin whales would have been followed by a switch from a migratory to resident behavior, sustained by the permanently highly productive ecosystem of the GoC (Mercado-Santana et al., 2017). The year-round foraging activity recorded for fin whales in the GoC (Jiménez López et al., 2019) is evidence of the low movement of whales from this region. Similarly, a suspension of the southward migration in humpback whales in Cape Columbine, South Africa, was observed in response to locally abundant prey (Best et al., 1995). Also, an exception to the seasonal migrations of humpback whales by year-round sightings as well as records of feeding and breeding behavior in a limited area of the Arabian sea has been reported (Pomilla et al., 2010). However, although such a residence pattern may explain the genetic structure detected between NP and GoC, the geographic distribution of the shared, closely related haplotypes, suggests the existence of low and asymmetric gene flow from the GoC to the NP. This migration pattern and the distribution of haplotypes between the two regions fits a continent-island migration model, where migrants coming from a small population have no noticeable effect on the gene frequencies of the large one (Nagylaki, 1992;Bürger and Akerman, 2011). Regardless, our migration rate estimates should be taken with caution, since they are derived from a single mtDNA marker. Further studies are needed in order to confirm patterns of asymmetrical migration in NP. Nevertheless, even if no gene flow was detected between NP and GoC, one haplotype detected in the GoC clearly has a genetic affinity with NP haplotypes, consistent with rare or even exceptional migration events from the NP to the GoC, as previously proposed by Rivera-León et al. (2019).
One of the most important contributions of this study is that our findings challenge the validity of the pygmy fin whale sub-species, B. p. patachonica (Burmeister, 1865), which has previously been called into question (Archer et al., 2019). Pygmy fin whales are described as being relatively short (total length approximately < 20 m), dark in color with black baleen, and occurring from approximately 55 • S to the Equator along the South American Pacific coast (Clarke, 2004). We could not positively assign any of the SH samples to either B. p. quoyi or B. p. patachonica, as we did not have total body length measurements, or observed significantly darker individuals, nor a whale with black baleen among the over 50 fin whales sighted around the area. Additionally, aside from the type specimen residing in the Museo Argentino de Ciencias Naturales Bernardino Rivadavia (MACN) in Buenos Aires, Argentina, there are no known specimens that have previously been positively identified as pygmy fin whales for morphological comparison. However, given the time scale and geographic scope of our SH samples, we would have expected to have sampled at least some pygmy fin whales off Chile. Our samples were obtained during the austral summer season, which is when fin whales of the subspecies B. p. quoyi are expected to be southbound, migrating to sub-Antarctic and Antarctic waters (Mackintosh, 1942). Moreover, five of the six individuals that were satellite tagged by Sepúlveda et al. (2018), which were also included in our study, remained in lower latitudes during the summer months (sampled period), within the latitudinal distribution range that is proposed for pygmy fin whales (Clarke, 2004). Additionally, Cooke (2018) reviewed information on body size of fin whales collected in modern whaling activity off the northern (∼20 • S), central (30-38 • S), and southern (44 • S) coasts of Chile. The author mentions that these whales were smaller than those of higher latitudes, potentially meeting the small size criterion pygmy fin whales of Clarke (Clarke, 2004).
Although significant control region differentiation was obtained in all other comparisons, none was detected between SEP fin whales and those from the rest of the Southern Hemisphere. Thus, it is evident that the differentiation that might exist between regions within the SH is not at the same level as that between subspecies within the Northern Hemisphere, between Hemispheres, or even between the Gulf of California and the greater North Pacific suggesting the existence of a single evolutionary unit in SH. Similarly, low genetic differentiation and high gene flow among regions within the SH has been reported for humpback whales (Jackson et al., 2014). By contrast, Bryde's whales (Balaenoptera edeni) show strong genetic differentiation between the South Pacific and South Atlantic oceans, however this species is not distributed further south than approximately 40 • S (Pastene et al., 2015) and would be effectively separated by the South American continent.
Based on the absence of any detectable structure within the Southern Hemisphere and the insufficient support of the initial subspecies description by Clarke (2004) (see review in Archer et al., 2019), it appears that the pygmy fin whale, B. p. patachonica, may not be a valid subspecies. Alternatively, it is possible that pygmy fin whales are not present in the study area or we have not yet sampled it. Thus, until evidence that positively supports the validity of B. p. patachonica can be produced, such as genetic, ecological, behavioral, or morphological data (Crandall et al., 2000), we propose that all Southern Hemisphere fin whales, including those present in middle latitudes of southeastern Pacific belong to the earlier described subspecies B. p. quoyi (Tomilin, 1946). The smaller size of fin whales along the central Chilean coast may be because these individuals are juveniles of B. p. quoyi, and the different migration behavior they present could be related to energy saving, as proposed by Cooke (2018). The existence of a single taxa within the Southern Hemisphere is important for the assessment of fin whales, primarily relevant to the long-term management issues defining conservation priorities (Moritz, 1994) being crucial in conservation legislation (Robertson et al., 2014).

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: GenBank (accession: MW266395-MW266537).

ETHICS STATEMENT
The animal study was reviewed and approved by the Permits 1502/2013, 1803/2015, and 2982/2015 from Chilean Undersecretariat of Fisheries and Aquaculture and with approval and permits issued by the Mexican Wildlife Agency (Dirección General de Vida Silvestre, Subsecretaría de Gestión para la Protección Ambiental, Secretaría del Medio Ambiente y Recursos Naturales).