Population dynamics of oceanic manta rays (Mobula birostris) in the Raja Ampat Archipelago, West Papua, Indonesia, and the impacts of the El Niño–Southern Oscillation on their movement ecology

Our aim was to collect sightings data on oceanic manta rays (Mobula birostris) within the Raja Ampat Archipelago to better understand their population dynamics within the region. These data were compared with environmental variables to seek correlates that may explain any variations in observed sightings frequency. Combined, it is hoped this knowledge will be used to aid effective management of this species in the region.


Manta rays (previously genus Manta but recently subsumed into
Mobula) (Marshall, Compagno, & Bennett, 2009;White et al., 2018) are large planktivores found circumglobally in tropical and subtropical waters. In addition to targeted fisheries, manta rays are incidentally captured in a wide variety of fisheries and gear types (Croll et al., 2016). Their low reproductive rates and conservative life history characteristics make them highly susceptible to population declines, and fisheries are likely to be unsustainable even at low catch rates (Dulvy et al., 2014). Although data on population trends are limited, several studies have reported declines in manta and devil ray sightings and fisheries landings, suggesting that fisheries are negatively impacting populations in many regions (Croll et al., 2016;Lewis et al., 2015;Ward-Paige, Davis, & Worm, 2013). Population declines have led to the implementation of several international, national and local management strategies for mobulid rays in the last decade (Croll et al., 2016;Lawson et al., 2017).
Oceanic manta ray (Mobula birostris) populations appear to be regionally distinct, and there is evidence for fisheries-induced declines in several isolated subpopulations (Lewis et al., 2015;Moazzam, 2018;Stewart, Beale, et al., 2016). In the closely related reef manta ray (Mobula alfredi), local declines have also been recorded in response to fisheries (Croll et al., 2016;Lawson et al., 2017;Marshall et al., 2011). In Indonesia, high-effort gill net fisheries in the Lembeh Strait captured 1,424 manta rays (Mobula spp.) in a 10-month period (Cochrane, 1997), leading to an apparent local extirpation (D. Djalal & A. Doali, personal communication). Indonesia was reported to be the third largest exporter worldwide (Heinrichs, O'Malley, Medd, & Hilton, 2011) for Mobula gill rakers up until the nationwide ban of fishing manta rays in January 2014 (KEPMEN-KP 4/2014), with some communities continuing to illegally harvest manta gills to date (E. Setyawan, personal communication). The majority of these reports are based on fisheries-dependent data, sightings per unit effort trends and anecdotal reports; there is a need for more quantitative data on abundance and trends in manta populations.

K E Y W O R D S
citizen science, climate change, conservation, demographics, marine megafauna, mobulid rays to affect population dynamics (Lo-Yat et al., 2011;Pearce & Phillips, 1988). Mass mortalities of coral reefs have occurred due to long periods of high-temperature exposure with far-reaching impacts on reef-reliant species viability (Glynn, 1999;Graham et al., 2007;Lo-Yat et al., 2011;Quiñones et al., 2010;Wallace et al., 2006). As ENSO events appear to be increasing in frequency and strength (National Oceanic & Atmospheric Administration, 2017), understanding their impacts on manta ray species will be important for developing effective management strategies.
Within the Raja Ampat Archipelago and along the coastline of West Papua, there are a number of cleaning stations that are visited seasonally by both reef and oceanic manta rays (Setyawan et al., 2018). This area is renowned for both its high marine biodiversity and its strong ocean currents (Mangubhai et al., 2012). Known as the Indonesian Throughflow, warm western Pacific Ocean water pushes down into the Ceram and Halmahera Seas and through Indonesia to the Indian Ocean (Gordon, 2005;Gordon, Susanto, & Vranes, 2003). This generally oligotrophic water is regularly mixed vertically by south-east monsoon ocean upwelling throughout the region, leading to local productivity hotspots during Autumn months (Tang, Kester, Ni, Kawamura, & Hong, 2002;Wyrtki, 1962). The region is a heavily trafficked tourism destination, with a number of live-aboard vessels and resorts providing access to cleaning stations and other coral reef habitats where manta rays are frequently encountered.
In this study, we describe the demographic characteristics of the oceanic manta ray population, sighting trends across time, the relationship of environmental covariates to sightings, and perform a mark-recapture analysis. In particular, we emphasize the apparent importance of the El Niño-Southern Oscillation in driving sightings trends and possible distributional changes in the regional oceanic manta ray population.

| Sighting records
The majority (75.4 %) of photo-IDs were collected by tourists and submitted to the project's photo-ID database. All data were collected from cleaning stations. We used the unique ventral spot patterns found on manta rays to identify and catalogue individual oceanic manta rays, a technique applied to both species of manta rays (Marshall & Pierce, 2012). Starting in 2011, data were collected annually by the authors during the tourism season, running from September to June. The strong monsoon winds during the months of July and August prevented trips from running and therefore precluded the collection of data (no data were submitted from other sources during this period either). Photo-ID images collected prior to August 2011 were added to the catalogue only when date and location were verifiable by the observer. These images were mostly from April 2011 but include sightings back to May 2004.
All submitted photographs and videos were dated, and all usable images were subsequently compared against the existing ID catalogue to determine whether the individual was a new capture or recapture. A new sighting record was then added to the database to reflect this. A sighting was defined as a capture record on a specific site on a specific day; a recapture was defined as a sighting of a previously identified individual on a different day, or on a different site at least one-hour post-initial sighting. Similar to Stevick, Palsbøll, Smith, Bravington, and Hammond (2001), very low-quality images were discarded to reduce identification errors. Secondary verification of IDs was performed manually by an outside expert and a third time using the Manta Trust's "IDTheManta" matching software by trained research staff (The Manta Trust, 2018).

| Physical characteristics
Additional information recorded for each individual included colour morph, sex, estimated size, sexual maturity (identified via pregnancy or presence of mating scars for females and calcified claspers for males), and external injuries as per Marshall and Bennett (2010) and Stevens (2016).

| Sighting seasonality and survey effort
We pooled sightings data by month and split them into equal periods of six months. We classified August through January as the Autumn season and February through July as the Spring season. This split maintains equal sightings opportunities for both seasons while separating the two major influencing weather patterns in the area: the very windy south-east monsoon in Autumn and the rainier north-west monsoon in Spring. We recorded survey effort from Autumn 2011 through Spring 2016. The effort value is reported as the number of divers who dived on the primary oceanic manta ray cleaning site Magic Mountain per month in the Misool region, southern Raja Ampat. Effort data prior to Autumn 2011, or from other regions, were not available and so not included in the analysis. While this is an imperfect metric of effort, 77% of sightings were collected at Magic Mountain, and precise records of diver visitation to the site were recorded by Misool Eco Resort, making it the best available proxy for survey effort.  (Bamston, Chelliah, & Goldenberg, 1997). We averaged the MEI across the seasonal sighting periods, from February to July for the Spring season and August to January for the Autumn season for use in our analysis. To examine a possible mechanism for increased sightings in the region, we also explored the relationship between sightings and sea surface temperature (SST) and chlorophyll-a (Chl-a). We used the Xtractomatic R package (Mendelssohn, Bessey, & Foley, 2018) to download satellite-derived seasonal averages of SST and Chl-a for the study region. We used Bayesian linear models to explore relationships between the number of manta sightings per season, and MEI, SST and Chl-a covariates.

| Climate conditions
We considered all seasons together, only Autumn seasons, and only Spring seasons independently to determine if covariate effects were most relevant in a certain season. We custom coded the linear models in R using the software Just Another Gibbs Sampler (JAGS) (Plummer, 2004).

| Mark-recapture modelling
To generate estimates of abundance, survival, sighting probability and recruitment to the population, we used a POPAN mark-recapture model in the R package RMark (Cooch & White, 2011;Laake, 2013). The model estimates four main parameters: survival probability (Phi), entry probability (pent), sighting probability (p) and superpopulation size (N). We considered time-varying and fixed values for each parameter and possible covariate relationships of ENSO and sex. We tested models with every possible combination of time-varying/fixed and covariate relationships for every parameter, which resulted in a suite of 125 candidate models (Appendix S1: Table S1.1). We used Akaike's information criterion for small sample sizes (AIC c ) to compare candidate model fits and select the best-fit model (Burnham and Anderson, 2004), and we report the top 19 models in the results (where AICc weight ≥ 0.002). Given the paucity of photo-ID contributions prior to 2010, we only included sightings from the Spring 2010 season onwards in the mark-recapture analysis. Analysis code is provided in Data S1.

| Sighting records
We identified a total of 588 individual oceanic manta rays within the 856 total photo-ID sightings dated between 17/05/2004 and 13/06/2016 (Table 1). Only 12 sightings were submitted from prior to Spring 2010. Records were submitted from 18 dive sites within Raja Ampat and West Papua.
Due to the nature of collecting data from volunteer sources, time stamps were often inaccurate or missing from submitted data. Only 36% (n = 306) of sightings were recorded with accompanying time of day, 213 of which were from Misool Eco Resort operations. Sighting times were heavily influenced by the trip schedule of the resort.
Sightings were most frequent between 12:00 and 13:00, accounting for 26% of records, while the two-hour period between 11:30 and 13:30 accounted for 43.4% of all sighting records.

| Sighting regions
The Misool region accounted for 84.8% of all oceanic manta ray sightings (n = 726), forming the major part of the database. The

| Population data
It was possible to determine the sex of 565 (96.1%) of the identified oceanic manta rays ( pectoral fin similar to reef manta rays (Marshall & Bennett, 2010).
F I G U R E 2 Sighting frequencies of oceanic manta rays (Mobula birostris) (Marshall et al., 2009;White et al., 2018)  A total of 110 females had mating scars, and almost all mating scars were ventral (only two were recorded solely from the dorsal side).

| Monthly sightings variation
Mean monthly total sightings were 13.2 (±2.5 SE), with a distinct peak during February (mean 26.0 ± 16.0 SE) and again during April-May (mean 24.0 ± 9.9 SE). Overall, February had the highest mean total sightings recorded, while April had the highest mean recapture rates of 11.4% (±6.53 SE). Peak Autumn months' recapture rates oc-

| Climate conditions and sightings
Between

| Mark-recapture modelling
The top four POPAN mark-recapture models were within two Akaike information criterion (AICc) points of one another and had similar AICc weights (range 0.143 to 0.316) ( Table 2) and are therefore all reasonable fits to the sightings data. The parameter estimates of all four top models were similar, and we report the best-fit model results here and include the next three model results in the Supporting Information (Appendix S1: Tables S1.2, S1.3 & S1.4). The The most parsimonious model (  shown that females are significantly more likely to be observed at cleaning stations than males (Couturier et al., 2014;Stevens, 2016).

| D ISCUSS I ON
Sex-specific site selection factors, including proximity to food sources, birthing grounds or reproduction opportunities (Marshall & Bennett, 2010), may explain the female bias observed as our sampling coverage is entirely from cleaning stations. The cause of the sudden increase in male sex ratios in 2013 could be linked to an increase in reproductive activity drawing males to cleaning stations, which may act as a lek (Stevens, 2016), but could also be explained by relatively low sample sizes. This period follows a La Niña event and likely an increase in local zooplankton densities correlated with SST drops (Beaugrand & Reid, 2003), both of which have been linked with increased reproductive activity of reef manta rays (Stevens, 2016). Sexual maturity rates may be higher than reported, as we often lacked photographs of claspers or female wingtips to confirm maturity status, and eight of 29 visibly pregnant manta rays did not have mating scars, suggesting that our measures of sexual maturity underestimate actual levels.
Significant linear relationships were found between increased oceanic manta ray sighting frequency and an increase in the MEI and a reduction in SST. While ENSO is perhaps best known for its impacts in the Eastern Pacific, the El Niño-La Nina cycle also has a major influence on prevailing winds and precipitation in Indonesia and elsewhere in the western Pacific (Tang et al., 2002;Zhang, Sumi, & Kimoto, 1996). We found an exponential relationship between the number of sightings of oceanic manta rays and the MEI. The fact that almost all of the newly identified individuals in 2015 and 2016 (the peak ENSO period) were apparently mature adults suggests that this was not a recruitment event in the traditional sense of newborns first recruiting to the adult population. Instead, it appears to be a pool of individuals that were previously not observed and who entered the study region for the first time in our study during the strong El Niño event.
Previous study by Stewart, Beale, et al. (2016)  The 2015-2016 increase in manta sightings highlights an inherent limitation in applying mark-recapture methodologies to migratory species with very large ranges. Mark-recapture models assume that all individuals in the population have an equal probability of capture (Laake, 2013). However, in cases where survey effort covers only a portion of the population's entire range, this assumption is not met.
The possibility of a range expansion of the regional oceanic manta ray population during strong El Niño events is supported by our covariate analyses, which also reveal a potential driver for this expansion. Over the study period, Autumn ocean and atmospheric conditions are most heavily impacted by ENSO variability, whereas Spring conditions remain relatively stable (Figure 3). During El Niño TA B L E 2 Model selection weighting (Akaike information criterion-AICc) for the POPAN mark-recapture models (n = 19)  (Chowdhury, Chu, & Schroeder, 2007;Gordon, 2005;Rejeki, Munasik, & Kunarso, 2017). It is likely that this enables south monsoon winds to push cooler waters from the Ceram Sea into the study area ( Figure 3) and increases vertical mixing of the water column, which together cause a shallowing of the thermocline (Tang et al., 2002;Wyrtki, 1962). SST change has been shown to be associated with changes in the composition and distribution of plankton (Edwards & Richardson, 2004;Rutherford, D'Hondt, & Prell, 1999), with knock-on impacts on the distributions of various large marine planktivores (Cotton, Sims, Fanshawe, & Chadwick, 2005;Simmonds & Isaac, 2007;Wilson, Polovina, Stewart, & Meekan, 2006). While sightings counts were higher in Spring, it is likely the Autumn conditions drive these peaks by increasing vertical mixing and perhaps leading to increased food availability that promotes extended periods of improved foraging and consequently reproduction opportunities (Clark, 2010;Stevens, 2016).
Sightings began to increase in Spring 2015 (Figure 2), but it was not until the cool water mass pushed beyond Magic Mountain (approximately 2.25ºS latitude) that sightings increased exponentially.
Oceanic manta rays, as filter feeders, are likely tied closely to local primary and secondary productivity (Stevens, 2016;Stewart, Hoyos-Padilla, Kumli, & Rubin, 2016). However, surface primary productivity measured by satellite-derived Chl-a did not vary substantially across the study period despite large swings in SST and MEI (Appendix S1: Figure S1.1), and there was not a clear relationship between sightings and Chl-a values. A number of studies suggest that oceanic manta rays target mesopelagic and vertically migrating prey (Burgess et al., 2016;Rohner et al., 2017;Stewart, Hoyos-Padilla, et al., 2016;. Consequently, surface Chl-a values may not represent prey availability and distribution in the region, which may explain the lack of a relationship between Chl-a and sightings. Another possible explanation is that oceanic manta rays may have a thermal envelope with an upper bound of approximately 29ºC. Oceanic manta rays are seen less frequently than reef manta rays in Raja Ampat, and the increase in sightings during anomalously cold conditions could be due to typical water temperatures in the region being above oceanic manta rays' thermal optimum. However, manta sightings remained high throughout Spring 2016 even as water temperatures returned to their typical Spring highs. This suggests that the range expansion may have been driven by food availability or a different, unknown factor, rather than temperature preference.
The mark-recapture analysis identified MEI as the most important factor influencing the probability of recruitment to the population.
Additionally, years with high MEI values also had the highest sighting probabilities for manta rays in our mark-recapture models. The apparent immigration into the study area during the El Nino would mean there was not equal probability of capture over time, and the total abundance estimate from our mark-recapture analysis should be interpreted as a minimum population size, as a large number of individuals may be inaccessible to sampling within Raja Ampat.
The superpopulation estimate of oceanic manta rays in Raja Ampat appears to be of similar size to populations in the Eastern Pacific: a photo-ID study of oceanic manta rays in Pacific Mexico has identified 715 individuals since 1978 (Rubin, 2017). The fact that Raja Ampat and the Ceram Sea, located in a comparatively oligotrophic region (Longhurst, Sathyendranath, Platt, & Caverhill, 1995), can support a population similar in magnitude to the rich waters of the Eastern tropical Pacific suggests the west Papua area may contain abundant food resources that could be important to other species with tropical distributions (Lehodey et al., 1997). Modelled estimates of survival are lower than expected for these long-lived species and lower than reported estimates from reef manta ray populations (Couturier et al., 2014;Deakos, Baker, & Bejder, 2011).
It is likely that this is a result of our incomplete sampling of the entire population distribution and quite likely reflects emigration as opposed to elevated mortality. This theory is supported by satellite tag tracks (Stewart, Beale, et al., 2016) (Blaber, 2006;Dharmadi & Sumadhihargs, 2008;Dharmadi & White, 2003).
Survey effort was concentrated in the Misool region; this may explain why little regional crossover was observed during the study.
We expect as more photo-IDs are submitted by tourists across Raja Ampat, there will be an increase in recaptures between regions.

| Conclusion and implications
Oceanic manta ray distributions appear to be impacted by ENSO-re-

ACK N OWLED G EM ENTS
This study was made possible by the government and communities of Raja Ampat, Indonesia, and the Indonesian Ministry of Marine Affairs and Fisheries. We thank all those involved in the collection and contribution of ID photographs, the support and contributions of Conservation International Indonesia and the Manta Trust. We also thank the in-kind support of Misool Eco Resort and Misool Baseftin during the study period which made collecting these data possible.

DATA ACCE SS I B I LIT Y
The R code and further data supporting the results are provided in the Data S1.