Quantifying the Effects of Diver Interactions on Manta Ray Behavior at Their Aggregation Sites

Manta rays (Mobula birostris, Mobula. cf. birostris, and Mobula alfredi), the largest mobulid rays, are subjected to exploitation and overfishing in certain parts of the world. Tourism has been supported as a sustainable alternative for the conservation of the species, and a potential source of economic spillover to local populations. Nevertheless, the effects of tourism over these highly social animals remains unknown. Manta rays aggregate at three sites in Mexico: Oceanic manta rays (M. birostris) in The Revillagigedo Archipelago and Banderas Bay in the Pacific. Caribbean manta rays (M. cf. birostris) around Isla Contoy National Park in the Caribbean. We analyzed the behavior of manta rays using video data collected by local researchers and tourism operators to determine how diver behaviors and techniques (SCUBA and free diving) affect them. Diver activities were grouped into passive and active categories. We described 16 behaviors and grouped them into four behavioral states: Directional, erratic, attraction and evasion to divers. We modeled the sequence of behaviors exhibited by manta rays via first order Markov chains. Our models accounted for passive and active diver behavior when modeling the changes in manta behavior. Manta rays in Banderas Bay and Revillagigedo displayed a higher frequency of erratic behaviors than at Isla Contoy, while Banderas Bay manta rays transitioned to evasion behaviors more often. Manta rays responded similarly in both sites to active divers. At freediving sites, manta rays from Isla Contoy displayed evasion less frequently than at Banderas Bay. Changes in manta ray behavior were similar for both sites, but mantas in Banderas Bay transitioned to evasion more with active divers. The increased food availability for Isla Contoy manta rays could be the reason for the reduced response toward divers in this site. The existence of additional stressors such as both traffic in Banderas Bay could be causing the mantas in this site to respond more frequently to active divers. This study, the first of its kind in oceanic and Caribbean manta rays, highlights that regulations and the use of best practices are vital for achieving longer and less disturbing encounters for both manta rays and divers.


