Contrasting Habitat Use and Population Dynamics of Reef Manta Rays Within the Nusa Penida Marine Protected Area, Indonesia

: Manta rays (Mobula spp.) are highly valued in nature-based tourism globally. In Indonesia, although manta rays are protected, critical information is lacking on their habitat use, population dynamics and movements. We investigate the population structure and residency patterns of reef manta rays (Mobula alfredi) in the Nusa Penida Marine Protected Area (MPA). From photo-identification data logged by citizen scientists and trained observers (mantamatcher.org), we identified 624 reef manta rays from 5,913 sightings (January 2012–April 2018) based on their unique ventral coloration patterns. Year-round records were collected from two shallow (<20 m) reefs – Manta Bay (MB; n = 3,029 sightings) and Manta Point (MP; n = 3,058) – that are used frequently by tourism operators. Maximum likelihood techniques and a Markov movement analysis were used to model residency patterns and movement between these sites within the MPA. Manta rays at MB were predominantly male (64%, n = 261 individuals), with immature males (14%, n = 59) being sighted most frequently (39%, n = 1,170). In contrast, few immature individuals were sighted at MP (6%, n = 28), and they were sighted on few occasions (2%, n = 45), while mature female manta rays comprised 26% (n = 127) of the MP community and were the most frequently sighted (48%, n = 1,413). Lagged identification rates indicated high site fidelity at each location. However, 44% (n = 278) of individuals moved between the two sites and cumulative discovery curves showed a continued recruitment of individuals over disturbance to this threatened species by tourism, regulations for the number of boats and interactions, especially during key reproductive times should be considered. Further, strict ﬁshing regulation in the area is recommended as ﬁshing gear entanglement was identiﬁed as a threat to this population.

Manta rays (Mobula spp.) are highly valued in nature-based tourism globally. In Indonesia, although manta rays are protected, critical information is lacking on their habitat use, population dynamics and movements. We investigate the population structure and residency patterns of reef manta rays (Mobula alfredi) in the Nusa Penida Marine Protected Area (MPA). From photo-identification data logged by citizen scientists and trained observers (mantamatcher.org), we identified 624 reef manta rays from 5,913 sightings (January 2012-April 2018) based on their unique ventral coloration patterns. Year-round records were collected from two shallow (<20 m) reefs -Manta Bay (MB; n = 3,029 sightings) and Manta Point (MP; n = 3,058) -that are used frequently by tourism operators. Maximum likelihood techniques and a Markov movement analysis were used to model residency patterns and movement between these sites within the MPA. Manta rays at MB were predominantly male (64%, n = 261 individuals), with immature males (14%, n = 59) being sighted most frequently (39%, n = 1,170). In contrast, few immature individuals were sighted at MP (6%, n = 28), and they were sighted on few occasions (2%, n = 45), while mature female manta rays comprised 26% (n = 127) of the MP community and were the most frequently sighted (48%, n = 1,413). Lagged identification rates indicated high site fidelity at each location. However, 44% (n = 278) of individuals moved between the two sites and cumulative discovery curves showed a continued recruitment of individuals over the 6 years of the study. In addition, the behaviors displayed by the manta rays differed markedly between the two sites: MB appears to be a foraging ground, especially for juveniles, and potentially a nursery, while MP is used mainly for cleaning and courtship, indicating a social and reproductive site. Reproductive behavior coincided with the peak annual sightings in May. To prevent

