Reef manta rays forage on tidally driven, high density zooplankton patches in Hanifaru Bay, Maldives

Study site

The Maldives has a large resident population of M. alfredi which undertake biannual migrations linked to the changing monsoons (seasons) within the archipelago (Anderson, Adam & Goes, 2011; Fig. 1A). During the Southwest Monsoon, or Hulhangu (April–November), numerous M. alfredi frequent foraging aggregation sites on the eastern side of the nation’s atolls (Harris et al., 2020). One site, Hanifaru Bay, is situated on the eastern edge of Baa Atoll, and attracts large feeding aggregations of this species annually (Stevens, 2016; Harris et al., 2020). Hanifaru Bay is a small reef inlet (700 m long by 200 m wide) which forms part of a core marine protected area within the Baa Atoll UNESCO Biosphere Reserve (5°17′N, 73°15′E; Fig. 1B). The shallow (maximum depth 22 m) inlet is periodically inundated with zooplankton-rich water. Motorised boat activity and SCUBA diving are prohibited in Hanifaru Bay due to the high numbers of manta rays and other megafauna that access the inlet (Murray et al., 2020).

Data collection

Fieldwork was conducted in the lead up to the new Moon in August 2017, under Ministry of Fisheries Permit No. (OTHR)30-D/PRIV/2017/280, and Ministry of Environment Permit No’s. EPA/2017/RP-01 & EPA/2016/PSR-M02. This time of year was chosen because strong lunar tides appear to overcome the force of the prevailing monsoonal current, drawing plankton-rich water from outside the atoll edge into Hanifaru Bay (Harris & Stevens, 2021). The currents form a back eddy, trapping and concentrating plankton in this shallow reef inlet, resulting in M. alfredi foraging opportunities, which peak during spring and high tides (Stevens, 2016; Harris & Stevens, 2021). Sampling was conducted during daylight hours and across the tidal cycle from 13–21 August 2017. Zooplankton was collected by two people using a 200 µm-mesh net with a 50 cm diameter mouth. The net was towed by hand at the surface for a ∼50 m transect between two coral features at the eastern end of Hanifaru Bay (Fig. 2A), at a speed of ∼22 m per minute. A flowmeter was fitted to the plankton net to allow calculation of the volume of water sampled. Flowmeter calibration was performed prior to the field trip in a swimming pool of known length to establish an accurate measurement of distance per flowmeter revolution. Samples were kept on ice and fixed with 10% buffered formalin solution at the end of each day.