INTRODUCTION
Wildlife tourism and ecotourism are potential tools to alleviate the economy of local communities while supporting conservation of their natural resources (Stronza et al., 2019). In order to accomplish this, proper management and scientific monitoring programs need to be established on site, especially when dealing with endangered species or animals susceptible to sensorial disturbances and habitat degradation (Higham, 2007;Buckley, 2011).
Tourism can be especially beneficial for the management of marine wildlife traditionally subjected to extractive exploitation such as sharks and rays (Mazzoldi et al., 2019). As such, elasmobranch wildlife tourism has been growing rapidly in the last several decades, being a sustainable alternative to extractive fishing and as an important source of income (Zemah-Shamir et al., 2019). Manta rays (oceanic, Mobula birostris, Caribbean, Mobula. cf. birostris and reef, Mobula alfredi) are a group of elasmobranchs in which tourism management programs have been particularly promoted due to their vulnerable status, high exploitation and slow growth and reproductive rates (Murray et al., 2020).
Manta ray tourism has been extensively developed in some parts of the world with revenues that range from 2.4 million US dollars in Hawaii to up to 14 million US dollars in the Revillagigedo Islands, Mexico (Osada, 2010;Ruiz-Sakamoto, 2015). More than 10 countries are engaged in manta ray tourism generating an estimated 140 million US dollars globally in 2013 (O'Malley et al., 2013). Most tourism-related studies at manta ray aggregation sites focus on conservation efforts or on the economic impacts of the newly established tourism industries (Atkins, 2011;Venables et al., 2016). Only a few studies have tried to measure the effects of tourism on reef manta rays' behavior (Murray et al., 2020), but due to the generally less accessible aggregation sites of oceanic and Caribbean manta rays, there are no studies to date that have attempted to measure the impacts of tourism on this species. Manta rays are protected by international conservation efforts and legislation such as Appendix II of the Convention on International Trade in Endangered Species of Wild Fauna and Flora (CITES, 2013). In Mexico, local laws such as NOM-059-SEMARNAT-2010 prohibit their capture, possession, extraction and commercialization. Despite these efforts, oceanic manta rays were re-classified as "endangered" on the IUCN's Red list of Threatened Species in 2020 (Marshall et al., 2020). This fact highlights the necessity to evaluate and improve our management and conservation strategies, as well as our understanding of the ecology and behavioral responses of these species in its different habitats.
Ongoing research programs in countries with several manta ray aggregation sites such as Mexico (Graham et al., 2012;Stewart et al., 2016a) have been studying the ecological aspects of the species while following best practices manuals used in other areas. The (i) development of local best practices manuals and protocols for tourists, and (ii) minimizing potential adverse effects by tourism service providers are necessary to improve our research and conservation efforts on the species (Lawrence et al., 2016).
Behavior catalogs or "ethograms" are lists used in animal behavior studies that describe the most common behaviors displayed by animals under certain conditions (Altman, 1974). Opportunistic video footage has proven to be an effective method to capture and analyze animal behavior (Ghaskadbi et al., 2016). This method has been used previously to collect data from species such as bull sharks Carcharhinus leucas (Pasos-Acuña, 2018) and killer whales Orcinus orca (Pagel et al., 2016). Behavior catalogs can be used as a baseline to improve the quality of encounters and to develop best practice manuals (Lawrence et al., 2016).
The baseline produced by behavior catalogs can be further developed into models by creating categorical behavioral states using the entries of the catalog. Each category includes a behavior or a series of behaviors so that they can be measured as a sequence over a time period. Probabilities can then be calculated for the transition of these sequences or "states" by assuming the sequences of states are realizations of a Markov chain (Lusseau, 2003). Markov chains model the changes in a sequence of events (states) by estimating the probability of a state to remain constant or to change to a different state using transition probability matrices (Lusseau, 2003;Peters et al., 2013). These models are useful tools to examine changes in animal behavior across time and at different time scales to evaluate the possible causes of the observed behavioral changes (Lusseau, 2004).
Given the interest in the development of sustainable ecotourism activities with threatened species (Lawrence et al., 2016) and the paucity of studies regarding the negative effects of poor practices by tourism operators on oceanic manta rays (Ruiz-Sakamoto, 2015) the goal of this study was to describe and compare the behavior of oceanic and Caribbean manta rays in the presence of divers at their main aggregation sites in Mexico, by creating behavior catalogs and developing first order Markov models using the observed behavioral states. This study will provide decision makers in Mexico with baseline information about the behavior and the effects of tourism on manta rays to improve management and conservation efforts for this emblematic species.

Study Sites
This study was carried on using data from the three main manta ray aggregation sites in Mexico (Figure 1): (1) The Revillagigedo Archipelago, a pelagic group of islands located 600 km west of the Mexican Pacific coast.
(2) Banderas Bay, a highly productive coastal upwelling area between the states of Jalisco and Nayarit in the Mexican Pacific (Stewart et al., 2016a). Research projects and tourism activities amongst the three sites have been producing basic knowledge about the seasonality FIGURE 1 | Study sites. The general location of each study site is denoted by red rectangles. Each site is represented by a letter as follows: (A) Revillagigedo archipelago (B) Banderas Bay (C) Isla Contoy National Park. The number of recorded interactions on each site is showed next to a symbol that represents the diving gear used (snorkeling mask for freediving and diving tank for SCUBA diving). and activities of oceanic and Caribbean manta rays (Rubin, 2010;Graham et al., 2012;Stewart et al., 2016a).
The Revillagigedo Archipelago is located in the Mexican Pacific (18 • 20 and 19 • 20 N; 110 • 45 and 114 • 50 W) 400 km off Baja California Sur. Due to its high degree of endemism, ecological and economic importance, the archipelago was protected as a National Biosphere Reserve in 1994, and several years after listed as a World Heritage Site in 2016 (UNESCO, 2016), and finally it became the Revillagigedo National Park in 2017, the largest marine protected area in North America. Oceanic manta ray aggregations are common from October to January, perhaps due to feeding or cleaning areas around the islands (Hull et al., 2006). Oceanic manta ray encounters at Revillagigedo are frequent and expected, and there are several dive sites close to cleaning stations with an extremely high probability of interacting with oceanic manta rays while they visit cleaning stations (Rubin, 2010). This situation along with the high abundance of charismatic fauna such as sharks and dolphins in the area has prompted the development of a strong live aboard tourism industry in the site, with around 70 tourism cruises traveling toward the archipelago every season (Hull et al., 2006). It is important to note that divers in this site interact with oceanic manta rays moving along the water column in relatively deep sites were the almost never come in contact with the bottom.
The Banderas Bay study site is located in the coastal area between the states of Jalisco and Nayarit in the Tropical Eastern Pacific (20 • 07 and 21 • 08 N; 105 • 10 and 105 • 45 W). Its proximity to a 2,000 m deep submarine canyon contributes to strong local upwelling events (Plata et al., 2006). Oceanic manta rays are present at this site year round with a major visitation period from January/February to June (Stewart et al., 2016a). Ongoing studies suggest that oceanic manta rays may use the deep areas of the site as a feeding grounds during their migratory route (Stewart et al., 2019). Around 1.5 million tourists arrive at Banderas Bay each year (Medina and Arnaiz, 2017), many of which partake in marine tourism activities such as whale watching, snorkeling, SCUBA diving and recreational fishing. These activities along with the need to transport local workers to the main commercial areas in Puerto Vallarta, and the presence of artisanal fisheries cause high levels of boat traffic in many parts of the bay (Everitt et al., 2008). SCUBA divers' interaction with oceanic manta rays happens opportunistically during shallow (less than 30 feet deep) recreational dives around rocky reefs, while free divers interact with oceanic manta rays primarily during field research activities with mantas that are swimming but not feeding near the surface of the water. While free divers attempt to comply with the proper codes of conduct, the surveying nature of their interactions forces them to swim after the manta rays and display behaviors that would be common to uninformed free divers in order to gather invaluable data for oceanic manta ray study and conservation. Isla Contoy National Park (21 • 27 and 21 • 32 ; 86 • 46 and 86 • 47 W) is located in the Mexican Caribbean off the coast of the state of Quintana Roo. The water surrounding the Island benefit from seasonal upwelling events which sustain a manta ray population that exhibits peaks in abundance from July to September (Hacohen-Domené et al., 2017). These events also cause an important Whale Shark (Rhincodon typus) aggregation subjected to substantial tourism pressure (Hacohen-Domené, 2015) of an estimated 60,000 tourists yearly in nearby areas such as Isla Holbox (Ziegler et al., 2016). Although the tourism industry focuses on the whale shark seasonal aggregations, the cooccurrence of whale sharks and manta rays is a well-known event that is even advertised by some tourism service providers (pers. obs., KF). Recent phylogenetic studies confirm that these manta rays belong to a third putative species (M. cf. birostris, Hinojosa-Alvarez et al., 2016;Hosegood et al., 2020). However, there are no studies that show significant behavior differences in these mantas. They form similar aggregations to the mantas observed in the Mexican Pacific, and the nature of their interaction with divers is practically identical to our other free diving site (Banderas Bay). Therefore, we considered that our behavioral analysis is as equally applicable for the Caribbean manta rays as it is for their Pacific counterparts. When we do not refer to one species specifically, we are referring to both of them simply as "manta rays" hereafter. Free diver interaction in this site is very similar to the observed for Banderas Bay, with scientist and researchers that attempt to comply with the codes of conduct when interacting with manta rays, but they also have to swim after the Caribbean manta rays in order to gather scientific data vital for the protection of the species.

Data Collection
Our data was collected from opportunistic video footage from the three study areas. Most videos were collected from the data bases of each site research team. These data bases consist of videos taken during monitoring trips and photo identification attempts of the manta rays, as well as videos collected from volunteers and tourists who shared their footage of interactions during recreational dives and casual encounters with manta rays. We gathered additional footage by visiting local dive operators in La Paz, Baja California Sur (for the Revillagigedo Archipelago), Puerto Vallarta (for Banderas Bay SCUBA) and Isla Mujeres (for Isla Contoy). All videos were taken using underwater video cameras. Videos from the Revillagigedo were provided by researchers, who collected them from tourists during liveaboard trips to the archipelago. SCUBA diving footage was collected exclusively from tourists and dive operators. Most of the videos collected from free divers in Banderas Bay and Isla Contoy were recorded during photo identification attempts in well-known monitoring sites, all of these were taken by researchers or by volunteers under the supervision of researchers. Researchers in both freediving sites follow the general recommendations of best practices manuals for interacting with manta rays. Since videos are taken in order to record the unique spot pattern of individual manta rays, researchers and trained volunteers in both sites start recording as soon as they jump in to the water. We had access to both, successful and unsuccessful photo ID attempts, as well as videos in which the free divers did not attempt to take a photo ID due to the speed or position of the manta ray. These videos were categorized according to the free diver behavior displayed during the interaction. Some of which, such as swimming toward a manta ray are very similar to what an uninformed and eager tourist would do when trying to interact with the animals.
As a preliminary analysis, we recorded the most frequent behaviors of manta rays as well as their response to the divers' behavior. This information was used to create a behavioral unit catalog for the diver-manta ray interactions for both free and SCUBA divers. The catalog was constructed following the design proposed by Murray et al. (2020) for reef manta rays, which describes the behavior of divers according to their level of disturbance toward the manta rays, and the behaviors of manta rays according to their swimming patterns ( Table 1).
We conducted manta ray individual focal sampling of behaviors following Altman (1974). In each video, we recorded the sequence, duration and frequency of the behavior units of the focal manta ray in 10 s intervals. This first analysis by intervals was used for the initial construction of behavior catalogs. Afterward, we constructed a time series by measuring the duration of each manta ray behavior by second in a sequential manner until the end of each manta-diver interaction for the three study sites, which we used for the more detailed analyses described below. The presence of additional manta rays was not considered for the analysis of the focal manta ray behavior. Diver behavior was determined by observing all divers in frame interacting with a manta ray, including the diver behind the camera when they interacted with the manta too. We were unable to determine the exact number of divers interacting with a manta ray from the available footage, therefore we did not include this factor in subsequent analyses. We considered interactions to end when the focal manta ray left the frame for more than 5 s. A random set of videos was analyzed by the first author and the researchers from each site in order to check for observer bias. Once an agreement was reached about the interpretation and meaning of each behavior until only minor differences were found when observing the same video, the first author proceed to systematically analyze the remaining video data.
Additionally, we recorded the ventral markings of each individual manta ray when possible and compared them between one another. These markings are unique to each manta ray (Marshall et al., 2011) and along with other markings such as missing cephalic lobes or scars were used to make sure that different manta rays were analyzed during each interaction and to avoid pseudo replication. Videos from the same diving trip, site and date provided by and individual source (tourist or researcher) in which a manta ray could have potentially being recorded twice were discarded, keeping the longest video with the highest resolution. Manta rays from separate sets of videos in which the individual could not be identified were considered to be different (unidentified) manta rays for the purposes of this study. Videos where a manta ray's movements or a diver's position could not be determined were discarded. Note that we lack control footage of mantas on their own. This was decided given the difficulty to collect videos of manta rays undisturbed by a diver cameraman. Our results should be interpreted with this caveat in mind. However, our statistical design takes this into account by presenting the results in direct changes in probability. We compare this changes assuming that divers are always interacting with manta rays, using variables such as diver behavior and site as predictors. This is further detailed in the Bayesian model construction section bellow.
Our sample contained a total of 1,159 videos, 685 of which were high quality enough to be further analyzed: 244 from Revillagigedo, 343 from Banderas Bay for oceanic mantas and 98 from Isla Contoy for Caribbean mantas. We recorded 745 interactions in total from these videos: 293 from Revillagigedo (SCUBA), 79 from Banderas Bay SCUBA, 261 from Banderas Bay freediving and 112 from Isla Contoy National Park (freediving). The mean duration of interactions (± standard error) as well as the number of individually identified manta rays per site are shown in Table 2.
TABLE 1 | Specific behaviors unit description and behavior states of mantas in our study sites.

Manta behavioral unit Description
Behavior states Site Directional swimming Manta swims without changing its general course. Pectoral fins do not remain still more than 5 s Directional All Non-directional swimming Manta swims constantly changing its course or turning slowly without sudden movements. Pectoral fins do not remain still more than 5 s Erratic All Directional cruising Manta cruises without moving its pectoral fins for more than 5 s and without changing general course.

Directional All
Non-directional cruising Manta cruises without moving its pectoral fins for more than 5 s and changes its course or turns slowly without sudden movements.

Directional Isla Contoy
Listed under "sites" are the study areas were each behavior was observed. Researchers, collected from researcher data bases) duration in minutes, mean duration + standard error in seconds of recorded interactions between manta rays and divers, number of individual manta rays identified by their ventral markings and % of identified individuals that appeared in two or more and three or more videos.
Frontiers in Marine Science | www.frontiersin.org Personality based analysis were out of the scope of this study. Less than 50% of the manta rays could be identified, an even lower proportion appeared in more than one video of interactions ( Table 2). Recent studies suggest personalities could affect the response of certain individuals, even in aquatic animals, and that these "bolder" individuals could have an effect over the evolution of their species (Wolf and Weissing, 2012). The nature of our data makes it impossible to apply the required experimental design for these analysis (Yang et al., 2020).
We observed 16 behaviors in total by manta rays across all sites (Table 1). While some behavioral units were present at every site, others like the feeding behaviors (which were only observed in Isla Contoy) and "swimming toward bubbles" (possible only in sites with SCUBA divers) were unique to some sites. Videos were considered for analysis only if manta ray behaviors were clearly distinguishable along with their change in course and movements.
To simplify the analysis, we grouped our 16 behavior units into four broad behavioral states. These states were named after the technical but complete terms "Taxis" and "Kinesis" described as defined by Kendeigh (1961) Taxis includes behaviors caused in direct response to a defined stimulus. In this study, we defined divers as the stimuli causing these responses to manta rays. Kinetic behaviors are defined as movements displayed without a clear source of stimulus. In this study, kinetic behaviors encompass the different movements displayed by manta rays which were not caused by divers. In order to simplify the interpretation of our results, we present a detailed explanation of each behavioral state along with alternative ecologically relevant terms used for each state hereafter: (a) Directional movements/Kinetic directional. In this state, manta rays kept a clear directional swim without changing their general course or doing maneuvers such as turns and somersaults. Hereafter referred as directional. (b) Erratic movements/Kinetic non-directional. This category grouped kinetic behaviors in which the manta rays did not keep a clear swimming path or in which they twitched, turned or stopped swimming due to unknown stimuli. This category did not include seemingly negative, aggressive or evasive movements. Hereafter referred as erratic. (c) Attraction toward divers/Positive taxis: Positive taxis behaviors are those in which the studied individual moves toward the source of the stimuli (Kendeigh, 1961). This category included only two behaviors, where the manta ray is clearly attracted to divers: swimming toward bubbles and observing a diver. Hereafter referred as attraction. (d) Evasion toward divers/Negative taxis: Negative taxis behaviors are those in which the focal animal moves away from a defined source of stimuli (Kendeigh, 1961). In this study, this category included several behaviors in which the manta ray either swims away or changes its swimming pattern in response to divers' activities.
Hereafter referred as evasion.
We observed four diver activities toward manta rays from the video analysis: passive observation, following, diving underneath and obstruction of movement (Table 3). We grouped these categories of diver activities according to their position and movements with respect of the manta rays: Passive and Active. Changes in diver behavior during an interaction were unusual. If any of the divers displayed the active behavior at any point during the video, the general diver behavior was considered to be active for the whole interaction.

