ORIGINAL RESEARCH article

Front. Mar. Sci., 30 April 2026

Sec. Marine Megafauna

Volume 13 - 2026 | https://doi.org/10.3389/fmars.2026.1808805

Integrating behavioral movement and environmental preferences to map critical habitat of whale sharks using long-term satellite tracking in the Indo-Pacific Ocean

  • 1. Department of Geography, Faculty of Mathematics and Natural Sciences, Universitas Indonesia, Depok, Indonesia

  • 2. Focal Species Conservation Program, Ocean and Science Department, Konservasi Indonesia, South Jakarta, Jakarta, Indonesia

  • 3. Department of Oceanography, Faculty of Fisheries and Marine Science, Universitas Diponegoro, Semarang, Indonesia

  • 4. Center for Coastal Rehabilitation and Disaster Mitigation Studies, Universitas Diponegoro, Semarang, Indonesia

  • 5. Elasmobranch Institute Indonesia, Denpasar, Bali, Indonesia

  • 6. Conservation International Asia-Pacific, Auckland, New Zealand

  • 7. Department of Biology, Faculty of Mathematics and Natural Sciences, Universitas Indonesia, Depok, Indonesia

Abstract

Understanding the spatial ecology of wide-ranging marine megafauna is essential for identifying critical habitats and designing effective conservation strategies. Whale sharks undertake extensive movements far beyond well-studied aggregation sites, yet the spatial structure and environmental drivers of these movements remain poorly resolved. Here we integrate long-term satellite tracking data from 2015 to 2025 for 70 whale sharks tagged at four major aggregation sites in the Indonesian archipelago, including Cenderawasih Bay, Kaimana, Saleh Bay, and the Gulf of Tomini, with state-space modeling and MaxEnt habitat suitability analyses to provide a basin-scale assessment of whale shark movement ecology and critical habitats in the central Indo-Pacific. Behavioral states were classified into foraging and migratory movements and then modeled across regions, seasons, sexes, and life stages using oceanographic and seafloor geomorphic predictors. Whale shark habitat selections were strongly structured by behavior and region. Aggregation sites were dominated by foraging behavior and characterized by shallow productive habitats with predictable prey availability, whereas non-aggregation regions functioned primarily as migratory corridors shaped by mesoscale oceanography and seafloor geomorphic features such as canyons and escarpments. Across demographic groups, year-round suitable habitat was limited and largely confined to a few aggregation sites, particularly Cenderawasih Bay and Saleh Bay, highlighting their role as irreplaceable functional habitats. Seasonal non-aggregation habitats in the Flores Sea, Ceram Sea, Timor Sea, and Banda-Arafura transition zone provided important but transient opportunities for foraging and migration. Distinct sex and life-stage specific habitat preferences indicate the need for demographic-specific conservation approaches. Overall, our findings demonstrate that effective whale shark conservation requires combining site-based protection of persistent aggregation habitats with connectivity focused and transboundary management that accounts for seasonal movements across national waters and Areas Beyond National Jurisdiction.

1 Introduction

The whale shark (Rhincodon typus), the largest fish in the world, is a highly migratory species found throughout tropical and warm-temperate waters. The species is currently listed as Endangered on the IUCN Red List following population declines exceeding 50% over the past three generations (). Although conservation measures have expanded globally, including protection of key aggregation sites and national legal protections, significant knowledge gaps remain in identifying the broader habitats used by whale sharks outside these focal areas. This gap is particularly concerning because whale sharks undertake extensive movements that frequently extend beyond managed areas, exposing them to cumulative anthropogenic threats such as fisheries interactions, vessel strikes, and unregulated tourism activities (; ). Improving spatial understanding of whale shark habitat use is therefore essential for guiding effective conservation and management actions.

Spatial ecology provides a critical framework for addressing this challenge by linking species movement patterns with environmental variability to identify where and when important habitats occur (). For highly mobile marine megafauna such as whale sharks, understanding spatial patterns is particularly important because their wide-ranging movements across ocean basins (; ). Mapping species distribution, movement, and habitat preferences can therefore inform a range of conservation applications, including identifying ecological hotspots (), guiding marine spatial planning () and marine protected area design (), assessing vessel collision risk (), supporting sustainable marine tourism (), and improving fisheries management ().

Hosting over 60% of the global whale shark population () and encompassing the broadest spatial distribution (), the Indo-Pacific subpopulation is a critical region for advancing understanding of whale shark spatial ecology. However, existing habitat suitability studies remain geographically limited or rarely account for behavioral states and demographic groups (; ; ; ; ), constraining our understanding of life history specific critical habitats. Long term satellite tracking, coupled with state space modeling approaches that infer behavioral states from movement persistence index, enables the development of robust habitat suitability models for migration and feeding that capture variation across time and space (). This habitat suitability provides strong conservation value for supporting the design of MPAs and transboundary management strategies that align with the ecological needs of whale sharks across behaviors, sexes, and life stages ().

We integrate more than ten years of satellite tracking data from 70 whale sharks across four major aggregation sites in Indonesia, capturing movements across the Indonesian Archipelagic Waters and into the Pacific and Indian Oceans. This represents one of the longest and most extensive whale shark satellite tracking studies in the world. We (1) classify behavioral states using state-space models to distinguish migratory and foraging movements and their environmental associations; (2) develop habitat models to map critical habitats and quantify environmental drivers of behavioral movement and it’s habitat suitability across regions, seasons, sexes, and life stages; and (3) highlight conservation implications for area-based and transboundary management. This integrative approach provides the world first a basin-wide assessment of whale shark movement ecology and critical habitats, delivering actionable insights for the conservation of this globally threatened species.

2 Materials and methods

2.1 Study site and environmental setting

The study area spans nine locations across the Indo-Pacific, comprising four whale shark aggregation sites in Indonesia, namely Saleh Bay (SB), the Gulf of Tomini (GO), Kaimana - Raja Ampat (KA), and Cenderawasih Bay (CB), which represent coastal habitats, and five non-aggregation regions representing archipelagic and oceanic habitats, including the Southeastern Indian Ocean (IO), the Indonesian Archipelagic Seas (INDAS), and the Northern and Southern Pacific Oceans (NPO and SPO). These areas encompass 18 marine zones ranging from semi-enclosed seas to major ocean basins and cover 13 national jurisdictions plus adjacent high seas (Figure 1; Supplementary Figure 1).

