In contrast to muscle tissue, the FA profiles of sub-dermal tissue smaller than 40 mg DW became highly variable, indicating a minimum reliable tissue quantity of 40 mg DW (= 184 mg WW), which is more than three times the minimal requirement for muscle. This is potentially due to the difference in PL concentration between the two types of tissue of the lipid profile. Combined with the lower lipid content, the lower relative PL contribution in the sub-dermal tissue may explain the comparatively larger minimum sub-dermal tissue quantity, as the ST, which are found in higher abundance in the sub-dermal tissue, do not contribute to the FA pool. This larger minimum tissue quantity required for sub-dermal tissue compared to muscle may limit the applicability of many aforementioned non-lethal biopsy methods. For example, the biopsy method yielding the second highest tissue volume provided only 80– 172 mg WW of sub-dermal tissue (Daly and Smale, 2013), which is not sufficient for reliable FA analysis. Only the Reeb and Best's dart head (Reeb and Best, 2006) which retained an average of 0.35 cm3 of sub-dermal tissue when trialed by Jaime-Rivera et al. (2013), obtained potentially suitable tissue quantities. Furthermore, biopsies from small elasmobranchs are unlikely to yield sufficient tissue, as the thickness of the sub-dermis is greatly reduced. For example, sub-dermal tissue layers in atlantic sharpnose shark, scalloped hammerhead and dusky smooth-hound sharks ranged 0.02–0.16 cm (Motta, 1977), compared to white sharks averaging 1.1 cm (Jaime-Rivera et al., 2013) and whale sharks exceeding 2 cm (Rohner et al., 2013).

### Degradation

The consistently low levels of FFA in muscle and sub-dermal tissue throughout the degradation trial contrasts with findings across marine taxa, which highlight large increases in FFA from enzymatic hydrolysis of several non-polar lipid classes (Fernández-Reiriz et al., 1992; Kaneniwa et al., 2000; Losada et al., 2005). The difference between our findings and the pervasive trends in previous studies may be attributable to species- and taxa-specific enzymatic processes. Rudy et al. (2016) and Kaneniwa et al. (2000) hypothesized that total lipid content drove the species-specific differences in the level of observed lipid class and FA degradation amongst teleost species, with the "fatty" fish most susceptible. Compared with the six teleosts assessed in Rudy et al. (2016), white sharks were orders of magnitude leaner, with muscle containing 0.6% lipid WW and sub-dermal tissue 0.4% lipid WW (vs. 10.3–2.9% WW in teleosts). The low lipid content may explain the lack of discernable lipid class degradation across both tissues and the comparative stability in FA profiles within the sub-dermal tissue. Given the aim of determining the functional limitations of using elasmobranch specimens not immediately frozen, for example from fisheries bycatch and shark mitigation measures, our results indicate that lipid classes from muscle and subdermal tissues are not convoluted by degradation within a 4 day period.

The lipid-poor sub-dermal tissue also showed no discernible shift in FA profile or level of dispersion through exposure to ambient temperature for 4 days. However, the FA profiles derived from muscle tissue immediately changed, with a decrease in dispersion observed after 24 h, potentially compromising the ability to distinguish between individual samples. This advocates for exploring the use of sub-dermal tissue over muscle in situations when samples have been left at ambient temperature, and should be the subject of controlled feeding trails to assess the capacity for sub-dermal tissue to reflect diet. Our earlier findings, however, highlights that such FA profiles based on sub-dermal layers cannot be directly compared to FA profiles from muscle and that this discrepancy should be accounted for.