Bayesian Model Construction
We modeled the sequence of manta ray behaviors over time via a first-order Markov chain at each site. The evolution of the behaviors is governed by a transition probability matrix, , with entries: Where S t can take on the values of each of the four manta behavioral states, and γ ij represents the probability of transitioning from state (S) i to state j at time t (in 1 s units). The probability that the manta ray is first observed in a behavior is modeled by the initial state distribution, δ, with entries In order to understand the effect of the diver's behavior on the manta rays, we allowed for the entries of to be functions of the diver's behavior through the use of a multinomial logit link function considering the transition probabilities of manta ray behavioral states as response variables and diver behaviors as binomial categorical predictor, with w t = 0 for passive divers at time t and w t = 1 for active divers.
Because the divers were only in passive or active mode during the recordings, we obtain two transition probability matrices that describe the sequence of behaviors exhibited by the manta rays when the divers are active A or when the divers are passive P . We fit the model in a Bayesian framework in Stan (Carpenter et al., 2017) and specified prior distributions for β 0 and β 1 as follows: Listed under "cumulative" are the broad groups for diver activities: Passive and Active.
We obtained 1,000 posterior draws from each of 4 chains, for a total of 4,000 posterior draws for each parameter. We calculated 95% credible intervals for each entry of A and P following McElreath (2015) by using the BayestestR package (Makowski et al., 2019). Note that since attraction manta ray behaviors were only observed in sites with SCUBA diving, models for freediving sites have fewer transition parameters (9 in free diving vs. 16 in SCUBA). This process generated a set of transition matrices for each site, separated by diver behavior and diving method. For oceanic manta rays: Revillagigedo active SCUBA divers, Revillagigedo passive SCUBA divers, Banderas Bay active SCUBA divers, Banderas Bay passive SCUBA divers, Banderas bay active free divers and Banderas Bay passive free divers. And for Caribbean manta rays: Isla Contoy active free divers and Isla Contoy passive free divers.
We calculated effect indices ( I ) in order to measure the paired difference between the effect of diver behaviors (passive or active in each site), comparable sites (SCUBA: Banderas Bay and Revillagigedo, freediving: Banderas Bay and Isla Contoy) and diving method (SCUBA or free diving in Banderas Bay). Each index had a range from −1 to 1, with positive values denoting a major effect in transition probabilities from the reference variable. This value shows which variable is more likely to change the behavior state of a manta ray for entries of the matrix representing state transitions (for example, transitioning from directional state at t time to erratic at t +1 ), or to make the mantas remain in those states for entries denoting the probability to remain in the same behavior state (the probability of transitioning from directional at t time toward directional state at t +1 ).
The indices were calculated as: A diver effect index estimated by subtracting the transition probability values of passive divers from those of active divers for each site, resulting in positive index values when manta rays displayed a higher transition probability with passive divers and vice versa.
A site effect index estimated for SCUBA sites by subtracting the values of Banderas Bay with passive and active divers, respectively, from their equivalent in Revillagigedo resulting in positive values when manta rays displayed a higher transition probability in Banderas Bay. Given the low number of interactions with passive free divers, we estimated a site effect index for free diving sites using only the active diver data set. This index compared the behavior of oceanic manta rays with active free divers against Caribbean manta rays with active free divers by subtracting the values of Banderas Bay with from those of Isla Contoy, resulting in positive index values when Banderas Bay manta rays displayed a higher transition probability.
A diving gear effect index, obtained by subtracting the values of Banderas Bay free diving with passive and active values, respectively, from their equivalent in Banderas Bay with SCUBA diving, with positive index values when manta rays displayed a higher transition probability with free divers.
Differences in diver effect indices were considered significant if the credible interval did not cross 0 (McElreath, 2015).

