Ocean Heat Content Reveals Secrets of Fish Migrations

For centuries, the mechanisms surrounding spatially complex animal migrations have intrigued scientists and the public. We present a new methodology using ocean heat content (OHC), a habitat metric that is normally a fundamental part of hurricane intensity forecasting, to estimate movements and migration of satellite-tagged marine fishes. Previous satellite-tagging research of fishes using archival depth, temperature and light data for geolocations have been too coarse to resolve detailed ocean habitat utilization. We combined tag data with OHC estimated from ocean circulation and transport models in an optimization framework that substantially improved geolocation accuracy over SST-based tracks. The OHC-based movement track provided the first quantitative evidence that many of the tagged highly migratory fishes displayed affinities for ocean fronts and eddies. The OHC method provides a new quantitative tool for studying dynamic use of ocean habitats, migration processes and responses to environmental changes by fishes, and further, improves ocean animal tracking and extends satellite-based animal tracking data for other potential physical, ecological, and fisheries applications.


Introduction
Movements and seasonal migration ecology of vertebrate animals are largely driven by dynamic environmental features [1][2][3][4]. A strong relationship between oceanic temperature patterns and the migratory forays of many large pelagic tunas, billfishes, sharks, and tarpon has long been recognized [2,5,6]. Fronts and eddies are well known dynamic ocean habitats for fishes [7][8][9], but knowledge of their utilization has been based mostly on qualitative evidence from ecological theory or extractive fisheries data [7,8,10]. However, significant advances in the geolocation accuracy of electronic tags over the past two decades as a result of improved GPS and Argos positioning used on sea turtles [11,12], sea birds [13,14], and marine mammals [15,16] have provided direct quantitative evidence that many marine vertebrates use ocean fronts and eddies as preferred habitats. In these studies modeled ocean current vectors and satellite remote sensing data (e.g., sea surface height (SSH, altimetry data); sea surface temperature (SST, infrared and microwave); and, sea surface chlorophyll (SSC, ocean color imagery) were used to distinguish these important ocean features. However, due to the relatively short surface intervals of most fishes, it has been difficult to fully utilize the advanced features of geolocation tags to produce accurate positions. A few recent studies have used these Argos geolocation tags on sharks [17,18], striped marlin (Kajikia audax) [19] and Atlantic tarpon (Megalops atlanticus) [18], fishes which spend time at or near the surface. Nevertheless, due to their infrequent surfacing, geolocation data lack the temporal and spatial resolution required to study utilization of mesoscale fronts and eddies. As a default, some studies have used the lightlevel derived geolocations from electronic tags and ocean models to explore the large-scale ocean feature utilizations by marine predator fishes [20,21].
Over the past decade, thousands of pop-up satellite archival tags (PSATs) with high frequency temperature-depth-light sensors have been deployed on coastal and pelagic fishes [2,6,22,23]. PSATs provide previously unavailable concurrent observations of fish movements and hydrographic conditions over relatively large ocean areas. Geolocation estimates between deployment and pop-up positions are derived from light, time and depth data recorded by the PSAT [24]. Inherent restrictions limit the accuracy of light-based longitudinal locations to ± 1 degree [24], and improvement of these estimates has been possible by re-processing the lightbased positions using a Kalman filter (KF) and an unscented Kalman filter (UKF), then incorporating SST (KF-SST, UKF-SST) [25]. Despite advanced KF-SST and UKF-SST methods, the homogenous distribution of SST (Fig 1a) in tropical regions greatly denigrates usefulness of SST in refining movement tracks and in explicitly delineating fish utilization of dynamical ocean features. Among the variables (currents, SSH, SST, SSC) used to derive ocean features, only SST is directly recorded by the satellite tags. In this paper, we introduce a new method to estimate movements, and habitat association with ocean fronts and eddies, of electronically monitored marine animals, while avoiding the deficiencies inherent in SST-based geolocation estimates.
Ocean heat content (OHC), a fundamental measure of the heat stored in the upper surface layers, has been used for more than four decades by oceanographers and meteorologists to predict hurricane intensity [26,27]. The minimum SST for hurricane formation is 26°C [27]. Expressed in kilo joules cm -2 , OHC is estimated by integrating temperature from the depth of the 26°C isotherm to the ocean's surface. OHC is typically obtained from satellite remote sensed SST data and vertical thermal structure measurements. For hurricane intensity forecasting, OHC is routinely estimated by an Atlantic regional temperature and salinity (SMARTS) climatology model that systematically merges daily thermodynamic fields from satellite remote sensing (i.e., radar altimetry and SST) with ocean observing system data [28,29].
While evaluating the usefulness of OHC estimated from PSATs as a supplemental data source for SMARTS, we discovered that time-series OHC maps appeared to be superior for revealing pelagic fish habitats as compared to maps based solely on SST. During the tropical cyclone season (June 1 -November 1), SST maps of the Caribbean Sea and Gulf of Mexico regions of the central western Atlantic are generally featureless during most days (Fig 1a). In contrast, OHC maps revealed unseen structure in the upper thermocline of the water column due to the integrated heat content (Fig 1b). Thus, the goal of this research was to introduce and develop a new methodology using OHC to refine movement tracks of PSAT-monitored large pelagic fishes and to use examples to demonstrate how the refined movements tracks help us to explore the quantitative relationships between fish migration paths and dynamic oceanic features (fronts and eddies).
previously published [30][31][32][33][34][35] and unpublished studies. Light intensity data were processed using the global positioning software WC-AMP (Wildlife Computers, Redmond, WA, USA) to provide initial geolocations [36]. The light data provided sunrise and sunset times used to generate initial geolocation estimates. Then an SST-corrected Kalman filter [25,37] (KF-SST or UKF-SST) was applied to refine light-derived geolocations. To relocate the points that were either on land or in shallow waters, we modified the estimates by also filtering depths based on 2 minute grid ETOP02 bathymetry data and the daily maximum depth from the tag [38].
Profiles of depth and temperature (PDT) along the migration tracks of individual fishes were mapped based on the continuous PSAT data binned at 3-h time intervals (Fig 2a). The midpoint of temperature minimum and maximum at each depth from PDT data was used to generate the vertical temperature profile of each tagged fish. The missing data bins were linearly interpolated with the nearest neighbor bins for both vertical and horizontal dimensions. However, in our analysis when a tag had, on average, more than three gaps of five days each, it was eliminated from further analysis. These eliminations comprised about 10% of the tags originally available to us. For fish track analysis, OHC values (Fig 2b) were calculated by integrating the thermal energy from the depths of the 26°C isotherm (D26) to the ocean surface according to the following equation [26]: where, c p is the specific heat constant of water, ρ is sea water density, and T z is water temperature at depth z. For fishes inhabiting cooler waters, the OHC concept still holds, and simply requires adjusting the lower thermal limit downwards such as 20°C for tropical tunas and billfishes, or 10°C for temperate tunas (e.g., bluefin). We generated OHC maps (OHC M ) based on the global hybrid coordinate ocean model (HYCOM) by the US Navy Oceanographic Office (NAVO), using the real-time ocean data assimilation method [39]. We used the standard output, which has a spatial resolution of 0.08 deg (~7 km) and a daily time step. Global HYCOM has been shown to accurately resolve mesoscale processes such as meandering currents, fronts, and oceanic eddies [39]. We used the above equation and HYCOM temperature data to estimate the OHC, then we re-gridded the OHC values using linear interpolation methods to 2x2 min (0.033 x 0.033 degree) spatial resolution maps (OHC M ) to match the ETOP02 [40] bathymetric grid.
We determined that migration tracks could be further optimized using OHC maps and a practical genetic optimization algorithm (GA) [41]. Our GA-OHC filtering method was similar to methods described previously [42]; however, we substituted SST maps with OHC maps. Here, we defined OHC T as the OHC value estimated from the tag at time t, and OHC M as the OHC value derived from the HYCOM model at position P t at time t from either a KF-SST or UKF-SST derived track. In theory, if the position P t , OHC T and OHC M are estimated without error, then OHC T and OHC M should have the same value (i.e., ΔOHC = OHC M -OHC T = 0). If the position P t is not accurate, ΔOHC will not equal to zero. Thus, we needed to find a new position P 0 t where ΔOHC will be zero or minimized, a simple task if a single position was considered. However, an individual fish movement track over several months consists of hundreds to thousands of auto-correlated positions. Thus, the minimization procedure had to satisfy all track positions. To accomplish this goal, we applied a GA to minimize the sum of the square root of ΔOHC for all track positions. The optimization search method is schematically illustrated in Fig 3. At each time step, a search area is defined by a circle of a radius R, a variable ranging from 4-50 km that decreases as a function of generation time. R was also constrained by either land or shallow bathymetry (i.e., where water depths were less than the maximum depth of fish at the time t). The current position P t was indicated by the dot at the center of the circle, and P t-1 was the fish's position at previous time step, while P t+1 was the fish's position at the subsequent time step. The mid-point between P t-1 and P t+1 is indicated by the center position (CP). In the simulations, the new position (P 0 t ) takes a particular pixel according to one of nine scenarios: 1. The pixel with minimum absolute difference between the tag OHC (OHC T ) and the model 2. The pixel with shortest distance to the center of the circle (P t ) from the subset of minimum ΔOHC (a subset pixels where ΔOHC is > (ΔOHC min + 0.1ΔOHC min ).
3. The pixel with longest distance to the center of the circle (P t ) from a subset of minimum ΔOHC.
4. The pixel with shortest distance to the CP of previous and next time step from a subset of minimum ΔOHC.

5.
The pixel with the smallest OHC and shortest distance to the center of the circle P t .
6. The pixel with the smallest OHC and longest distance to the center of the circle P t .
7. The pixel with largest OHC and shortest distance to the center of the circle P t .
8. The pixel with largest OHC and longest distance to the center of the circle P t .
9. The pixel with shortest geometric average distance to previous (P t-1 ) and next (P t+1 ) position, that is, minð The nine scenarios spread pixels over the entire range within the circle to overcome potential local minimizations during a given simulation. The basic idea was to use a GA to optimize the  selection of the scenarios so as to not only minimize the sum square root of the error for an entire track, but also keep the distance between each position within the limit of known movement rates. At the beginning of GA track simulation for a particular fish, 25 sequences (vectors) of movement rules were created for a movement track of N positions, where each vector will have a "gene code" of length N that is randomly generated, i.e., the movement rule for each time step was randomly selected from one of the 9 scenarios described above. At the end of each generation (i.e., completion of N movements), a selection function (SF) was calculated for each of the 25 tracks as the root-mean-square (RMS) error of OHC between the observed OHC T and predicted OHC M for all N positions as For the next generation, the three tracks (populations P1, P2, P3) with the lowest SFs were selected as parents and these retained the same genes (same movement rules) as their parents (P1, P2, P3). For the other 22 populations, the first 11 populations had the genes crossed (i.e., mixing movement rules) between P1 and P2 at a random crossing point, and then for the remaining 11 populations the genes were crossed between P1 and P3. Finally, for these 22 populations, 10% of their genes were allowed to mutate (i.e., change in movement rules) at random. Simulations ran until the SF reached a constant level within a specified tolerance (<0.01) [41,42].
To increase speed and efficiency of the simulations, the OHC M data were generated ahead of the tracking GA procedure and saved into a data base for all tags. Additional OHC M data will be added to the database as they become available with time. In general, each simulation of a 100-day track took about 15 minutes to complete using currently available middle-tier specification desktop computers (i.e., Intel(R) Core™ i7 CPU @2.93 GHz).
To evaluate the accuracy of the GA-OHC method, a sinusoidal track extended over known OHC M maps in the Gulf of Mexico (GoM) was generated as the "true" unbiased track, and the OHC M values at the track positions were taken as the unbiased OHC T . Then, an altered track was generated by adding systematic sinusoidal biases of 0.5 degrees in latitude and random normal errors N(0,σ) on both the longitudes and latitudes of the "true" track. Finally, the GA-OHC method was applied to the altered track and unbiased OHC T to generate an estimated track, then compared with the true track. The "true" track was 45 days long and was simulated bi-monthly using 2007 OHC M data (Table 1) to evaluate the effects of different ocean conditions, i.e., seasonal changes in OHC. To evaluate the effect of OHC variation errors, an OHC average track was estimated by averaging all longitude and latitude values within a circle of a 0.5 degree radius at each point location where the OHC M values were within OHC T ±1 range. The mean and standard deviation of the differences between the average OHC track and the "true" track were calculated. Similarly, to repeat the validation simulations in a different geographic area, a 90 day blue marlin track from the South Atlantic Ocean was used as the "true" unbiased track, and the OHC M values at those track positions were taken as the unbiased OHC T . Actual simulations were performed quarterly using 2005 OHC M data.
For the purposes of demonstration and application of the new method, the refined movements tracks of 124 pelagic fishes monitored with PSATs were used to explore the quantitative relationships between fish migration paths and oceanic features (fronts and eddies). We used the Belkin-O'Reilly algorithm [43] to compute the gradient of OHC M and classified areas where gradients greater than 50 kJ cm -2 per grid cell were determined as zones of ocean discontinuities (ZOD). Fish track positions were then compared with the ZOD and all positions within one grid cell distance from ZOD were classified as being associated with ocean fronts and eddies. To simplify the presentation, the proportion of time that each fish spent associated with ZODs was calculated and categorized into three levels of eddy and front association: low (0-30%); medium (30-60%); and high (>60%). All numerical and visualization procedures were computed using IDL (Interactive Data Language, http://exelisvis.com), the procedures are available directly from the first author upon request. Fish handling and tagging procedures were followed Prince et al. [44], and was conducted under permits from the National Marine

Results
The UKF-SST filtered tracks (Fig 4a) were further refined with the GA that strives to minimize the RMS error between the observed OHC T and predicted OHC M . The resultant GA-OHC filtered track (Fig 4b) revealed fish movements with much greater detail than those tracks produced by the UKF-SST method alone. The evolution of learning for the GA-OHC filter is shown by a decrease in the RMS as a function of generation time (Fig 4c). The strong similarity of the final OHC values is shown by comparing the OHC T to the OHC M (Fig 4d).
Another example is presented for a blue marlin track (Fig 5) in the South Atlantic Ocean. The GA-OHC filtered track (Fig 5a) revealed fish movements with much greater detail than those tracks produced by the UKF-SST method (blue line). The 95% confidence interval (grey shaded area) was calculated for each location using ±5 kJ cm -2 OHC variation around the estimated location. The confidence intervals are plotted as probability ellipses with 2 standard deviations of longitude, latitude as the major axes and angle of rotation determined by the covariance. The comparisons of longitude and latitude from the UKF-SST and GA-OHC track estimation methods are shown (Fig 5b and 5c). The final OHC M values (red crosses) are optimized very close to the values of OHC T (Fig 5d) as compared to the OHC values from the UKF-SST track (blue line). The validation results of the simulated sinusoidal track (Fig 6a) showed that the GA-OHC filter improved the precision for determining the "true" track. The mean difference between the estimated track and the "true" track was less than 0.001 degrees and less than 0.01 degrees standard deviation for longitude and latitude during all periods of the year. The results were confirmed through validation of tagged fish that showed similar results in accuracy (Fig 6b). The effect of OHC variation (Table 1) was demonstrated by comparing the OHC average track to the "true" track. The mean differences in both longitude and latitude for all validation simulations were less than 0.04 degree, which is consistent with the 2-minute grid resolution (0.033 degree) of the model. For the GoM track and all periods of the year, the mean of the standard deviation (msd) ranged from 0.145 to 0.211, and 0.154 to 0.221, with OHC variation of ±1 and ±5 kJ cm -2 , respectively. For the South Atlantic Ocean track, the msd ranged from 0.144 to 0.185, and 0.162 to 0.203, with OHC variation of ±1 and ±5 kJ cm -2 , respectively. The OHC msd values (Table 1) within the circle of the 0.5 degree radius from each GA-OHC position ranged from 12.74 to 25.81 kJ cm -2 , and 6.84 to 8.57 kJ cm -2 , for the GoM and South Atlantic Ocean, respectively. This suggests the presence of a year-round OHC gradient for both geographic areas. On average, the OHC gradient is much steeper for the GoM than the South Atlantic Ocean.
The GA-OHC filtered tracks revealed detailed movements of fishes that were not apparent from track estimates using UKF-SST alone. For example, the detailed movements of one yellowfin tuna (Thunnus albacares) around an eddy and the Loop Current (LC) in the GoM are shown (Fig 7). From March 27 to April 9, 2012, this tuna was moving among weak fronts off the Mississippi River (Fig 7a) before reaching an eddy centered in the GoM on April 10 th . The eddy structure is discernible on the OHC M and PSAT temperature-depth profile (Fig 7c). The   fish left the eddy on April 23 (Fig 7b). For the next two months it moved between weak fronts on the continental slope. On June 30 th , it arrived on the western side of a newly formed eddy that was shed from the LC (Fig 7d), then proceeded around the boundary to the east side of the eddy by July 13 (Fig 7e).
The 124 pelagic fishes monitored with PSATs, including yellowfin tuna, blue marlin (Makaira nigricans), white marlin (Tetrapturus albidus), and sailfish (Istiophorus platypterus), exhibited a strong preference for fronts and eddies, the dominant features of ocean discontinuities (Table 2). From 2009 to 2013 the area of "ocean discontinuities" accounted for about 18% of total available habitat (range of 10.6% to 27.1%, mean of 17.9%, standard deviation of 3.8%) in the GoM and Atlantic Ocean (50°W to 100°W longitude, 0°to 50°N latitude). More than 50% of these tagged fish spent greater than 60% of their time associated with ocean discontinuities. To further highlight these associations, we provide here some individual illustrative tracks. The yellowfin tuna previously described showed clear preferences for the edges of the eddy (Fig 8a) and the LC (Fig 7e). Another yellowfin tuna exhibited a high affinity for the edge of the same eddy located in the center of the GoM (Fig 8b). Remarkably, this particular tuna swam around the periphery of the eddy many times over 20 days (March 22 to April 11), rarely crossing over the eddy. One bluefin tuna was observed using the eastern edge of the LC while in the GoM, then the western edge of the Gulf Stream as it migrated north in the western North Atlantic Ocean (Fig 8c). Front and eddy edge use were also observed for billfishes. Blue marlin appeared to prefer the edges of cold-and warm-core eddies, and also the LC (Fig 8d). Similarly, the edge of the Gulf Stream front was used by white marlin off Cape Hatteras (Fig  8e). Finally, a sailfish monitored in the Florida Keys used the edge of the Florida Current, then transitioned to the edge of the LC while moving north-northwest in the GoM (Fig 8f).

Discussion
Our validation results (Fig 6 and Table 1) clearly demonstrated that the GA-OHC is effective in elucidating the "true" track from biased location data. The GA technique was very efficient in addressing the complexities that have hindered calculus-based methods [41]. No previous fish tracking studies [47,48] have incorporated GA methods to match tag SST with satellite SST data for the purpose of improving location estimates from the light level data. Teo et al. [47] reported a mean error ranging from 0.73 to 1.41 degrees and standard deviation ranging from 0.54 to 1.28 degrees, while no error estimates were provided by Domeier et al. [48]. Nielsen et al. [37] and Lam et al. [25] included SST data into the KF and UKF state-space model to estimate the most probable tracks. They reported a standard deviation of 0.5 degree in longitude, and 1.2 degree in latitude. Our validation results are not directly comparable to earlier studies as regards accuracy basically because we lacked Argos positions from double-tagged fish for direct comparisons; however, these results did give a relative sense of the method's precision. Our GA-OHC method had mean differences less than 0.04 degrees between the "true" and estimated tracks for both longitude and latitude, and msd of 0.15-0.22 degrees (~15-25 km) for both longitude and latitude when considering OHC errors of ±1 and ±5 kJ cm -2 (Table 1). This represents a large reduction in position error, which is much smaller than mesoscale eddies (<100 km [49,50]). While all of these methods do not necessarily require the presence of a spatial gradient in the data per se (i.e., SST or OHC), it greatly facilitates model fitting. In tropical and sub-tropical ocean regions SST data are more or less uniform throughout large areas (Fig 1a), while the OHC data revealed the ocean with mesoscale features such as fronts and eddies (Fig 1b). Our study showed that OHC gradients are present year-round in the GoM and South Atlantic Ocean (Table 1) with much higher OHC gradients in the GoM. When the objective is to explore the relationship between fish migration paths and ocean features, an exact geolocation (i.e., errors in SST or OHC maps) may be less important than when the fish are in the proximity of these features. Use of an OHC gradient of 50 kJ cm -2 to discriminate ocean features should produce robust estimates of fish associations with those dynamic ocean features.
Besides SST and OHC, other variables such as SSH, SSC and ocean current vectors, could be used to differentiate mesoscale ocean features in oceanographic and ecological studies [10,13,43,51]. However, we did not use them for fish tracking purposes because these variables were not recorded or estimated by any satellite fish tags under present technologies. The OHC is a unique metric because it can be estimated from ocean circulation models and satellite tags, despite not being measured or recorded directly by the sensors presently aboard either satellites or the electronic fish tags. By combining the power of the genetic algorithm and OHC, the GA-OHC filtering method could provide, with a little refinement, a powerful tool for improving the accuracy of movement tracks for satellite tagged fishes.
Through the novel application of OHC in this study, a clear connection emerged between the behaviors of fish and spatiotemporal dynamics of oceanographic features. The GA-OHC filtered tracks suggested that fishes use oceanographic habitat features as ecological road maps to guide migrations. It was not clear whether they were responding directly to the physical features, or were searching for prey, where the prey also may have responded to the physical features, or both. Ocean fronts and eddies support high biological productivity at all trophic levels [7,52], and are well known as "hotspots" of marine fish habitats [4,53]; however, most of this knowledge has been garnered from ichthyoplankton surveys [54], recreational and commercial fisheries catch data, or ecological theory [7,8]. Compared to UKF-SST tracks (Fig 4a), the GA-OHC track methodology (Fig 4b) greatly improved the level of detail elucidating fish movements in a dynamic ocean environment (Figs 7 and 8) over those previously presented for large pelagic fishes monitored with electronic tags [2]. We believe the reason GA-OHC generated more realistic movement tracks and revealed previously unknown levels of front and eddy use during fish migrations was not only due to higher track resolution, but also because OHC integrates the vertical thermal properties of oceanographic features in a way that is similar to what fishes do when they naturally swim through those features and use ambient information to discriminate zones of ocean discontinuity.
These detailed movement tracks present tremendous implications for the use of satellite tagging data to study fish behavior in relation to habitat features in the seemingly monotonous seascape of the vast ocean. For example, the GoM is the only known spawning ground for the bluefin tuna in the western North Atlantic Ocean. The bluefin tuna track (Fig 8c) showed the fish was located at the eastern edge of the LC, an area known for spawning activity based on high concentrations of newly hatched larvae [54] and adult bluefin tuna catches [55]. The idea that ocean fronts and eddies are "hotspots" of marine fish habitats is not necessarily a new Table 2. Front and eddy utilization by satellite-tracked fish. Summary of proportion of total time at liberty associated with fronts and eddies by individual PSAT-tagged yellowfin tuna, blue marlin, white marlin, and sailfish based on GA-OHC filtered tracks. Individual tracks were grouped into high, moderate and low utilization.  concept [4,9,53] per se, what is novel is the explicit empirical association of fishes with these prominent oceanographic features. The inherent physical characteristics of fronts and eddies are known to concentrate, retain, and enhance the productivity of living resources [7,8,56]. While fronts and eddies may only account for a small fraction of the vast ocean, they play a major role in shaping the seascape of marine ecosystems as key zones of intense predator-prey interactions across every trophic level. This study demonstrates that OHC T estimated from PSAT data are precise and can be useful as another integrated data source to refine OHC maps and improve the prediction of oceanic states in models. Development of a satellite tag capable of GPS quality location and instantaneous transmissions of position and multi-sensor data could be highly beneficial to the SMARTS system, resolving smaller-scale oceanic variability in the OHC maps that could improve hurricane intensity forecasting [28]. Improving OHC measurements at the boundaries of fronts and eddies is essential to the generation of accurate and precise OHC maps. Because large pelagic fishes have an affinity for these areas, using a network of satellite-tagged fish to monitor OHC values and augment other platforms is highly advantageous and cost-effective, as compared to the more traditional use of XBTs, floats, drifters, and gliders that are very difficult to keep in these boundary zones. Fish could greatly assist the resolution of structures around oceanic fronts and eddies at much smaller scales where important air-sea interactions occur under strong winds.
Given the potential availability of this new GA-OHC method, re-analyses of existing satellite tagging data will almost certainly shed new light on movement and behavior, and afford scientists a powerful tool to develop and investigate new hypotheses. Besides facilitating improvements in track prediction, OHC maps may also be applicable to other aspects of marine ecosystem research. These might include, for example, replacement of SST maps with OHC as a primary ecosystem indicator variable to help define, identify, and model essential habitats and improve estimates of the spatial abundance distribution for fishery abundance forecasts [57], conservation and management policy [58,59], reproductive ecology [60], and climate changes [61].