Each net tow was accompanied by an in-water observer recording manta ray activity in vicinity of the tow. This included: (1) manta ray abundance; (2) behaviour (Feeding, Non-feeding—when manta rays were present but not feeding, and Absent—when manta rays were not present); and (3) most common feeding strategy employed (as described in (Stevens, 2016; Fig. 2B).

Temperature and salinity data were collected at 1 s intervals from a CTD unit deployed at the site for the study duration (except for ∼24 hrs from 17–18 August for battery exchange). Temperature and salinity ranges were relatively small throughout the study (28.6–29.2 °C and 34.3–35.0 ppt respectively). These variables were excluded from the models as their inclusion resulted in missing values. Local tide data were obtained from a government representative from the Ministry of Environment.

Sample processing

Zooplankton samples were processed in the CSIRO Plankton Ecology Laboratory in Brisbane, Australia. Flowmeter readings and the area of the net mouth were used to estimate the volume of water filtered.

Zooplankton biomass

Zooplankton samples were split into two halves using a Folsom splitter (Harris et al, 2000). The first half was used to determine dry mass, with each sample oven-dried at 70 °C for 24 hrs prior to weighing. Zooplankton dry mass (hereafter referred to as biomass) per unit volume of filter-seawater for each tow was calculated by dividing the dry mass of the sample (mg) by the volume of filtered water (m³): ${\displaystyle Biomass\left(mg{m}^{-3}\right)=Drymass\left(mg\right)/Volumeofwaterfiltered\left(m^3\right).}$

Zooplankton identification

The second half of the sample was used to examine size structure and community composition via a 2400 dpi ZooScan system and microscopy. The Hydroptic v3 ZooScan (EPSON Perfection V700 Flatbed) is a high resolution, waterproof scanner that digitises particles for size and biovolume measurements (Gorsky et al., 2010). An aliquot of each sample was prepared using a Stemple pipette of known volume and placed on the scanning tray. To avoid overlap, particles were manually separated. Once separated, the sample was scanned and particles were extracted into vignettes for categorisation into broad taxonomic groups (24 groups) using Plankton ID software (Version 1.2.6) and manual validation (Gorsky et al., 2010). Objects classified as sand, fibre, detritus, bubbles and shadows were excluded from further analysis (as per (Rohner et al., 2015). For visualisation, taxa that comprised <5% of the total abundance were grouped as “other”, and these included cnidaria, polychaetes, echinoderm larvae, bryozoan larvae, fish larvae, salps, and various classes of arthropods.

To investigate which species were responsible for the majority of the biomass at the site when overall biomass values in the water were high, samples were analysed taxonomically via microscopy. A subsample was prepared using a Stempel pipette, and organisms were identified and counted in a Bogorov tray using a microscope. Dominant members were identified to genus or species with assistance from trained plankton taxonomists at the CSIRO Plankton Ecology Laboratory (Eriksen et al., 2019).

Zooplankton size structure

A size distribution of the sample particles, known as a Normalised Biomass Size Spectra, was produced to analyse the size structure of the zooplankton community (Vandromme et al., 2012). Spherical biovolume was calculated from the size measurements obtained from ZooScan. Each particle was assigned to one of 50 logarithmic size categories based on its spherical biovolume. The sum of the spherical biovolume of the particles in each size class (mm³) was standardised by the fraction of sample scanned and the volume of water filtered (m³), and normalised by dividing this value by the width of the size class measured in biovolume (mm³). Both axes of the Normalised Biomass Size Spectra use a logarithmic scale.

Drivers of zooplankton biomass and manta numbers

To investigate potential drivers of zooplankton abundance and M. alfredi visits to Hanifaru Bay, we constructed generalised linear models (GLMs) using R (R Core Team, 2019). Separate analyses were conducted for two response variables: (i) Zooplankton biomass (mg m⁻³), with a Gamma error structure and log-link function; and (ii) Manta ray abundance (number of M. alfredi observed during zooplankton sampling), with a negative binomial error structure and log-link function (Poisson error structure was overdispersed). We visually inspected diagnostic plots to assess assumptions of homogeneity of variance and normality. Independent variables in both models were Tide (hours from high tide) and Behaviour (Feeding, Non-feeding and Absent). To account for the circular nature of Tide (∼12-hr cycle), the variable was transformed using a truncated Fourier series (a harmonic function of sines and cosines). This ensures that the cyclical nature of this variable is captured, while guaranteeing that the response values predicted at the extremes of the variable range are the same (i.e., the same prediction for Zooplankton biomass or Count at times of 0 and 24 h). For the Manta ray abundance model, the Behaviour variable was reduced to two categories—namely Feeding and Non-feeding, and Zooplankton biomass (mg m⁻³) was included as an independent variable. Models were plotted on the response scale using the package “visreg” in R (Breheny & Burchett, 2017).

Critical thresholds for feeding behaviour and strategy

We assessed whether there might be a critical threshold for M. alfredi feeding—i.e., a level above which the likelihood of feeding increases dramatically. We thus used a GLM with a binomial error structure to analyse manta ray behavioural response (Non-Feeding = 0, Feeding = 1) in relation to zooplankton biomass (mg m⁻³) as a predictor. The critical density threshold was taken as the zooplankton biomass at which the proportion of feeding was 0.5. A theoretical prey density threshold was plotted for comparison, based on findings by Armstrong et al. (2016). Their study assumed general morphometrics (average disc width of 3.5 m, mouth opening of 0.3 m, and mass of 100 kg) and swim speeds (2 knots when feeding) for M. alfredi. These assumptions are also applicable in the current study.

Feeding samples were categorised into either Solo feeding (Straight, Surface and Somersault) or Group feeding (Piggy-back and Chain) based on the most common strategy observed in the manta rays (Stevens, 2016). A GLM with a binomial error structure was used to analyse manta ray feeding strategy response (Solo = 0, Group = 1) in relation to zooplankton biomass (mg m⁻³) as a predictor. The critical density threshold was taken as the zooplankton biomass at which the proportion of Group feeding was 0.5.

Zooplankton community analysis

To determine how different the zooplankton communities were for the M. alfredi behaviours (Feeding or Non-Feeding), non-metric multidimensional scaling was used based on abundance counts of the different taxonomic groups from the Zooscan analysis. The Bray Curtis distance measure was used because it is unaffected by joint absences of taxonomic groups in samples. To account for abundance of certain taxa, data were transformed using a root transformation. To test for differences in community composition between M. alfredi behaviours (Feeding and Non-Feeding), we performed an adonis analysis, a multivariate analysis of variance. Both the adonis and non-metric multidimensional scaling were conducted using the “vegan” package in R (Oksanen et al., 2007).

Zooplankton community composition

There were differences in the zooplankton community composition between Feeding and Non-feeding samples when analysed using non-metric multidimensional scaling on the Zooscan taxonomic counts (Fig. 6A). The 95% confidence ellipses for Feeding and Non-feeding were not overlapping, implying that they were significantly different zooplankton community compositions, and this was confirmed by the adonis analysis (F = 9.42, df = 1,53, p = 0.001). Crustaceans (such as copepods) were more associated with Feeding samples, compared to gelatinous taxa (such as chaetognaths and eggs), which were more associated with Non-feeding samples.

Calanoid copepods comprised 66.3% of Feeding samples compared to 46.7% of Non-feeding samples (Figs. 6C and 6D respectively). Chaetognaths were 5.3% of Feeding samples, and 11.9% of Non-feeding samples. Fish eggs were less than 2% of Feeding samples, and 13.0% of Non-feeding samples. Based on microscopy, Undinula vulgaris (juveniles and adults) was the dominant calanoid copepod species in both Feeding and Non-feeding samples (25.0% and 30.7% respectively, Fig. 6B).

Zooplankton size structure

Analysis of the size structure of zooplankton from Hanifaru Bay revealed that the biovolume of zooplankton increased in the majority of size categories when M. alfredi was feeding (Fig. 7). The biovolume of zooplankton was significantly higher across particle size categories during M. alfredi Feeding events than Non-feeding events (Mean total standardised biovolume: Feeding = 288.4 and Non-feeding = 172.1; t = −2.66, df = 51.38, p = 0.01). Feeding and Non-feeding samples had similar biovolumes of small and large particles, but Feeding had significantly more moderate-sized particles (from 10^−1.2 to 10^0.5 mm³).

Summary

Zooplankton concentrations influence the number of M. alfredi present and their observed behaviour in Hanifaru Bay. Rapid changes in zooplankton are observed across the tidal cycle, and M. alfredi feed when biomass reaches a critical density which is higher than predicted to meet their theoretical metabolic requirements. Mobula alfredi foraging occurs when the zooplankton community is dominated by calanoid copepods, and feeding is less likely when there are greater numbers of gelatinous taxa (such as chaetognaths or eggs). Taxonomic analysis reveals that the large-bodied copepod, Undinula vulgaris, dominates the zooplankton environment at Hanifaru Bay, suggesting oceanic incursions may play an important role in bringing zooplankton to this small reef inlet.

Tidal influence on zooplankton density and manta ray foraging

Manta rays feed when zooplankton biomass is high, which is typically observed on the high to ebbing tide at Hanifaru Bay. Oceanographic investigations in Hanifaru Bay suggest tidal currents draw zooplankton into the shallow reef systems of the atoll, where they become trapped inside due to a back-eddy mechanism created by the unique shape of the reef system and the combination of the lunar and monsoon currents (Hosegood pers comms). Tides are known to influence the distribution and abundance of zooplankton around island inlets in the Great Barrier Reef (Alldredge & Hamner, 1980), and have been shown to influence manta ray feeding behaviour at aggregation sites in Indonesia (Dewar et al., 2008), eastern Australia (Armstrong et al., 2016), and the Chagos Archipelago (Harris et al., 2021). Therefore, short-term in situ observations of zooplankton concentrations in relation to tidal cycles and manta ray behaviour can help inform when M. alfredi is likely to be observed in Hanifaru Bay.

Animal movements and productivity hotspots

Large planktivores seeking to exploit ephemeral food sources in surface waters are likely to respond to currents and water movements that concentrate zooplankton. Cetorhinus maximus forage along thermal fronts (Sims & Quayle, 1998), R. typus targets regions of upwelling (Ryan et al., 2017), and surface foraging in M. alfredi is often tidally driven (Dewar et al., 2008; Armstrong et al., 2016). In conjunction with responding to physical oceanographic cues, animals that seek patchily distributed prey sources are also likely to congregate in areas where they have previously encountered energetically rewarding prey abundances, resulting in larger numbers of animals in reliable foraging regions. For example, M. alfredi predictably switches to the down-current side of the atolls in the Maldives in response to monsoonal winds and primary productivity (Harris et al., 2020). Area-restricted search theory predicts that animals will remain localised in areas where they have a higher probability of encountering prey (Bailey et al., 2019), and this perhaps explains why some M. alfredi individuals remain in Hanifaru Bay when not feeding. This location also has two cleaning stations used by M. alfredi (Stevens, 2016), and it is hypothesised manta rays will frequent cleaning stations in close proximity to foraging opportunities (Armstrong et al., 2021). Manta rays in Palmyra Atoll used area-restricted searching when adjacent to ledges or channels with high plankton concentrations, but their movements were more random at larger spatial scales (Papastamatiou, De Salles & McCauley, 2012). Area-restricted searching has also been observed in two dolphin species (Tursiops truncatus and Delphinus delphis) in areas of high prey availability, and where they have had previous successful foraging experience, suggesting memory plays a role in their movement ecology (Bailey et al., 2019). For M. alfredi, the apparent preference for returning to the same cleaning stations over time (Armstrong et al., 2021), suggests they may form a cognitive map of shallow reef environments, and this is likely the case for known productivity hotspots as well.

High critical feeding threshold for manta rays at Hanifaru Bay

The critical prey density threshold for M. alfredi feeding in Hanifaru Bay (53.7 mg m⁻³) is more than four times higher than that in east Australia where M. alfredi feeds (11.2 mg m⁻³; Armstrong et al., 2016), and in east Africa where R. typus feeds (12.4 mg m⁻³; Rohner et al., 2015). It is also double the theoretical prey density threshold calculated to meet the metabolic requirements for M. alfredi (25.2 mg m⁻³; Armstrong et al., 2016), which may explain why this site hosts such a large feeding aggregation of this species. However, these large planktivorous elasmobranchs are assumed to feed in the mesopelagic layer (Couturier et al., 2013; Burgess et al., 2016), so an understanding of the prey densities available at these depths is required to gauge the relative importance of aggregations sites such as Hanifaru Bay in meeting these species’ daily energetic requirements. Sampling zooplankton at depth remains a logistical challenge for researchers, but with technological advances, such as satellite tags equipped with accelerometer data loggers, and unmanned video submersibles (Stewart et al., 2018), these inferences can be better investigated.

Manta ray feeding strategies

In the current study, plasticity in M. alfredi feeding strategies in response to changes in prey biomass in Hanifaru Bay is not supported. This contrasts with work previously conducted in Hanifaru Bay that found manta rays were significantly more likely to employ group feeding strategies as prey density increased (Stevens, 2016). The previous work was based on a qualitative visual index for prey density, with data obtained over a long time period (>5 years) and included aggregations upwards of 150 animals. However, we did see that when zooplankton biomass values were very high, over 200 mg m⁻³, that only Group feeding strategies were used, and no Solo feeding was seen. But in either scenario, it is uncertain whether the observations are due to true cooperative feeding strategies, or simply that coordinate movements reduce collisions with other manta rays (Stevens, 2016).

Zooplankton composition and size

Differences in the composition of the zooplankton community were observed between M. alfredi Feeding and Non-feeding events, and M. alfredi was observed feeding when the overall biovolume of zooplankton was greater. Calanoid copepods dominate the zooplankton community for manta rays at Hanifaru Bay during the time of the study. This is similar to findings for M. birostris at Isla de la Plata, Ecuador, where zooplankton samples collected during feeding and non-feeding events were dominated by calanoid copepods, but also cyclopoid copepods (Burgess, 2017). The diet of mobulids (M. japonica, M. thurstoni, M. tarapacana, and M. birostris) in the Philippines is more diverse, with euphausiid krill dominating, but also records from stomach content analysis of squid, fish and copepods (Rohner et al., 2017). Investigations into foraging aggregations of C. maximus off Plymouth in the UK have also revealed the dominance of copepods during feeding events, with threefold increases in this taxa reported (Sims, 1999). In contrast, the food environment at R. typus feeding aggregations appear to be more diverse, with R. typus feeding on sergestid shrimps, fish spawn, calanoid copepods, and small bait fish (Motta et al., 2010; Fox et al., 2013; Rohner et al., 2015). These findings suggest that the food environment for large planktivores feeding in predominately tropical or sub-tropical waters is varied, and is likely driven by pulses in productivity, rather than particular taxa.

Undinula vulgaris is the most prominent species observed in zooplankton samples from Hanifaru Bay. It is a key species in tropical areas due to its large size and tendency to swarm in high numbers, making it a good food resource for planktivorous fishes (Alvarez-Cadena, Suárez-Morales & Gasca, 1998). This species has been observed at numerous large planktivore feeding aggregation sites, including those visited by M. alfredi in eastern Australia (Couturier et al., 2013; Armstrong et al., 2016), R. typus in the Gulf of Tadjoura, Djibouti (Boldrocchi et al., 2018), both M. birostris and R. typus in the Gulf of California (Notarbartolo-di Sciara, 1988; Lavaniegos, Heckel & De Guevara, 2012), and both M. alfredi and R. typus in the Philippines (Canencia & Metillo, 2013; Yap-Dejeto et al., 2018). Undinula vulgaris is considered an indicator of the influence of neritic-oceanic waters in reef environments, and its local distribution can suggest oceanic water sources (Alvarez-Cadena, Suárez-Morales & Gasca, 1998). Further investigation into the ecology of U. vulgaris in tropical environments may aid our understanding of how vital swarms of this species are for supporting large tropical planktivores, and whether their distribution and abundance is likely to be impacted by a rapidly changing climate.

Caveats to zooplankton sampling

There were three main limitations in the sampling design of the current study. First, as boat engines are prohibited in the area, towing of the plankton was done by hand and this resulted in tow speeds slower than recommended (22 m vs 60 m per minute) for minimising capture-avoidance by larger zooplankton (Harris et al, 2000). However, the use of a net with a relatively large mouth area (50 cm diameter) would help mitigate this. Further, the use of slower towing speeds (Davies & Beckley, 2010) and a relatively large mesh size for tropical waters (200 µm; Eriksen et al., 2019) should reduce the pressure wave from the towed net, and thus reduce concerns about net avoidance by zooplankton. Second, zooplankton sampling was limited to surface waters. This could influence the findings, particularly in relation to investigating the zooplankton dynamics when manta rays are employing different feeding strategies. For example, somersault feeding and cyclone feeding strategies do not necessarily occur in surface waters. This issue could be overcome with the use of drop nets, vertical free-fall nets that sample on the way down with weighted rings for propulsion (Eriksen et al., 2019), to provide coverage throughout the water column. Lastly, it is likely our relatively short sampling duration failed to detect group feeding dynamics, and our results could suffer from small sample size. More work needs to be done to assess whether the presence of higher zooplankton biomass is positively correlated with group feeding events. In particular, higher sampling replication of the eight different feeding strategies is required to tease apart how these relate to zooplankton dynamics. One of the key issues for sampling design is that the magnitude and exact timing of the biggest feeding events each season is somewhat variable (Harris & Stevens, 2021). This highlights the need for studying multiple feeding events in a consistent way to distil the general pattern. To further examine the zooplankton and manta ray foraging dynamics at this site, we suggest a future sample size at least double the current study, and attempts should be made to sample larger feeding aggregations (50+ manta rays).