Muscle segments stored at −20◦C showed significant FA profile shifts in both assessment periods, highlighting concerns regarding the capacity to accurately use archived samples. Results in this study suggest that although there may be some level of FA degradation, the time frame at which this occurs and processes involved remains unclear. It is also plausible that the difference in the 3–6 months group is not driven by the time spent in the freezer, but by the influence of unassessed biotic factors (e.g., individual's state of maturity, sex, season of capture). The comparison of FA profiles from archived samples stored for 1– 16 years did not provide further clarification and showed no clear differences in FA profiles. Furthermore, neither trial's FA profiles decreased in dispersion, a pattern characteristic of FA degradation in the ambient temperature trial. Regardless of the degradation that might be occurring through long-term storage, differences between locations (NSW vs. QLD) remained, further suggesting that frozen samples may retain viable and indicative FA signatures.

The shift in the relative proportions of individual FAs of the muscle tissue illustrates the complex nature of FA degradation at both 20 and −20◦C. Our study found that SFA, driven primarily by 18:0, can remain constant during some time periods, but also decreased drastically through other periods. The MUFA, unchanged at 20◦C, demonstrated some resistance to degradation, with no shifts in either individual MUFA, or the PMUFA. Unexpectedly, they showed variable patterns of alteration in the early month of storage, suggesting that they are prone to degradation at −20◦C, consistent with findings across other taxa (**Table 5**, e.g., teleosts in Rudy et al., 2016 and octopus in Gullian-Klanian et al., 2017). PUFA are more reactive owing to their numerous double-bonds and are especially prone to degradation (Refsgaard et al., 1998; Paola and Isabel, 2015; Rudy et al., 2016; Gullian-Klanian et al., 2017). However, shifts in relative levels of PUFA of white sharks, including key dietary indicators 22:6ω3 (DHA) and 20:5ω3 (EPA), were only distinguishable in the ambient temperature trial, and not in either the short- or long-term −20◦C analysis (with the exception of 18:2ω6). Additionally, the polyene index (EPA+DHA/16:0), a well-established metric for tissue degradation, thought to be ubiquitous across taxa (Jeong et al., 1990; Paola and Isabel, 2015),showed no decrease across any trials (**Table 5**). The present study shows that white shark muscle PUFA might not show the stark degradation seen in the muscle tissue of other species. Given the relative importance of PUFA, as essential FAs and key dietary markers, these findings suggest that elasmobranch samples may retain these key FAs throughout extensive storage at −20◦C.

#### CONCLUSION

Our findings indicate that muscle and sub-dermal tissue contain distinct FA profiles and differing individual FAs, many of which are key trophic indicators. As such, these tissues are not directly comparable. They may, however, present complementary trophic information reflecting differing time frames, providing the opportunity to garner additional information from non-lethal biopsies. The minimum tissue amount for sub-dermal tissue was 40 mg DW (184 mg WW), whereas muscle samples as small as 12 mg DW (equating to 49 mg WW) retained consistent FA profiles. This makes FA analysis an ideal tool for elucidating trophic ecology of rare or endangered elasmobranchs for which lethal sampling is inappropriate. Degradation of muscle tissue occurred within the first 24 h at ambient temperature, unlike sub-dermal tissue, which revealed no discernible degradation across 4 days. As such, the use of deceased organisms, from shark mitigation strategies, by-catch, or beach strandings should be undertaken with caution, ensuring that preservation occurs within 24 h. Muscle tissue appears to retain viable and indicative FA signatures across long periods of frozen storage (up to 16 years), advocating for the use of archived samples, especially in cases where sampling opportunities are rare or opportunistic. Overall, lipid class and FA analysis can be reliably assessed from small tissue quantities derived from minimally invasive, nonlethal biopsies, deceased elasmobranchs preserved within 24 h and archived samples, proving a robust toolset for elucidating the trophic ecology of rare and endangered wildlife.

#### REFERENCES

#### AUTHOR CONTRIBUTIONS

LM, CH, CB, HP, PN, JW, and BB: conceived and designed the experiments; BB, JW, CB, and CH: provided tissue samples; LM: performed the experiments, analyzed the data with the help of CH, CB, PN, and HP; LM wrote the manuscript with the advice of CH, CB, HP, PN, JW, and BB; all authors provided editorial advice.

#### FUNDING

This work was supported by a research grant from the Save Our Seas Foundation (grant number RPF14/553).

#### ACKNOWLEDGMENTS

This work was supported by the Save Our Seas Foundation. We thank the Queensland Department of Agriculture and Fisheries Shark Control Program as part of the QLD large shark tagging research program and the QLD white shark tracking program for providing tissue samples. We also thank the New South Wales Department of Primary Industries Shark Meshing Program and the South Australian Research and Development Institute for collecting and storing tissue material. The white shark tourism operators; Andrew Fox and the Fox Shark Research Foundation, Calypso Star Charters and Adventure Bay Charters staff and crew provided ongoing support and in-kind contributions, facilitating sample collection within South Australia.

into the feeding ecology and ecophysiology of a complex top predator. PLoS ONE 9:7877. doi: 10.1371/journal.pone.0097877

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2017 Meyer, Pethybridge, Nichols, Beckmann, Bruce, Werry and Huveneers. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# iDNA at Sea: Recovery of Whale Shark (Rhincodon typus) Mitochondrial DNA Sequences from the Whale Shark Copepod (Pandarus rhincodonicus) Confirms Global Population Structure

Mark Meekan1 , Christopher M. Austin2, 3, 4, 5, Mun H. Tan2 , Nu-Wei V. Wei 5 , Adam Miller 2 , Simon J. Pierce6 , David Rowat 7 , Guy Stevens 8 , Tim K. Davies 9 , Alessandro Ponzo10 and Han Ming Gan2, 3, 4 \*

*1 Australian Institute of Marine Science, Indian Ocean Marine Research Centre (MO96), Crawley, WA, Australia, 2 Centre for Integrative Ecology, School of Life and Environmental Sciences, Deakin University, Geelong, VIC, Australia, 3 School of Science, Monash University Malaysia, Petaling Jaya, Malaysia, 4 Genomics Facility, Tropical Medicine and Biology Platform, Monash University Malaysia, Petaling Jaya, Malaysia, 5 Research Institute for the Environment and Livelihoods, Charles Darwin University, Darwin, NT, Australia, 6 Marine Megafauna Foundation, Truckee, CA, United States, 7 Marine Conservation Society, Seychelles, Seychelles, 8 The Manta Trust, Dorchester, United Kingdom, 9 MRAG Ltd., London, United Kingdom, 10 Large Marine Vertebrates Research Institute Philippines, Bohol, Philippines*

The whale shark (*Rhincodon typus*) is an iconic and endangered species with a broad distribution spanning warm-temperate and tropical oceans. Effective conservation management of the species requires an understanding of the degree of genetic connectivity among populations, which is hampered by the need for sampling that involves invasive techniques. Here, the feasibility of minimally-invasive sampling was explored by isolating and sequencing whale shark DNA from a commensal or possibly parasitic copepod, *Pandarus rhincodonicus* that occurs on the skin of the host. We successfully recovered mitochondrial control region DNA sequences (∼1,000 bp) of the host via DNA extraction and polymerase chain reaction from whole copepod specimens. DNA sequences obtained from multiple copepods collected from the same shark exhibited 100% sequence similarity, suggesting a persistent association of copepods with individual hosts. Newly-generated mitochondrial haplotypes of whale shark hosts derived from the copepods were included in an analysis of the genetic structure of the global population of whale sharks (644 sequences; 136 haplotypes). Our results supported those of previous studies and suggested limited genetic structuring across most of the species range, but the presence of a genetically unique and potentially isolated population in the Atlantic Ocean. Furthermore, we recovered the mitogenome and nuclear ribosomal genes of a whale shark using a shotgun sequencing approach on copepod tissue. The recovered mitogenome is the third mitogenome reported for

#### Edited by:

*Rob Harcourt, Macquarie University, Australia*

#### Reviewed by:

*Celine Frere, University of the Sunshine Coast, Australia Gail Schofield, Deakin University, Australia*

\*Correspondence:

*Han Ming Gan<EMAIL_ADDRESS>
#### Specialty section:

*This article was submitted to Marine Megafauna, a section of the journal Frontiers in Marine Science*

Received: *08 October 2017* Accepted: *07 December 2017* Published: *18 December 2017*

#### Citation:

*Meekan M, Austin CM, Tan MH, Wei N-WV, Miller A, Pierce SJ, Rowat D, Stevens G, Davies TK, Ponzo A and Gan HM (2017) iDNA at Sea: Recovery of Whale Shark (Rhincodon typus) Mitochondrial DNA Sequences from the Whale Shark Copepod (Pandarus rhincodonicus) Confirms Global Population Structure. Front. Mar. Sci. 4:420. doi: 10.3389/fmars.2017.00420*

**42**

the species and the first from the Mozambique population. Our invertebrate DNA (iDNA) approach could be used to better understand the population structure of whale sharks, particularly in the Atlantic Ocean, and also for genetic analyses of other elasmobranchs parasitized by pandarid copepods.

Keywords: eDNA, sharks, minimally-invasive sampling, copepod, control region, population genetics, parasite, commensal

#### INTRODUCTION

The term environmental DNA (eDNA) refers to DNA extracted from cells that are not collected directly from a target organism, but are obtained from the environments they inhabit such as oceans, river water, soil, and air (Ficetola et al., 2008; Fonseca et al., 2010; Andersen et al., 2012; Sigsgaard et al., 2016). Invertebrate-derived DNA (iDNA) is an offshoot of this approach that involves the extraction of genetic material of animals via the flesh-eating or haematophagous invertebrates that parasitise them (Schnell et al., 2015; Schubert et al., 2015; Lee et al., 2016). To date, most iDNA studies have focused on terrestrial vertebrates and have extracted host DNA from insects, ticks, or leeches. There have been no analogous studies in marine environments, despite the potential usefulness of the approach for sampling large marine megafauna such as cetaceans, sirenians, pinnipeds, marine reptiles (sea turtles), elasmobranchs, and teleosts. For such taxa, iDNA sampling offers both ethical and practical advantages as it is minimally-invasive and thus preferable to the direct collection of blood or tissue, particularly where target species are rare or endangered.

The whale shark (Rhincodon typus), similar to many large marine vertebrates, poses significant challenges for researchers trying to understand their biology and ecology and for managers attempting to develop conservation strategies (Graham and Roberts, 2007; Rowat et al., 2009). Tagging, genetic and modeling studies suggest that individuals can disperse widely (Sequeira et al., 2013), although there is evidence of isolation between populations in the Indo-Pacific and Atlantic oceans (Vignaud et al., 2014). However, the scale at which population structure can be discerned is dependent on sample sizes, with low numbers reducing analytical power to reject the null hypothesis of a panmictic population (Castro et al., 2007). For whale sharks, low sample size tends to be a result of the rarity of the species and the difficulties in obtaining biopsies (usually of skin tissue) for genetic analyses, which requires appropriate permits and ethical approvals for invasive sampling of a species that is categorized as "Endangered" by the International Union for the Conservation of Nature and the use of trained divers and technicians. In this situation there are many benefits to the use of an iDNA approach, provided that external invertebrate parasites or commensals that harbor the intact DNA of the whale shark host can be successfully identified. Although eDNA techniques have recently been used to obtain whale shark DNA (Sigsgaard et al., 2016), iDNA still offers a major advantage because it enables haplotypes to be linked to individual whale sharks for the analysis of the genetic structures of populations.

Copepods (Crustacea, Maxillopoda, Copepoda) are an excellent candidate species for iDNA studies on marine vertebrates. Free-living forms are a dominant element of marine zooplankton and may even exceed insects in terrestrial environments in terms of sheer abundance of individuals. Less recognized is that approximately half of the nearly 30,000 described species live in parasitic or commensal associations with a diverse range of taxa, including fish and mammals (Boxshall, 2005). The order Siphonostomatoida contains largely parasitic copepods that feed on the blood, epidermal tissue, or mucus of many marine teleost fishes, sharks, and rays. Approximately 550 genera, representing nearly 40 families, are placed in the order and include economically important species such as sea lice (Brachiura) that parasitise farmed fish (Gunn and Pitt, 2012). Eleven siphonostomatoid families have been reported as symbionts of a diverse range of elasmobranchs (Dippenaar, 2009) and one family, the Pandaridae, is composed of 23 genera that include species that are ectoparasites or commensals of large sharks (Cressey and Boyle, 1978; Walter and Boxshall, 2017). Within the Pandaridae, the genus Pandarus currently has 17 recognized species of which P. rhincodonicus (Norman et al., 2000) is noteworthy as it is appears to be associated exclusively with the whale shark, where it is predominately found on the leading and trailing edges of fins and on the lips. Although it is thought to be a commensal that feeds off bacteria and other microorganisms on the skin of the shark (Norman et al., 2000), dietary studies of the species are incomplete.

Here, we demonstrate that P. rhincodonicus sampled from sites across the Indian Ocean contain sufficient DNA from their whale shark host to routinely recover mitochondrial sequences that allow analyses of the genetic structure of host populations, and the recovery of the complete whale shark mitogenome and ribosomal nuclear genes. The implications of these findings are discussed in relation to the use of copepods for iDNA studies of marine teleosts and sharks and the ecological status of Pandarus spp. as commensals.

### MATERIALS AND METHODS

#### Sampling and Sequencing

A total of 45 copepods were collected from 31 individual whale sharks sampled in Baa Atoll in the Maldives, at Ningaloo Reef, Western Australia, off Tofu Beach, Mozambique and off the coast of Mahe in the Seychelles as part of an earlier study using approved procedures under the Western Australian Department of Parks and Wildlife (WADPW) and University of Tasmania (UTAS) ethics permits (WADPW: SF009814 and SF009227; UTAS: 2255 and 2307) (Vignaud et al., 2014) (**Supplementary Table 1**). Copepods were scraped off the edges of the fins or lips using a plastic knife by a snorkeler who swam alongside an unrestrained animal, collected in an aquarium net and brought back to the vessel where they were preserved in 100% ethanol. DNA extraction was performed on whole copepod specimens using the DNeasy Blood and Tissue kit (Qiagen, Halden, Germany). Of the 45 copepods we collected, 44 were selected for the amplification of whale shark the mtDNA control region, using primers WSCR1-F and WSCR2- R, and following the PCR protocol as described by Castro et al. (2007). PCR products were purified using the Viogene Gel/PCR DNA Isolation System Kit (Viogene Biotek Corp, Taitung, Taiwan) and sequenced on an ABI 3130xl Genetic Analyzer (Applied Biosystems, Foster City, CA) located at Charles Darwin University.

The remaining copepod, isolate 366.1, was used for partial genome sequencing. Briefly, 1 µg of genomic DNA (gDNA) was sheared to 300 bp using a Covaris Focused Ultrasonicator (Covaris, Woburn, MA), and subsequently processed with Truseq DNA library prep kit (Illumina, San Diego, CA) according to the manufacturer's instructions. This was followed by nextgeneration sequencing on a fraction of a HiSeq 2000 lane (Illumina, San Diego, CA), with a run setting of 2 × 100 bp, for the initial goal of recovering the copepod mitogenome (Austin et al., 2016) and microsatellite loci (unpublished).

#### Mitochondrial Control Region Analysis

Authenticity of the DNA sequences of the mitochondrial control region of whale sharks were validated by BLAST searches, and sequences were aligned with previously generated sequences (Castro et al., 2007; Vignaud et al., 2014) using MAFFT version 7.310 with the option "–adjustdirection" activated (Yamada et al., 2016). The 5′ and 3′ ends of the alignment were manually trimmed using MEGA 6 (Tamura et al., 2013) so that each aligned control region had the same flanking sequence (**Supplementary Figure 1**). The trimmed alignment was subsequently de-gapped using Seqret (Rice et al., 2000) and clustered with cd-hit-est at a 100% sequence identity and 100% sequence length coverage cut-off as implemented using the –c 1.0 and –a 1.0 setting (Li and Godzik, 2006).

#### Analysis of Population Differentiation

The mitochondrial control regions of whale sharks generated using iDNA were combined with 613 sequences published by previous studies of whale shark populations (Castro et al., 2007; Ramírez-Macías et al., 2007; Schmidt et al., 2010; Vignaud et al., 2014). Identical iDNA sequences from the same whale shark host were removed prior to genetic analysis, to eliminate the chance of biasing estimates of population genetic differentiation (**Supplementary Table 1**). Arlequin suite version <IP_ADDRESS> (Excoffier and Lischer, 2010) was used to estimate global and pairwise φST. A median-joining haplotype network was constructed using PopArt version 1.7 (Leigh and Bryant, 2015) and to simplify network representation, only the common haplotypes (a haplotype with more than one sequence representation) were used to generate the network graph.

### Recovery of Whale Shark Whole Mitochondrial Genome and Nuclear Ribosomal Genes

Assembly used a baiting and iterative mapping approach as implemented in MITObim version 1.8 (Hahn et al., 2013) from the complete mitogenome of a whale shark sampled off Taiwan (Accession Number: NC\_023455) as the reference sequence (Alam et al., 2014). In-silico circularization and annotation of the assembled mitogenome followed (Iwasaki et al., 2013; Gan et al., 2014). To improve alignment accuracy and specificity, mitogenome coverage was calculated by mapping the raw paired-end reads to the assembled mitogenome sequence using Bowtie2 version 2.3.2 (Langmead and Salzberg, 2012) with the setting "—score-min L,0.2,−0.2." BRIG was used to visualize the mitogenome annotation and read mapping coverage (Alikhan et al., 2011). Similar Bowtie2 mapping setting was used to map the raw reads to a whale shark assembly contig containing the nuclear ribosomal genes (18S rRNA and 28S rRNA) and the read alignment was subsequently visualized using Integrative Genomics Viewer (Thorvaldsdottir et al., 2013).

### RESULTS

### Mitochondrial Control Region Analysis

A ∼1,000 bp DNA sequence fragment of the whale shark mitochondrial control region was successfully generated from the extracted gDNA from 44 copepods representing 31 sharks. Haplotype clustering confirmed a 100% sequence match for DNA sequences amplified from copepods collected from the same whale shark host (**Supplementary Data 1**). An initial alignment of the 31 newly-generated iDNA control region sequences with the 613 sequences provided by previous studies (Castro et al., 2007; Schmidt et al., 2010; Vignaud et al., 2014; Walter and Boxshall, 2017) yielded a 1,910 bp product (**Supplementary Data 2**). Trimming using 5′ and 3′ -ends produced a final alignment length of 750 bp (**Supplementary Data 3**), from which a total of 163 unique haplotypes were evident (**Supplementary Data 4**). The Ningaloo Reef population had the highest number of sampled control region sequences (n = 163) and also had the greatest number of haplotypes (h = 36) (**Figure 1A**). Interestingly, the observed haplotype diversity for the Isla Holbox population was 50% lower (h = 11) than that of Djibouti population (Vignaud et al., 2014), despite similar sample sizes (**Figures 1A,B**).

Of the 31 deduplicated iDNA sequences reported in this study, five were found to represent replicate individuals from the Ningaloo population previously sequenced by Vignaud et al. (2014) using skin samples (**Supplementary Table 1** and **Supplementary Data 4**). Despite the relatively low number of sequences reported by our study that represented sequences from new whale sharks (N = 27), our study contributed a 39 and 22% increase in the number of individuals sampled from the Mozambique and Seychelles populations, respectively (**Figure 1A**). In addition, four novel haplotypes for control region sequences of whale sharks were identified from copepod samples.

### Population Structure of the Whale Shark

Global φST was weak but significantly different from zero, indicating limited gene flow among some sampling locations. This pattern was driven by a single population, with all pairwise φST calculations associated with the Isla Holbox population differing significantly from zero (average φST of 0.2), whereas all estimates of φST did not (**Figure 2**). This implied that the population at Isla Holbox was genetically isolated from an otherwise panmictic population spread across the Indo-Pacific Ocean. Three common haplotypes occurred at similar frequency across the sampling range (**Figure 3**). A few private haplotypes occurred in the Ningaloo Reef, Philippines and Seychelles populations. Consistent with φST estimates, the population at Isla Holbox was the most genetically differentiated, as it lacked one of the three common ancestral haplotypes (H108, 110, 125, 150, and 152), and had a high frequency of others that were rare (H67, H68, and H76).

### Recovery of the Whale Shark Mitogenome and Nuclear Genes from Shotgun Sequencing of Copepod Tissue

The complete mitochondrial genome of the whale shark was successfully assembled from partial genome sequencing of an iDNA sample from a copepod (isolate 366.1; **Figure 4**). Mapping of the 14 million paired-end reads to the reconstructed whale shark mitogenome gave a mapping rate of 0.006% (900 reads) representing approximately 5× mitogenome coverage (**Figure 4**). The mitogenomic composition of the whale shark were extracted was very similar (>99% identity) to those reconstructed from sharks from Taiwanese waters (Accession codes: KC633221 and NC\_023445). A comparison of the 13 genes coding for mitochondrial proteins indicated one or two nucleotide mismatches, mostly associated with non-synonymous mutations found in the atp6, cox1, cytb, nad2, nad4, and nad5 genes. In addition, we also observed reads mapping to the whale shark 18S and 28S rRNA genes (**Supplementary Figure 2** and **Supplementary Data 5**), indicating the presence of whale shark nuclear DNA in the copepod extracted gDNA, a phenomenon that will be useful for future population genetic studies of these sharks based on nuclear markers.

### Data Deposition

Mitochondrial control regions were submitted to NCBI under the accession numbers MF872682-MF872725. Raw reads from the shotgun sequencing of the copepod were submitted to Sequence Read Archive under the run number SRR4111090 and the reconstructed whale shark mitogenome has been assigned the accession number MF872726.

### DISCUSSION

Our study is the first demonstration of the use of iDNA sampling of a marine invertebrate in order to obtain mitochondrial genetic information from an elasmobranch host. We successfully amplified mitochondrial DNA fragments of whale sharks in the size range of ∼1,000 bp from a copepod, one of the largest host DNA fragments to be recovered by an iDNA study (Schnell et al., 2012, 2015; Lee et al., 2016; Pérez-Flores et al., 2016; Rodgers et al., 2017). This suggests the presence of largely intact whale shark DNA in or on copepods, although the exact source could not be identified as copepods were analyzed whole. Targeted efforts to amplify DNA fragments of whale sharks from the intestine and epidermal layers of P. rhincodonicus will be useful to resolve this issue and clarify the role of the copepod as a commensal or parasite. We found iDNA sequences from multiple copepods sampled from the same whale shark host to be identical, suggesting copepods have a persistent, long-term association with their host shark, which contrasts with the more generalist and mobile host associations of other invertebrate ectoparasites in terrestrial and aquatic ecosystems, such as leeches, blow flies, and mosquitoes (Calvignac-Spencer et al., 2013). However, these findings should be treated with some degree of caution and a multi-locus approach using microsatellite

markers (Olson et al., 2012) might offer a more powerful means to validate these findings. As reported by Vignaud et al. (2014), our estimates of genetic structure across the Indian, Pacific, and Atlantic oceans indicated the presence of two distinct populations, one in the Indo-Pacific and the other in the Atlantic Ocean. Isla Holbox was the only location sampled in the Atlantic Ocean and the minimally-invasive iDNA approach could potentially be used to improve the coverage of sampling in this region.

We also demonstrated that complete mitochondrial genome and nuclear ribosomal RNA sequences of whale sharks could be obtained using iDNA sampling. Although not the first for the species (Alam et al., 2014; Chen et al., 2016), the sequences generated by our study represent the first for the Mozambique population of whale sharks. The successful recovery of the mitogenome and nuclear ribosomal RNA from the shotgun sequencing of P. rhincodonicus was probably due to the high copy number of mitochondrial and nuclear ribosomal genes in somatic cells, coupled with high sequencing depth (∼14 million paired-end reads). Efficiency of the recovery of whale shark mitogenomes from copepod tissue samples might be improved by sequencing gDNA extracted from the part of the copepod body that contains higher concentrations of whale shark tissue. In addition, given the intactness of whale shark gDNA (>1,000 bp) that was present in P. rhincodonicus, long range PCR of additional mitogenome regions followed by high throughput amplicon sequencing could also be explored to dramatically increase the coverage of the whale shark mitogenome for population genetic and phylogeographic applications (Deiner et al., 2017; Pavlova et al., 2017). Given that mtDNA only represents a fraction of the evolutionary history of a species due to its strict maternal inheritance and small genome size (16 kb), studies of the population genetics of whale sharks may benefit from the combination of both mtDNA and nuclear data

transcriptional orientation of each protein coding gene is indicated by the red arrows.

(Godinho et al., 2008; Schmidt et al., 2009; Vignaud et al., 2014). The presence of nuclear ribosomal RNA reads in the copepod shotgun sequencing library suggests that it is possible to amplify both nuclear DNA and mtDNA of whale sharks from the copepod and the recent publication of the first draft genome for the whale shark (Read et al., 2017) will greatly facilitate the design of a targeted next-generation sequencing panel to enhance future studies of whale shark population genetics.

The Pandaridae comprises 23 genera and more than 100 species, many of which are ectoparasites or symbionts of large sharks (Walter and Boxshall, 2017), and there are nearly 30 other families of copepods that include species that are parasites of fishes (Boxshall, 2005). These offer an opportunity to extend iDNA methods to other taxa of large marine vertebrates. Conversely, the presence of gDNA of whale sharks in the extracted DNA of P. rhincodonicus will need to be considered as a possible contaminant in future whole genome sequencing or population genomic projects that target this copepod.

#### AUTHOR CONTRIBUTIONS

MM and CA: Original concept, sample collection, analyses, interpretation of results, writing, proof-reading. MT: Bioinformatics analysis. N-WW: Control region sequencing. AM: interpretation of results, writing, proof-reading. SP, DR, GS, TD, and AP: sample collection proof-reading. HG: analyses, writing, interpretation of results, proof-reading.

#### FUNDING

Funding for this study was provided by the SeaWorld Research and Rescue Foundation, the Save our Seas Foundation, the Monash University Malaysia Tropical Medicine and Biology Platform, Monash University Malaysia School of Science, Quadrant Energy, Australian Institute of Marine Science, The Western Australian Department of Environment and Conservation, The Four Seasons Resorts Maldives. Funding for SP came from Aqua-Firma, the Shark Foundation, the Save Our Seas Foundation and two private trusts.

#### REFERENCES

### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars. 2017.00420/full#supplementary-material

Supplementary Table 1 | Sample collection information.

Supplementary Figure 1 | Visualization of aligned control region sequences (one representative from each study). Arrows indicate regions upstream and downstream of the 5′ and 3′ ends, respectively, that will be trimmed off for subsequent haplotype clustering.

Supplementary Figure 2 | Bowtie2 alignment of Illumina reads generated from copepod gDNA to the whale shark genomic region containing 18S and 28S nuclear rRNA genes. Red and blue arrow bars indicate forward and reverse pair-end reads, respectively.

Supplementary Data 1 | Cd-hit clustering of trimmed whale shark control region sequences that were amplified from copepods.

Supplementary Data 2 | Alignment of whale shark control region sequences.

Supplementary Data 3 | Trimmed alignment of whale shark control region sequences based on the conserved region depicted in Supplementary Figure 1.

Supplementary Data 4 | Cd-hit Clustering of trimmed whale shark control region sequences.

Supplementary Data 5 | Paired-end reads mapped to the whale shark genomic region containing the 18S and 28S ribosomal RNA genes.

ocean populations of the whale shark, Rhincodon typus. PLoS ONE 4:e4988. doi: 10.1371/journal.pone.0004988

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The reviewer GS declared a shared affiliation, with no collaboration,with several of the authors, CA, MT, AM, and HG, to the handling Editor.

Copyright © 2017 Meekan, Austin, Tan, Wei, Miller, Pierce, Rowat, Stevens, Davies, Ponzo and Gan. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Route Fidelity during Marine Megafauna Migration

Travis W. Horton1 \*, Nan Hauser 2 , Alexandre N. Zerbini 3, 4, 5, Malcolm P. Francis 6 , Michael L. Domeier 7 , Artur Andriolo5, 8, Daniel P. Costa9 , Patrick W. Robinson9 , Clinton A. J. Duffy 10, Nicole Nasby-Lucas 7 , Richard N. Holdaway 11 and Phillip J. Clapham3

1 Department of Geological Sciences, University of Canterbury, Christchurch, New Zealand, 2 Center for Cetacean Research and Conservation, Avarua, Cook Islands, 3 National Marine Mammal Laboratory, Alaska Fisheries Science Center, National Marine Fisheries Service, NOAA, Seattle, WA, United States, 4 Cascadia Research Collective, Olympia, WA, United States, 5 Instituto Aqualie, Juiz de Fora, Brazil, 6 National Institute of Water and Atmospheric Research, Wellington, New Zealand, 7 Marine Conservation Science Institute, Fallbrook, CA, United States, 8 Laboratório de Ecologia Comportamental e Bioacústica, ICB, Universidade Federal de Juiz de Fora, Juiz de Fora, Brazil, 9 Ecology and Evolutionary Biology, University of California, Santa Cruz, Santa Cruz, CA, United States, 10 Department of Conservation, Auckland, New Zealand, 11 Palaecol Research Ltd., Christchurch, New Zealand

#### *Edited by:*

Lars Bejder, Murdoch University, Australia

#### *Reviewed by:*

Phil Bouchet, Centre for Marine Futures, School of Animal Biology, University of Western Australia, Australia Nuno Queiroz, University of Porto, Portugal

*\*Correspondence:*

Travis W. Horton<EMAIL_ADDRESS>
#### *Specialty section:*

This article was submitted to Marine Megafauna, a section of the journal Frontiers in Marine Science

*Received:* 25 September 2017 *Accepted:* 11 December 2017 *Published:* 21 December 2017

#### *Citation:*

Horton TW, Hauser N, Zerbini AN, Francis MP, Domeier ML, Andriolo A, Costa DP, Robinson PW, Duffy CAJ, Nasby-Lucas N, Holdaway RN and Clapham PJ (2017) Route Fidelity during Marine Megafauna Migration. Front. Mar. Sci. 4:422. doi: 10.3389/fmars.2017.00422 The conservation and protection of marine megafauna require robust knowledge of where and when animals are located. Yet, our ability to predict animal distributions in space and time remains limited due to difficulties associated with studying elusive animals with large home ranges. The widespread deployment of satellite telemetry technology creates unprecedented opportunities to remotely monitor animal movements and to analyse the spatial and temporal trajectories of these movements from a variety of geophysical perspectives. Reproducible patterns in movement trajectories can help elucidate the potential mechanisms by which marine megafauna navigate across vast expanses of open-ocean. Here, we present an empirical analysis of humpback whale (Megaptera novaeangliae), great white shark (Carcharodon carcharias), and northern elephant seal (Mirounga angustirostris) satellite telemetry-derived route fidelity movements in magnetic and gravitational coordinates. Our analyses demonstrate that: (1) humpback whales, great white sharks and northern elephant seals are capable of performing route fidelity movements across millions of square kilometers of open ocean with a spatial accuracy of better than 150 km despite temporal separations as long as 7 years between individual movements; (2) route fidelity movements include significant (p < 0.05) periodicities that are comparable in duration to the lunar cycles and semi-cycles; (3) latitude and bedrock-dependent gravitational cues are stronger predictors of route fidelity movements than spherical magnetic coordinate cues when analyzed with respect to the temporally dependent moon illumination cycle. We further show that both route fidelity and non-route fidelity movement trajectories, for all three species, describe overlapping in-phase or antiphase sinusoids when individual movements are normalized to the gravitational acceleration present at migratory departure sites. Although these empirical results provide an inductive basis for the development of testable hypotheses and future research questions, they cannot be taken as evidence for causal relations between marine megafauna movement decisions and geophysical cues. Experiments on model organisms with known sensitivities to gravity and magnetism, complemented by further empirical observation of free-ranging animals, are required to fully explore how animals use discrete, or potentially integrated, geophysical cues for orientation and navigation purposes.

Keywords: navigation, gravity, moon, humpback whale, great white shark, elephant seal, tracking, g-space

#### INTRODUCTION

During some of the most spectacular yet least understood events in nature, millions of animals migrate between widely separated habitats without getting lost. In fact, animals from at least three different phyla are able to relocate previously inhabited sites, including chordates, arthropods, and molluscs (Switzer, 1993). With respect to marine megafauna, humpback whales (Megaptera novaeangliae), great white sharks (Carcharodon carcharias), and northern elephant seals (Mirounga angustirostris) maintain interannual site fidelity to specific seasonal habitats (Oliver et al., 1998; Wedekin et al., 2010; Anderson et al., 2011) despite undertaking long-distance migrations that span thousands of kilometers of open ocean (Le Boeuf et al., 2000; Nasby-Lucas and Domeier, 2012; Garrigue et al., 2015). Direct observation of the same whales, sharks, and seals in the same areas, year after year, demonstrates that all three species have well-developed navigational abilities that enable high levels of spatiotemporal movement accuracy and precision. Despite our awareness of these remarkable movements, a mechanistic understanding of how marine megafauna navigate remains elusive.

One of the main reasons why we do not yet understand the mechanics of marine megafauna navigation is the fact that we have not yet identified the coordinate space in which navigation occurs. Several fundamental questions remain unanswered for most migratory species, including: What exogenous cues are used for orientation purposes? What reference frame(s) and reference datum(a) are applied during the transduction and neurological processing of these cues? To what extent are endogenous cues integrated with exogenous cues during navigation?

In this study, we explore these knowledge gaps through empirical analyses of humpback whale, great white shark, and northern elephant seal route fidelity movements from both magnetic and gravitational geophysical perspectives. We focus our analyses on magnetic inclination and vertical gravitational acceleration cues (see Nomenclature for a glossary of terms) as experimental studies have suggested both can serve as exogenous sources of orientation information during animal movement (Wiltschko and Wiltschko, 1972; Larkin and Keeton, 1978; Keeton, 1979; Kanevskyi et al., 1985; Light et al., 1993; Lohmann and Lohmann, 1994; DeVries et al., 2004; Putman et al., 2011; Blaser et al., 2013, 2014). By focussing our analyses on route fidelity movements, our research is directly relevant to positional orientation during migration, one of the greatest unknowns in animal navigation science (Gould, 2004). By considering a diverse suite of species and scientific disciplines, our research represents a direct response to growing calls for more integrative research on animal migration and navigation (Bowlin et al., 2010; Hays et al., 2016).

Humpback whales, great white sharks and northern elephant seals are ideally suited to empirical analysis of animal movement due to the fact that telemetry datasets for all three species include remarkable examples of asynchronous migratory route fidelity. Route fidelity is similar to site fidelity in that both refer to the repeated utilization of migratory destinations at distinctly different times. In contrast to site fidelity, route fidelity refers to the repeated utilization of well-defined migration routes by either the same individual during multiple independent migrations (i.e., intra-individual route fidelity) or different individuals migrating separately (i.e., inter-individual route fidelity). In this study, we present and analyse 22 humpback whale, great white shark, and northern elephant seal route fidelity movements, including multiple examples of intra- and inter-individual route fidelity. Our dataset includes humpback whales in the South Atlantic and South Pacific Oceans, great white sharks in the North and South Pacific Oceans and northern elephant seals in the North Pacific Ocean.

#### MATERIALS AND METHODS

#### Satellite Telemetry

Humpback whales were tracked using published methods (Garrigue et al., 2015). In brief, a carbon-fiber pole (Heide-Jørgensen et al., 2003) or a modified pneumatic line-thrower (Heide-Jørgensen et al., 2001) pressurized to ∼10 bar with compressed air was used to implant transdermal location-only SPOT radio-frequency platform transmitting terminal (PTT) satellite tags (Wildlife Computers, Redmond, WA, U.S.A.) into the upper flank of each whale near the base of the dorsal fin. Transmitters were duty cycled to maximize battery life as described in Zerbini et al. (2006) and Hauser et al. (2010). In all references to PTT tag numbers in the current study, the two digits to the right of the decimal point correspond with the abbreviated Julian calendar year in which the satellite-monitored movement was initiated.

Great white sharks were tracked using similar technology deployed by different methods (Domeier et al., 2012; Francis et al., 2012). In brief, SPOT5 PTT tags were affixed using 3– 4 small plastic or stainless steel bolts in the apex of the dorsal fin of each temporarily restrained individual. Unlike marine mammals, great white sharks do not have to surface to breathe. Thus, position data from this species are much more sporadic than they are for whales and seals. Location estimates for white sharks were received only when animals spent enough time at the surface, with the dorsal fin above the waterline, allowing the ARGOS (Argos, CLS Group) satellite array to receive three or more consecutive transmissions from the tag. Surface swimming behavior varied among individual sharks, resulting in varying resolution migration data for each shark.

Northern elephant seals were also tracked using PTT satellite tags (Robinson et al., 2012) which were attached to the head of each seal during onshore residence in Año Nuevo State Park (ca. 37.11◦N; 122.33◦W), California, U.S.A. Elephant seal locations were estimated using the ARGOS satellite system based on transmissions received when the PTT was out of the water.

Route fidelity movements were identified from larger satellite telemetry datasets by visual inspection of PTT track maps using GIS software. Movement tracks, or portions of tracks, that visually overlapped for extended periods (i.e., coincident animal location symbols or track lines for multiple days in geographic coordinates) when viewed at a scale of 1:10,000,000 were identified as route fidelity movements. At this scale, an 8 point map symbol is approximately 10 km wide. The three South Atlantic humpback whale route fidelity movements were identified from inspection of a total of 12 long-distance migration tracks (i.e., 25%). The four South Pacific humpback whale route fidelity movements were identified from inspection of a total of 13 intra-tropical movement tracks (i.e., 31%). The two South Pacific great white sharks demonstrating both intra- and interindividual route fidelity were identified from inspection of longdistance PTT tracks of three different sharks (i.e., 67%). The intra-individual route fidelity movements of two North Pacific great white sharks were identified from inspection of the PTT tracks of 11 individual sharks (i.e., 18%). The four northern elephant seal tracks were identified from inspection of a total of 74 long-distance migration tracks (i.e., 5%).

Raw animal location estimates downloaded from the ARGOS system were processed using a 20 km/h velocity filter and combined to determine single average daily locations for each individual. Transmissions were not received by the ARGOS system from all PTT tags on all days. Average daily locations were only determined for those calendar dates on which velocityfiltered locations were received. Gaps in the tracking datasets were not filled by interpolation due to the fact that we do not know the coordinate space in which navigation was performed (Horton et al., 2014). All animal tracking research reported here was carried out in accordance with animal ethics consents given to the authors by their home institutions and/or relevant government agency.

#### Astronomical and Geophysical Variables

We determined multiple astronomical, magnetic, and gravitational cues present at the whale, shark and seal locations recorded by the PTT animal tracking devices (see Nomenclature for a glossary of terms). Astronomical cues, including moon illumination, were calculated using published astronomical algorithms (Meeus, 1991). Moon illumination is a unitless time-dependent quantity ranging between 0 (i.e., new moon) and 1 (i.e., full moon) across the average 29.53-day-long synodic (i.e., moon illumination) cycle. We used it as a direct proxy for time to facilitate comparative analyses of multiple asynchronous individual telemetry tracks in a single panel.

With respect to magnetic cues, we determined main field plus rock anomaly field magnetic coordinates, including magnetic inclination, field intensity and declination, for all animal locations using the Enhanced Magnetic Model (Maus, 2010). Of the seven different magnetic variables used to define the position of Earth's magnetic field from a geocentric perspective in both spherical (F, I, D) and Cartesian (X, Y, Z, H) coordinate spaces (for definitions see Nomenclature; Horton et al., 2014), magnetic inclination (I) is the most widely associated with animal orientation and navigation (e.g., Wiltschko and Wiltschko, 1972; Light et al., 1993). Extensive experimental research has demonstrated that the movement behaviors of diverse species change in response to changes in magnetic inclination. These results have been variably interpreted as evidence for utilization of magnetic inclination as: (1) a navigational compass that facilitates directional orientation (Wiltschko and Wiltschko, 1972), or, (2) part of a metaphorical bi-coordinate "magnetic map" that facilitates positional orientation (Putman et al., 2011). Regardless of how magnetic inclination is used during animal navigation and orientation, compelling experimental evidence demonstrates that it is part of the system by which many migratory species find and follow specific movement trajectories.

Thus, a key question emerges: Are route fidelity movements compatible with orientation relative to magnetic inclination? We explored this question by plotting magnetic inclination vs. moon illumination for the route fidelity movements reported above. Considering the fact that these highly directional movements were performed at distinctly different times in distinctly different locations, strong correlations and systematic nonrandom patterns in the data would be suggestive of a potential spatiotemporal orientation behavior informed by cues associated with magnetic inclination.

With respect to gravity, we determined local gravitational accelerations associated with both latitude and bed-rock density using the International Gravity Formula (Götze, 2014) and the International Gravimetric Bureau's 2 × 2 arc-min (i.e., ∼3.7 × ∼3.7 km) World Gravity Map (Balmino et al., 2012), respectively. Latitude (gL) and bed-rock (gB) vertical gravitational accelerations are reported in Gals (i.e., cm sec−2 ).

Of the multiple factors that determine gravity (i.e., g) at a given place and time, latitude has the largest effect. At 20◦ north or south latitude, the theoretical gL is 978.637 Gal. In contrast, gL is 981.070 Gal at 50◦ north or south latitude, equating to a 2.433 Gal range in gL over this 30◦ range. Importantly, gL is a trigonometric function of latitude with the most rapidly changing gL values occurring in the middle latitudes and the most gradually changing gL values occurring near the geographic equator and poles. For example, an animal that migrates from the equator to 30◦ south or north latitude would experience a change in gL of 1.292 Gal, roughly half the range in gL experienced by migrating between 20◦ and 50◦ latitude. To facilitate understanding, a 30,000 kg humpback whale requires 900 N (i.e., 202 lbs) less buoyancy force to remain effortlessly afloat in a tropical habitat, where gL is ∼979 Gal, relative to a high-latitude habitat, where gL is ∼982 Gal.

The shape and density of Earth also affects local g. However, the difference in geologically imparted Bouguer gravity anomalies from one location to the next is on the order of 10 to 100 mGal, roughly 1 to 10% as large as the latitudinal effects on g across the geographic range of the tracks we report. Since both latitudinal and geological effects on g are temporally independent, we additively combined these two gravitational variables into a single spatially dependent coordinate, gL+gB, that is determined solely by an animal's location.

The whales, sharks, and seals we studied experienced a wide range of gravitational conditions during their individual migrations. Latitude dependent gravity values (gL) ranged between 978.287 Gal (i.e., cm/sec2 ) and 981.737 Gal. Local bedrock-dependent gravity values (gB) ranged between +0.025 Gal and +0.587 Gal. Using these spatially dependent gravitational variables, we constructed bivariate plots of all three species' movements through gravitational coordinates, or g-space. We define g-space as the bicoordinate system that includes the sum of gL and gB as the spatially dependent variable and moon illumination as the temporally dependent variable.

As an animal moves, its gL value will change as its latitudinal position changes, its gB value will change as the density of the Earth's crust over which it is swimming changes, and the magnitude of the tidal gravity vector will change in unison with lunar phase. In contrast, for an animal that stops its migratory movement and begins residence in a single area, only moon illumination would continue to change, and the gL and gB values it experiences would remain constant. Pure east-west movements would experience changing gB values and moon illuminations while the latitudinal component (gL) remains constant.

#### Data Analysis Methods

The satellite telemetry, astronomical, and geophysical data generated by our research were analyzed using a variety of statistical tools. Piecewise linear regression breakpoint analysis (Muggeo, 2008) performed on latitude/longitude time-series plots was used to estimate the date on which open ocean migration was initiated at the individual scale. Latitude timeseries plots were used on all individuals with the exception of the North Pacific great white sharks, where longitude timeseries were used due to the predominantly meridional nature of these tracks. Spectral analysis using the Lomb-Scargle algorithm (Scargle, 1982; Press et al., 1992; Hammer et al., 2001) was used to determine significant (α = 0.05) periodicities present in latitude time-series plots at the individual scale. Sinusoidal regression (Press et al., 1992; Hammer et al., 2001) was performed on geophysical (i.e., magnetic and gravitational) cue vs. moon illumination plots in order to determine the proportion of data variance explained (i.e., coefficient of determination; R 2 ) by a sine-function model of the data.

## RESULTS

### Route Fidelity Movements

Inspection of the humpback whale, great white shark, and northern elephant seal tracking datasets revealed multiple examples of migratory route fidelity (**Table 1**). Migration tracks exhibiting route fidelity were followed by: humpback whales tagged off the coasts of Brazil and Rarotonga, Cook Islands; great white sharks tagged off southern New Zealand and Guadalupe Island, Mexico; northern elephant seals tagged in Año Nuevo State Park, California, U.S.A. Our satellite tracking results demonstrate that these diverse marine megafauna demonstrated a remarkable ability to find and follow near-identical paths despite the fact that none of these individuals was ever within 100 km of the same geographic location at the same Julian calendar date and time.

### South Atlantic Humpback Whale Route Fidelity Movements

Three humpback whales tagged off Brazil (PTT identification numbers: 87760.08; 87761.08; 87769.08) swam near-identical south then southeast directed asynchronous migration paths away from Abrolhos Bank off the eastern coast of Brazil (ca. −39◦W; −19.8◦ S; **Figure 1A**; Zerbini et al., 2006) and toward higher-latitude feeding grounds in the South Atlantic Ocean (ca. −50 to −60◦ S; **Figure 1A**; Zerbini et al., 2011). All three whales followed approximately the same ∼1,550 km southerly path away from Abrolhos Bank, during the first 12 (PTT 88760.08), 15 (PTT 87769.08), and 16 (PTT 87761.08) days of their temporally distinct migrations (**Figure 1A**). These three whales migrated through a migratory corridor that was at most ∼100 km wide despite swimming across a vast expanse of open-ocean at different times. Average swimming speeds during these southdirected movements ranged between 4.3 and 5.0 km/h (±1.1 to 1.4 km/h, ±SD).

Following these initial movements, two of the whales (PTT 87760.08 and PTT 87769.08) turned to a southeast heading until both tags stopped transmitting another 17 and 24 days later, respectively (N.B. PTT 87761.08 stopped transmitting ca. −38◦W; −34◦ S). During these southeast-directed movements, 88760.08 and 87769.08 swam an additional 1,991 and 2,581 km, respectively, through a <85 km wide open-ocean corridor despite migrating ∼15 days apart (**Figure 1A**). The average swimming speeds during these southeast-directed movements were 4.8 and 4.4 km/h (±0.9 and ±1.5 km/h, ±SD) for whales 87760.08 and 87769.08, respectively. For comparison, we include two additional migration tracks of humpback whales that migrated through the same corridor, but along distinctly different geographic coordinate trajectories, in 2005 (PTT 10946.05, **Figure 1A**) and 2012 (PTT 111871.12, **Figure 1A**).

### South Pacific Humpback Whale Route Fidelity Movements

Four humpback whales tagged off Rarotonga, Cook Islands (PTT identification numbers: 37282.07; 120946.14; 120947.14; 121195.14) swam similar asynchronous northwest directed migration paths away from Rarotonga, in the southern Cook Islands (ca. −159◦W; −21◦ S; **Figure 1B**), to Tutuila, American Samoa (ca. −171◦W; −14◦ S; **Figure 1B**). All four of these route fidelity movements began in mid-September, albeit 7 years apart (one in 2007, three in 2014). All four whales followed similar ∼1,300 km long northwest directed paths to Tutuila, American Samoa, during the first 12 to 14 days of their intratropical movements despite migrating at distinctly different times TABLE 1 | Platform transmitting terminal (PTT) tag deployment information.

\*Average Velocity was determined during the period of continuous and directed open-ocean movement. Individuals not known to perform route fidelity movements (i.e., non-route fidelity tracks) are indicated by a double asterix (\*\*).

(**Figure 1B**). The Rarotonga-Tutuila migratory corridor used by these whales was at most 150 km wide in the area ∼300– 400 km northwest of Rarotonga and 50–100 km wide elsewhere (**Figure 1B**). The open-ocean swimming speeds for the four intra-tropical movements ranged between 3.9 and 5.0 km/h.

#### South Pacific Great White Shark Route Fidelity Movements

Two great white sharks departed the waters off southwest New Zealand in early to mid-July, 2014 and 2015 (PTT 55614), and early to mid-September, 2014 (PTT 55615). Although the two sharks followed distinctly different paths between New Zealand and the southeast Great Barrier Reef, Australia, shark 55614 followed a near-identical route between New Zealand and Australia in both 2014 and 2015 (**Figure 2A**). Additionally, sharks 55614 and 55615 both followed a similar return migration path from coastal waters off Byron Bay, Australia, to coastal waters off southwest New Zealand (**Figure 2A**) in late November, 2014 (PTT 55615), and early December, 2014 (PTT 55614).

The satellite tracking results acquired for these two sharks include examples of both intra-individual (i.e., PTT 55614, New Zealand to Australia in 2014 and 2015) and inter-individual (i.e., PTT 55614 and PTT 55615, Australia to New Zealand in 2014) route fidelity. Like both the South Atlantic and South Pacific humpback whale datasets, satellite tracking demonstrates that these two sharks can find and follow <150 km wide migratory corridors during temporally separated open-ocean movements that are well in excess of 1,000 km distance during multiple weeks of continuous swimming. For comparison, we include a third great white shark (PTT 55612) that performed a round-trip migration between Stewart Island, New Zealand, and the Loyalty Islands, New Caledonia (ca. 167◦E; −21◦ S; **Figure 2A**), during 2013.

#### North Pacific Great White Shark Route Fidelity Movements

Two great white sharks tagged off Guadalupe Island, Baja, Mexico (PTT identification numbers: 19787; 20720) swam similar to near-identical asynchronous west-southwest directed migration paths away from Guadalupe (ca. −118◦E; 29◦N) to the central northeast Pacific shared offshore foraging area, or SOFA (ca. −134◦W; 23◦N; **Figures 2B,C**; Domeier et al., 2012). These movements were repeated each year from 2009 to 2012, by shark 19787 (**Figure 2B**), and 2009 to 2011, by shark

FIGURE 2 | Route fidelity and non-route fidelity movements of satellite-monitored great white sharks, including: (A) great white sharks in the South Pacific Ocean, (B,C) great white sharks in the North Pacific Ocean. Symbol sizes and color hue correspond with average daily track velocity (km/h) as shown in the legends. Unique platform transmitting terminal (PTT) tag numbers to the left of the PTT number decimal point correspond with individual sharks and are represented by symbol shape; digits to the right of PTT number decimal point correspond with the abbreviated year in which the whale was tagged and are represented by symbol color family (i.e., red, yellow, brown, etc.). Addition (PTT 55612.13) and black/gray diamond (PTT 55615.14) symbols in (A) correspond with long-distance non-route fidelity migrations. 20720 (**Figure 2C**). The tracking data indicate that shark 19787 consistently departed Guadalupe several weeks earlier than shark 20720. Specifically, shark 19787 departed Guadalupe between:

Whereas, shark 20720 departed Guadalupe between:

Piecewise linear regression breakpoint analysis of longitude vs. time profiles for these seven separate migrations suggests that they likely began on, or about, January 15, 2009, January 24, 2010, February 4, 2011 and January 9, 2012 for shark 19787, and March 22, 2009, April 1, 2010 and March 18, 2011 for shark 20720.

Despite the relatively low temporal resolution of the Mexico great white shark dataset, these data reinforce the observation that great white sharks are capable of migrating through welldefined <150 km wide migratory corridors with a high degree of inter-annual route fidelity (**Figures 2B,C**). Four years in a row, shark 19787 followed the same west/southwest-directed path from Guadalupe to the SOFA (**Figure 2B**). Shark 20720 followed a similar but consistently more southerly path than 19787 every year between 2009 and 2011 (**Figure 2C**).

#### North Pacific Northern Elephant Seal Route Fidelity Movements

Four northern elephant seals tagged at Año Nuevo State Park, California, U.S.A. (PTT identification numbers: 39454.04, 13365.04, 39455.05; 62036.11) swam similar to near-identical asynchronous northwest directed migrations between the central California coast and the central North Pacific Ocean/Gulf of Alaska (**Figure 3**). In contrast to previous research documenting intra-individual route fidelity in male northern elephant seal tracks (Le Boeuf et al., 2000), the seal tracks we analyzed demonstrate inter-individual route fidelity achieved by three females and one male. These seals departed the California coast at Año Nuevo State Park between June 9, 2004 and February 21, 2011 (**Table 1**). Despite migrating across a >6 year period, these four seals followed the same ∼150 km wide migration route during their individual 30–38 day, and 2,500–4,000 km long, northwest-directed migrations (**Figure 3**).

Like the humpback whale and great white shark datasets we analyzed, satellite tracking demonstrates that these four seals were able to find and follow a well-defined migratory corridor during asynchronous long-distance open-ocean movements spanning several weeks of continuous swimming. Although we focus our analysis on the initial route fidelity stage of these elephant seal migrations, these same seals dispersed across ∼3 million km2 of the North Pacific Ocean during subsequent non-route fidelity stages of their individual migrations. For comparison, we include two additional migration tracks of northern elephant seals that followed distinctly different

FIGURE 3 | Route fidelity and non-route fidelity movements of satellite-monitored northern elephant seals. Symbol sizes and color hue correspond with average daily track velocity (km/h) as shown in the legends. Unique platform transmitting terminal (PTT) tag numbers to the left of the PTT number decimal point correspond with individual seals and are represented by symbol shape; digits to the right of PTT number decimal point correspond with the abbreviated year in which the whale was tagged and are represented by symbol color family (i.e., green, yellow, black, etc.). Yellow multiplication (PTT 133774.14) and addition (PTT 120078.14) symbols correspond with long-distance non-route fidelity migrations.

geographic coordinate trajectories in 2014 (PTT 120078.14 and 133774.14, **Figure 3**).

#### Timing of Route Fidelity Movements

Satellite-monitored PTT tracking of humpback whales, great white sharks, and northern elephant seals in the Atlantic and Pacific reveals that all three species have navigational systems capable of reproducing near-identical movements across vast expanses of featureless ocean (**Figures 1**–**3**). At a minimum, these remarkable examples of navigational reproducibility require both precise and accurate spatial orientation. The fact that all three of the sharks that were tracked for more than 1 year (i.e., PTT identification numbers: 55614, 19787, 20720) also showed a migratory fidelity to specific times of the year further suggests there is a temporal component to open-ocean navigation and migration.

Indeed, one of the distinctive features of the route fidelity movements we report is the fact that they occurred at different times. The asynchrony of these movements is important: little navigational information can be gained from analysis of two or more marine migrants that are swimming together as the coordinate space trajectories followed by these animals would be indistinguishable from any geophysical or environmental perspective. Thus, it is the combined effects of the movement asychrony and the temporal dependence of orientation cues available from the environment that make route fidelity movements potentially novel indicators of navigational behavior. We focused our temporal analysis of the route fidelity movements described above on whether or not a systematic temporal pacing was present.

Spectral analysis using the Lomb-Scargle periodogram algorithm (Scargle, 1982) performed on individual route fidelity latitude-time datasets demonstrates that significant (p < 0.05) periodicities are present in all of the route fidelity movements with 15 or more average daily PTT locations. The fact that significant periodicities were not detected in the route fidelity tracks with <15 locations in total does not necessarily mean that these movements lacked temporal pacing. Rather, the relatively small number of PTT locations in these spatially and temporally shorter tracks may simply preclude periodicity detection using the Lomb-Scargle algorithm.

Significant periodicities were detected in three South Atlantic humpback whale, two South Pacific white shark and four North Pacific elephant seal route fidelity movements (**Figure 4**). These nine tracks represent the nine longest distance and duration annual movements in our dataset. None of the shorter duration intra-tropical Rarotonga humpback whale tracks, nor any of the Guadalupe Island white shark tracks, included a detectable significant periodicity. Spectral analysis of the nine longest duration tracks detected significant periodicities with average peak powers of 15 days (SD = 3 days; n = 3; **Figure 4**) and 27 days (SD = 2 days; n = 7; **Figure 4**). These significant periodicities are not unlike the period (range = 29.3–29.8 days) and semiperiod (range = 14.6–14.9 days) of the lunar illumination cycle (i.e., synodic cycle). Periodicities in these ranges are not entirely unexpected given the growing number of studies demonstrating that lunar illumination is strongly correlated with a variety of organismal behaviors including animal movement (e.g., Larkin and Keeton, 1978; Grau et al., 1981; Baird et al., 2003; Tsukamoto et al., 2003; DeVries et al., 2004; Fraser, 2006; Pinet et al., 2011; Erisman et al., 2012; Schmidt et al., 2012; Sudo et al., 2014).

The tracks we studied further demonstrate that synodic periods are also present in the time gaps separating individual movements (Figure S1). For example, the three 2008 South Atlantic whales tracks that individually include significant semisynodic periodicities (**Figure 4**), also demonstrate an intertrack semi-synodic separation that is maintained throughout the southward migrations of these whales (Figure S1A). This ∼15 day separation persists despite a reduction in movement velocity between ∼30◦ S and ∼40◦ S latitude in all three tracks (**Figure 1A**). Two other South Atlantic humpbacks that followed distinctly different geographic coordinate paths (i.e., non-route fidelity movements; **Figure 1A**) across the same southeastdirected migratory corridor in 2005 (PTT 10946.05) and 2012 (PTT 111871.12) performed their migrations 36 and 52 lunar synodic cycles before and after the 2008 whales, respectively (Figure S1A).

Despite the fact that no significant periodicities were detected within any of the individual Rarotonga humpback whale route fidelity movements, these four intra-tropical movements occurred approximately 0.5, 86.0, or 86.5, synodic cycles apart (Figure S1A). Similar time gaps, ranging between 0.5 and 70.5 synodic cycles, separate both the route fidelity and non-route fidelity movements of the white sharks and elephant seals we tracked (Figures S1B–D) suggesting there is a strong temporal component to marine megafauna movement.

Time-series analyses of the humpback whale, great white shark, and northern elephant seal movements we report demonstrate that these movements include significant and systematic temporal periodicities that are similar in duration to the synodic and semi-synodic cycles. However, the synodic cycle is just one of the many quasi-monthly lunar cycles (i.e., sidereal, anomalistic, tropical, etc.) modulated by the relative position of the moon. The empirical results we report demonstrate that humpback whales, great white sharks and northern elephant seals are capable of performing long distance movements that include both spatial and temporal fidelity to well-defined migratory trajectories.

#### Movements in Magnetic Coordinates

Our analyses demonstrate that magnetic inclination is a strong predictor of route fidelity movements (**Table 2**), consistent with the hypothesis that magnetic inclination informs navigation. At the population level, route fidelity movements appear as overlapping in-phase or antiphase sinusoids when magnetic inclination values are plotted against moon illumination (**Figure 5**). These unexpected systematic and symmetrical correlations require temporal pacing: even individuals that followed the same track at the same speed would not be expected to show this pattern unless their movements were initiated at similar, or antithetical, times in the synodic cycle. The systematic nature of the movement trajectories shown in **Figure 5** is reinforced by the fact that the overlapping in-phase/antiphase TABLE 2 | Significant (p < 0.05) sinusoidal regression correlation coefficients for route fidelity movements in geophysical coordinates relative to moon illumination.

(n) Corresponds with the number of average daily locations included in the open-ocean route fidelity movement analyzed. Gray shading highlights the largest correlation coefficient determined for each route fidelity movement. Bold-face and italicized font indicates the second largest correlation coefficient determined for each route fidelity movement.

pattern persists despite temporal separations between individual tracks as long as 7 years. Although these magnetic inclination vs. moon illumination trajectories are highly unexpected, they are consistent with the synodic and semi-synodic periodicities revealed by spectral analysis.

Perhaps even more surprising, however, is the finding that the route fidelity movements of individual species appear to overlap despite extreme geographic separations between the different populations studied. For example, both the South Atlantic and South Pacific humpback whale route fidelity tracks, and both the South Pacific and North Pacific great white shark route fidelity tracks, describe near continuous magnetic inclination vs. moon illumination trajectories (**Figures 5A–D**). Six of the seven humpback whale movements we report departed coastal breeding grounds, in areas with magnetic inclination values ca. −38 to −39◦ , within a few days of full or new moon (**Figure 5A**). Similarly, most of the South Pacific and North Pacific great white sharks arrived at, or departed from, coastal locations in Australia or Mexico, in areas with magnetic inclination values ca. ±53 to ±55◦ , within a few days of full or new moon.

Sinusoidal regression with respect to the moon illumination cycle reveals that magnetic inclination is the strongest magnetic

FIGURE 5 | Magnetic inclination vs. moon illumination plots. (A,B) Humpback whales in the South Atlantic and South Pacific Oceans, (C,D) great white sharks in the South Pacific and North Pacific Oceans, (E,F) northern elephant seals in the North Pacific Ocean. Symbols as in Figures 1–3. In (B,D,F), average daily location symbols have been removed, and the tracks of PTT numbers 87769.08, 120946.14, 55614.15, 55615.15, 19787.10, 19787.11, and 39455.05 are plotted against the upper reverse moon illumination axis. These seven tracks are thus mirrored across the 0.5 (i.e., 50%) moon illumination value. (B,D,F) are included in order to show the symmetrical in-phase or antiphase distribution of the data with respect to the synodic cycle. These 7 mirror-image tracks are plotted as lighter colored lines. All other tracks (n = 15) are plotted against the lower moon illumination axis and are shown as darker colored lines. Magnetic inclination values in (C,D) are plotted as absolute values due to the fact that the two white shark populations reside in opposite magnetic hemispheres. Magnetic inclination axis values span the same range in all panels to facilitate comparisons between the three species.

spherical coordinate predictor of individual and populationlevel route fidelity movements (**Table 2**). However, these same movements also exhibit significant sinusoidal correlations with respect to both magnetic field intensity and magnetic declination in the majority of tracks analyzed (**Table 2**). Thus, it is not possible to conclusively rule out potential roles for any of the spherical magnetic coordinates in whale, shark or seal navigation based on these analyses alone.

However, magnetic inclination explains a higher proportion of the sinusoidal regression model variance at the individual track level (median = 96%; average = 93%; range = 83–99%; SD = 4.2%; n = 19) than magnetic field intensity (median = 92%; average = 84%; range = 42–99%; SD = 16.4%; n = 19) or magnetic declination (median = 79%; average = 80%; range = 52–99%; SD = 14.2%; n = 17). In contrast, when all of the route fidelity tracks from all populations of all species are combined, magnetic field intensity is the strongest magnetic coordinate predictor, explaining 33% of the model variance (**Table 2**), followed by magnetic inclination (28%) and magnetic declination (21%).

#### Movements in Gravitational Coordinates

G-space plots of marine megafauna route-fidelity migrations reveal that these long-distance open-ocean movements describe highly symmetrical trajectories (**Figure 6**) not unlike the magnetic inclination trajectories reported above. Several aspects of these results are notable.

First, all three South Atlantic humpback whales departed the southeast corner of Abrolhos Bank (ca. gL+gB = 978.7 to 978.8 Gal) at antithetical positions of the synodic cycle, 2–3 days prior to full or new moon (**Figure 6A**). These whales maintain this symmetry, about a quarter moon (i.e., 0.5 or 50% moon illumination) mirror plane, throughout their migrations toward higher-latitude feeding grounds (**Figure 6B**). For example, 88760.08 and 87769.08 both pass the gravitational half-way point in their migrations (ca. gL+gB = 980.42 Gal; −33.1◦W; −37.4◦ S) when the moon is either full or new (**Figure 6A**). This spatiotemporal symmetry is a direct consequence of the fact that all three whales systematically reduced their swimming speeds by ∼1.5–2.0 km/h as they approached the gravitational half-way point of their southward migrations (**Figure 6A**).

Second, similar symmetrical g-space trajectories are present in the four South Pacific humpback whale route fidelity tracks (**Figure 6A**). Despite migrating at distinctly different calendar dates across a 7 year period, all four whales departed Rarotonga, Cook Islands (ca. −159.8◦W; −21.2◦ S) at antithetical positions of the moon illumination cycle within hours of full or new moon. All four whales arrived off the southeast coast of Tutuila, American Samoa (ca. −170.8◦W; −14.4◦ S), near the subsequent full or new moon (**Figure 6A**). As recognized in the South Atlantic dataset, these movements are symmetrical across the gravitational half-way point (ca. 979.01 Gal; −167◦W; −16.3◦ S in geographic coordinates; **Figure 6B**). This symmetry exists despite each whale: (1) following slightly different routes during the first week of their movements to Tutuila (**Figure 1B**); and (2) swimming at different speeds during different stages of their movements (**Figure 1B**).

Third, although the great white shark tracks we analyzed predominantly demonstrate intra-individual route fidelity, these tracks have many of the same characteristics as the humpback whale movements when plotted in gravitational coordinates (**Figure 6C**). For example, during its near-identical 2014 and 2015 northward migrations to Australia, white shark 55614 departed southwest New Zealand within 24 h of the July, 2014, full moon and the August, 2015, new moon. During these geographically coincident northward migrations, 55614 passed the gL+gB midpoint value of 980.42 Gal within 24 h of both the July, 2014, new moon and August, 2015, full moon (**Figure 6D**), at a geographic position located ∼350 km southeast of Sydney, Australia. During its intervening 2014 southward migration, 55614 passed the same 980.42 Gal gravitational midpoint within 24 h of the December, 2014, new moon (**Figure 6C**) at a geographic position located ∼500 km east of where it passed the same g-space midpoint during its northward passage (**Figure 2A**).

Despite following distinctly different northward and southward geographic coordinate routes during its seasonal migrations between New Zealand and Australia (**Figure 2A**), the g-space trajectories of white shark 55614's movements are all but indistinguishable from the g-space trajectories followed by three South Atlantic humpback whales (**Figure 6A**) and at least one other south Pacific white shark (PTT 55615; **Figures 6C,D**). Our analyses further show that South Pacific white sharks 55614 and 55615 followed diametrically opposed g-space trajectories during their 2014 return migrations from Australia despite performing these migrations approximately 15 days apart (**Figures 6C,D**, Figure S1C).

Fourth, the seven examples of intra-individual route fidelity performed by the two North Pacific great white sharks tagged off Guadalupe Island provide further empirical evidence of the observed pattern of symmetrically distributed route fidelity movements in gravitational coordinates (**Figure 6C**). For example, every year between 2009 and 2012, white shark 19787 swam a near-identical route from Guadalupe Island to the SOFA (**Figure 2B**), passing the gL+gB midpoint (ca. 979.47 Gal) of its migration route near the: (1) January, 2009, new moon; (2) January, 2010, full moon; (3) February, 2011, full moon; (4) January, 2012, new moon (**Figures 6C,D**). Similarly, great white shark 20720 passed the same gL+gB midpoint (ca. 979.47 Gal) within 48 h of both the March new moon and the April new moon during its 2009 and 2011 migrations to the SOFA, respectively (**Figures 6C,D**).

Fifth, g-Space plots of the near-identical routes followed by four North Pacific northern elephant seals further reinforce the movement patterns described above. All four seals followed gspace trajectories that were highly symmetrical across the gL+gB gravitational mid-point between Año Nuevo and the Aleutian Trench (ca. 981.1 Gal; **Figure 6E**) despite departing the coast at Año Nuevo State Park at distinctly different positions in the moon illumination cycle. These same g-space trajectories are also highly symmetrical across a mirror plane projected through the 50% moon illumination position in the synodic cycle (**Figure 6F**), not unlike the pattern present in both the humpback whale and white shark trajectories. These mirror-image g-space trajectories

FIGURE 6 | Latitude and bedrock dependent gravity vs. moon illumination "g-space" plots. (A,B) Humpback whales in the South Atlantic and South Pacific Oceans, (C,D) great white sharks in the South Pacific and North Pacific Oceans, (E,F) northern elephant seals in the North Pacific Ocean. Symbols as in Figures 1-3. In (B,D,F), average daily location symbols have been removed, and in (B,F) the same tracks are plotted against the upper reverse moon illumination axis as in Figures 5B,F. In (D), PTT numbers 19787.09, 19787.12, 20720.09, 20720.10, and 20720.11 are plotted against the upper reverse moon illumination axis and PTT numbers 19787.10 and 19787.11 are plotted against the lower moon illumination axis. All mirror-image tracks plotted against the upper moon illumination axis are displayed as lighter color lines. gL+gB corresponds with the sum of the latitude (gL ) and bedrock (gB) dependent vertical gravitational accelerations at each individual location (see Nomenclature).

are present due to the fact that all four seals passed the latitude and bed-rock dependent 981.1 Gal position within 6-h of a quarter-moon (**Figures 6E,F**) despite these movements being separated by as much as 7 years in time.

As part of our analysis of the recurrent pattern of symmetrical g-space trajectories observed in all three species, we tested the gravitational datasets for significant correlations using sinusoidal regression. Like we found for the spherical magnetic coordinate movement trajectories, sinusoidal regression of the gravitational trajectories demonstrates that latitude and bedrock dependent gravity is a strong predictor of the individual route fidelity movements (median = 98%; average = 97%; range = 83–99%; SD = 4.1%; n = 20; **Table 2**). The highly symmetrical in-phase or antiphase g-space trajectories we report reinforce the strong inter-annual pacing of route fidelity movements with respect to the moon illumination cycle (**Figure 6**).

Our results demonstrate that the spatially dependent gravitational cue, gL+gB, is a stronger predictor than spherical magnetic coordinates at the individual animal, population, species and inter-species levels (**Table 2**). This unexpected finding is supported by the facts that: (1) gravity is the strongest predictor of 15 of the 22 route fidelity movements we analyzed, whereas magnetic inclination is only the strongest predictor for 4 of the individual movements (**Table 2**); (2) gravity is the strongest route fidelity movement predictor for 3 of the 5 marine megafauna populations we studied (**Table 2**); (3) gravity is the strongest predictor of the route fidelity movements of all three species (**Table 2**); (4) gravity cues explain 69% of the sinusoidal regression model variance when all of the route fidelity data are concatenated, whereas magnetic inclination only explains 28% of the variance (**Table 2**; **Figure 7**).

### DISCUSSION

#### Geophysical Navigation

None of the significant correlations we report demonstrate causality between geophysical orientation cues available from the environment and navigational decisions. Such causal relations can only be established by experimental testing under controlled conditions. However, the strong and systematically patterned correlations we report can be used as an empirical data-based platform from which hypotheses can be proposed.

For example, magnetic coordinate projections of an animal's repeated utilization of a well-defined geographic coordinate migration route will describe distinctly different magnetic coordinate trajectories due to changes in Earth's main magnetic field through time. Yet, the opposite scenario is also true. The geographic coordinate paths followed between migratory destinations might instead systematically shift through time in response to magnetic secular variation in situations where magnetic cues are the primary source of orientation information. Thus, the extent to which an individual migrant uses magnetic cues for navigational purposes can be further explored through multi-annual tracking of long-lived individuals. For large and elusive species that are difficult to study under controlled conditions, long-term repeat tracking represents an important opportunity to better understand their navigation at sea.

The widespread deployment of PTT tags on animal migrants facilitates longitudinal studies, and long-term tracking studies are viable for some species (e.g., Berthold et al., 2004; Alerstam et al., 2006). Based on the empirical results presented in the current study, we hypothesize that the magnetic coordinate trajectories followed by individual migrants over multiple migratory cycles will describe overlapping and symmetrically distributed paths when plotted against the moon illumination cycle despite changes in the raw geocentric magnetic coordinate values caused by secular variation of the magnetic field. If true, this hypothesis would provide evidence in support of temporally modulated navigation with respect to Earth's magnetic field. Given the relatively small changes in magnetic field conditions from 1 year to the next, tracking studies that span a decade or more will produce the strongest results.

With respect to gravity, the possibility that animal navigation is informed by cues derived from Earth's spatially and temporally dynamic gravitational field was suggested at least 40 years ago (Larkin and Keeton, 1978). However, the possibility that animal orientation is informed by gravitational cues remains largely untested. The empirical results we present suggest further experiments, like those performed by Blaser et al. (2013, 2014), Fisahn et al. (2015), and Cresci et al. (2017) will improve our understanding of the role, if any, gravity plays in animal orientation. We particularly encourage experimental tests on model organisms (e.g., zebrafish, Danio rerio; honey bees, Apis mellifera; eels, Anguilla spp.) that are sensitive to both magnetic (Kirschvink, 1981; Tesch et al., 1992; Osipova et al., 2016) and gravitational cues (Korall and Martin, 1987; Moorman and Shorr, 2008; Cresci et al., 2017). Integrated analysis of telemetry datasets will further help identify the ways in which geophysical cues are used for navigational purposes.

Long-distance inter-hemispheric migrants are also attractive targets for future studies. The presence or absence of symmetrical magnetic or gravitational coordinate movement trajectories across either of the spatially distinct magnetic or geographic equators has the potential to provide significant insight into how navigation is performed. The variable separation between the magnetic and gravitational equators at different longitudes presents additional opportunities to determine if emergent patterns in magnetic or gravitational trajectories are reproducible between tracks of related species that separately migrate across the Pacific and Atlantic Oceans.

Tests for non-random and reproducibly patterned movement trajectories, such as the symmetrical geometry of the g-space trajectories we report, can also be more deeply explored for a variety of geophysical and astronomical orientation cues in migratory domains with distinctly different cue distributions. For example, animals that migrate through middle latitude positions from higher or lower latitude habitats are particularly attractive due to the eccentric and anomalous geometries of magnetic, gravitational and astronomical cues through both space and time. Specific targets for future research in this area might include comparisons between populations that migrate across the South Atlantic magnetic anomaly vs. other ocean basins (Figures S3–S8).

Our analyses suggest that gravity might play an important role in long-distance animal navigation, either in concert with magnetic cues or in isolation. Confirming and elucidating this role requires integrated experimental testing and animal tracking. We hypothesize that spatially dependent gravity cues, when plotted against temporally dependent gravity cues, such as tidal gravity, will describe reproducible movement trajectories at the individual level. We further hypothesize that similarly reproducible movement trajectories will be less pronounced, if not absent, when the same individual's movements are analyzed with respect to other geophysical and astronomical cue/coordinate systems.

#### Route Fidelity vs. Non-route Fidelity

Although this study specifically addresses route fidelity movements, it is important to consider whether or not non-route fidelity movements describe geophysical coordinate trajectories that are similar to the route fidelity movement trajectories we report.

Perhaps one of the most unexpected results animal tracking studies have revealed is the extreme diversity in movement trajectories followed by animal migrants when plotted in geographic coordinates (e.g., Block et al., 2011). Few individuals follow the same path. Our research reinforces this observation: of the animals we tracked, we found only 13% (i.e., 15 out of 113 individuals) achieve intra- or inter-individual route fidelity.

However, our dataset may underestimate the prevalence of route fidelity due to the fact that the vast majority of the individuals we tagged were tracked for only one migratory cycle or less. It is also possible that our use of two-dimensional geographic coordinate space projections, when classifying movements as either route fidelity or non-route fidelity, is a flawed approach. Perhaps geographically distinct individual movements become more alike when they are viewed from a different perspective? In an effort to explore this possibility, we compared route-fidelity to non-route fidelity movement projections in geophysical coordinates.

Non-route fidelity humpback whale geophysical coordinate trajectories (Figures S2A,B) have a similar overall sinusoidal shape and pattern as the route fidelity tracks from the same population (**Figures 6A,B**, **7A,B**). These similarities include multiple segments of both magnetic and gravitational coordinate non-route fidelity movement trajectories that overlap the route fidelity movements for extended periods (Figures S2A,B). However, there are also periods, spanning several days of animal movement, during which the non-route fidelity geophysical coordinate trajectories do not overlap the route fidelity trajectories (Figures S2A,B).

The non-route fidelity great white shark trajectories show a similarly ambiguous pattern. These tracks describe largely coincident magnetic and gravitational trajectories during the first-half of the northward movements away from New Zealand for both sharks, yet, distinctly different trajectories in the secondhalf of the northward movements (Figures S2C,D). Similar variations are also apparent during the southward movement of 55612.13 (Figures S2C,D).

Similarly equivocal patterns are also present in the nonroute fidelity northern elephant seal trajectories. In this case, neither of the non-route fidelity tracks overlap the magnetic coordinate route fidelity trajectories (Figure S2E). However, in marked contrast to the magnetic trajectories, both of the nonroute fidelity tracks extensively overlap with the gravitational route fidelity trajectories (Figure S2F).

The varying degrees of similarity between non-route fidelity geophysical coordinate trajectories and route fidelity trajectories suggests that the navigational system(s) being utilized are calibrated at the individual level.

To explore this possibility, we normalized both magnetic inclination and gL+gB values to the values present in core areas inhabited by each individual immediately prior to the onset of long-distance movement. For the northern elephant seals, magnetic inclination and gL+gB values were normalized to the values present at Año Nuevo State Park, California, U.S.A. For the South Pacific great white sharks and the South Atlantic

bivariate spatiotemporal coordinate spaces. Lines in (A) colored as in Figure 5. Lines in (B) colored as in Figure 6.

FIGURE 8 | Site-normalized South Atlantic humpback whale magnetic inclination and gravitational trajectories. In this figure geocentric magnetic inclination and latitude and bedrock dependent gL+gB values have been normalized to the values present at the locations inhabited immediately prior to the onset of long-distance open-ocean migration at the individual level. Site-normalized magnetic inclination (A) and gravitational acceleration (B) values are displayed as percentages and are plotted against the moon illumination cycle as in Figures 5–7. Route fidelity movements are shown as light blue lines and non-route fidelity movements are symbolized as indicated in the figure legends. Non-route fidelity track symbol colors correspond with day of the year and exceed 365 days when the tracked movement spanned Julian calendar years.

FIGURE 9 | Site-normalized South Pacific great white shark magnetic inclination and gravitational trajectories. In this figure geocentric magnetic inclination and latitude and bedrock dependent gL+gB values have been normalized to the values present at the locations inhabited immediately prior to the onset of long-distance open-ocean migration at the individual level. Site-normalized magnetic inclination (A) and gravitational acceleration (B) values are displayed as percentages and are plotted against the moon illumination cycle as in Figures 5–7. Route fidelity movements are shown as dark red lines and non-route fidelity movements are symbolized as indicated in the figure legends. Non-route fidelity track symbol colors correspond with day of the year and exceed 365 days when the tracked movement spanned Julian calendar years.

humpback whales, magnetic inclination and gL+gB values were normalized to the values present at each individual's last known continental shelf location prior to the onset of open-ocean migration. The South Pacific humpback whale and North Pacific great white shark tracks were not included in this analysis due to uncertainties regarding departure site locations for these individuals.

When compared to site-normalized route fidelity migrations, the site-normalized non-route fidelity geophysical coordinate movement trajectories exhibit several distinctive features (**Figures 8**–**10**). First, site-normalized non-route fidelity magnetic inclination trajectories generally do not overlap with route fidelity magnetic inclination trajectories. Second, site-normalized magnetic inclination trajectories exhibit limited symmetry across the 50% moon illumination mirror plane (**Figures 8A**, **9A**, **10A**). Third, both site-normalized route fidelity and site-normalized non-route fidelity gravitational coordinate trajectories predominantly plot as either overlapping in-phase or anitphase sinusoids for all three species (**Figures 8B**, **9B**, **10B**). Fourth, in comparison to the magnetic inclination trajectories, site-normalized gravitational trajectories exhibit a more pronounced symmetry across both the gL+gB midpoint values in the migratory domain utilized by each species and the 50% moon illumination mirror plane (**Figures 8B**, **9B**, **10B**).

The systematic nature of the highly symmetrical site-normalized trajectories shown in **Figures 8**–**10** suggest there is a temporally modulated triggering and/or pacing to long-distance migratory movement behavior. Considering the high specific gravity and extreme crystallographic symmetry of magnetic biominerals, we hypothesize that exogenous magnetic and gravitational cues are integrated components of a spatiotemporal orientation system that is calibrated at the individual level.

#### Future Directions

Knowledge gaps in our understanding of how animals navigate limit our ability to assess and anticipate the sensitivity of migrating animals to perturbations resulting from environmental change, anthropogenic activities, and predator-prey distribution (Hays et al., 2016). Empirical approaches, such as animal telemetry, can be used to develop and assess data-based models of animal movement, with the strongest models creating opportunities to inform conservation and management decision making with respect to both space and time. At a fundamental level, telemetry studies remain uniquely powerful ways to inform marine conservation through databased demonstrations of when and where marine megafauna are located across vast expanses of open-ocean (Shillinger et al., 2008; Gredzens et al., 2014; Maxwell et al., 2016; Dawson et al., 2017). Technological advances in animal telemetry, combined with advances in remote sensing (e.g., Figures S3–S8), have created unprecedented opportunities

FIGURE 10 | Site-normalized North Pacific northern elephant seal magnetic inclination and gravitational trajectories. In this figure geocentric magnetic inclination and latitude and bedrock dependent gL+gB values have been normalized to the values present at the locations inhabited immediately prior to the onset of long-distance open-ocean migration at the individual level. Site-normalized magnetic inclination (A) and gravitational acceleration (B) values are displayed as percentages and are plotted against the moon illumination cycle as in Figures 5–7. Route fidelity movements are shown as dark green lines and non-route fidelity movements are symbolized as indicated in the figure legends. Non-route fidelity track symbol colors correspond with day of the year and exceed 365 days when the tracked movement spanned Julian calendar years.

to retrospectively extract and analyse the geophysical and oceanographic conditions experienced by individual migrants during long-distance migration. Integrated analysis of the data produced by satellite telemetry and remote sensing tools strengthens our understanding of movement behavior at the individual scale.

The data and analyses we report provide a platform for future research that specifically targets both the spatial and temporal aspects of long-distance animal migration. The recently proposed "chord and clock" model of animal navigation (Horton et al., 2014) explicitly includes exogenous temporal cues, consistent with the highly correlated movements we report here. The integrated spatiotemporal perspective that defines the "chord and clock" model provides a novel parallel to the widely accepted, yet purely spatial, positional and directional orientation frameworks proposed by Griffin (i.e., Type III "true" navigation; Griffin, 1952) and Kramer (i.e., "map and compass" navigation; Kramer, 1961).

Future research on the biogeophysics of animal navigation will facilitate the development of a mechanistic understanding of how animals find their way. Given the speed at which oceanic environments are currently changing, technologicallydriven data-based tests of the observation that geophysical cues are strong predictors of the open-ocean movements of diverse marine megafauna are urgently required.

#### ETHICS STATEMENT

All animal tracking research reported here was carried out in accordance with animal ethics consents given to the authors by their home institutions and/or relevant government agency. Northern elephant seal tracking was performed in accordance with approvals granted by the University of California at Santa Cruz Institutional Animal Care and Use Committee and under National Marine Fisheries Service permits #786-1463 and #87-143. Humpback whale tracking was performed in accordance with approvals granted by the Brazilian

#### REFERENCES

Environmental Agency (IBAMA), permit #009/02/CMA/IBAMA and process #02001.000085/02-27 and by the Cook Islands Government, Office of the Prime Minister, to Cook Islands Whale Research. Great white shark tracking was performed in accordance with approvals granted by Secretaria de Medio Ambiente y Recursos Naturales (SEMARNAT), Comision Nacional de Areas Naturales Protegidas (CONANP), and the New Zealand Department of Conservation.

#### AUTHOR CONTRIBUTIONS

TH and RH conceived of the study. TH designed the current study, performed the magnetic, gravitational and astronomical analyses, prepared all of the figures and wrote the manuscript. NH, AZ, MF, MD, AA, DC, PR, CD, NN-L, and PC provided the humpback whale, great white shark and northern elephant seal PTT tracking data used in the current study. All authors contributed to the revision of the initial manuscript.

#### ACKNOWLEDGMENTS

Shell Brasil funded humpback whale tagging in Brazil. Great white shark tracking in New Zealand was funded by: New Zealand Ministry of Business, Innovation and Employment; New Zealand National Institute of Water and Atmospheric Research; New Zealand Department of Conservation. Northern elephant seal tracking was funded in part by U.S. Office of Naval Research grant N00014-08-1-1195 to DC and conducted under NMFS Permit 14636.

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars. 2017.00422/full#supplementary-material

with pigeon homing–a GPS tracking study. J. Exp. Biol. 217, 4057–4067. doi: 10.1242/jeb.108670

whales Megaptera novaeangliae in the Southwest Atlantic Ocean. Mar. Ecol. Prog. Ser. 313, 295–304. doi: 10.3354/meps313295

Zerbini, A. N., Andriolo, A. R., Heide-Jørgensen, M. P., Moreira, S. C., Pizzorno, J. L., Maia, Y. G., et al. (2011). Migration and summer destinations of humpback whales (Megaptera novaeangliae) in the western South Atlantic Ocean. J. Cetacean Res. Manag. 3, 113–118.

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2017 Horton, Hauser, Zerbini, Francis, Domeier, Andriolo, Costa, Robinson, Duffy, Nasby-Lucas, Holdaway and Clapham. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

### NOMENCLATURE

**Anomalistic month—**the period required for the moon to return to perigee following the preceding perigee (average = 29.6 days) **Apogee (lunar)—**the point in its orbit at which the moon is farthest from Earth

**Bouguer gravity anomaly (spherical)—**location dependent difference in gravitational acceleration between normal gravity and observed gravity caused by variations in the shape and density of Earth (for a more explicit definition, see Balmino et al., 2012)

**D (magnetic declination)—**location dependent angle in the horizontal plane between Earth's magnetic field and geographic north expressed positive to the east

**Declination (lunar)—**angle between the moon and the celestial equator

**F (magnetic field)—**location dependent magnetic field flux density (informally: intensity) of Earth's magnetic field expressed, in Standard International Units, as nanotesla (nT); radial distance coordinate {ρ} in spherical coordinate space {ρ, θ, ϕ} descriptions of Earth's magnetic field

**Gal—**unit of gravitational acceleration (1 Gal is equivalent to 1 cm s−2 )

**gB—**spherical Bouguer gravity anomaly

**gL—**latitude-dependent gravitational acceleration

**H (magnetic)—**horizontal component of Earth's magnetic field (F)

**I (magnetic inclination)—**location dependent angle in the vertical plane between Earth's magnetic field and the horizontal (expressed positive downwards)

**Perigee—**the point in its orbit at which the moon is closest to Earth

**PTT—**platform transmitting terminal

**Sidereal month—**the period required for the moon to complete one full orbit relative to a fixed star's position (average = 27.3 days)

**SPOT—**smart position or temperature transmitting tag

**Synodic—**the period required for the moon to complete one full illumination/phase cycle (average = 29.5 days)

**Tropical month—**the period required for the moon to complete one full orbit relative to the ecliptic (average = 27.3 days)

**X (magnetic)—**geographic north-south component of Earth's magnetic field (F) in the horizontal plane

**Y (magnetic)—**geographic east-west component of Earth's magnetic field (F) in the horizontal plane

**Z (magnetic)—**vertical component of Earth's magnetic field (F)

# Thermal Imaging and Biometrical Thermography of Humpback Whales

Travis W. Horton1 \*, Alice Oline2 , Nan Hauser 3 , Tasnuva Ming Khan4 , Amelie Laute3 , Alyssa Stoller 3 , Katherine Tison3 and Peyman Zawar-Reza5

*1 Department of Geological Sciences, University of Canterbury, Christchurch, New Zealand, 2 Environment Program, Colorado College, Colorado Springs, CO, United States, 3 Center for Cetacean Research and Conservation, Rarotonga, Cook Islands, 4 Cornell University, Ithaca, NY, United States, 5 Department of Geography, University of Canterbury, Christchurch, New Zealand*

#### Edited by:

*Lars Bejder, Murdoch University, Australia*

#### Reviewed by:

*Nuno Queiroz, University of Porto, Portugal Phil Bouchet, Centre for Marine Futures, School of Animal Biology, University of Western Australia, Australia*

> \*Correspondence: *Travis W. Horton<EMAIL_ADDRESS>
#### Specialty section:

*This article was submitted to Marine Megafauna, a section of the journal Frontiers in Marine Science*

Received: *28 September 2017* Accepted: *12 December 2017* Published: *21 December 2017*

#### Citation:

*Horton TW, Oline A, Hauser N, Khan TM, Laute A, Stoller A, Tison K and Zawar-Reza P (2017) Thermal Imaging and Biometrical Thermography of Humpback Whales. Front. Mar. Sci. 4:424. doi: 10.3389/fmars.2017.00424* Determining species' distributions through time and space remains a primary challenge in cetacean science and conservation. For example, many whales migrate thousands of kilometers every year between remote seasonal habitats along migratory corridors that cross major shipping lanes and intensively harvested fisheries, creating a dynamic spatial and temporal context that conservation decisions must take into account. Technological advances enabling automated whale detection have the potential to dramatically improve our knowledge of when and where whales are located, presenting opportunities to help minimize adverse human-whale interactions. Using thermographic data we show that near-horizontal (i.e., high zenith angle) infrared images of humpback whale (*Megaptera novaeangliae*) blows, dorsal fins, flukes and rostrums record similar magnitude brightness temperature anomalies relative to the adjacent ocean surface. Our results demonstrate that these anomalies are similar in both low latitude and high latitude environments despite a ∼16◦C difference in ocean surface temperature between study areas. We show that these similarities occur in both environments due to emissivity effects associated with oblique target imaging, rather than differences in cetacean thermoregulation. The consistent and reproducible brightness temperature anomalies we report provide important quantitative constraints that will help facilitate the development of transient temperature anomaly detection algorithms in diverse marine environments. Thermographic videography coupled with laser range finding further enables calculation of whale blow velocity, demonstrating that biometrical measurements are possible for near-horizontal datasets that otherwise suffer from emissivity effects. The thermographic research we present creates a platform for the delivery of three important contributions to cetacean conservation: (1) non-invasive species-level identifications based on whale blow shapes and velocities recorded by infrared videography; (2) reduced ship-strike rates through automated thermographic cetacean detection systems deployed in high traffic areas; (3) monitoring the spatial and temporal distributions of endangered animals in remote habitats.

Keywords: humpback whale, infrared, biometrical thermography, emissivity, blow velocity

## INTRODUCTION

Since the International Whaling Commission (IWC) banned commercial whaling in 1986, many baleen whale species have shown signs of recovery (Thomas et al., 2016). However, all species listed as Least Concern under the International Union for Conservation of Nature and Natural Resources (IUCN) Red List of Threatened Species framework also include threatened subpopulations classified as Vulnerable, Endangered, or Critically Endangered (Thomas et al., 2016). The conservation status, recovery and health of whale populations is very much site and context specific: modern human threats to whales, including ship strikes and entanglement in fishing gear, are not evenly distributed with respect to both space and time (Thomas et al., 2016).

Despite the spatially and temporally dynamic challenges associated with cetacean conservation and protection in the postwhaling era, progress has been made. For example, revisions to shipping lane positions, vessel traffic management plans and mandatory maximum vessel speeds along the eastern coast of North America correlate with significant reductions in North Atlantic right whale (Eubalaena glacialis) deaths due to ship strikes (Laist et al., 2014; Thomas et al., 2016). However, the risks associated with ship strikes remain high elsewhere. Necropsies performed on stranded whales demonstrate that at least one humpback whale (Megaptera novaeangliae), one fin whale (Balaenoptera physalus), and two blue whales (Balaenoptera musculus) are killed by ship strikes off the California coast every year (Redfern et al., 2013). Similar analyses suggest that, on average, one Bryde's whale (Balaenoptera edeni) is killed every year by ship strikes in New Zealand's Hauraki Gulf (Constantine et al., 2015). True cetacean mortality rates due to human activities at sea are almost certainly higher (Kraus, 1990), however, and the annual loss of even a single individual can be significant for smaller populations of long-lived species with low recruitment rates (Laist et al., 2001).

In an effort to reduce the risks ships pose to large whales, the IWC has developed a 3-year (2017–2020) Strategic Plan that seeks to increase the development and use of whale avoidance technologies (Cates et al., 2017). Acoustic and infrared automated cetacean detection systems are attractive and emerging tools for enhanced cetacean conservation (Zitterbart et al., 2013; Nowacek et al., 2016). The ability to detect whale blows, as far away as 5 km using around-the-clock 360◦ infrared scanners outfitted with rigorous detection algorithms (Zitterbart et al., 2013), will benefit many, including marine mammal observers onboard large vessels and land-based scientists studying whale movement behavior (e.g., Perryman et al., 1999) and human-whale interactions along rapidly changing coastlines (e.g., Graber, 2011). Infrared thermography can also facilitate the non-invasive collection and monitoring of fundamental biometrical information, including thermal physiology, injury diagnoses and population surveys (McCafferty, 2007).

Infrared cetacean detection systems also create opportunities for conservation biologists and cetacean ecologists to document the spatial and temporal distribution of animals utilizing remote or inaccessible environments. For example, the Oceania subpopulation of humpback whales, the only migratory humpback whales in danger of going extinct (Childerhouse et al., 2008), seasonally inhabit ∼10 million km2 of the tropical South Pacific Ocean. Yet, only a handful of scientists, spread across an area of ocean the size of China, actively study these whales. Automated detection systems have the potential to dramatically improve our knowledge of when and where these endangered whales are utilizing highly understudied breeding/calving ground habitats.

However, thermal imaging also has several important limitations. Infrared imaging systems are not inexpensive, particularly so for current high sensitivity models with cryogenically cooled detectors or large focal lengths capable of long-range applications. Infrared detectors also require a direct line of site to the target, yet they can also lose functionality through interaction with sea-spray. The data streams generated by infrared imaging systems are large, creating challenges with data handling, analysis and signal processing. Thermal cameras are also highly inaccurate when imaging scenes from nearhorizontal positions due to emissivity effects (Masuda et al., 1988; Cuyler et al., 1992; see Nomenclature).

The effects of emissivity on the brightness temperatures recorded by a thermal imaging device are extremely relevant to cetacean thermography. For example, as a whale exhales, its breath pushes sea water present in the near-surface water column, or nasal depression, or both, into the overlying atmosphere. From observation points at or near sea-level, this spouting of water droplets immediately and drastically changes the angle at which the whale's blow is being measured by the thermal camera. For example, a 2 m high whale blow will be measured perpendicularly (i.e., measured at a 0◦ zenith angle) from an observation point located 100 m distant and 2 m above sea level. In contrast, the adjacent flat ocean's surface will be measured sub-horizontally at an 89◦ zenith angle. Similar to blows, emergent dorsal fins, flukes or rostrums will also be measured at a relatively low zenith angles in relation to the adjacent ocean's surface. These rapid changes in the angle at which the object is being imaged will have large effects on the surface brightness temperatures estimated by the thermal imaging device due to the effect zenith angle of observed radiation has on sea water emissivity (Masuda et al., 1988).

The research we present was driven by three primary objectives, all aligned to the IWC's strategic goal of developing large whale avoidance technologies. We sought to: (1) quantify infrared image brightness temperature and brightness temperature anomaly (BTA) values for humpback whale blows, dorsal fins, flukes, and rostrums in both tropical breeding/calving ground and sub-polar feeding ground habitats; (2) calculate humpback whale blow height and blow velocity through coupling of infrared videography with laser-range finding; (3) evaluate the effects of emissivity on thermal imaging data collected from high zenith angle (i.e., oblique to target) positions. Achievement of these objectives creates a platform from which a variety of cetacean conservation tools can be further developed and delivered.

#### METHODS

Thermal images of humpback whale surfacing features were collected using a Forward Looking Infrared camera (FLIR A615, FLIR Systems, Inc.,) and analyzed using FLIR Tools+ software (FLIR Systems, Inc.,). The FLIR A615 we used had a focal length of 24.6 mm, 25◦ × 19◦ field of view, F-number of 1.0, infrared resolution of 480 × 640 pixels and a detector pixel pitch of 0.017 mm pixel−1 . The camera's detector comprised an uncooled Vanadium Oxide (VoX) long–wavelength (i.e., 7.5– 14µm) microbolometer (see Nomenclature) with a thermal sensitivity of <0.05◦C. Infrared images were captured every 0.04 s (i.e., 25 Hz) but frame rates as high as 200 Hz can be achieved with the A615's high-speed windowing option. The A615 was powered by a small 12-volt battery externally strapped to the camera's casing. The A615 was also connected to a FZ-G1 ToughPad tablet computer (Panasonic Corporation) via a high-speed USB cable. A GoPro Hero4 camera (GoPro, Inc.,) was affixed to the top of the A615 for contemporaneous visible wavelength image collection. This study was carried out in accordance with the recommendations of the Cook Islands Government. The protocol was approved by the Office of the Prime Minister, Cook Islands Government.

In Rarotonga, infrared and visible wavelength images were collected either ∼2 m above the ocean surface while onboard a Cook Islands Whale Research vessel, or from shore-based positions ∼5–10 m above sea level on the island's northwest coast (**Figure 1**). In Sitka Sound, all images were recorded ∼4 m above the ocean surface while onboard a commercial whale watching cruise arranged by the Sitka Sound Science Center as part of the annual Sitka Whale Fest (e.g., **Figures 1D,E**). Despite these variable imaging heights, our entire dataset was collected at >85◦ zenith angles (i.e., <5 ◦ off horizontal) due to the range in distances at which whales were imaged. A Nikon Forestry Pro laser rangefinder was used to determine whale distances whenever possible. All measurements were made during Beaufort wind force scale numbers 2–4 and similarly ranked World Meteorological Organization (WMO) Sea State codes.

Brightness temperatures were extracted from individual thermal images using the line measurement tool available in FLIR Tools+. Two lines for temperature data extraction were drawn across each image: the first line was drawn vertically through the background scene immediately adjacent to the targeted whale feature (i.e., blow, dorsal fin, fluke, rostrum, **Figure 1**), and the second line was drawn vertically such that it passed through the maximum brightness temperature included within the targeted whale feature. Thermal benchmarks included within each image, such as the steep thermal gradient across the ocean–atmosphere boundary, were used to align the pixels included in each line's thermal profile (**Figure 2**). Once aligned, the brightness temperatures recorded by each line were subtracted from each other in order to calculate BTA-values at the individual pixel scale for each whale feature analyzed (**Figure 2**).

Because the A615's pixel pitch and focal length were known, independent measurement of whale distances by laser rangefinding allowed us to estimate blow height from thermal image pixel measurements by combining the optical lens equation,

$$
\left(\frac{1}{\frac{blow}{distance}\{m\}}\right) + \left(\frac{1}{\frac{image}{distance}\{m\}}\right) = \frac{1}{\frac{focal}{length}}\tag{1}
$$

with the magnification equation,

$$
\left(\frac{image}{height} \text{(m)}\right) = \left(\frac{image}{distance} \text{(m)}\right)
$$

$$
\left(\frac{blue}{distance} \text{(m)}\right)
$$

and the camera's pixel pitch,

$$pixel\,pitch\left(\frac{m}{pixel}\right) = \frac{image}{image} \, (m) \tag{3}$$

Equation (1) can be rearranged to,

$$\frac{1}{\frac{image}{distance} \text{(m)}} = \left(\frac{1}{\frac{focal}{length} \text{(m)}}\right) - \left(\frac{1}{\frac{blue}{distance} \text{(m)}}\right) \tag{4}$$

Equation (2) can be arranged to,

$$
\begin{split}
\langle \text{llow}\_{\text{left}} \rangle\_{\text{(m)}} &= \left( \frac{1}{\left( \begin{matrix} 1 \\ \text{diameter} \end{matrix} \right)} \right) \times \left[ \begin{pmatrix} \text{image}\_{\text{(m)}} \\ \text{height} \end{pmatrix} \right] \\
&\times \left( \begin{matrix} \text{blow} \\ \text{distance} \end{matrix} \right) \end{split}
\tag{5}
$$

and equation (3) can be rearranged to,

$$\underset{height}{image}(m) = \text{pixel pitch}\left(\frac{m}{pixel}\right) \times \underset{height}{image}(pixel). \tag{6}$$

Substituting equations (4) and (6) into equation (5) gives,

$$\begin{aligned} \mathop{\text{bigow}}\limits\_{height} \{ m \} &= \left[ \left( \frac{1}{\underset{length}{\text{(}m\text{)}}} \right) - \left( \frac{1}{\underset{distance}{\text{(}m\text{)}}} \right) \right] \\ &\times \left[ \left( \underset{pitch}{\text{(}pixel} \left( \frac{m}{pixel} \right)} \times \underset{height}{\text{(pixel)}} \right) \right] \\ &\times \left( \underset{distance}{\text{(}ldistance} \left( m \text{)} \right) \right] \end{aligned} \tag{7}$$

parentheses (i.e., 0◦–10◦C) correspond with sub-polar thermal image brightness temperatures shown in (E).

(G) images of a nostril and adjacent rostrum at ∼10 m distance in tropical waters; visible (H) and infrared (I) images of a blow, rostrum and dorsal fin at 40 m distance in tropical waters; visible (J) and infrared (K) images of a footprint at ∼30 m distance and 50 s following fluke in tropical waters. Temperature scale numbers in

which simplifies to,

$${}^{block}\_{
left}\left(m\right) = \left[\frac{\begin{pmatrix} \text{pixel}\left(\frac{m}{\text{pixel}}\right) \times \text{image}\left(\text{pixel}\right) \times \text{distance}\left(m\right) \\ \text{left}}{\text{length}}\right)\right]$$