Quantitative Results
Passive observation was the most frequent diver behavior in Revillagigedo, while following was the most frequent in the other sites. Obstruction was the least frequent diver behavior (Table 4).
Behaviors units from the kinetic category (directional swimming, non-directional swimming and course change) were the most frequent in all sites, with the exception of Banderas Bay with active divers, where the defensive movements behavior unit was the most frequently observed ( Table 5).
For SCUBA sites, Revillagigedo oceanic manta rays did not display the avoidance somersault nor the defensive movement behaviors when divers were in passive observation, nor did they display cleaning movements when being chased. While oceanic manta rays in Banderas Bay did not display the non-directional cruising, avoidance somersault nor the defensive movement behaviors with divers in passive observation.
For free diving sites, oceanic manta rays in Banderas Bay, directional swimming behavior was the most frequent behavior, followed by defensive movements and directional cruising regardless of diver behavior. On the other hand, Caribbean manta rays in Isla Contoy National Park did not display the non-directional cruising, avoidance somersault nor the defensive movement behaviors with passive free divers, and were the only mantas to display feeding (Supplementary Figure 1) behaviors.

Markov Models
The following results are presented as a series of box plots blocks in order to ease interpretation. Each box represents the transition probabilities of each effect index (diver effect for Figure 2, diving gear effect for Figure 3 and site effect for Figures 4, 5). Each block contains four boxplots (of the transition probabilities) for SCUBA sites and three boxplots for freediving sites. The initial state is shown on the top of each block, while the state transitioned to is marked by the x-axis. The y-axis represents the value of each index, which as stated previously denotes which variable had a stronger effect over the change of behavior of the manta rays. Positive values denote a stronger effect of the reference value for each index. A higher value for the diver index (regardless of if it has positive or negative values) denotes a stronger effect of the variable. For example, a diver effect over the probability of a manta to transition from directional toward attraction behavior of 0.15 in one site and 0.50 in another would mean that passive divers (the reference value) had a higher effect over this change   Definitions of each behavior unit can be found in Table 1. (D) Isla Contoy freediving. Credibility intervals are shown as numbers below and above each box. The y-axis shows the diver effect index over the transition parameters using Passive divers as a reference value. Positive values relate to a higher effect of Passive diver behavior and negative values relate to a higher effect of Active diver behaviors. Each block contains four transition probabilities for SCUBA sites and three for the freediving sites; the initial state is shown on the top of each block, while the state transitioned to marks the x-axis. Transitions were the credible interval did not cross 0 are highlighted in blue for positive values and red for negative values. in behavior in both sites, but the effect was stronger in site 2. The credible interval, which is a dispersion parameter akin to the confidence interval (McElreath, 2015) is shown as decimal numerals above and below each box. Our analysis does not have a "significance" threshold per se, but the credible interval can be used to get a notion of the level of significance of each transition probability. With this in mind, we highlighted the boxes in which the credible interval did not cross 0. Transitions not highlighted could be likely to occur because of reasons other than diver behavior on each site.