Figure 1

).

The aggregation sites reflect Indonesia’s ecological diversity. Saleh Bay is influenced by Indian Ocean dynamics (), the Gulf of Tomini combines coastal and deep basin characteristics (), Kaimana lies along the shallow Arafura shelf (), and Cenderawasih Bay is connected to the South Pacific and Bismarck Seas (). These sites provide relatively stable and productive habitats that support feeding aggregations of juvenile whale sharks (; ; ). In contrast, non-aggregation regions are environmentally heterogeneous, characterized by deeper bathymetry, canyon and escarpment systems, stronger seasonal oceanographic variability, and dynamic currents that structure migratory corridors and episodic foraging opportunities (Supplementary Figures 39).

2.2 Data collection, processing, and analysis

A conceptual workflow summarizing the main analytical steps used in this study, including data cleaning, behavioral state classification, environmental variable selection, and habitat modeling, is presented in Supplementary Figure 10. The workflow illustrates how satellite tracking data were integrated with environmental predictors to identify patterns of habitat preference and spatial distribution. Detailed methodological settings and parameterization are provided in the Supplementary Methods.

2.2.1 Satellite tracking

Between 2015 and 2025, a total of 70 whale sharks were satellite tagged, comprising 64 large juvenile males (LJM), one adult male (AM), and five large juvenile females (LJF). The dataset was strongly dominated by LJM individuals (91%). Tagging was conducted at four major aggregation sites in Indonesia: Saleh Bay (21 LJM, one AM, two LJF), Cenderawasih Bay (34 LJM), Kaimana (seven LJM, one LJF), Raja Ampat (one LJF), and the Gulf of Tomini (three LJM).

All sharks were equipped with fin-mounted SPLASH satellite tags (Wildlife Computers, USA) programmed to transmit Argos location data upon surfacing. Tagging procedures followed previously established protocols (; ), including morphometric measurements and sex determination based on clasper morphology. Argos locations were classified into standard accuracy classes, retaining classes 3, 2, 1, 0, A, and B, while class Z records were excluded.

To account for location errors and temporal autocorrelation, Argos tracks were analyzed using a correlated random walk state-space model (crw-SSM) implemented in the R package aniMotum (). The model estimated regularized whale shark locations using a 72-hour time step (), selected based on model performance (Akaike Information Criterion) and biological relevance.

Behavioral states were inferred using a move persistence model, where the gamma (γ) parameter (ranging from 0 to 1) was used to distinguish area-restricted search (foraging; close to 0) from directional movement (migration; close to 1; ). Differences in behavioral states across regions and seasons were evaluated using chi-square and Fisher’s exact tests.

Environmental variables associated with whale shark locations were extracted using ArcGIS Pro spatial analysis tools. Environmental differences among regions, seasons, behavioral states, sexes, and life stages were evaluated using Kruskal–Wallis tests followed by Dunn’s post hoc tests, and visualized using ridge plots in R ().

2.2.2 Environmental predictors

A suite of oceanographic and seafloor geomorphic variables was compiled to characterize environmental conditions associated with whale shark movement and habitat suitability (Table 1). These variables were used as predictors and were selected based on their ecological relevance to prey availability, ocean circulation, and habitat structure.

Table 1

ParameterHypothesisSpatial resolutionTemporal resolutionUnitData source
Oceanography
Sea surface temperature (SST)Influences the distribution of marine megafauna based on species-specific thermal preferences and tolerance limits ().4 kmSeasonal composites; monthly data (2002-2024)°C()
Sea surface chlorophyll-a (SSC)Proxy for primary productivity, indicating plankton availability as a food source ().4 kmmg.m-3
Sea surface height (SSH)Variations in sea surface elevation driven by oceanographic processes may affect prey availability and thermal structure ().27.75 kmSeasonal composites; monthly data (2015-2023)mCopernicus Marine Service*
https://doi.org/10.48670/moi-00148
Eddy kinetic energy (EKE)Enhances local productivity and redistributes prey, influencing megafauna aggregation patterns ().27.75 kmcm2/s2
Surface current velocityAffects spatial distribution of marine animals through plankton transport, energy flux, and navigation efficiency of migratory species ().9.21 kmm/sCopernicus Marine Service**
https://doi.org/10.48670/moi-00021
Geomorphology
Bathymetry (depth)Physical constraint influencing vertical distribution based on species-specific depth preferences and tolerance ().500 mStatic (2023)m()
Seafloor slopeRepresents seabed gradients that can induce localized upwelling and prey concentration in productive zones ().Vector data converted into raster distance layers (predictor variables) using Euclidean distance at a 4-km spatial resolution.N.Adegrees
Distance to ridgeGeomorphic complexity associated with current flow modification and enhanced upwelling, increasing local productivity and providing navigational cues for migration ().km()
Distance to seamountSeamounts enhance vertical mixing and prey aggregation, potentially acting as ecological hotspots ().km
Distance to escarpmentSteep seabed features influence circulation patterns and prey distribution, contributing to habitat heterogeneity ().km
Distance to canyonSubmarine canyons facilitate nutrient transport and prey concentration, potentially shaping movement and foraging behavior ().km

Environmental and seafloor geomorphic variables used in habitat suitability modeling.

*Copernicus Marine Service, SEALEVEL_GLO_PHY_L4_REP_OBSERVATIONS_008_047.

**Copernicus Marine Service, GLOBAL_MULTIYEAR_PHY_001_030

Environmental datasets were harmonized to consistent spatial grids suitable for regional analyses. Interpolation procedures and spatial resolution settings used to generate predictor layers are described in Supplementary Methods. To minimize multicollinearity among predictors, correlation analyses were conducted prior to model construction, and highly correlated variables were excluded following standard species distribution modeling practices. The final set of predictors used in each model is provided in Supplementary Table 2.