INTRODUCTION
Many marine megafauna species, such as sea turtles (Bowen and Karl, 2007), cetaceans (Baird et al., 2008), large teleost fish (Robichaud and Rose, 2001;Rooker et al., 2008), and elasmobranchs (Hueter et al., 2005), range widely, but within these large areas show high fidelity to defined areas for reproduction and feeding (Chapman et al., 2015). Further, a specific demographic within a population might depend on limited habitats during different stages in their development (e.g., coastal nursery areas versus open ocean foraging grounds) (Hueter et al., 2005;Chapman et al., 2015). While Marine Protected Areas (MPAs) are often created to protect critical habitats of threatened species, these areas are sometimes designated without in-depth knowledge of how the habitats are used by the target species. An understanding of habitat use is needed to achieve conservation goals, which can also require that tourism and other threats are managed to avoid disruption to the biology and behavior of the threatened species in question (e.g., Hueter et al., 2005).
Manta rays are large, pelagic, filter-feeding mobulid rays (Mobula spp.) found in equatorial and tropical waters (e.g., Couturier et al., 2012), that are threatened with extinction. Several known aggregation areas exist within the Coral Triangle region, particularly within Indonesia (Marshall et al., 2009;Couturier et al., 2012;Germanov and Marshall, 2014). Manta rays are long-lived, slow to mature, and exhibit low fecundity (Marshall and Bennett, 2010;Stevens, 2016;Rambahiniarison et al., 2018). Owing in part to these conservative life history traits and to overfishing across their range (Dulvy et al., 2014), both Mobula alfredi (Krefft, 1868) and the larger M. birostris (Walbaum, 1792) are currently listed as vulnerable to extinction on the IUCN Red List of Endangered Species (Marshall et al., 2011). Concerns over the population status of manta rays led to their listing on the Convention on the International Trade of Endangered Species (CITES) and the Convention on the Conservation of Migratory Species (CMS), which aim to curb their consumptive over-exploitation by fisheries (Ward-Paige et al., 2013). Until recently, Indonesia was ranked third for the highest annual mobulid catches, including manta rays (Heinrichs et al., 2011). In 2014, Indonesia protected both manta ray species within their entire exclusive economic zone (an area of over 6 million km 2 ) through the Ministerial Decree of Marine Affairs andFisheries No. 4/2014 (Ministry of Marine Affairs andFisheries, 2014), in the hopes of slowing population declines. Currently, limited information is available for Indonesian manta ray populations to assist in tracking the current status of these slow-growing species.
Manta ray protection within Indonesia came largely in response to the growth of the manta ray tourism industry. Indonesia is ranked second in the world for manta ray tourism, with an estimated value of USD $15 million per year (O'Malley et al., 2013). Prior to their national protection, manta rays were first protected within three sanctuaries east of Bali: The Nusa Penida MPA (200 km 2 ), West Manggarai including Komodo National Park (7,000 km 2 ), and Raja Ampat (11,655 km 2 ). Together, these sites encompass the bulk of the manta ray tourism industry in Indonesia (O'Malley et al., 2013). In addition to CITES and CMS regulations in Indonesia, evidence-based conservation strategies are needed to protect the remaining manta ray populations from the potential threats associated with unregulated tourism (reviewed by Tyne et al., 2014;Trave et al., 2017;Stewart et al., 2018a). Globally, as awareness grows of the chronic stresses to manta rays associated with an increase of tourism, measures aimed at curbing these impacts are underway or planned in some popular manta ray dive sites (e.g., Raja Ampat, Indonesia, Kasmidi and Gunadharma, 2017;and West Hawaii, United States, Division of Boating and Ocean Recreation, 2016). Other hazards to manta rays (reviewed by Stewart et al., 2018a) include entanglement in fishing lines (Croll et al., 2016) and boat mooring lines.
The Nusa Penida MPA, located 18 km south-east of Bali Island (Figure 1), encompasses the habitats for M. alfredi (hereafter manta rays). The waters of the Nusa Penida MPA, established in 2010 (District Fisheries and Marine Agency, 2010), are very productive (Ayers et al., 2014), with complex oceanographic conditions owing to their proximity to a major channel of the Indonesian Through Flow (ITF) current (Murray and Arief, 1988;Nyegaard, 2018). The islands of Ceningan, Lembongan, and Penida, located within the MPA, are home to 48,000 people and have been a tourist destination since the 1990s. In 2015, the MPA received an estimated 200,000 tourists per year, with an annual estimated increase of 7.7% (Barr et al., 2017). Manta ray watching is a major drawcard for tourist visits to the MPA, with annual estimates ranging from 4,200 (Barr et al., 2017) to 10,440 people (Aquatic Alliance, personal communication) or ∼ 63,500 manta dives and ∼ 15,000 manta swims (O'Malley et al., 2013). As manta rays are present in the MPA year-round (Germanov and Marshall, 2014), this manta ray population is subjected to chronic pressure from tourism. Currently, the number of boats allowed to enter manta ray habitats is unregulated, as is general boat traffic in the area (Barr et al., 2017). Further, codes of conduct for manta ray interactions are voluntary and are not strictly enforced, leading to concerns over low compliance (Barr et al., 2017). Moreover, basic information about the local manta ray population structure, habitat use, and movement patterns is lacking in the region, yet is required to refine spatial management within the MPA and promote ecologically sustainable tourism (Trave et al., 2017).
Here, we used a citizen science approach accessing over six years of publically contributed photographs and individualbased sighting records of manta rays within the Nusa Penida MPA, which were logged on the online database 'Manta Matcher' Wildbook 1 , to investigate population dynamics and habitat use. The aim of this work is to inform evidence-based management plans to ensure the long-term sustainability of manta ray tourism within the Nusa Penida MPA for both operators and manta rays. This was done by comparing the habitat use and population dynamics at two, high-intensity tourism sites, located only 12 km apart, to better understand the fine-scale habitat use of the manta population in the region. Given the proximity of the two sites, we predicted that they would share common aspects of population demography and behavior. With this study, we provide in-depth information on: (1) population structure, (2) injury rates, (3) fine scale habitat use, (4) local movement, and (5) behavior, and provide specific recommendations for effective manta ray management within the Nusa Penida MPA. 1 www.mantamatcher.org

Description of Study Sites
The Nusa Penida MPA is located in a unique oceanographic region, with a shallow (200 m) continental shelf adjacent to deep water basins of up to 4,200 m depth. The three islands located within the Nusa Penida MPA (Figure 1, created using QGIS Development Team, 2016, v 2.18) are in the direct path of the ITF current due to their proximity to the Lombok Strait (Murray and Arief, 1988). This major current separates two significant marine and terrestrial biogeographical realms (i.e., the 'Wallace Line, ' Barber et al., 2000). The current flows largely from north to south between the larger islands of Bali and Lombok, with its strength varying seasonally (i.e., the north-west vs. south-east windward seasons, Mayer and Damm, 2012). The ITF brings warmer water (Tillinger, 2011) and nutrients (Ayers et al., 2014) from the Pacific to the Indian Ocean. Deep water upwelling (Ningsih et al., 2013) south of Bali further contributes nutrient rich waters to the region, creating productive foraging grounds within the Nusa Penida MPA for sunfish (Thys et al., 2016;Tito and Susilo, 2017;Nyegaard, 2018), and potentially other megafauna like manta rays.
Although 11 sites are dived frequently (number of days > 20 per year) in the Nusa Penida MPA, the majority of manta ray sightings (∼ 93%), based on year-round dive logs for 2016 (n = 809) and 2017 (n = 753), by one dive operator (Supplementary Figure 1), were at two sites, Manta Bay (MB) and Manta Point (MP) (Figure 1). Thus, these two sites were selected for detailed investigation. Both are shallow bays (7-25 m depth) fringed by steep rocky cliffs and deep water, separated by a straight-line distance of approximately 12.4 km. A range of habitats used by manta rays, including cleaning stations, foraging aggregation areas, and reproductive grounds are found at the two sites (Germanov and Marshall, 2014).

Data Collection
The online database 'Manta Matcher' Wildbook (see text footnote 1) ) was used to collate data (date, time, location, and identifying photographs of manta rays) from the public and researchers (Germanov and Marshall, 2014) (Supplementary Figure 2). Citizen science was an important part of this process and has been demonstrated to provide basic information when data are limited (e.g., Jaine et al., 2012;Rohner et al., 2013;Couturier et al., 2014;Germanov and Marshall, 2014;Robinson et al., 2016;Norman et al., 2017;Pierce et al., 2018) and engage local communities in marine stewardship (Miller-Rushing et al., 2012). Approximately 40 observers were trained to take photos of manta ray ventral spot patterns, record sex and maturity data, estimate size, and to note relevant behaviors (discussed below), as well as the maximum number of boats present at the sites. Tourists were regularly informed about 'Manta Matcher' through dive operator briefings and weekly educational presentations starting in 2012, resulting submissions from 276 different emails. A more traditional sampling approach with structured surveys (e.g., Deakos et al., 2011;Couturier et al., 2014) is required for a traditional mark-recapture approach to population analyses to estimate abundance and survival.
Photographs of the ventral surfaces of individual manta rays (Deakos et al., 2011;Kitchen-Wheeler et al., 2011;Marshall et al., 2011;Couturier et al., 2014) were logged year-round by visiting tourists to the Nusa Penida MPA and members of the local dive community (i.e., residents and dive professionals). The earliest available photographs and sighting records date back to June 2004, while year-round survey effort by trained observers started in 2012 (Supplementary Table 1). Only records collected between January 2012 and April 2018, when greater than 400 sightings were logged in any 1 year, were used for the analyses. Further, data from 2012, 2015, and 2018 were excluded from seasonal analyses as year-round survey effort was lower, or not available, during these years. For the purposes of this study, we assign a survey day as the unit of survey effort, as the number of dives per day and hours of diving were not recorded on a regular basis.
All sighting records were validated manually for accuracy, completeness, and adequate photo quality by trained observers. Sightings lacking information on location or date and those with an indistinct ventral image of the manta ray were excluded. Photographs were manually matched to an ID catalog with the assistance of built-in (Town et al., 2013;see Germanov and Marshall, 2014) and/or external software ('MantaUtil, ' Winstanley, 2016) by trained observers and errorchecked by a regional research manager, who assigned new individuals if no match was found. The sex and maturity status of manta rays were assessed based on the absence (female) or presence (male) of claspers, clasper size and/or clasper scarring for males, and pregnancy bulge and/or presence of pectoral fin mating scars for females (Marshall and Bennett, 2010;Marshall et al., 2011). Immature females were not assigned throughout the study because of the difficulty in definitively classifying females lacking maturity indicators in the absence of accurate size measurements. Manta ray sizes, based on disk width, were estimated relative to an entity of known size (i.e., swimmer or dive equipment) and placed into 0.25 m size classes. Only size estimates recorded by two observers that spanned the entire study period were considered (n = 260) in an effort to reduce observer bias. Injuries to the manta rays, including cephalic and pectoral fin truncations or disfigurement, and fishing line cuts and entanglements, were also noted (Supplementary Figure 3).
Manta ray behavior was assessed from sighting data based on mouth and cephalic fin positioning, the presence and absence of cleaner fish, and records of social interactions, including mating trains as described by Marshall and Bennett (2010), Deakos (2012), Jaine et al. (2012), and Stevens et al. (2018). Based on these previous studies, behavior was categorized into four mutually exclusive categories: foraging, cleaning, cruising, and courtship (Supplementary Figure 4). Two additional nonexclusive behavior categories were included -foraging/cleaning and courtship/cleaning to account for instances when an individual was observed performing two behaviors within the same encounter (maximum 60-min dive time). 'Foraging' manta rays were defined as individuals with open mouths and unrolled cephalic fins near the surface of the water column, or those who were completing full vertical rotations ('barrel-rolling'). 'Cleaning' manta rays were defined as individuals hovering over reef patches ('cleaning stations') with cleaner fish species (Chaetodon kleinii, Thalassoma lunare, Labroides bicolor, and Labroides dimidiatus) in close proximity. 'Cruising' manta rays were those sighted with closed mouths and rolled up cephalic fins. 'Courtship' behavior was recorded when several male manta rays pursue a mature female (in the absence of foraging behavior), or when two or more individuals gave acrobatic displays as described by Stevens et al. (2018).

