We collated whale shark identification data from 25 whale shark sites in the Atlantic, Indian and Pacific Oceans (see Figure 1; Table 1). A summary of the data used, and the temporal scale covered at each site is presented in Table 1. Data were collected by photographing the left flank of the whale sharks, behind the gill slits and above the pectoral fin (Arzoumanian et al., 2005). Each image was assigned or matched to an individual whale shark, the latter referring to an individual already within a database. Every time a whale shark was encountered, or ‘identified, this was recorded on a log or spreadsheet with all encounter information within each study site.

Whale shark identification images were collected by researchers and citizen scientists, and collated and curated by in-country researchers (e.g. Figure 2). Identification data were uploaded onto Wildbook for Whale Sharks (WWS) for catalogue matching and validation, with the exception of Djibouti, Seychelles, Peru, Indonesia, Hawai'i, Bahia de La Paz (Mexico), Bahia de Los Angeles (Mexico), Thaa Atoll (Maldives) and South Ari Atoll (Maldives), who used localised catalogues to compare and identify individual animals assisted by program I³S Classic (Van Tienhoven et al., 2007). Both I³S and WWS use machine-automation to match spot patterns against a reference library – user-built by the former, and from global submissions by the latter. Here, we used a single sighting per shark per day across all sites.

We calculated Lagged Identification Rate (LIR) for each site independently, and produced site-specific estimates of abundance, residency and mortality or permanent emigration from the resulting models.

We first prepared the data for input into the program SOCPROG 2.9 (Whitehead, 2009), with each site comprising three columns: date, study site and whale shark id. Data was then inputted into SOCPROG, and we used the ‘Movement’ module in the program to estimate the LIR. We set the sampling period to ‘day’, to establish that each individual whale shark was recorded on a per day basis.

Eight models exploring different population scenarios, as pre-set in the program, including closed and open population with varying combinations of immigration, re-immigration and mortality, were then fitted to each dataset (Table 2). For example, Model H tests that the population is open and that there is emigration, re-immigration and mortality of individual animals within this population between sampling periods (days as set above), and produces estimates of abundance (per day), mean residency inside and outside the study site, and morality rate. Abundance is an estimate of the average number of individuals per sampling occasion (day). Mortality or permanent emigration refers to the probability of an animal dying or leaving the population permanently (Whitehead, 2019) and is used to calculate apparent survival (φ), which is calculated as 1 minus the mortality estimate. The models cannot distinguish between mortality or permanent emigration (Whitehead, 2019). Residency (in days) is estimated as the mean time spent within the study site (residency in), and the mean time spent outside the study site before returning to the site (residency out). Estimates of the latter are imprecise however (H. Whitehead, pers. comm.), which leads to our improved estimates of residency below.

We used the quasi-Akaike information criterion (QAIC) for model selection when there was data over-dispersion as determined by SOCPROG, and the Akaike information criterion (AIC) otherwise (Burnham & Anderson, 2002; Whitehead, 2007). All models were run simultaneously, and we selected the models following Whitehead (2019), where a ΔAIC or ΔQAIC value of ≤ 2 indicates substantial support for the model, whilst a value of 4–7 indicates considerably less support, and a value of > 10 indicates essentially no support (Burnham & Anderson, 2002). For the purpose of this study (comparable results and new residency metric) we used Model H for all sites (see Results). The best-fit model was then bootstrapped for 100 repetitions through random resampling with replacement to produce parameter confidence intervals and standard errors (Buckland & Garthwaite, 1991). The LIR was assessed from the plots (Whitehead, 2019), where a falling LIR value indicates individuals leave the site (emigration and mortality), and increases over time indicate individuals return to the site (re-immigration).

The model estimates of residency described above are influenced by sighting and resighting of individual animals at the study site, and we wanted to test model outputs by comparing them to a known residency value. This was possible because at Oslob, Philippines, (OSL) whale sharks have been monitored daily, year-round since 2012, allowing for near complete coverage. Daily photo-ID surveys (see Araujo et al., 2014 for details) at the site began on 31 Mar 2012 and concluded on 31 Dec 2019. We assume all sharks were counted at the site, given the effort (daily 3 x 1 hr surveys) and the relatively small size of the search area (~480 x 170 m²), but we acknowledge the possibility of missed sharks. Of all global 25 sites, this is the most complete coverage of a site with 2,791 days of sampling effort. Given its accessibility from shore, this was not replicable at any of the other sites.