2.2.3 Habitat suitability modelling

Habitat suitability was modeled using a maximum entropy [MaxEnt; 3.4.1; ()] with presence-only whale shark occurrence data derived from satellite tracking. To capture ecological variability across demographic and behavioral groups, occurrence data were stratified by region, behavioral state, sex, life stage, and season, resulting in 73 habitat possible models (Supplementary Table 2). This stratified approach was applied to explicitly characterize distinct environmental preferences associated with each behavior and sex-life stage across regions under differing seasonal conditions, rather than assuming a single homogeneous response.

Prior to modeling, tracking data were spatially filtered to reduce spatial autocorrelation and sampling bias. Background points representing available environmental conditions were generated within the accessible area of tracked sharks. Detailed procedures for spatial filtering, background selection, and model parameterization are provided in Supplementary Methods. Model performance was evaluated using the area under the receiver operating characteristic curve (AUC). Habitat suitability maps were converted to binary predictions using a 10-percentile training presence threshold, and the resulting binary outputs were used to estimate seasonal habitat suitability frequency.

3 Results

A total of 10,958 Argos locations were obtained (Supplementary Figure 1) from whale sharks tagged at Saleh Bay (5,174), Cenderawasih Bay (3,807), the Gulf of Tomini (1,337), and Kaimana (488). Tracking durations ranged from 23 to 990 days (mean ± SD = 434 ± 263 days), with the longest duration from Saleh Bay (472 ± 296 days) and Cenderawasih Bay (433 ± 244 days; Supplementary Table 1). The dataset was dominated by large juvenile males (94% of locations), with smaller contributions from adult males (1%) and large juvenile females (5%).

Most individuals from all aggregation sites moved into one or more non-aggregation regions (62-100%), with no cross visitation among aggregation sites. In the Indonesian Archipelagic Seas, whale sharks originated primarily from Saleh Bay (14 individuals, ~45%), followed by Cenderawasih Bay (7, ~23%), Kaimana (7, ~23%), the Gulf of Tomini (2, ~6%), and Raja Ampat (1, ~3%). Use of the Arafura and Timor Seas was dominated by individuals from Kaimana (7, ~58%), with additional contributions from Cenderawasih Bay (2, ~17%), Saleh Bay (2, ~17%), and Raja Ampat (1, ~8%). More site-specific patterns were observed in oceanic regions, where sharks from Cenderawasih Bay accounted for all movements into the North and South Pacific Oceans, while individuals from Saleh Bay were the primary users of the southeastern Indian Ocean (Supplementary Figure 1).

3.1 Behavioral movements

From 10,958 satellite-transmitted locations, invalid quality class “Z” records and temporal duplicates were removed, retaining 6,170 locations (56%) from 64 individuals. Of the 70 tagged whale sharks, six individuals were excluded from the crw-SSM analysis due to insufficient tracking data. The mean γ value was 0.55 (median = 0.58, SD ± 0.24) and was used as the threshold for behavioral classification. Locations with γ > 0.55 were classified as transiting behaviors associated with migratory movements (52.45%, 3,236 of 6,170), while γ < 0.55 indicated area-restricted search (ARS; 48%, 2,934 of 6,170) associated with foraging activities (Figure 1). Swimming speed differed significantly between foraging (0.32 ± 0.56 m/s) and migratory behavior (0.96 ± 0.97 m/s; Wilcoxon rank-sum test, W = 2,437,787, p < 0.0001).

Spatial patterns of γ values revealed clear functional differences in whale shark habitat use. Aggregation sites were characterized by lower γ values (mean 0.49 ± 0.24) and were dominated by foraging behavior (59%), whereas non-aggregation regions showed higher γ values (0.66 ± 0.21) with predominantly migratory movements (71%). Kaimana exhibited the lowest mean γ (0.39 ± 0.22) and the highest proportion of foraging (84%; Figure 1, Supplementary Figure 11), while all aggregation sites showed mean γ < 0.55, including Cenderawasih Bay (0.49 ± 0.26; 58% foraging), Saleh Bay (0.50 ± 0.23; 58%), and Tomini Bay (0.51 ± 0.25; 62%). Among non-aggregation regions, only the South Pacific Ocean had a mean γ below 0.55 (0.54 ± 0.21; 58% foraging), whereas the highest migratory dominance occurred in the Southeastern Indian Ocean (0.78 ± 0.18; 74%), North Pacific Ocean (0.69 ± 0.17; 85%), Indonesian Archipelago Seas (0.64 ± 0.20; 68%), and Arafura-Timor Seas (0.62 ± 0.20; 61%).

Seasonal effects were significant across most regions Supplementary Tables 46), including Cenderawasih Bay (χ² = 176.57, df = 3, p < 0.0001), Saleh Bay (χ² = 17.84, df = 3, p < 0.001), and all non-aggregation regions (p < 0.001), with temporary seasonal shifts toward foraging observed in the Arafura-Timor Seas, Indonesian Archipelago Seas, and South Pacific Ocean, while Kaimana showed no significant seasonal variation (p = 0.22) despite persistent foraging dominance.

3.2 Habitat suitability

Overall, the performance of the 73 MaxEnt models demonstrated a high level of reliability, with a mean Area Under the Curve (AUC) value of 0.82 (SD = 0.08; Supplementary Table 2). The details of habitat suitability for each aggregation site and non-aggregation region are described below.

3.2.1 Aggregation sites

Across aggregation sites, only four out of 11 predictors, including sea surface chlorophyll-a (SSC), sea surface temperature (SST), bathymetry, and slope were consistently available at the spatial scale of aggregation areas to characterize whale shark behavioral movement (Figure 2A). Seafloor geomorphic complexity such as escarpments and canyons was only present in the Gulf of Tomini. These predictors consistently differentiated whale shark habitats across aggregation sites, behavioral states, demographic groups, and seasons (Supplementary Tables 68).

Figure 2

3.2.1.1 Saleh Bay