Analyses of Residency and Movement
We use the term 'site fidelity' to denote the return of an individual to a site of previous residency after a periodic absence greater or equal to the residency period, following Chapman et al. (2015). However, a caveat of solely using this terminology is that some individuals might regularly use both study areas within the Nusa Penida MPA, moving frequently between the two. Further, periods between sightings varied greatly between individuals, precluding us from assigning a standard residency and absence period. In these uncertain scenarios, where individuals could be moving between areas across the large home range, we use the term 'site affinity' to describe our results (Couturier et al., 2011;Jaine et al., 2014).
Frontiers in Marine Science | www.frontiersin.org A modified maximum likelihood approach was used to compare manta ray re-sighting data against residency models, implemented in the program SOCPROG 2.8 (Whitehead, 2009). These statistical models were previously used for manta rays in Hawaii by Deakos et al. (2011) and in several more recent studies on whale sharks Rhincodon typus (Fox et al., 2013;Robinson et al., 2016;McCoy et al., 2018) and cetaceans (e.g., Chabanne et al., 2017). This approach determines the spatial and temporal distribution of sampling effort using the re-identification data itself, making this approach suitable for opportunistic and presence-only sighting data (Fox et al., 2013). Residency patterns for individuals within the study areas were investigated by using the 'Movement Analyses' module of SOCPROG 2.8 to calculate the lagged identification rate (LIR), the probability of re-sighting an individual manta ray after a variable lag time (Whitehead, 2009). These empirical LIR data, based on a per day time lag, were then compared to a series of model scenarios encompassing both closed and open populations. Open population models were developed by incorporating varying situations with emigration, immigration, re-immigration, and/or mortality ( Table 1). The lowest value from the quasi-Akaike information criterion (QAIC), used to account for over-dispersion of the data (Whitehead, 2007), was used to select which model best approximated the residency characteristics of manta rays at each study site. The sole exception was for female manta rays at MP, where AIC was used to select a best-fit model, as the data were not overdispersed. The QAIC difference (or QAIC) between the bestfit model and any other indicates how well the data support the less favored model as follows: QAIC < 2 = substantial support; QAIC 4-7 = considerably less support; and QAIC > 10 = essentially no support (Whitehead, 2007). The LIR analyses were extended to a 'within/between' analysis to estimate the probability that an individual manta ray, identified first within one site, will be re-sighted at the other site after a specified time lag (in days). This analysis effectively provides a significance test for population-level mixing between sites. Model fits were bootstrapped 1,000 times to generate standard errors (SE).
The transition probabilities between sites were calculated using a parameterized Markov movement model (Whitehead, 2009), in which an individual has a certain probability of moving from one area to another at the specified time lag (1 year; Table 2). This model includes movements to and from both MB and MP and also accounted for movements to a third, hypothetical area (i.e., an area 'outside' of MB or MP). Optimized values of transition probabilities were bootstrapped 1,000 times to generate SEs and the maximum number of evaluations was set to 100,000. Mortality, including permanent emigration from all areas within the MPA, was considered in the model.