To improve the residency values obtained from the LIR models described above, we used the empirical daily photo-ID data from Oslob. For example, a shark first identified on 31 Mar 2012, sighted on 299 additional days, and last identified on 31 Dec 2019, would have a residency index of 0.12 (300/2831) or mean 43.8 days per year (d.yr^-1). The value (2831) is the number of days between 31 Mar 2012 and 31 Dec 2019. We calculated this for all individual whale sharks at Oslob, except those first identified in 2019 to reflect inter-annual returnees at the site and not add bias to the probability of recapture.

With the empirical residency values from Oslob, we then calculated the ‘LIR residency proportion’ (∝_LIR) for all sites, as a2/a3, where a2 is the modelled residency time in, and a3 is the modelled time out of the study site from Model H in Table 2, to understand residency within study sites. We used the following linear equation to derive an adjusted residency time, in d.yr^-1, for all sites: residency = (∝_LIR * 124.6)/1.97, where 124.6 is the modelled residency time in for OSL (a2) and 1.97 is the ∝_LIR for OSL. We used a linear model to understand the relationship between adjusted residency values and LIR values from 1 day to ~1 year.

All statistical analyses including Chi-squared tests and linear regressions were performed using program R version 4.1.2 (R Core Team, 2021). We identified 8,976 individual whale sharks at 25 global sites, across 11,858 combined sampling days (Table 1). A small number of individuals (270) were identified at more than one site. These movements are reported elsewhere in the literature (summarised in Supplementary Table 3), but briefly, most of these took place within close geographic localities (e.g. within the Gulf of California or along Africa’s eastern coast, Ramírez-Macías et al., 2012b; Andrzejaczek et al., 2016 respectively). The number of individual whale sharks identified at each site ranged from 35 (Thaa Atoll, Maldives) to 1,313 (Yucatan, Mexico) per site, with a mean of 359.0 ± 300.1 S.D. sharks per site (Table 1). Total effort, in sampling days, varied greatly across sites, with an overall mean of 474.3 ± 610.0 d total sampling days, ranging from 68 d (Qatar) to 2,791 d (Oslob, Philippines). The first identification was from Donsol in Apr 1998, and the most recent from Hawai'i in Sep 2020, with an overall mean across sites between first and last identification of 11.1 ± 6.1 yr [range 4.1 (EAS) – 21.2 yr (Donsol, Philippines)].

On average, 5.0 ± 5.4 individual whale sharks were identified per sampling day (range 1.1 – 22.7) across all sites. A significant male bias was detected in the overall sex ratio of sightings, with 72.7% of individuals across all sites being male (Chi-squared test, χ² = 9.93, p < 0.005). The exceptions were Galapagos (1.2% male), St Helena Island (53% male), and Al Lith, Saudi Arabia (50% male). Most sites were juvenile dominated (assuming maturity at 8-9 m, Norman & Stevens 2007), with an overall mean size of 5.7 ± 1.5 m (range 3.6 – 10.8 m; Table 1).

Model H, that tested for an open population where emigration, re-immigration and mortality of individual animals occurs between sampling occasions, was the best-fit model in 22 of our sites (ΔAIC or ΔQAIC < 2), and within the group of best models (ΔAIC or ΔQAIC < 4) in the remaining 3 sites (Supplementary Table 2). For Peru, the best-fit models with ΔAIC or ΔQAIC < 2 were Models C and E, where emigration and mortality, and emigration and re-immigration, occur respectively between sampling occasions – all which are also tested within Model H. For Thailand and Belize, the best-fit model was Model G which is structurally identical but parameterised differently than Model H. To compare sites using the same model, we therefore continued with model H for all analyses.