Habitat suitability models for Saleh Bay demonstrated moderate predictive accuracy (0.75 ± 0.05) and were strongly structured by dynamic oceanographic conditions. SSC peaked during the southeast monsoon from June to August (e.g., foraging large juvenile females: 0.75 mg/m³; migrating large juvenile females: 0.65 mg/m³; Supplementary Table 6) and consistently emerged as the primary driver of foraging habitat selection across demographic groups (juvenile females: 30 ± 12%, adult males: 40 ± 10%, juvenile males: 27 ± 2%; Figure 3). This pattern promoted fine scale habitat partitioning along depth and productivity gradients among sex and life stage classes (Figure 4). SST decreased significantly during the southeast monsoon (28.47-28.57 °C; p < 0.0001). Migratory sharks occupied shallower bathymetry (-53.09 m [-31 to -95.25 m]) than foraging sharks (-79.48 m [-16.41 to -111.81 m]), and most environmental differences were significant across seasons (p < 0.001; Supplementary Table 7). Overall, 91% (77-97%) of Saleh Bay was identified as suitable habitat, dominated by extensive migration corridors used by large juvenile males (up to 96% spatial coverage, with 41% suitable year-round; Supplementary Table 12).

Figure 3

Figure 4

3.2.1.2 Cenderawasih Bay

In Cenderawasih Bay, habitat suitability showed moderate predictive accuracy (0.81 ± 0.03) and reflected moderately productive, persistently warm within shelf slope environments. SSC differed significantly between behaviors, with lower values during foraging (0.71 [0.45-1.06 mg/m³]) than migration (1.17 [0.64-2.44 mg/m³]; Z = -5.04, p < 0.0001), and emerged as the dominant predictor for both foraging (47 ± 17%) and migratory habitats (43 ± 7%). Foraging sharks occupied shallower bathymetry (-82 [-48 to -124 m]) than migrating sharks (-143 [-108 to -208 m]; Z = 15.47, p < 0.0001). SST remained uniformly warm across seasons (30.86 [30.43-31.10 °C]), while seasonal variability was strongest in SSC, peaking during seasonal transition I from March to May (p < 0.001; Supplementary Table 8). Suitable habitats covered much of the bay (74-78%), but year-round suitability was limited (12-18%), indicating strong seasonal structuring of habitat use (Figures 5A–H).

Figure 5

3.2.1.3 Kaimana

Model performance in Kaimana was moderate (0.73 ± 0.03). Habitat use was centered on shallow, highly productive coastal waters with limited seasonal variability. Foraging large juvenile males occupied high productivity environments (SSC = 1.12 mg/m³), significantly exceeding those in other aggregation sites (p < 0.0001). Bathymetry was the shallowest among sites (-27 [-22 to -33 m]; p < 0.05-0.0001) with low to moderate slopes (0.85 [0.7-1.0°]), while SST remained uniformly warm (30.81 [30.59–31.03 °C]). Seasonal variation in SSC was not significant (p > 0.05), contrasting with patterns in Cenderawasih Bay and Saleh Bay. SSC consistently emerged as a major predictor in MaxEnt models (28 ± 21%). Spatially, suitable habitats covered most of the area (84%), with substantial seasonal overlap (63%) concentrated in Triton Bay, supporting Kaimana’s role as a stable seasonal feeding aggregation dominated by foraging large juvenile males (Figures 5I, J).

3.2.1.4 Gulf of Tomini

In the Gulf of Tomini, habitat suitability showed moderate accuracy (0.82 ± 0.06) and was centered on deep habitats shaped by complex seafloor geomorphology. Large juvenile males occupied the deepest environments among aggregation sites, particularly during migration (-893 [-249 to -1,563 m]; p < 0.0001). Migrating sharks occurred significantly closer to canyons than foraging sharks (8 and 30 km, respectively; Z = 4.91, p < 0.0001). Distance to canyon emerged as the dominant MaxEnt predictor for both migration (52%) and foraging (44 ± 3%). SSC remained low (0.18 [0.15-0.23 mg/m³]) and SST showed minimal seasonal variability (30.48 [30.42-30.53 °C]). Slopes were the steepest among aggregation sites (>2°; p < 0.0001), indicating that habitat differentiation was driven primarily by bathymetry and geomorphic features rather than productivity. Suitable habitat was limited for migration (10%) but extensive for foraging (66%), with minimal seasonal overlap (0.56%), reflecting strong seasonal shifts in feeding areas from coastal margins during the northwest monsoon from December to February (Figure 5L) to canyon associated habitats during seasonal transition II from September to November (Figure 5M; Supplementary Figure 3C).

3.2.2 Non-aggregation regions

At non-aggregation regions, whale shark movements were structured by a broader set of dynamic oceanographic variables, mesoscale processes, and complex seafloor geomorphology, with clear contrasts among regions, behaviors, sex life stages, and seasons (Figure 2B; Supplementary Tables 911).

3.2.2.1 Indonesian Archipelagic Seas