Statistical Analyses
The correlation between the annual number of survey days and sightings was estimated using Pearson's product-moment correlation, cor.test function of the R statistical software (R Core Team, 2018). Differences in the number of individuals and sightings between the sexes at MB and MP were tested

Sightings and Survey Effort
Data from 6,087 sightings (excluding daily duplicate data for the same individuals; n = 163) of manta rays were collected within the Nusa Penida MPA from 25 June 2004 to 9 April 2018. Sightings came from 965 unique dates, with ∼ 97% of the records (5,913 sightings, excluding daily duplicates) logged between January 2012 and April 2018 and from our two main focal sites of MB and MP (Supplementary Table 1). Data prior to 2012 allowed 136 individual manta rays to be identified from 174 sightings within the MPA. After 2012, as education about manta rays amongst tour operators in the area improved, the number of manta ray identification photos and sighting information available increased (see Supplementary Figure 2). The average annual sighting rate from 2012 until 2017 (i.e., the full survey years) was 917 sightings/year ± 125 (± 1 SE), while 412 sightings were recorded for the first quarter of 2018. The number of annual manta ray sightings varied across the study period (Figure 2A and Supplementary Table 1), with sightings being positively correlated with survey effort (r = 0.92, p < 0.01) for all study years. As the days where no manta rays were sighted (i.e., absence data) are not logged with 'Manta Matcher, ' we used survey dives logged by trained observers as a proxy for survey effort (Supplementary Figure 2). The number of annual manta ray sightings submitted to 'Manta Matcher' varied between MB and MP. In the earlier study years (i.e., 2012-2014), more sightings were recorded (daily mean and total) at MB than MP, while there were more sightings at MP in the later years (i.e., 2015-2017) (Figure 2A and Supplementary  Table 1). These observations reflect the gradual increase in reported surveys at MP and the steep decrease of reported surveys at MB after 2014, influenced by a shift in site usage by tourism operators (i.e., more dive operators using MP and more snorkel operators using MB; EG pers. obs.).
Logged presence/absence data indicated variable survey effort by trained observers between years and months (Supplementary Figure 5), with a lack of complete year-round effort in 2012 and 2015. Further, annual weather patterns, such as the north-west monsoon season (December-January), that result in unfavorable weather conditions and limit access to the study sites, regularly contributed to a seasonal reduction in survey effort. Nevertheless, there was an approximately twofold increase in manta ray sightings for MP in the month of May ( Figure 2B). In contrast, sightings at MB were only marginally higher from February to April and in July than other months of the year.
A total of 624 individual manta rays were identified from sightings within the Nusa Penida MPA between 2012 and 2018. Of these individuals, 407 were sighted at MB and 494 were sighted at MP, with 277 sighted at both sites. A discovery curve shows a steep rise in newly identified manta rays within the Nusa Penida MPA until approximately 400 individuals, but the curve kept increasing to 624 individuals, though at a slower rate, until 2,260 days, with no evidence of an asymptote (Supplementary Figure 6). When treated separately, neither of the discovery curves for the two study sites show signs of approaching an asymptote. However, a brief period of slower discovery (i.e., a plateau) of new individuals was observed at MB after the identification of approximately 280 individuals at 800 days, until 1,600 days or between April 2014 and December 2015 (i.e., years 2 and 3 of data collection), after which the numbers of individuals identified increased to 407 after 2,260 days.

Population Structure
Overall, more male (n = 344; 55%) than female (n = 243; 39%) manta rays were identified in the Nusa Penida MPA (χ 2 1 = 17.378, p < 0.001; Figure 3), while the sex was unknown for 37 individuals (6%). There was a significant association between sex and site for individuals (χ 2 1 = 13.864, p < 0.001), and between sex and site for sightings (χ 2 1 = 1080.2, p < 0.001). Of the entire population, at least 53% (n = 129) of the females were considered to be sexually mature, based on their pregnant appearance or the presence of mating scars. Among the males, 81% (n = 277) were considered sexually mature based on the size and shape of their claspers, with the remaining 19% (n = 66) being immature. In addition, 123 manta rays (∼ 20% of the total identified) were classed as immature at some point throughout the study, with 57 known to reach maturity. The population structure of individuals and the demographics of sightings differed between MB and MP (Figure 3). At MB, males made up 64% (n = 261) of identified individuals, whereas only 29% (n = 120) were female (χ 2 1 = 52.18, p < 0.001). Males were also more frequently encountered (n = 2,388; 80%) than females (n = 564; 19%) at MB (χ 2 1 = 1127, p < 0.001) and immature males had the highest proportion of the total sightings (n = 1,170; 39%). At MB, three males and one female were estimated to be 1.5 m in width on their first recording and another thirteen individuals were estimated to be ≤ 2 m, while 40 individuals were estimated to be ≥ 3 m. At MP, the majority of the individuals were mature males (n = 240; 49%) and females (n = 127; 26%), with overall more males (n = 268; 54%) than females (n = 212; 43%; χ 2 1 = 6.53, p = 0.011). However, females were more frequently encountered (n = 1,759; 60%) than males (n = 1,142; 39%; χ 2 1 = 131.23, p ≤ 0.001). Immature males were rarely encountered at MP, comprising just 2% of sightings (n = 45), a nearly 20-fold difference to MB. Five manta rays were estimated to be ≤ 2 m in width on their first recording with the smallest being 1.75 m, while 130 were estimated to be ≥ 3 m. The maturity status could not be determined (unknown) for 66 and 85 females (on 16 and 17% of total sightings) at MB and MP, respectively.

Injury Rates
Notably, 87 (∼ 14%) of identified manta rays had either cephalic fin, pectoral fin, and/or fishing line injuries, including 11 that were observed to be pregnant during the study period (three  of the individuals were sighted pregnant two times with a span of at least 1 year in between). A total of 48 manta rays (7.7%) were observed entangled in fishing line or with fishing line marks (Supplementary Figure 3) and 49 (7.9%) had permanent injuries, either cephalic (n = 29; 4.6%) and/or pectoral fin (n = 21; 3.3%) truncations or disfigurements.