Revillagigedo SCUBA
Oceanic manta ray transition probabilities were similar regarding diver behavior, with diver index values ranging from −0.15 to 0.15 (Figure 2A). Our models showed that oceanic manta rays interacting with active divers remained in directional behaviors (91% of posterior draws < 0) and erratic behaviors (92% of draws < 0), and changed their behavior from evasion to directional (72% of draws < 0) more often than with passive divers. Oceanic manta rays with passive divers on the other hand, did not remain in evasion behavior for long, changing to other behavior states such as attraction and erratic. Oceanic mantas also changed their behavior from attraction to directional more often with passive than with active divers (91% of draws > 0). Note that the probability of remaining in the evasion state was very similar between diver behaviors, with 53 % of the posterior draws being < 0 (active divers) and 43% of the posterior draws being > 0 (passive divers).

Banderas Bay
For SCUBA divers, oceanic manta ray behaviors did not show differences in the transition probability diver index greater than 0.01 between diver behaviors. Oceanic manta rays had a slightly higher probability of remaining in erratic behaviors with active SCUBA divers, while they changed from erratic to directional behaviors more often with passive SCUBA divers. Oceanic manta rays had a higher probability to remain in evasion behaviors with active divers (99% of posterior draws < 0). With passive divers on the other hand, oceanic manta rays were more likely to change their behavior from evasion to directional behaviors (97% of draws > 0) ( Figure 2B).
In contrast, only the evasion behavior showed differences according to free diver behavior ( Figure 2C). Oceanic manta rays interacting with passive free divers changed their behavior from evasion to directional (93% of draws) and to erratic (99% of draws) behaviors more often than with active free divers. Furthermore, oceanic manta rays with active free divers had a much higher probability of remaining in evasion behaviors (98% of the posterior draws).
For the comparison between free divers and SCUBA divers in Banderas Bay, since transition probability matrices for SCUBA divers had one additional state (attraction) we removed the transitions to and from attraction behaviors in order to compare FIGURE 5 | Site effect of each manta transition probability parameter with Divers in Active behaviors for Freediving sites. Credibility intervals are shown as numbers below and above each box. The y-axis shows the site effect index over the transition parameters using Banderas Bay freediving as a reference value. Positive values relate to a higher effect in Banderas Bay and negative values relate to a higher effect in Isla Contoy. Each block contains four transition probabilities; the initial state is shown on the top of each block, while the state transitioned to marks the x-axis. Transitions were the credible interval did not cross 0 are highlighted in blue for positive values and red for negative values. the two cases. We made a separate comparison for passive and active divers, shown in Figures 3A,B, respectively.
Oceanic manta ray behavioral changes with passive SCUBA and free divers were very similar ( Figure 3A). Credible intervals were wider for evasion due to the low number of observations of this state. In contrast, oceanic manta rays interacting with active free divers showed a higher probability to stay in evasion behaviors, and to change their behavior from the other states to evasion (Figure 3B). Regarding active SCUBA divers, oceanic mantas had a higher probability to change their behavior to directional and erratic states, even when transitioning from evasion.

Banderas Bay vs. Revillagigedo SCUBA
Oceanic manta rays showed different changes in behavior amongst the two sites. With passive divers (Figure 4A), oceanic manta rays in Banderas Bay changed their behavior from erratic to directional (99% of draws > 0) and evasion behaviors (98% of draws > 0) more often, while in Revillagigedo they showed higher probability of remaining in erratic (99% of draws < 0), evasion (98% of draws < 0) and attraction behaviors (87% of draws < 0).
For directional and erratic behaviors, oceanic mantas in both sites had a very similar response toward active divers, the site index showing a comparatively small range (from −0.1 to 0.1) with narrow credible intervals ( Figure 4B). Oceanic mantas did not show differences when changing their behavior from evasion to the other states between the two sites, the wider credible intervals were probably caused by a the relatively small number of instances of evasion behavior observed with SCUBA divers, with differences in the transition probability site index of up to 0.15 in some posterior draws. However, oceanic mantas did show a higher probability to change their behavior from directional, erratic and attraction behavior states to directional states with active divers in Banderas Bay.