In the Indonesian Archipelagic Seas, habitat suitability models showed moderate performance (0.86 ± 0.08) and revealed strong behavioral and sex-based habitat segregation across bathymetry, geomorphic features, and oceanographic conditions. Large juvenile males exhibited clear behavioral segregation (p < 0.05-0.0001; Supplementary Table 10), with foraging occurring in shallower (-1,270 [-718 to 1,829 m]; Supplementary Table 9), more productive waters (0.63 [0.29-0.87 mg/m3]) with weaker currents (0.19 [0.10-0.30 m/s]), while migration shifted to deeper habitats (-2500 [-2,108 to 2,887 m]) closer to escarpments (15 [10–20 km]), characterized by stronger currents (0.24 [0.18-0.31m/s]) and higher mesoscale activity (540 [194–789 cm2/s2). Bathymetry was the main predictor of foraging habitat for large juvenile males (19 ± 10%). Migration habitats consistently follow deep slope corridors, particularly through the Flores Sea (Figures 6A–D). Sex specific differences were also evident during migration (p < 0.05-0.0001), with females occupying shallower (-1,360 [-803 to -2,595 m]), more productive waters (0.51 [0.21-0.80 mg/m3]) farther from canyons (34 [32–36 km]) and escarpments (33 [22–52 km]) than males (p < 0.05; Figures 6A-D, I-L; Supplementary Figure 3A). Foraging habitat of large juvenile females represented the largest critical habitat (74%; Supplementary Table 12) but was seasonally constrained, with limited year-round suitability (<1%), highlighting strong control by dynamic oceanographic variability.

Figure 6

3.2.2.2 Southeastern Indian Ocean

Habitat suitability models for the Southeastern Indian Ocean showed moderate performance (0.81 ± 0.06). Suitable habitats for both migratory and foraging behaviors were predominantly offshore (Figure 7), with seasonal shifts toward the southern Java coast during the southeast monsoon and seasonal transition II (Figures 7C, D, G). Migrating large juvenile males (-3,835 [-1,897 to -5,044 m]) and adult males (-2,722 [-1,673 to -3,770 m]) occupied extremely deep habitats across seasons, significantly deeper than other non-aggregation regions (p < 0.05-0.0001), with no bathymetric difference between life-stages (p > 0.05). Bathymetry consistently emerged as a key predictor of migration habitat suitability for adult males (24 ± 5%). Migrating adult male utilized habitats closer to canyons (45 [11–79 km], Z = -3.99, p < 0.01; Figure 7G; Supplementary Figure 3A), with gentler slopes (>2 [1.8-2.3°], Z = 3.63, p < 0.05) and stronger currents (0.27 m/s), while SSC remained generally low, peaking only during seasonal transition II for large juveniles (3.91 mg/m3). Estimated suitable habitat covered 47-85% of the region, dominated by extensive offshore migration corridors, whereas foraging habitats were limited and coastal (Figure 7E). Year-round suitable habitat was minimal (0.27%), indicating predominantly seasonal habitat use.

Figure 7

3.2.2.3 Arafura and Timor Seas

Maximum Entropy models for the Arafura and Timor Seas showed moderate performance (0.86 ± 0.07). Whale shark habitat use was concentrated in shallow areas near canyon and escarpment features, with enhanced productivity and SST cooling during the southeast monsoon (26.9 °C). The Banda-Arafura transition zone was consistently identified as highly suitable for both sexes and behaviors (Figure 8). SSC contributed to foraging habitat selection of large juvenile males (11 ± 6%), peaking during seasonal transition II (5.62 mg/m³), while eddy kinetic energy indicated strong seasonal mixing. Migrating large juvenile females occupied significantly deeper habitats (2,740 [-2,635 to -2,968 m]) than males (-986 [-143 to 1,668 m]; Z = -6.80, p < 0.0001) and occurred closer to canyons (12 [2–27 km], Z = -9.06, p < 0.0001) and escarpments (12 [10-15km], Z = -7.20, p < 0.0001), particularly during the southeast monsoon, when canyon contribution peaked (58%), and during seasonal transition II (24%). Estimate models further confirmed that suitable habitats were strongly structured by canyon features (Figures 8J, K; Supplementary Figure 3A). Suitable habitat coverage ranged from 4% to 79%, dominated by migratory habitats of large juvenile males, while year-round suitability was limited (0.15-3%), indicating strong seasonal forcing.

Figure 8

3.2.2.4 South Pacific Ocean

In the South Pacific Ocean, habitat suitability showed moderate accuracy (0.86 ± 0.10) and was strongly structured by mesoscale dynamics, strong currents, seasonal productivity, and slope related features. Foraging habitats dominated estimated suitability (86%) relative to migration (57%) and shifted seasonally from productive coastal waters along northern Papua during the northwest monsoon to offshore areas during the southeast monsoon (Figures 9E–G), with no year-round suitable habitat identified (seasonal overlap ≤10%). Foraging large juvenile males experienced significantly higher eddy kinetic energy than migrating individuals (526 vs 429 cm²/s²; p < 0.01), particularly during the northwest monsoon and seasonal transition I (p < 0.05 to 0.0001; Supplementary Table 11). EKE contributed more to foraging (25 ± 32%) than migratory models (9 ± 5%). Foraging sharks also occurred closer to escarpments (19 [7–43 km]), occupied shallower depths (-1,191 m [-506 to -2,814 m]), and experienced steeper slopes (2.74° [1.59-3.44]) than migrating sharks (all p < 0.0001), while current velocity remained high but did not differ between behaviors. SSC was elevated during foraging in the northwest monsoon and seasonal transition I (1.19 and 1.18 mg/m³) and declined sharply during the southeast monsoon (0.17 mg/m3). These patterns indicate strong seasonal shifts in slope associated foraging habitats driven by dynamic mesoscale processes.

Figure 9

3.2.2.5 North Pacific Ocean

Among all regions, the North Pacific Ocean showed the highest model performance (0.88 ± 0.05). Habitat use was concentrated in deep, oligotrophic offshore waters structured by strong currents and mesoscale dynamics, with limited behavioral differentiation. Migration habitats dominated estimated suitability (86%) relative to foraging (57%), with minimal year-round habitat (1%) and limited seasonal overlap (9-24%), indicating persistent migration corridors aligned with the North Equatorial Countercurrents (Figures 9A-D, 10A; Supplementary Figures 6, 7). Estimated foraging habitats during seasonal transitions I and II were largely distinct, with minimal overlap (0.11%; Figures 9F, H), despite both seasons targeting complex seafloor geomorphic features in the region. Migrating large juvenile males occupied some of the deepest habitats among non-aggregation regions (-3,907 [-3,168 to -4,599 m]; p < 0.0001) under high current velocities (0.29 m/s). SSC remained uniformly low, while eddy kinetic energy was elevated across behaviors and seasons, contributing strongly to habitat suitability for migration (20 ± 14%) and foraging (16 ± 14%). Sharks generally occurred far from major geomorphic features, such as canyons (330 [128–496 km]), escarpments (65 [12–132 km]), seamounts (120 [45–198 km]), indicating predominant use of dynamic pelagic corridors, although occasional associations with geomorphic features were observed near Micronesia and the Marshall Islands (Figure 1).

Figure 10

3.3 Distribution of critical habitat based on sex and life stage

3.3.1 Large juvenile males

Across 10.9 million km² of estimated migration habitat, the North Pacific Ocean dominated suitability (38%; Figure 10A; Supplementary Figure 12), while Indonesia held the largest national share (39%), followed by Micronesia (20%), high seas (15%), and Australia (11%). Year-round migration habitat for large juvenile males was extremely limited (123,790 km²; 1.14%) and concentrated mainly in the North Pacific (70%), particularly within the Marshall Islands, Indonesia, Palau, and Micronesia, with Indonesia contributing the largest share (39%). Tri-seasonal migration habitat covered 940,784 km² (9%), largely within the North Pacific and distributed across Indonesia (35%), high seas (21%), and Micronesia (19%).

Foraging habitat suitability spanned 5.7 million km², with the Southeastern Indian Ocean contributing the largest share (23%; Figure 10B; Supplementary Figure 13) but only seasonally. Indonesia encompassed the greatest share of suitable foraging habitat (39%), followed by Papua New Guinea (20%) and Australia (18%). Year-round foraging habitat was highly restricted (52,789 km²; 0.93%), concentrated in the Timor Sea (70%), primarily within Indonesian waters, while tri-seasonal foraging habitat (182,766 km²) was led by the Timor Sea, Banda Sea, and Arafura Sea.

3.3.2 Adult males

Estimated migration habitat covered 2.2 million km², with almost all suitable area (99%) located in the Southeastern Indian Ocean (Figure 11A; Supplementary Figure 14). The largest proportion of suitable habitat occurred in high-seas waters (42%), followed by Indonesia (27%) and Australia (11%). No year-round migration habitat was identified in this region. Instead, suitability was limited to bi-seasonal habitats totaling 596,647 km², predominantly within high-seas areas (45%), indicating strong seasonal dependence of migration space use in the Southeastern Indian Ocean.

Figure 11

Foraging habitat was highly restricted, covering 1,894 km² and occurring only in Saleh Bay (Flores Sea). This habitat was predominantly bi-seasonal (64%), with smaller proportions of single-season (23%) and tri-seasonal (13%) suitability (Figure 11B; Supplementary Figure 15).

3.3.3 Large juvenile females

Migration habitat suitability spanned approximately 1.2 million km², with the Banda Sea emerging as the primary region, accounting for 39% of the estimated area (Figure 12A; Supplementary Figure 16). Suitable habitats were predominately concentrated within Indonesian jurisdiction (98%), beside Timor-Leste and Australia. Year-round migration habitat was extremely limited (12,707 km²) and occurred exclusively in Indonesian waters, mainly within the Ceram Sea (33%) and Halmahera Sea (31%). In contrast, tri-seasonal migration habitats were more extensive and mostly distributed across the Molucca Sea (34% of 206,918 km²) and the Banda Sea (29%).

Figure 12

Foraging habitat suitability covered about 1.5 million km², concentrated largely in the Java Sea and Banda Sea (Figure 12B; Supplementary Figure 17). Nearly all suitable foraging areas were located in Indonesian waters (99.5%), with only a small fraction extending into Timor-Leste. Persistent habitat suitability was minimal, with only 0.02% year-round and 0.07% tri-seasonal, both primarily restricted to Saleh Bay.

4 Discussion

This study presents a basin-scale assessment of whale shark habitat use in the central Indo-Pacific by combining long-term satellite tracking, behavioral state, and habitat suitability modeling. Three key findings emerge. First, habitat use is strongly structured by behavioral state and region, with aggregation sites serving as predictable foraging habitats, while non-aggregation regions function as migratory corridors and opportunistic foraging grounds that were shaped by mesoscale oceanography and seafloor geomorphology. Second, year-round suitable habitat is limited and largely concentrated in a few aggregation sites, particularly Cenderawasih Bay and Saleh Bay. Third, habitat preferences differ across sexes and life stages. These findings highlight the importance of protecting key aggregation habitats while strengthening connectivity-focused management to support fisheries management, sustainable marine tourism, and vessel strike mitigation.

4.1 Methodological advances, performance, limitations, and improvements

The stratified modeling approach used in this study enabled the identification of region-specific environmental drivers of whale shark movements while accounting for differences in behavior and demographic groups (). These results demonstrate the value of segregating habitat models by behavior and demographic group to accurately capture life stage-specific habitat use and improve ecological niche for conservation planning ().

The modeling framework accounted for predictor collinearity, spatial autocorrelation, and sampling bias (; ; ), resulting in moderate to high model performance across most models. However, 21 of the 73 models were based on small sample sizes (<10 locations; Supplementary Table 2). Although several of these models still achieved moderate predictive performance, consistent with the robustness of MaxEnt in data-limited contexts (), their output should be interpreted with caution. For models derived from very limited occurrences, habitat suitability patterns should be considered preliminary indications rather than definitive ecological relationships, and further data collection will be required to validate these patterns.

Nevertheless, habitat suitability estimates for adult males and large juvenile females should be interpreted with caution, as the dataset is overwhelmingly dominated by juvenile males (~91%), with only a few tracked individuals representing females and adults. In particular, models for adult males and large juvenile females in Saleh Bay, the Indonesian Archipelagic Seas, the Arafura and Timor Seas, and the southeastern Indian Ocean were derived from only one to two individuals, limiting the robustness and generality of these habitat inferences. Consequently, the spatial patterns identified for these demographic groups should be considered preliminary and interpreted as indicative rather than definitive. Despite this limitation, extended tracking durations for adult males (up to 730 days) and large juvenile females (471 ± 371 days) still provide valuable insights into individual-level movement patterns and highlight priority areas for future research targeting underrepresented demographic groups.

Future research should prioritize expanding demographic representation, particularly adult males, females, and neonate whale sharks to at least >10 individuals per region to better capture behavioral variability (). Integrating prey-field or trophic data (), linking habitats with anthropogenic risk layers (; ), and coupling movement-based models with stranding and mortality data will further strengthen conservation prioritization, particularly in Indonesia, where whale shark stranding rates are among the highest globally ().

4.2 Functional segregation of whale shark habitats

Our results reveal clear functional segregation in whale shark habitat use between aggregation and non-aggregation regions and between foraging and migratory behaviors. Aggregation sites were dominated by area-restricted search behaviors, indicating localized feeding supported by predictable prey availability (; ; ) and fine-scale productivity associated with coastal inputs and lift-net fisheries (; ; ).

In contrast, non-aggregation regions were characterized primarily by migratory movements across broader pelagic areas. In these regions, whale sharks appear to track dynamic productivity features, such as upwelling zones, mesoscale eddies, and deep-water prey fields associated with canyon and escarpment systems. These findings support previous observations that whale sharks frequently move between predictable feeding hotspots and transient pelagic habitats (; ).

Most aggregation sites overlapped with lift-net fisheries, where sharks exploit concentrated prey resources (; ), likely representing an energetically efficient strategy in oligotrophic waters (; ; ) consistent with the energetic knife-edge theory (; ). Overall, aggregation sites function as stable feeding focal areas (; ), while non-aggregation regions serve as migratory corridors and opportunistic foraging habitats (; ; ). Differences in core habitat use among aggregation sites [Supplementary Figure 18 (, in prep)], despite connected broader home ranges, highlight their functional uniqueness and underscore the need to incorporate behavior-specific habitat roles into conservation and threat assessments.

4.3 From aggregation hubs to pelagic networks: structuring whale shark movement ecology

Recent observations from the Bird’s Head Seascape (Cenderawasih Bay and Kaimana) and Saleh Bay indicate that whale shark aggregations occur year-round, although individual residency is variable and generally short (; ). Residency metrics show clear differences among sites, with the lowest residency in Kaimana, intermediate levels in Saleh Bay, and the highest in Cenderawasih Bay. These contrasts likely reflect variation in aggregation-site carrying capacity (), driven by prey availability, competition, and human-mediated foraging associated with lift-net fisheries, consistent with satellite telemetry showing a higher proportion of individuals leaving Kaimana (Supplementary Figure 19).

Habitat structure and carrying capacity appear to be key drivers of whale shark movement ecology. Within the movement ecology framework (), transitions between aggregation and non-aggregation habitats arise from interactions between internal state, navigation capacity, motion capacity, and environmental drivers. Movements are influenced by prey availability () and habitat configuration (), while navigation responses to productive oceanographic features and geomorphic habitats further structure movement pathways (; ; ; ), with realized movement patterns reflecting trade-offs between energetic costs and prey encounter rates ().

In Kaimana, satellite tracking indicates strong seasonal redistribution from the aggregation site to surrounding pelagic habitats (Supplementary Figure 19). Use of the aggregation area declines markedly during the southeast monsoon, while movements increase toward the Arafura, Ceram, and Banda Seas. These shifts likely reflect reduced local foraging opportunities and increased competition, combined with opportunities to exploit monsoon-driven productivity associated with canyon systems and Indonesian Throughflow currents (; ). Movements into the Banda and Timor Seas also coincide with higher foraging speeds (0.85 m/s) compared with other non-aggregation habitats (0.70 m/s; Supplementary Table 13), suggesting movement pathways structured along canyon networks (Figure 8).

In contrast, whale shark habitat in Cenderawasih Bay remains suitable year-round, supported by river-enhanced productivity during the northwest monsoon and localized feeding opportunities around lift-net fisheries (Figures 5A–H). Abundant anchovies associated with lift-net fisheries provide a consistent food source for whale sharks (). The bay’s large and sheltered environment supports the highest residency among sites, although 62% of individuals temporarily leave the study area concentrated only around the southern area of the bay where lift-net fisheries are operating (). Most movements remain moderate and restricted to nearby Pacific waters and the Ceram and Banda Seas (Supplementary Figure 19). These departures likely reflect individual foraging strategies rather than strong carrying-capacity constraints, as sharks exploit alternative prey resources in deeper offshore habitats. Movements are influenced by regional environmental cues, including coastal upwelling (), mesoscale eddies (), and major current systems () that facilitate connectivity and optimize prey encounters ().

In Saleh Bay, suitable habitat is also predicted year-round, but residency durations are shorter than in Cenderawasih Bay (). Movements away from the aggregation site occur mainly during the northwest monsoon from December to February and seasonal transition I from March to May, when productivity declines locally (Supplementary Table 8). During these periods, sharks move toward nearby productive regions such as the Flores Sea () and southeastern Indian Ocean (), where coastal upwelling, mesoscale eddies, and complex seafloor geomorphology provide alternative foraging opportunities.

Tagging data from the Gulf of Tomini remains limited, but satellite tracking suggests that departures occur primarily during seasonal transition I and the southeast monsoon from June to August, when habitat suitability within the gulf decreases. Sharks appear to move toward the Molucca Sea, where seasonal upwelling and complex seafloor geomorphology create favorable foraging conditions (; ).

4.4 Sexual and ontogeny niche partitioning

Our findings support and extend evidence of strong demographic segregation in whale sharks by sex and life stage. Globally, coastal aggregation sites are dominated by juvenile males (; ; ), a pattern mirrored across all aggregation sites in this study (; ; ), whereas adult-dominated aggregations are rare and typically offshore (; ). Sex-specific growth trajectories, with males maturing earlier at smaller sizes and females growing larger and more slowly (), likely explain juvenile male dominance in productive coastal systems and greater female reliance on deep pelagic prey (; ).

Consistent with this, we observed marked sex and life stage-specific environmental preferences. In Saleh Bay, adult males occupied deeper, lower-SSC habitats than juveniles, suggesting trophic differentiation (; ). In the Arafura Sea, large juvenile females preferentially used deep habitats near canyons and escarpments, with limited surface foraging signals, indicating likely deep feeding. Movements of satellite tracked female whale sharks from Ningaloo (Western Australia) to the Arafura Sea during the southeast monsoon () supports our predictions. These patterns highlight the importance of incorporating demographic heterogeneity into conservation planning (; ).

4.5 Conservation implications

A key conservation finding of this study is the extremely limited extent of year-round suitable whale shark habitat across demographic groups and behaviors. Only a few aggregation sites, most notably Cenderawasih Bay and Saleh Bay, function as persistent year-round habitats, underscoring their role as irreplaceable ecological anchors within an otherwise highly dynamic Indo-Pacific seascape. This supports prioritizing long-term protection of these sites and aligns with findings from Tanzania, where predictable multi-year residency enables effective site-based management ().

Beyond aggregation areas, several pelagic regions including the Flores Sea, Ceram Sea, Timor Sea, and the Banda-Arafura transition zone serve as seasonal migration and foraging habitats. These regions likely provide important opportunities for sharks to exploit transient productivity pulses and may enhance population resilience by buffering seasonal prey variability (; ). At broader scales, whale shark habitat use is structured by major oceanographic features including the Indonesian Throughflow, regional current systems, and interconnected canyon networks. These physical features create migration corridors and seasonal feeding opportunities across the Indo-Pacific.

These findings demonstrate that site-based protection alone is insufficient without connectivity-focused and transboundary approaches spanning national waters and Areas Beyond National Jurisdiction (). Effective conservation will require strengthening MPA networks around persistent habitats (), implementing dynamic or seasonal management in key corridors (), and integrating habitat information into fisheries () and shipping governance (). Positioned as a regional connectivity hub, Indonesia can play a strategic role in advancing coordinated Indo-Pacific conservation aligned with global biodiversity targets ().

Statements

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material. further inquiries can be directed to the corresponding authors.

Ethics statement

The research permit was approved by the Cenderawasih Bay National Park Authority under permits SIMAKSI SI.18/BBTNTC-2/TEK/2015, SI.46/BBTNTC-2/TEK/2015, and SI.05/BBTNTC-2/TEK/2016 issued to Abraham Sianipar. Additional research and tagging permits were granted to Conservation International Indonesia and Konservasi Indonesia by the Raja Ampat MPA Management Authority, the Kaimana MPA Management Authority, the Department of Marine and Fisheries of West Nusa Tenggara (for Saleh Bay), and the Department of Marine and Fisheries of Gorontalo in coordination with BPSPL Makassar (for the Gulf of Tomini). All field activities were conducted in accordance with the ethical standards for animal care and use of the Research and Innovation Agency of Indonesia (Approval No. 199/KE.02/SK/10/2023) and complied with relevant local legislation and institutional requirements.

Author contributions

MP: Conceptualization, Data curation, Formal analysis, Investigation, Software, Validation, Visualization, Writing – original draft, Writing – review & editing. AW: Conceptualization, Supervision, Writing – review & editing. AS: Data curation, Funding acquisition, Project administration, Writing – review & editing. AH: Data curation, Writing – review & editing. ES: Data curation, Investigation, Writing – review & editing. IS: Data curation, Writing – review & editing. RM: Data curation, Writing – review & editing. ME: Conceptualization, Funding acquisition, Investigation, Project administration, Supervision, Writing – review & editing. JS: Supervision, Writing – review & editing. MM: Supervision, Writing – review & editing.

Funding

The author(s) declared that financial support was received for this work and/or its publication. The research was generously funded by Conservation International and Konservasi Indonesia donors, including the Sunbridge Foundation, the MacArthur Foundation, the Wolcott Henry Foundation, Save the Blue Foundation, The Alchemy of Change Fund, Stellar Blue Fund, The Charles Engelhard Foundation, The Paine Family Trust, the David and Lucile Packard Foundation, MAC3 Impact Philanthropies, and Sea Sanctuaries Trust. Individual donors included Audrey and Shannon Wong, Dawn Arnall, Marie-Elizabeth Mali, Ray and Barbara Dalio, Katrine Bosley, Michael Light, Daniel Roozen, Sarah Argyropoulos, and Jill Warnick. Additional support was provided by Ant International, SEA Aquarium Singapore, Mowilex, Citizen Watch, Saison Technology, and guests of the MV True North expedition vessel. The funders were not involved in the study design, collection, analysis, interpretation of data, the writing of this article, or the decision to submit it for publication.

Acknowledgments

The authors express their sincere gratitude to the Government of Indonesia, particularly the Ministry of Marine Affairs and Fisheries, the National Research and Innovation Agency, the Raja Ampat MPA Management Authority, the Kaimana MPA Management Authority, the Department of Marine and Fisheries of West Nusa Tenggara, and the Department of Marine and Fisheries of Gorontalo, as well as BPSPL Makassar and Denpasar, for their invaluable support of the whale shark satellite telemetry research program in Indonesia. Their strong commitment to marine science and conservation was essential to the successful implementation of this pioneering research. We also acknowledge the crucial contributions of Cenderawasih Bay National Park officers, Iman Tilahunga, Fahri Amar, and Nugraha Maulana, and the fishing communities of Cenderawasih Bay, Raja Ampat, Kaimana, Saleh Bay, and Gorontalo, whose local knowledge, logistical assistance, and technical support were fundamental to the whale shark tagging operations.

Conflict of interest

The author(s) declared that this work was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

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

References

Summary

Keywords

behavioral movement, habitat suitability, Indo-Pacific, state-space models, whale shark

Citation

Putra MIH, Wirasatriya A, Sianipar A, Hasan A, Setyawan E, Syakurachman I, Mambrasar R, Erdmann M, Supriatna J and Manessa MDM (2026) Integrating behavioral movement and environmental preferences to map critical habitat of whale sharks using long-term satellite tracking in the Indo-Pacific Ocean. Front. Mar. Sci. 13:1808805. doi: 10.3389/fmars.2026.1808805

Received

11 February 2026

Revised

16 March 2026

Accepted

26 March 2026

Published

30 April 2026

Volume

13 - 2026

Edited by

Xuelei Zhang, Ministry of Natural Resources, China

Reviewed by

Gonzalo Mucientes Sandoval, Spanish National Research Council (CSIC), Spain

Jun Liu, Sichuan University, China

Updates

Copyright

*Correspondence: Mochamad Iqbal Herwata Putra, ; Masita Dwi Mandini Manessa,

†Present address: Mark Erdmann, Re: Wild, Austin, TX, United States

Disclaimer

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.

Outline

Figures

Cite article

Copy to clipboard


Export citation file


Share article

Article metrics