Residency and Movement
Of the 624 uniquely identified manta rays, the majority were encountered more than once (82%), with 181 individuals (29%) sighted more than 10 times, and 30 individuals (5%) sighted 31-99 times (Figure 4). The mean number of re-sightings per individual was 9.5 (± 0.46). The majority (n = 514; 82%) of these individuals were re-sighted across multiple years and 18% (n = 110) were sighted in only 1 year. The longest time between re-sightings of an individual was 13.8 years (including pre-2012 data). Overall, of the ten most re-sighted individuals, nine were male and one was a female. The most re-sighted individual was a male that was seen 99 times, but at MB only, between February 2012 and October 2017. This individual was classified as a juvenile at the first encounter, in February 2012, but matured over the monitoring period and was reclassified as an adult in July 2016. Similarly, the other nine most re-sighted individuals at MB were all male, with six out of 10 classified as juveniles upon first encounter (Supplementary Figure 7A). The most re-sighted female was seen 56 times (54 times at MP and 2 at MB) between June 2015 and April 2018. Her maturity status was unknown until March 2017, at which point she was classified as mature as she appeared to be pregnant. Similarly, the other nine most re-sighted individuals at MP were females, and all were seen pregnant at least once during the study period (Supplementary Figure 7B). The average period between re-sightings across all individuals (i.e., lag) was 118 days (± 3), ranging from 1 day to 5.3 years. For the most sighted male and female, greater than 90% of re-sightings had a lag period of less than 2 months, with the longest lag periods of these individuals being 234 and 239 days, respectively. We found that 277 (44%) individuals moved between the two study sites, and tested whether the two sites were fully mixed (a single population) using LIR analysis ( Figure 5A). The two LIR curves did not converge during the study period, indicating that, although there is interchange between MB and MP, most individuals have a high affinity for one site over the other (i.e., communities were not mixed), at least when considering time lags of 1,500 days (∼ 4.1 years) or less. This notion is supported by the Markov movement model, which indicated high site affinity rates for MB (64.7%) and MP (73.7%) ( Table 2A). A higher percentage of individuals moved from MB to MP (35.3 ± 9.9%) than vice versa (16.4 ± 6.5%) ( Table 2A). Further, a higher percentage of individuals move away from MP (9.9 ± 15.7%) to 'outside' areas than from MB (0.0 ± 1.6%), indicating that the MP community is more mobile.
Results of the 'within/between' analysis showed that the sites should be treated separately for detailed analysis. Based on the QAIC (Table 1A), the best fitting model for both MB and MP was Model H, which described a pattern of emigration, reimmigration and mortality. Models C and D, where emigration and mortality are included, were also supported (to a lesser degree) for MB (based on QAIC < 2). Although a similar number of sightings were recorded for both areas, re-sighting rates are initially (up to two-year time lag) much higher at MB than MP, indicating that manta rays tended to reside longer in MB than MP ( Figure 5B). Model H scenarios provide information on the proportion of time individuals spend within (residence time in) and out (residence time out) of an area, as well as mortality. Individual manta rays stayed approximately twice as long (2.0 days ± 187.7) at MB than at MP (0.9 days ± 0.2). Further, individuals sighted at MP spent an average of 2.9 (± 0.6) days outside the site, being absent nearly twice as long as those initially sighted at MB (1.6 ± 3.1 days). However, estimations of residence time within and out of MB were variable, suggesting that sex and/or demographic-specific (i.e., males vs. females, immature individuals vs. adults) or high individual differences in residence patterns exist. Mortality (which includes permanent emigration) was negligible for both sites (∼ 0).
To further explore sex specific differences in LIRs, we ran separate analyses for each site and sex (Tables 1B,C). Based on the QAIC (except for females at MP, which  was based on AIC), the best fitting models were Models E (closed population with emigration and re-immigration) and H (emigration, re-immigration and mortality) for males and Model G (emigration, re-immigration and mortality) and F (emigration and re-immigration) for females at MB and MP, respectively (Tables 1B,C). We further explored maturity influenced differences in LIR for males (immature females were not assigned) at MB (Supplementary Table 2) using sighting records (immature = 487; mature = 707) of individuals (immature = 59; mature = 127) whose maturity status remained the same throughout the duration of the study (January 2012 to April 2018). The best fitting models (based on QAIC) were Models C and D (emigration/mortality) for immature males and Model H for mature males. Site and sex-specific LIR indicated that, on average, female manta rays that use MB, and to a lesser extent MP, are two times more likely to be resighted at the site of first sighting during the first two years after their initial identification than males (Figures 5C,D and Tables 2B,C). However, immature males were twice as likely to be re-sighted at MB than mature males (Supplementary  Figure 8). The Model H scenario indicates the residence time in for mature males was 32.3 days (± 47.0), while Model D estimated the mean residence time for immature males to be ∼ 23-fold longer (745.3 days ± 169.9). Markov movement models indicated differential movement patterns with males, moving from MP to MB (24.5 ± 2.8% SE) at a higher rate than females (8.1 ± 2.2% SE). Further, a modestly higher percentage of females are staying within MP than males (85.8% vs. 75.5%). Notably, the model maximum number of iterations was exceeded, suggesting that the estimated SEs are inaccurate and larger than expected (H. Whitehead, Dalhousie University, personal communication). Nevertheless, taken together these results show a lower probability of movement between the study sites (i.e., higher residence) for female manta rays than mature males; although, immature males display high residency to MB.

Behavior
A range of manta ray behaviors (Supplementary Figure 4) were observed (n = 3,692) at both study sites, including foraging, cleaning, courtship, cruising and a mix of cleaning and courtship, and cleaning and foraging (Figure 6). However, the frequency of observed behaviors differed between the two sites, and a significant association between behavior and site was present (χ 2 3 = 2303.7, p < 0.001). From the 1,211 sightings with recorded behavior at MB, manta rays were foraging in 79% of sightings (surface ram feeding and near surface barrel-rolling), with cleaning and courtship behaviors observed in 7% and < 1.5% of sightings, respectively. In contrast, at MP (2,481 sightings with recorded behavior), manta rays were foraging in < 4% of sightings. The principal observed behaviors at MP were cleaning (59%), courtship (7%), and courtship combined with cleaning (8%) during single encounter sessions. Courtship behavior was observed year-round at MP but peaked in April-May ( Figure 6B). More cruising behavior was observed at MP (22%) than MB (13%).