Isla Contoy National Park Freediving
Caribbean manta rays displayed small changes in their transition probabilities from directional and erratic behaviors according to diver behavior, although they had a higher probability to remain in the directional and erratic behaviors with passive divers (Figure 2D). In contrast, Caribbean manta rays showed differences when changing their behavior from evasion toward other states. Caribbean manta rays with passive free divers showed an increased probability changing their behavior from evasion to directional (88% of draws > 0) and erratic (98 % of draws > 0). Caribbean manta rays tended to remain in the evasion behavior state with active divers (96% of draws < 0).
When comparing the transition probabilities from oceanic and Caribbean manta rays from Banderas Bay and Isla Contoy National Park, respectively, the results were similar to the comparison between active SCUBA divers. The range of the site effect index was relatively small (from −0.15 to 0.15) and credible intervals were relatively wide except when directional behaviors were the initial behaviors due to the relatively larger sample size of directional behaviors (Figure 5). Caribbean manta rays stayed in directional behaviors and changed their behavior to the directional state (from erratic 99% of draws < 0, and from evasion 89% of draws) than those from Banderas Bay, which in turn displayed higher probabilities of remaining in erratic (99% of draws > 0) and evasion behaviors (99% of draws), as well as a higher chance to change their behavior to evasion.

DISCUSSION
Tourism can impact targeted populations of marine wildlife negatively by disturbing their normal activities, stressing the animals and even causing physical damage (Higham, 2007). Nevertheless, if regulations and management plans address the possible effects caused by it, tourism can be a powerful tool for the management and sustainable use of animal resources (Buckley, 2011). Long living, late maturing and highly exploited marine species could be specially benefited from informed regulations over tourism activities and the use of tools such as best practices manuals, developed by the application of scientific knowledge and research (Lawrence et al., 2016). During the last decades, manta ray tourism has risen as a promising alternative to consumptive an extraction activities of these species (O'Malley et al., 2013). Sites such as Mozambique in which surveys have evaluated economic spillover and perceived value of manta ray tourism have shown an increase awareness and interest in the protection of the species by local residents (Venables et al., 2016). Despite the promotion given to manta ray tourism and the efforts of several organizations for the protection and study of manta rays (Murray et al., 2020), some of their populations are still declining due to exploitation, illegal fishing, slow growth rate and the lack of ecological knowledge about these species (Dulvy et al., 2014). The continual decline of the less studied oceanic manta rays recently caused the species to be moved from "vulnerable" to "endangered" by the IUCN's Red list of Threatened Species (Marshall et al., 2020). There are important gaps in our knowledge about all species of manta rays (particularly the recently described Caribbean manta), and is up to us to use the means we have at hand to keep trying to fill this gaps in order to improve our management of this species.
With this in mind, we analyzed a larger data set consisting of footage of manta rays collected from research surveys and tourism dives. We coupled technology with Bayesian models allow us to describe, analyze and compare the behavioral response of the manta rays toward divers for the first time. Our results confirm the beneficial effect (such as longer interactions and less diver avoidance) of following a calmed or passive approach over diver-animal interactions and shed light over the different factors that can lead to encounters in which the diver can interact with a manta ray without causing it to quickly leave the area. This kind of encounters would in turn benefit local tourism service providers by improving tourist visitation turnover and facilitating promotion by local research and conservation organizations.
Our results should be interpreted with the caveat that we only measured manta ray behaviors and responses during relatively short encounters with divers. While our results can be used for comparisons to other manta ray studies, they do not reflect how the animals would behave in the absence of divers, during interactions with other animals or with other anthropogenic stimuli such as boat traffic.

Divers
Although the passive and active behaviors were displayed in a similar manner by divers across all sites, site specific conditions affected the frequency of particular behaviors in some of them.
The predictability and high probability of encounters with oceanic manta rays in Revillagigedo highly contributes to the low incidence of following behavior by the divers to the oceanic manta rays (less than 30% of the encounters). In contrast, Banderas Bay SCUBA diving encounters with oceanic manta rays are opportunistic and infrequent. Banderas Bay's SCUBA diving activities focus on observation of coral and rocky reefs, and tourists are not typically briefed on how to act around oceanic manta rays. As a result, when an oceanic manta rays are seen during a dive, divers generally swim actively toward them to prolong the encounter (around 50% of our recorded encounters). Active free divers in Bahia Banderas and Isla Contoy National Park acted in the same manner, and with the same goals (scientific data gathering) in both study sites. These active behaviors were very similar to what would be expected from eager and uninformed free divers during tourism trips. The few instances of passive free divers were very valuable for this study, given the observed reduction in transitions toward negative behaviors when divers remained passive, even with freediving gear.
It is important to note that the diver behaviors observed in this study do not encompass all the possible behaviors that a diver may display during an interaction with manta rays. Divers in other aggregation sites have been observed swimming toward manta rays from the front (Murray et al., 2020), sitting on the bottom during natural aggregations (Perryman et al., 2021) or during aggregations caused by using underwater lighting (Osada, 2010). Note that these examples as well as the majority of research about tourism with manta rays has been made with reef manta rays, who tend to form larger aggregations during longer encounters than oceanic manta rays.
Several factors may favor the use of certain diving techniques over others. SCUBA diving may improve interactions in deep sites with good visibility such as Revillagigedo (Hull et al., 2006) or some parts of Hawaii (Osada, 2010). However, SCUBA interactions should not alter the functionality of aggregation sites by, for example, dispersing cleaning fish or damaging coral reefs due to poor buoyancy control (Murray et al., 2020). In these situations, and with the proper site-specific codes of conduct, free diving has the potential to be a viable option, especially considering that passive free divers cause similarly low levels of disturbance in manta rays according to our results. Free diving also has the added benefit of requiring less time for diver set up, thus increasing the potential duration of diving trips and tourist visitation turnover. Nevertheless, free diving is a viable option only in sites were the manta rays can be encountered swimming in the surface.

Effect of Environmental and Site-Specific Conditions
The comparison between sites with the same diving method (SCUBA and freediving) under similar tourism conditions (passive or active divers) highlighted some differences that could be attributed to the environmental conditions that we did not measured (such as visibility, temperature and food availability). As mammals, we struggle to recognize or understand the responses of marine animals to their environment (Bshary et al., 2002), but we are capable of identifying general behavior patterns, isolating the highest possible number of unknown variables in order to determine the responses of animals toward specific stimuli (Altman, 1974). This process can become overly complex when we consider highly social and mobile animals which cannot be subjected to traditional experimental designs such as manta rays.
Manta rays can respond to visual and environmental cues (Ari and Correia, 2008) and are known to display specific behaviors under certain conditions such as social feeding (Perryman et al., 2021) and reproductive aggregations (Stevens et al., 2018). These social and environmental variables can increase or decrease the tolerance of manta ray toward tourism according to the level of trade-off between disturbance and benefit that certain areas may provide (Barr and Abelson, 2019). The effect of other environmental factors such as visibility can be included in future studies in order to determine if evasive behaviors are more prominent under certain conditions. Turbid waters for example, can cause an increased risk in predation, while temperature, which can affect the metabolism and level of activity in elasmobranchs (Carrier et al., 2012) can affect the degree in which mantas react to the presence of divers.
The goal of this study was to measure the effect of the different diver behaviors over manta rays by using interaction data gathered from different aggregation sites. We observed clear differences that can be attributed to a number of environmental factors both known (such as visibility and habitat use) and unknown, but we also observed equally clear differences in the responses of manta rays toward specific diver behaviors. By considering these factors, future studies may build upon the bases provided by this study in order to better understand the reasons driving the behavioral responses of manta rays.
These site differences should be taken in consideration when developing tourism management plans for this as well as other manta ray aggregation sites under similar conditions. We cannot control the natural variability of environments in which mantas occur, but we can change, hopefully for the better, the way in which divers engage with manta rays.

Manta Rays
Directional behaviors are simple displacement toward or away an unknown stimulus (Altman, 1974) while taxis (positive and negative) involves clearly defined triggers (in this case, divers).
However, the interpretation of erratic behavior frequencies has to be viewed with caution. Studies have shown that environmental factors may indirectly affect the decision process of manta rays to remain in cerain sites (Barr and Abelson, 2019), it is possible that a similar trade-off effect could be influencing the frequency of erratic behaviors. At SCUBA diving sites the relatively high frequency of erratic behaviors may be linked to an adaptive response to remain longer in the same area. It is possible that the frequency of erratic behaviors observed at SCUBA sites could be due to potential foraging opportunities or the presence of nearby cleaning stations (Stewart et al., 2016b(Stewart et al., , 2019. Due to the nature of our data we cannot determine the exact length of an encounter. However, since mobile marine animals in general tend to get away from tourists when they are disturbed (Lawrence et al., 2016), we can infer that a higher transition probability toward evasion would most likely shorter encounters. Likewise, transitions toward attraction and erratic behaviors could increase the length of encounters since mantas in these behaviors are not moving away from divers nor swimming toward a specific destination.
We did not analyze the distance, speed or direction of active divers while swimming after the manta rays, nevertheless negative manta ray behaviors were uncommon or absent in all sites while divers were passive. This supports the results observed in reef manta rays by Murray et al. (2020) who, regardless of the fact that they measured the effect of different free diver approaching methods in much denser reef manta ray aggregations, also discourage following behaviors, finding similar disturbances to the observed in this study.
We observed a similar response from Caribbean and oceanic manta rays when analyzing these sites separately. Both species showed a higher probability to display evasive behaviors when interacting with active divers and a higher probability to switch from evasive toward erratic movements with passive divers. When comparing the responses of both species toward active divers directly via the site index, oceanic manta rays displayed a higher probability to stay evasive and to transition to evasive behaviors, while Caribbean manta rays displayed a higher probability of staying in directional behaviors and transitioning to them from evasive, erratic and positive behaviors. These differences may be attributed to the fact that Caribbean mantas were feeding during the encounters. Recent studies have shown that food availability affects the behavior and decision-making process of manta rays (Barr and Abelson, 2019), thus increasing the tolerance of feeding Caribbean mantas to free divers following them. Further studies can be developed in order to determine if these observations are caused by behavioral differences between the two species of manta rays. Additional studies with similar methodologies can also be applied in other Caribbean and oceanic manta ray aggregation sites. These studies should aim to compare the behavior of oceanic and Caribbean manta rays in similar habitats, preferably in absence of divers.

Revillagigedo
Although oceanic manta rays at Revillagigedo showed and increased probability of abandoning negative behaviors with passive divers, they displayed a similar proportion of attraction behaviors regardless of diver behavior. As mentioned before, Revillagigedo oceanic manta ray encounters occur in relatively deep sites with nearby cleaning stations (Hull et al., 2006) and in proximity to deep areas where feeding has been observed (Stewart et al., 2016b). Nevertheless, tourism encounters in this study were recorded in sites where divers were always suspended in the water column (and therefore free to swim after the oceanic manta rays or remain passive) and where oceanic manta rays were neither feeding nor being cleaned. Manta rays in the site have been subjected to tourism for decades, which along with the potentially long term memory of the manta rays suggested by other studies (Ari and Correia, 2008;Barr and Abelson, 2019), may be causing some degree of general habituation to divers in Revillagigedo. The remote nature of the archipelago along with the aforementioned factors are likely to contribute to the presence of attraction behaviors and the reduce frequency of evasion behaviors. However, it is important to note that oceanic manta rays had a higher probability to remain in directional behaviors which could imply a reduced level of engagement with active divers.

Banderas Bay
In contrast, oceanic manta rays in Banderas Bay tend to end the encounter quickly or transition to directional behaviors. These animals are exposed to constant, heavy marine traffic and noise disturbance due to their proximity to the important tourism and commercial town of Puerto Vallarta (Everitt et al., 2008). Additionally, while feeding may occur during the period that oceanic manta rays remain in the bay (Stewart et al., 2019) feeding events were never recorded for oceanic manta rays in Banderas Bay in this study (Supplementary Figure 1). Recent studies suggest that the site is used for resting or "basking" purposes (Fonseca-Ponce et al., 2021) by the manta rays, which could be the reason why they seem to react more abruptly toward active divers. The constant presence of stressors and the lack of a favorable resource that undermines the possible disturbance caused during interactions with divers (such as food and cleaning stations) may be contributing to the display of evasion behaviors, and the higher transition probabilities toward evasion behaviors.
Our model results suggest that oceanic manta rays in Banderas Bay stayed longer in the evasion behavior category in presence of active divers. On the other hand, while divers remained passive, the oceanic manta rays showed a higher probability of returning to directional, attraction and erratic behaviors. Both results demonstrate the importance of advising divers and tourism service providers not to swim after oceanic manta rays if they encounter them. Based on our results, swimming after oceanic manta rays is not only likely to prompt an evasive response in the animals, it also appears likely to shorten or end an encounter, which is counterproductive from a tourism perspective. Since encounters with oceanic manta rays are hard to predict in Banderas Bay, providing divers with preemptive briefings on best practices would most likely improve the opportunistic interactions. Instructing divers on how to behave during an oceanic manta ray encounter would not only increase the likelihood of the oceanic manta ray to remain in the area, but would also reduce the possible disturbance caused by the divers. In the long term, this could prompt tourist service providers to encourage best practices amongst each other, encouraged by the longer and more engaging nature of encounters between passive divers and oceanic manta rays.

Diving Gear Comparison
Since oceanic manta rays in Banderas Bay encountered by free divers and SCUBA divers presumably belong to the same population, this site provides a unique opportunity to evaluate whether manta rays respond differently to free divers and SCUBA divers. Oceanic manta rays appear to display evasion behaviors more often in the presence of free divers than in presence of SCUBA divers when the divers are active. However, oceanic manta rays showed a similar response to passive divers regardless of diving gear, which again, highlights the importance of following interaction guidelines. These observations should be taken in consideration for the development of tourism encounters with manta rays. If best practices manuals are developed and properly followed, encounters between manta rays and divers could have minimal behavioral impacts even in places where SCUBA diving is not possible.

Isla Contoy National Park
Caribbean manta rays in Isla Contoy National Park routinely feed at the surface in the study area, oceanic manta rays in Banderas Bay do not. Much of the data from Isla Contoy was obtained while Caribbean manta rays were clearly feeding even in the presence of free divers in the area. Interestingly our model results were nearly identical for both sites. This shows that evasion behaviors, while less common in Isla Contoy, tend to be displayed in a similar manner to that of the Banderas Bay manta rays.
Differences came up only when comparing the models for active divers from both sites. Oceanic manta rays in Banderas Bay are more likely to avoid active divers than those from Contoy Island. Although environmental factors such as marine traffic (for Banderas Bay) and the presence of food (for Isla Contoy) surely have an effect over manta ray behavior in both sites. Our results suggest that if evasion behaviors are being displayed by manta rays, they are more likely to keep avoiding active divers, but they would change their behavior (possibly resuming their previous activities) if divers are passive. This shows that the behaviors of divers do affect, at least partially, the behavior of manta rays in both sites.
Despite similar conditions in visibility and interaction with free divers, the site specific differences such as foraging opportunities, a higher presence of divers (e.g., due to whale shark tourism) in Isla Contoy National Park (Graham et al., 2012) and the aforementioned marine traffic in Banderas Bay as well as the possible effects speciation should be further analyzed. Diver behavior can have a lesser of higher effect over mantas according to site-specific characteristics that should be taken into account when designing tourism management plans for these and other manta aggregation sites.

CONCLUSION AND RECOMMENDATIONS
According to our findings, manta rays display different responses toward active and passive divers. Active diver behaviors caused manta rays to avoid divers or end encounters early with more frequency both in SCUBA and freediving encounters. The variations in duration and change of manta ray behaviors observed amongst sites may have been caused by a combination of factors. Site specific conditions such as proximity from the coast and feeding events, habitat use and possible habituation to divers are some likely causes. Therefore, we suggest that future studies, supported by additional technologies such as unmanned aerial vehicles or submarine cameras, should measure the base behavior as well as social interactions of undisturbed manta rays (i.e., in absence of divers). Additional studies are necessary to elucidate how the effects of diver behaviors intersect with (i) presence of other stressors (such as boat traffic and habitat degradation), (ii) environmental conditions, and (iii) occurrence of cleaning stations, to determine which sites are more favorable or not for ecotourism with manta rays. Another interesting line of research would be to address the individual differences in habituation to divers, and the effect of the number of divers present during interactions. Our study, the first of its kind in oceanic and Caribbean manta rays, underlines that regulations and the use of best practices are vital for achieving positive encounters for both manta rays and divers. By using and encouraging proper practices for the interaction with manta rays, tourism service providers can archive more lasting and memorable encounter and open the possibility for the development of sustainable tourism activities in areas were SCUBA or freediving is viable.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
MG-G: project conceptualization, data gathering, data structuring, coding and analysis, and manuscript writing and editing. MB-M: data review and manuscript reviewer. JS and JK: data contributor, data review, and manuscript reviewer. VL-B: data structuring, coding and analysis, statistical advisor, and manuscript reviewer. IF-P: aid in project conceptualization, Banderas Bay data contributor and data gathering, and manuscript reviewer. AZ-J: aid in project conceptualization, Banderas Bay data contributor and data gathering, design of study site Map, and manuscript reviewer. KF: Contoy Island National park data contributor and data gathering, and manuscript reviewer. All authors contributed to the article and approved the submitted version.

FUNDING
This project was part of a Master's Thesis at CIBNOR, and was supported by a CONACYT scholarship, number 715334.