$$-\left(\frac{\text{pixel}\left(\frac{m}{\text{pixel}}\right) \times \text{image}\left(\text{pixel}\right)}{\text{height}}\right).\tag{8}$$

Blow heights were estimated using Equation (8) every 0.04 s following blow initiation. Image pixel heights were measured using FLIR Tools+ and blow distances were measured by laser range-finding as described above.

#### RESULTS

In total, we determined BTA profiles for 174 humpback whale blows, 20 dorsal fins, 9 flukes, and 20 rostrums. An equivalent number of whale features were analyzed from each of the two study areas, with the exception of flukes, for which 6 were imaged in Alaska and only 3 were imaged in Rarotonga. Of the 87 blows analyzed in each study area, 32 Rarotonga blows and 16 Alaska blows were imaged at distances <150 m. Of these, only 10 blows from each study area were recorded in the 100–150 m range.

Average BTA profiles demonstrate that humpback whale blows, dorsal fins, flukes and rostrums appear as thermal anomalies of similar magnitude relative to adjacent ocean water (**Figure 3**). For example, 100–150 m distant blows in Rarotonga and Alaska appear as 20–30 pixel-wide thermal anomalies that are ∼3 ◦C warmer than the adjacent ocean (**Figure 3A**). Similarly, dorsal fins and flukes in both areas exhibited maximum BTA values ca. 3–4◦C (**Figures 3B,C**), whereas rostrums from both populations were ∼2–3◦C warmer than the adjacent ocean (**Figure 3D**). Ocean water temperatures were measured by perpendicular thermography and satellite observations in both study areas. These measurements indicate surface ocean water temperature was ∼24◦C in Rarotonga, and ∼8 ◦C in Sitka Sound, Alaska, at the time thermal images were recorded.

The shapes of the average dorsal fin, fluke and rostrum BTA profiles differ because these features were recorded across a large range of distances in each study area. Because the Rarotonga whales were generally imaged at closer ranges, the dorsal fin, fluke and rostrum BTA profiles are spread across a larger number of image pixels than the Sitka BTA profiles (**Figures 3B,D**). In other words, the Rarotonga whale features fill a larger portion of the 640 × 480 pixel thermal images because these images were recorded at closer distances. Despite these distance-related differences in BTA profile shape between the study areas, the maximum BTA values for humpback whale blows, dorsal fins, flukes and rostrums (indicated by arrows in **Figure 3D**) we recorded are not significantly different (p >> 0.05, two-tailed t-test, **Figure 3**).

Laser range-finding enabled quantification of the relationship between the pixel-length of individual blows and blow distance for the FLIR A615. As expected, blow pixel-lengths are larger for images recorded at closer range, and blow pixel-length decreases