Boats
The average annual number of boats (7.2) present at the two manta ray sites increased by 60% (Supplementary Figure 9 Figure 9A). The monthly average number boats at MP over all years was higher from April-October (10.9-14.0), than the other months (< 8.3), with peaks in May and September (Supplementary Figure 9B).
In contrast, the average number of boats at MB was relatively constant throughout the year (5.4-7.9), peaking in August.

DISCUSSION
We found the Nusa Penida MPA was an important habitat for a substantial number of manta rays (624 identifications in approximately 6 years), with manta rays present yearround and showing high site fidelity to two sites, ∼ 12 km apart. Individual re-sighting rates reported in the Nusa Penida MPA are higher than those reported for any other M. alfredi population in the world (Table 3), with 82% of the individuals re-sighted more than once. The manta rays display high site fidelity to the west coast of Penida Island and are rarely reported from other parts of the MPA. Our analyses highlighted that the two main study sites, MB and MP, are used for different purposes by different life stages of manta rays. MB is a foraging ground, used primarily by immature individuals, whereas MP is an important site for adult social and reproductive activity. Female and male manta rays also used the sites differently, with more males frequenting MB than MP. As a whole, the sex bias found in the Nusa Penida MPA (males: females = 1.4:1) is the largest reported bias toward males worldwide (Table 3). On the other hand, the proportion of females confirmed as sexually mature (53%) is greater than anywhere else in the world (Table 3). This clear segregation in habitat use by different manta ray demographics at almost adjacent sites is noteworthy for marine spatial planners and conservationists globally.

Manta Bay
Unlike other manta ray aggregation sites studied to date, MB was used mostly by males with a few resident females. The largest percentage of sightings at MB (38%) were of 59 immature males (one sighted 99 times during the study period). This is particularly noteworthy as few studies report such frequent sightings of juvenile manta rays (Stevens, 2016). Thus, the steady accumulation of new individuals being identified at this site may reflect recruitment of juveniles to an important foraging ground, or potentially a nursery, following birth. The smallest manta rays recorded in the duration of the study were approximately 1.5 m in disk width, which is within the size range of newborn M. alfredi (1.3-1.6 m; Stewart et al., 2018a). Future research looking into the size of manta rays using standardized measurement techniques (Costa et al., 2006;Deakos, 2010;Couturier et al., 2014) would further our understanding of how different age classes use these habitats. Continued recruitment of new, immature, individuals to MB suggests that females give birth nearby, which is a topic worthy of future research. The location of birthing grounds (i.e., pupping grounds) and (or existence of) nursery grounds of manta rays and other mobulids are major outstanding questions in mobulid ecology (Stewart et al., 2018a). The first nursery habitat described for M. birostris in Flower Garden Banks National Marine Sanctuary in the northwestern Gulf of Mexico (Stewart et al., 2018b), followed the most recent re-description of the term 'nursery' for elasmobranchs, proposed by Heupel et al. (2007), and hinges on: (1) more immature individuals are found in this area relative to others nearby; (2) immature individuals reside for extended periods in this area relative to others nearby; and (3) immature individuals show long-term habitat use of this area relative to others. Thus, nursery areas are expected to be areas that provide beneficial opportunities to the young, such as increased food availability (Heupel et al., 2007;McCauley et al., 2014), protection from predation, or areas for thermoregulation via basking behavior (Stevens, 2016;Stewart et al., 2018b). The high re-sighting rates of immature manta rays, the observation of foraging as the dominant behavior, and meeting the criteria proposed by Heupel et al. (2007) and others for a nursery area (McCauley et al., 2014;Stevens, 2016;Stewart et al., 2018b), provide strong evidence that MB forms part of a nursery habitat for the local M. alfredi population. The proposed nursery criteria are aligned with our study in that (1) more immature individuals are found at MB than the nearby MP; (2) much longer residency of immature individuals is observed at MB than at MP; and (3) manta rays identified as immature on first sighting show long-term habitat use of MB relative to MP.
In oligotrophic environments, immature and smaller sized manta rays (more often male because females reach a larger maximum size) possibly prefer to aggregate in shallow coastal areas that offer protection from predators (McCauley et al., 2014;Stevens, 2016). In the Nusa Penida MPA, smaller individuals might use MB as a prominent foraging ground because of reliable food availability. The semi-enclosed nature of the bay, located at the base of the ITF through the Badung Strait and neighboring channel between Ceningan and Penida Islands, and proximity to deep water to the south, could entrain plankton and thereby provide consistent shallow water foraging opportunities for these immature rays.