All sites followed a similar LIR pattern with a decline from 0 – 100 days following initial identification, and trailing off over time (Figure 3). For some sites, the LIR reached zero between 45 – 75 days (e.g. St Helena Island, Galapagos), indicating the whale sharks completely leave the area following initial identification or are affected by sampling – albeit the LIR increases at yearly intervals suggesting some individuals return seasonally (Figure 3). The LIR across all sites declined between 45.2 (Oslob) – 99.5 (Koh Tao, Thailand) % after ~ 1 year (mean 73.9%; Supplementary Table 4). Notably, the LIR at Mafia Island, Tanzania, Bahia de La Paz, South Ari Atoll and Oslob, displayed a steady, yet not dramatic, decline over time, highlighting longer mean residency periods following initial identification and increased associated probabilities of resighting individuals over time, as well as a lower mortality or permanent emigration (Figure 3). Contrastingly, the LIR at Koh Tao, Hawai'i and Galapagos declined rapidly to near zero, suggesting these animals likely do not reside within these study sites.

Modelled daily abundance ranged from 5.4 (Honduras) to 135.8 (Yucatan) individuals per day, with an overall average of 38.3 ± 36.1 individual whale sharks on any one day (Table 3). Modelled residency within study sites ranged from 1.0 (Galapagos) to 66.3 (Oslob) days, with a mean of 21.6 ± 20.3 days across all sites. Residency outside the study sites ranged from 2.2 (Peru) to 1140.4 days (Koh Tao) with a mean of 110.3 ± 243.1 days. Estimates of mortality or permanent emigration ranged from 0.06 (Yucatan) to 0.61 (Pintuyan, Philippines) δ.yr^-1, with an overall mean of 0.28 ± 0.13 across all sites – representing a mean apparent survival of 0.72 ± 0.13 (Table 3).