Manta Point
The manta rays sighted at MP were typically larger than those at MB, with mature females making up 40% of sightings. Thus, the recruitment of new individuals to MP is probably due mainly to immigration, not births. Notably, we found the highest percentage of mature females at this site ever reported (63%) globally ( Table 3). The most frequent behavior observed at MP was cleaning (67%). The removal of parasites and other external fouling and the promotion of wound healing is likely an important service provided by cleaner fish (Marshall, 2008;Stevens, 2016). Cleaning activity appears to be a daily ritual, with individuals sometimes spending hours at cleaning stations during the day (Dewar et al., 2008;Marshall, 2008;O'Shea et al., 2010).
Visits to reefs serving as cleaning stations may also provide an opportunity for social interactions in elasmobranchs, including mobulids (O'Shea et al., 2010;Oliver et al., 2011;Murie and Marshall, 2016;Stevens, 2016). Both cleaning and courtship behaviors were observed during 8% of sightings at MP. Courtship behavior is clearly initiated at predictable aggregation sites for mature individuals, with cleaning stations potentially acting as lekking sites (Stevens, 2016). Individuals might visit these aggregation sites during the breeding season in search of mates, rather than for cleaning (Deakos et al., 2011;Stevens, 2016). These suggestions are consistent with the peak in reported sightings at MP during the main reproductive months of April-May. The seasonal pattern of manta ray courtship behavior varies throughout the world (Marshall and Bennett, 2010;Deakos et al., 2011;Couturier et al., 2014;Stevens, 2016), with no single driver for initiating courtship identified to date. Reproductive behavior is likely to be linked to seasonal productivity in the water column, as fecundity has been linked to food availability (Stevens, 2016).

Site Fidelity and Productivity
Manta rays are very mobile species and range widely in both archipelago environments (Germanov and Marshall, 2014;Conservation International Indonesia, 2016) and along continental coastlines Jaine et al., 2014). However, recent studies (Couturier et al., 2018;Setyawan et al., 2018) show that they do not use the extent of their home range uniformly, and that individuals concentrate the majority of their activities (foraging, cleaning, and courtship) within specific critical habitats. While our study focused on sightings in two specific areas, it is evident that habitat usage within the Nusa Penida MPA is not uniform.
The year-round sightings of manta rays and their large population size in the Nusa Penida MPA are likely to be linked to the sustained productivity in the region (Surinati, 2009;Ayers et al., 2014;Tito and Susilo, 2017;Nyegaard, 2018). Food availability and sustained foraging opportunities are thought to be responsible for dictating the size of manta ray populations aggregating in a particular area (Deakos et al., 2011) and the frequency of their pregnancies (Stevens, 2016). In most of the manta ray populations studied worldwide, seasonal increases in manta ray abundance coincide with increased oceanic productivity (Dewar et al., 2008;Jaine et al., 2012) and in prey density (Armstrong et al., 2016;Rohner et al., 2017). For example, the largest estimated manta ray population globally is found in the Maldives (Kitchen-Wheeler et al., 2011;Stevens, 2016), where alternating monsoons result in year-round productivity (Anderson et al., 2011). The finescale oceanography of the Nusa Penida area is poorly known at present, though large differences in water temperature, indicating regional upwelling and the potential for nutrient mixing and increased nutrient availability to support primary productivity, were recorded over distances of about 30 km from the northern coast to the south-western coast during research on sunfish populations in the region (Tito and Susilo, 2017;Nyegaard, 2018). Studies that evaluate both the presence/movements of filter feeders and the concurrent biological and physical oceanographic characteristics Rohner et al., 2013) would enhance the understanding of manta ray ecology and the drivers of differential site use, increased seasonal abundance and reproductive behavior in this complex area.

Movement Patterns
Although the manta rays of Nusa Penida display some site preference, emigration and interchange between MB and MP does occur. Almost half of the individual manta rays recorded (n = 277; 44%) visited both locations during the study period. However, the interchange between the two sites was not symmetric and differed between the sexes. Higher transition probabilities for movement were observed from MB to MP (∼ 35%) than for the reverse movement (∼ 16%), raising the possibility that a longer-term study may demonstrate that the two sites, taken together, constitute a single population. The higher percentage of movement from MB to MP could be individuals leaving their nursery ground as they grow larger, to forage offshore and spend more time interacting with other mature individuals (McCauley et al., 2014;Stevens, 2016).
Further, the sex-based differences in habitat use for manta rays might be mainly linked to reproductive behavior (Deakos et al., 2011;Stevens, 2016). Similar to other pelagic species of marine megafauna including sharks, sea turtles and cetaceans (Hueter et al., 2005;Lee et al., 2007;Engelhaupt et al., 2009), female manta rays, as the sex with the greater parental investment (i.e., the 'limiting sex'), gain a greater choice of mates by residing in a popular aggregation area, assuming that ample foraging opportunities are available nearby. Males, in contrast, benefit from moving between aggregation areas in search of mates.
Consistent with this hypothesis, female manta rays from the Nusa Penida MPA are likely to reside in either of the sites for a longer time than males.
Moreover, the Nusa Penida MPA manta rays occasionally move to areas 'outside' MB and MP. Previous research (Germanov and Marshall, 2014;Conservation International Indonesia, 2016) shows that manta rays move long distances from the Nusa Penida MPA to locations such as the Gili Islands (80 km straight-line distance north-east), Sumbawa, and the Komodo National Park (up to 450 km straight-line distance east).

Implications for Management
The Nusa Penida MPA serves as year-round critical habitat for a substantial local manta ray population making it a high priority area for conservation and management. MB is used year-round for foraging, particularly by juvenile manta rays, while MP is an important area for cleaning and reproductive behaviors that peak in May each year. This has implications for MPA management, zoning and the types of activities to be allowed in these areas during distinct times of year. Management actions should address the following threats to the local manta ray populations: (1) disruption of manta ray behavior through habitat crowding and human disturbance from excessive tourism (Venables et al., 2016;Barr et al., 2017;Trave et al., 2017) and (2) entanglement in fishing gear and injury from continuing fisheries in the area (present study, Aquatic Alliance, personal communication).

Manta Ray-Focused Tourism Impacts
While tourism to the area is increasing, there are currently no regulations in place for manta ray interactions (Barr et al., 2017). A 2016 survey (Barr et al., 2017) directed at divers and snorkelers participating in manta ray watching activities, within the Nusa Penida MPA, reported participant conduct that is considered disruptive to the manta rays (e.g., closely following and touching manta rays was reported by 10% and 3.5% of participants, respectively). Disturbance stimuli to animals, such as loud boat engine noise, fast approaches by boats, in water chasing, touching, and crowding behavior by tourists (Norman, 1999;Bejder et al., 2006a,b;Anderson et al., 2010;Higham et al., 2016;Venables et al., 2016;Barr et al., 2017), is argued to be analogous to a predation risk (reviewed by Frid and Dill, 2002). Thus, tourism pressure at MB might disrupt foraging behaviors, reducing growth rates of immature manta rays and the fitness of mature individuals (Stewart et al., 2018a). Further, higher boat numbers at MP in May (Supplementary Figure 9) coincide with the seasonal increase in sightings at MP and in reproductive behavior. With the Nusa Penida MPA serving as an important reproductive ground, disruption of these important social behaviors is a concern. Based on the criteria to evaluate marine wildlife tourism practices outlined by Trave et al. (2017), we recommend that: (1) science based carrying capacity calculations of tourism operations be carried out to estimate the acceptable number of tour boats and diver interactions for the area (Ríos-Jara et al., 2013;Zelenka and Kacetl, 2014) potentially limiting the number of boats/divers/swimmers allowed at one time (Division of Boating and Ocean Recreation, 2016; Kasmidi and Gunadharma, 2017); (2) codes of conduct for diving and snorkeling with manta rays become mandatory (see Garrud, 2016;Venables et al., 2016), akin to regulations for whale shark interactions in Ningaloo Reef, Australia (Mau, 2008;Catlin and Jones, 2010); (3) a licensing system for tour operators with penalties for breaches be implemented (Mau, 2008;Division of Boating and Ocean Recreation, 2016); and (4) area-time closures be considered as a management option (Tyne et al., 2014;Setyawan et al., 2018) to protect the manta rays from disturbance during the peak time of mating, especially at MP.

Artisanal Fisheries
Small-scale fishing often takes place in MB and nearby coastlines (EG, pers. obs., Supplementary Figure 10). Over the course of the study, ∼ 14% (n = 87) of manta rays were observed either trailing hooks and lines or with cephalic and pectoral fin injuries and/or amputations, which occur when fishing lines or nets cut through the skin and cartilage skeleton. The true scale of the issue is under-represented as only animals that have survived the entanglement event are counted. Currently, the long-term impacts of these effects on individuals that survive entanglement are not known, nor whether foraging and swimming efficiency is impaired and if there are any impacts on fecundity. However, observations of pregnant individuals within the Nusa Penida MPA, with single cephalic fin amputations, suggest that manta rays retain their reproductive fitness even with these sub-lethal injuries. Further, the proportion of pregnant manta rays from the Nusa Penida MPA with injuries (0.11) was not significantly lower than the overall proportion of females with injuries (0.12; χ 2 1 = 0.676, p > 0.5), suggesting that injured manta rays are not substantially impaired reproductively. Nevertheless, it is unclear whether these injuries have an initial impact on the success of their pregnancies, as many related species (e.g., batoids) will abort their fetuses if their individual survival is threatened (entanglement, landing, predation; reviewed by Adams et al., 2018). Further, biomechanical modeling of manta ray pectoral fin movements indicate that the major thrust force comes from the distal portion of the pectoral fins (Liu et al., 2015), thus pectoral fin truncations could significantly impact manta ray swimming efficiency, energy consumption and ability to evade predation. While all fishing activities are officially prohibited in both MB and MP, increased enforcement of this regulation is necessary. As a precautionary measure, management could prohibit all fishing activities along the west coast of Nusa Penida, which would require strong compliance if it was to be effective in contributing to the recovery of manta rays in the Nusa Penida MPA.

CONCLUSION
Understanding localized habitat use, identifying nursery areas, as well as key times of year for reproduction, provides enhanced information for developing effective management plans for manta rays in this region related to specific manta ray behaviors. Context-dependent and adaptive management solutions are important for manta ray populations as well as for people, as manta ray watching tourism supports many livelihoods through activities that rely on healthy manta ray populations (i.e., snorkeling and SCUBA diving). With their 2014 protection in Indonesia, manta rays have entered the spotlight for conservation initiatives and serve as ideal flagship species to refocus future research on the overarching challenges and opportunities facing ocean health in the Coral Triangle (see also Germanov et al., 2018), the world's premier marine biodiversity hotspot.

ETHICS STATEMENT
This research was completed under Indonesian Research permit #458/FRP/E5/Dit. KI/XI/2015, #11K1/XI/2016; #86/KI/XI/2017. This study was carried out in accordance with the recommendations of the Animal Ethics Committee, Murdoch University. The protocol was approved (R2781/15) by the Animal Ethics Committee, Murdoch University. Publiccontributed photographic data were collected opportunistically by recreational divers and deposited by them into a public online repository (mantamatcher.org), which was developed explicitly to facilitate citizen science contributions.

AUTHOR CONTRIBUTIONS
EG and AM conceived the study. EG analyzed the data with LB, DC, and SP advising on the analysis. EG wrote the manuscript with contribution to drafting, critical review, and editorial input from LB, DC, DD, IH, NL, AM, SP, and MvK. All authors critically reviewed the manuscript with final proofs.

FUNDING
While no specific funding was received for this research, EG is supported by an Australian Postgraduate Award and Murdoch International Top Up (32608315), which supports her doctoral studies; Ocean Park Conservation Foundation (FH04_1516), which has funded a large portion of fieldwork and laboratory work for her doctoral studies; Foundation Fortuna, which supports the field work for her doctoral studies; Idea Wild, which has donated field sampling equipment; Mantahari Oceancare, which has funded manta ray research capacity building programs related to her doctoral studies; PADI Foundation (14668), which supports the field work for her doctoral studies; and private donors, which support the field work for her doctoral studies.

ACKNOWLEDGMENTS
This study would not have been possible without the generous support of Lembongan Marine Association, the local dive community and many citizen scientists to whom were are indebted. Further, we thank the Aquatic Alliance for spearheading this project and Marine Megafauna Foundation staff and volunteers for carrying it forward. Notably, we thank, H. Mitchell, P. Bassett, and L. Ellevog for their unwavering dedication to Indonesia's manta rays and enabling this study. Further, we appreciate the data processing assistance provided by S. Ecob, L. Auditore, R. Cooper, and E. Sinderson. We thank G. Winstanley who provided the use of 'MantaUtil' for streamlined data processing and facilitated 'batch' lodging of legacy data to 'MantaMatcher.' We thank H. Whitehead for the guidance he provided on the use of Markov movement models, M. Calver and R. Admiraal for statistical analysis assistance, and S. Venables for plotting using R. Bathymetry information used in the creating Figure 1 was available from: GEBCO_2014 Grid, version 20150318; www.gebco.net. We thank members of Faculty of Marine Sciences and Fisheries, Udayana University, Bali, Indonesia and C. O. Lafuente for assisting in creating Figure 1. We are grateful for the study permit from the Indonesian Ministry of Research (Ristek Dikti). This manuscript represents HIMB and SOEST contribution numbers 1755 and 10681, respectively. We are grateful for the reviewers' input, which much improved the manuscript.