We used a total of 42,357 unique identifications from 382 individual whale sharks at Oslob to adjust model estimates of residency given ~eight years of nearly continuous daily survey effort. Mean residency at Oslob was estimated at 124.6 ± 154.8 d.yr^-1 or an equivalent residency index of 0.34 ± 0.42. The ∝_LIR for all sites ranged from 0.01 (Hawai'i, Galapagos, Koh Tao) to 1.97 (Oslob), with an overall mean of 0.58 ± 0.49 (Table 3).

Adjusted residency estimates in d.yr^-1 varied greatly across all sites, ranging from 0.6 (Hawai'i) to 124.6 (Oslob) d.yr^-1, with a mean of 37.0 ± 32.3 d.yr^-1 (Table 3; Figure 4). Similarly, residency indices based on adjusted residency estimates ranged from 0.002 (Hawai'i) to 0.341 (Oslob). Declines in LIR over ~ 1 year were significantly correlated with increased adjusted residency values (coefficient = 178.11, adjusted r² = 0.73, p < 0.001).

In this study, we identified 8,976 individual whale sharks across 25 sites. Given the heavy bias observed towards juvenile males, this number is likely to represent only a small portion of the total population. However, based on our modelled average daily abundance of 38 whale sharks, and if we assumed independence between sites, ca. 950 individual whale sharks would be using these sites on any given day. Extrapolated to our adjusted residency of average 37 days per year, this would indicate that ca. 35,000 individual whale sharks use our 25 sites annually. Our average apparent survival of 0.72 suggests there is either high mortality or high emigration to other sites, not covered by our surveys. This means that there are many more whale sharks in the ocean than what we see, and although we cannot apply CMR models to our data, we can infer that global population size is significantly larger than the ~9,000 individuals identified in this study. Further work to estimate population size for the species (globally or within known subpopulations or spatial management units) is paramount to understand the species recovery potential (Akçakaya et al., 2018).

There is evidence of whale shark aggregations with more than 100 individuals on any given day (e.g. off the Yucatan Peninsula, Mexico, de la Parra Venegas et al., 2011), which is of special interest when dealing with an endangered, rare species and requires dedicated management. Our modelled daily whale shark abundance was highly variable across all sites (range 5 – 136 individuals per day), yet it does reflect the empirical knowledge and data from these sites. For example, our largest estimated abundance per day of any site was for Yucatan which is currently the largest known whale shark aggregation globally. Other sites with modelled n of more than 100 include offshore Qatar in the Arabian Gulf and St. Helena in the South Atlantic. The former site is associated with the sharks feeding on mackerel tuna Euthynnus affinis eggs (Robinson et al., 2013; Robinson et al., 2016), similar to the Yucatan site (de la Parra Venegas et al., 2011). The drivers for the latter are not yet fully understood, but St Helena is unique in hosting a mostly adult 1:1 male to female aggregation where courtship and attempted mating behaviours have been reported (Perry et al., 2020).

Apparent survival at these three sites ranged from 0.80 – 0.93, suggesting these animals display philopatric behaviour and return to the site over time (i.e. low permanent emigration or mortality). Interestingly these sites are known to be seasonal for the species, so seasonal emigration is likely (Acuña-Marrero et al., 2014; Perry et al., 2020). In St Helena, for example, average water temperatures fall outside the whale shark’s thermal range (<21°C) for part of the year, which may be why the species is not present. These aggregations are good candidates for dynamic protection, as their occurrence is seasonally predictable and in great numbers, yet not always within the same defined area (see Maxwell et al., 2015). This relatively new management approach might prove useful across different whale shark sites whether they utilise an area to feed or navigate through.

Our residency index of 0.22 for Mafia Island, Tanzania, falls within previously published values based on passive acoustic telemetry (0.15-0.39; Rohner et al., 2020) suggesting our adjusted residency metric based on sightings-derived presence-only data is valid. Small differences could be explained by the fact that the LIR was modelled for all 201 individuals identified at this site, and that 17% of individuals were only seen once, while the acoustic telemetry study estimated the residency index for 51 tagged sharks. Additionally, there could also be biases on which individual whale sharks were tagged for the passive acoustic study (i.e. more resident sharks are more likely to be encountered and therefore tagged). Similarly, our residency index of 0.08 (corresponding to ~28 days per year) falls within published residency indices for whale shark at Al. Lith, Saudi Arabia (0.05-0.26; Cochran et al., 2019). The relatively low residency estimates from both visual and acoustic data is further supported by satellite telemetry work from the area (Berumen et al., 2014) that shows tagged whale sharks seasonally shifting from coastal (Al Lith) to offshore habitat utilisation (Berumen et al., 2014).

Residency index at St Helena Island was estimated at 0.01-0.24 based on passive acoustic telemetry (Perry et al., 2020), with our adjusted residency index falling within that (0.10). Our adjusted metric therefore works for estimating residency, and can be presented in days per year – an intuitive metric more useful for management. Whale sharks might indeed be present for longer and not be detected by conventional photo-ID methods, and highlights how technical solutions (such as passive acoustic telemetry) can be employed to generate high-quality data (e.g. Cagua et al., 2015). Such methods are often expensive and can be invasive, making non-invasive methods like photo-ID valuable tools, especially in countries where tagging is not allowed (e.g. Maldives, R. Rees, pers. comm.).

The South Ari Marine Protected Area (SAMPA) in the Maldives has year-round sightings of whale sharks (Harvey-Carroll et al., 2021). SAMPA likely plays an important role as a developmental habitat for males given the extended residency of individuals there estimated within this study of mean 102 d.yr^-1, their small average size (5.4 m), and the male bias observed there (88%; Perry et al., 2018). It is likely that whale sharks occurrence at SAMPA is linked to feeding opportunities as observed in Oslob and Mafia Island (e.g. Rohner et al., 2015a), also sites with extended residency as per our results. The apparent survival at these three sites is higher than that observed at other sites (0.78 – 0.85), indicating whale sharks might reside for a considerable amount of time before permanently emigrating elsewhere – also consistent with developmental habitat theories (Heupel et al., 2007). Juvenile sharks at other coastal sites also reside for some time, yet not year-round as the aforementioned sites. For example, our adjusted residency values for whale sharks in Donsol and Pintuyan in the Philippines, Belize, and in Djibouti, still suggest they might spend a considerable amount of time at these sites (i.e. 60 – 70 d.yr^-1). These aggregations tend to peak during the boreal winter, starting in October (Djibouti) and ending in early June (Donsol). Boldrocchi et al. (2020) highlighted the foraging preference of whale sharks for copepods at DJI, similar to that observed in Pintuyan (Sanabria et al., 2019) and DON (R. Dungog, pers. comm., May 2018). Whale sharks in Djibouti also target swimming crab spawn (D. Rowat, pers. obs.). In Pintuyan, seasons can be highly variable (Araujo et al., 2017), yet the occurrence of whale sharks has been linked to zooplankton abundance and prey availability (Sanabria et al., 2019). In Belize, whale sharks aggregate in Apr-May each year to capitalise on reef snapper spawn (Heyman et al., 2001). These events tend to last for ~2 weeks, and whale sharks residing within the general area are also preying on thimble jellyfish (Linuche unguiculata), copepods and baitfish (R. Graham pers. obs.) and plausibly vertically migrating prey (Graham et al., 2006). Although all these sites have differences in seasonal use, they all show philopatry, with individuals returning at annual or inter-annual intervals, with one individual in Donsol having the longest reported philopatric behaviour herein of 21 years similar to that observed at Ningaloo Reef (Norman & Morgan, 2016a). The lower apparent survival of individual whale sharks in Pintuyan could be explained by the highly variable seasons at the site, with some years yielding few sightings, and 34% of individuals have only been sighted once (Table 1). These two facts point to a high permanent emigration rate from the site which would lower their estimated apparent survival over time (i.e. probability of return to the site). Within season, a mean estimated abundance of 77 individual sharks on any given day at DJI makes it a significant aggregation globally – contrasting with ~16.5 in Pintuyan, 38 in Donsol or ~7 in Belize. Donsol is the longest-running whale shark ecotourism endeavour in Asia (late 1990s) and hosts a considerable number of adults, juveniles (McCoy et al., 2018), and is globally important for the species’ reproductive ecology (Miranda et al., 2020). Basing conservation action on abundance alone might not make be the most effective for the species recovery – rather targeting areas of ecological importance such as nursery, mating and feeding grounds might prove more effective (Pierce et al., 2021).

Contrastingly to the sites above, whale sharks feeding predominantly on small fish in Honda Bay, Philippines, appear to reside for shorter time periods (mean ~8 d). They do however show periodicity at the site over time (Araujo et al., 2019b), suggesting that although feeding opportunities might be short during the low primary productivity months in the region (May-Nov, Cabrera et al., 2011), it is worth returning and capitalising on this feeding opportunity. Similarly, at Honduras, whale sharks feed on baitfish (Fox et al., 2013) and have adjusted residency of ~9 days. Whale sharks identified in Honduras move to and from Belize and Yucatan, indicating that, although some feeding opportunities exist off Utila, they travel to other more productive grounds dominated by other prey sources (e.g. snapper spawn, Graham and Roberts, 2006; sergestids and copepods, Motta et al., 2010; tunny spawn, de la Parra Venegas et al., 2011). In contrast, whale sharks in Madagascar are normally encountered in association with baitfish similar to that observed in Honda Bay (Diamant et al., 2018), but their residency time is around three times that observed in Honda Bay or Belize. This could be because Nosy Be is a big embayment whereas at the aforementioned sites whale sharks are likely feeding on frontal systems passing through (Ryan et al., 2017). Sharks at Nosy Be are also reportedly associated with other filter-feeding megafauna that target krill species, so a combination of both prey sources might provide the whale sharks with longer residency times (Diamant et al., 2018). The Authors showed that whale sharks spent some time at Nosy Be following satellite tag deployment, before moving west or southwest.

An unexpected result among the low-residency aggregations is Ningaloo Reef, where we estimated ~9 days per year – contrasting with previous studies at the site (33 d, Holmberg et al., 2009; ~44 d, Lester et al., 2020). There is evidence that whale sharks use a broad area along Western Australia beyond Ningaloo Reef as shown through satellite telemetry (Norman et al., 2016b), and that whale sharks are found on the Ningaloo Reef year-round (Norman et al., 2017b). Interestingly, apparent survival (0.82) was similar to that modelled by Lester et al. (2020) (~0.85) and Holmberg et al. (2009) (0.48 – 0.89). Both these studies included covariates to their CMR models: scarring by the former as described by Speed et al. (2008) (e.g. bites, abrasions, lacerations, etc.), and size by the latter, whereas we used a different approach altogether (LIR). Further modelling with the same dataset for Ningaloo Reef exploring both LIR and CMR modelling approaches could shed light to these discrepancies. For example, Harvey-Carroll et al. (2021) used scarring as a covariate for whale sharks in SAMPA and obtained similar apparent survival results through both the LIR and CMR approaches.

Some whale shark hotspots appear to be transitional waypoints rather than aggregations. Galapagos, Hawai'i and Koh Tao stand out with adjusted residency estimated values of <1 day. This is not surprising as whale sharks seen at these sites are normally cruising through the area and not commonly resighted within a short time-scale (<6 months; Acuña-Marrera et al., 2014; S Marcoux, pers. obs.; Magson et al. 2022). Satellite telemetry data from Galapagos and Tubbataha in the Philippines support this understanding (Hearn et al., 2016; Araujo et al., 2018). Little is known about the movement of whale sharks in the western South China Sea, including those encountered in KOH and neighbouring islands. To date, no major feeding ground has been identified, though reports of whale sharks off Cambodia’s offshore islands, and encounters off eastern Peninsular Malaysia exist (G Araujo unpub. data). Further investigation in this general area is paramount given the proximity and likely connectivity of whale sharks to China, where the largest single fishery for the species operated until recently (Li et al., 2012). Similarly, little is known about whale sharks in the Central Pacific, and further telemetry work in Hawai'i can help identify their feeding grounds, or connectivity across sites. It is likely that some of these whale sharks encountered along seamounts and offshore islands are foraging on frontal areas offshore (Ryan et al., 2017).

There are a number of limitations to using photo-ID data to understand population demographics. Even assuming no changes in the spot pattern of individuals, animals might still be misidentified and wrongly assigned as a new or other individual (Pierce et al., 2018). This is particularly difficult for growing photo-ID catalogues (hundreds of individuals) where visual photo-ID matching becomes infeasible and reliance on automated systems (i.e. I³S or WWS) is necessary. Another common caveat is that whale sharks might be present at the site yet not be encountered and identified during photo-ID survey efforts. Although some models can account for this (e.g. Torres et al., 2008), it is a limitation that other methods (e.g. passive acoustic monitoring, Janik et al., 2013) can account for on vocal marine species for example. In whale sharks, passive acoustic telemetry at one site showed that only few individuals (<5%) were missed by photo-ID during the survey season (Rohner et al., 2020). Detection probability is likely to vary among sites, however, depending on survey area, bathymetry, frequency, number of individuals and their behaviour (e.g. Cagua et al., 2015).

A limitation with the LIR movement models is the inability to deal with heterogeneity within the data, where population data is mostly split between resident (i.e. high resighting rate) and transient (i.e. seen once) individuals (H Whitehead, pers. comm., May 2020). Both behaviours are expected given that ideal whale shark feeding opportunities are seasonal (e.g. Heyman et al., 2001). Estimates of ‘residency time out’ are also imprecise within the LIR approach given the presence-only data used as input, and hence our approach herein to try and further understand residency patterns. In the present study, data for East Kalimantan and Peru followed these patterns and there were several best-fit models (Supplementary Table 2). The data heterogeneity issue has been previously highlighted in other taxa (e.g. marine turtles, Araujo et al., 2019a). These models are, however, indicative, not absolute, and their suitability is determined to a large degree by the ecological question being investigated. However, it is worth noting that given the difficulty of applying standardised methods across different locations (e.g. 25 described here), using methods like LIR to answer important ecological questions and future work to deal with heterogeneity should be encouraged. The ability of the LIR approach to use presence-only data (i.e. no absence or zero data) facilitates its applicability across different sites that might have different survey methods, like that described herein.

It is worth noting that our adjusted residency times in days per year are based on an aggregation of whale sharks whose residency behaviour has been noted to be different due to the provisioning activities (Araujo et al., 2014; Thomson et al., 2017; results herein). Given that the LIR uses the identification data itself as the unit of effort, and that the same methods were applied to all sites, the results obtained should not be biased. Whale sharks in Oslob display higher residency rates than at other sites, however, we used the ∝_LIR across all sites and the empirical data from Oslob to interpret the LIR residency values using presence-only data. Our extrapolation of this approach to all sites provided intuitive residency estimates relevant for management and allowed for validation of the results against both current and future acoustic data from each site.