Intraspecific differences in short-and long-term foraging strategies of reef manta ray (

Assessing the foraging ecology of a threatened species is necessary to understand their movement behaviour and habitat use patterns, which are essential for developing effective protection strategies. Here, the foraging ecology of reef manta rays ( Mobula alfredi ) in the Chagos Archipelago, a region encompassed by a vast no-take marine protected area (MPA), was investigated using stable isotope analysis of skin and muscle tissue. Enriched δ 13 C values suggest the population predominantly forages in nearshore environments. Skin δ 13 C values increased with increased rainfall, likely associated with the boosts in primary production and zooplankton biomass due to the coastal advection of seabird guano. Annual variations in δ 13 C values of skin and muscle were observed and are consistent with reduced nutrient transport associated with the effects of Indian Ocean Dipole oscillations, including a deepening of the thermocline, a suppression of cold-water upwelling, and reduced rainfall. Short-and long-term foraging strategies and locations were identified by applying hierarchical clustering, isotopic niche analysis, and Bayesian stable isotope mixing models to δ 13 C and δ 15 N of paired skin and muscle tissue samples. Two isotopically distinct groups of M. alfredi were identified, employing either local foraging strategies restricted to specific locations or wide-ranging strategies that likely mean they engage in regular migrations throughout the archipelago. Ninety-eight percent of M. alfredi were estimated to switch between strategies utilising and connecting multiple discrete nearshore habitats, emphasising their role in ecosystem functioning by facilitating the transport of nutrients across ecosystem boundaries. However, illegal, unreported, and unregulated fishing and lost or abandoned fishing gear commonly occur within the MPA. Locations of particular concern are Egmont Atoll as it is a highly active aggregation location and Peros Banhos Atoll where IUU frequently occurs and

Intraspecific differences in short-and long-term foraging strategies of reef manta ray (Mobula alfredi) in the Chagos Archipelago

Introduction
Developing effective protection and management strategies for threatened species requires an understanding of their movement ecology and habitat use patterns (Allen and Singh, 2016;Ogburn et al., 2017).These factors are intrinsically associated with foraging, which is fundamental to an animal's survival (Armstrong et al., 2016(Armstrong et al., , 2021;;Harris et al., 2020;Jeltsch et al., 2013;Stewart et al., 2018) and encompass a variety of elements, including dietary preferences, habitat suitability, energetic trade-offs and functional response (Brownscombe et al., 2022;Jeltsch et al., 2013;Jeschke et al., 2004).Species that predominantly depend on specific foraging habitats may have a limited range and are particularly at risk from localised threats and environmental change (Dulvy et al., 2014b;Pimm et al., 2014).
Reef manta rays (Mobula alfredi), a zooplanktivorous elasmobranch of the monogeneric Mobulidae family (mobulids) (Hosegood et al., 2020;Marshall et al., 2009;White et al., 2017), are widely distributed in the Indo-West Pacific in relatively small, seemingly fragmented regional populations that have limited home ranges (Couturier et al., 2012;Harris and Stevens, 2021;Marshall et al., 2022).Like many elasmobranch species, M. alfredi are under pressure from a multitude of anthropogenic threats, the most significant being target fisheries and bycatch [including from lost or abandoned (ghost) fishing gear], all of which are prevalent in the Indian Ocean (Croll et al., 2016;Curnick et al., 2021;Dulvy et al., 2014a;Filmalter et al., 2013;Lawson et al., 2017;Maufroy et al., 2017;O'Malley et al., 2017;Strike et al., 2022).Due to exploitation, M. alfredi populations have declined dramatically by > 90% in some regions (Lewis et al., 2015;Rohner et al., 2017).As population recovery is hindered by the species' conservative life-history traits (Ward-Paige, Davis and Worm, 2013;Dulvy, Fowler et al., 2014;Stevens, 2016), local extinctions may have already occurred, with more predicted in the next few years (Dulvy et al., 2014b;Marshall et al., 2022).While some protection exists, for example, local fishing bans, marine protected areas (MPAs), and international protection agreements and treaties [e.g., under the Convention on International Trade in Endangered Species (CITES, 2016) and the Convention on the Conservation of Migratory Species of Wild Animals (CMS, 2017)], illegal, unreported, and unregulated (IUU) fishing is a common issue (Collins et al., 2021b(Collins et al., , 2021a;;Fernando and Stewart, 2021;Hays et al., 2020).Therefore, it is essential to bridge current gaps in the knowledge of M. alfredi populations foraging ecology and concomitant movement patterns to highlight areas most likely to be at risk from fisheries and to inform and improve current regulation and enforcement strategies.Some of the current knowledge gaps can be addressed using stable isotope analysis of nitrogen (δ 15 N) and carbon (δ 13 C).Isotopic integration is a bottom-up process with geographically distinct biogeochemical and oceanographic processes influencing primary producers' δ 15 N and δ 13 C composition in aquatic environments (McMahon et al., 2013).By establishing an isotopic baseline for specific locations and applying stable isotope mixing models (SIMM), which incorporate the consumer's isotopic values and isotopic means of their prey (Parnell et al., 2013), foraging locations can be directly estimated.For example, Silva et al. (2019) incorporated mean isotopic values from bulk plankton samples collected over multiple years from locations in the Atlantic Ocean.These data allowed the identification of probable feeding areas of blue, fin, and sei whales and bridged knowledge gaps regarding their migration patterns (Silva et al., 2019).
Stable isotope analyses can also be used to assess changes over time by sequentially sampling the same individual (Kim et al., 2012;Malpica-Cruz et al., 2012), although this is challenging for highly mobile species in the wild (Silva et al., 2019).However, incorporation rates can vary by tissue type with more metabolically active tissue, such as liver, blood and potentially skin reflecting assimilated diet over weeks to months (Carlisle et al., 2012;Li et al., 2016a;Matley et al., 2015), while less active white muscle can be 300 -700 days (Kim et al., 2012;Logan and Lutcavage, 2010;Wyatt et al., 2019).Therefore, a multiple-tissue sampling approach can be used to investigate potential changes in foraging strategies and movement patterns over time (Curnick et al., 2019;MacNeil et al., 2005;Stewart et al., 2017).
The current study focuses on the M. alfredi population of the Chagos Archipelago, which occurs within a large no-take MPA that spans 640,000 km 2 (Hays et al., 2020;Sheppard et al., 2012).With little development, no commercial tourism, and low levels of pollution, anthropogenic pressures in the region are considered low (Hays et al., 2020).However, IUU fishing, ghost fishing gear, and fish aggregation devices pose a significant threat to marine life (Collins et al., 2021b;Curnick et al., 2021;Hays et al., 2020;Tickler et al., 2019).While the MPA is patrolled by an enforcement vessel, which has seized many threatened elasmobranch species (including manta rays) from illegal fishers (Collins et al., 2021b), it is estimated that only 10% of illegal fishing activity is detected (Jacoby et al., 2020;Price et al., 2010).Enforcement activity is prioritised in response to the perceived spatiotemporal risk of IUU, which is informed by scientific research in the region (Collins et al., 2021b;Jacoby et al., 2020).Therefore, it is essential that key areas where M. alfredi are at risk from IUU fishing and fishing gear entanglement are identified.Due to the remoteness and strict protection of the archipelago, M. alfredi research has been limited to satellite and acoustic telemetry to investigate broadscale movements and short-term fine-scale habitat use of this population (Andrzejaczek et al., 2020;Harris, 2019;Harris et al., 2021).Although these studies have provided valuable insight into the movement ecology of this M. alfredi population, they were limited to a small number of individuals, and their feeding ecology remained unstudied.
The current study aims to bridge current knowledge gaps in the foraging ecology of this M. alfredi population by: i) assessing demographic and spatiotemporal variation, ii) identifying short-and long-term foraging strategies, and iii) identifying distinct foraging habitats within the Chagos Archipelago.Concurrently, this research also aims to enhance our current understanding of the population's movement ecology, which can inform and develop improved protection and enforcement strategies in the region.These aims will be achieved using a multiple tissue sampling approach, hierarchical clustering, isotopic niche analysis, and Bayesian stable isotope mixing models to investigate intraspecific variation in M. alfredi δ 15 N and δ 13 C values and C:N ratios.

Study location
The Chagos Archipelago (Fig. 1) consists of seven discrete atolls and 58 low-lying islands (Hays et al., 2020;Sheppard et al., 2012).Within the archipelago, the four coral reef atolls (Egmont, Salomon, Peros Banhos, and Diego Garcia) which remain above sea level were identified as focal areas for M. alfredi based on observations, photographic identification records, acoustic and satellite telemetry (Andrzejaczek et al., 2020;Harris, 2019;Harris et al., 2021;Harris, unpublished data).Salomon and Peros Banhos Atolls are situated in the north of the archipelago, while Egmont and Diego Garcia are in the south, approximately 170 km and 200 km from the northern atolls, respectively.These four atolls vary in documented M. alfredi activity, with Egmont and Peros Banhos being the most and least frequented, respectively (Andrzejaczek et al., 2020;Harris, 2019;Harris et al., 2021;Harris, unpublished data).

Reef manta ray tissue sampling and processing
Fifty-five tissue biopsies were collected from 33 female (adults = 6, juvenile = 27) and 22 male (adults = 15, juvenile = 7) M. alfredi in Egmont (n = 46), Salomon (n = 8), and Peros Banhos (n = 1) between 2018 and 2022 (Table S1).Forty-eight contained both white muscle and skin, six samples contained only muscle, and one contained only skin.Two M. alfredi (CG-MA-0090 and CG-MA-0102) were sampled in two different years.Samples were collected from the posterior dorsal surface of individuals using a biopsy dart deployed by a modified Hawaiian hand sling while swimming behind the animal.Before sampling, the ventral side of each individual was photographed for identification (photo-ID), and their sex and visually estimated size class, which acts as a proxy for maturity (Table S2; Stevens, 2016), were recorded.Biopsy samples were kept on ice until returned to the laboratory (<2 h), where white muscle and skin were separated and then dried for 24 h at 60 • C. All samples were then shipped to the University of Plymouth for further processing.
Following recommendations in previous elasmobranch stable isotope studies (Hussey et al., 2012;Kim and Koch, 2012;Li et al., 2016b;Peel et al., 2019a), lipids and urea were extracted from all M. alfredi tissue.Lipids were removed by soaking in 2:1 chloroform-methanol solution in a 7 mL glass vial for 24 h, air-dried for 24 h and then oven dried at 60 • C for 24 h (Carlisle et al., 2017;Peel et al., 2019a).Urea was removed by soaking the samples in 1.5 mL of MilliQ (deionised water) for 24 h at room temperature then centrifuged for 3 min at 3000 rpm.The process was repeated three times, replacing the water each time (Marcus et al., 2017;Peel et al., 2019a).At the end of the third cycle, the samples were oven dried at 60 • C for a further 48 h (Marcus et al., 2017;Peel et al., 2019).Once dried, 800 μg and 1000 μg aliquots of skin and muscle, respectively, were weighed into tin capsules and sent to the University of Exeter for stable isotope analysis.The effect of lipid removal on skin and muscle tissue was examined (Supporting information Section 2, Fig. S1).

Plankton and fish sampling and processing
Plankton sample collection took place during six research expeditions between May 2018 and June 2022.A total of 84 plankton samples were opportunistically collected from the top 2 m of the water column from four atolls: Egmont (n = 43), Salomon (n = 18), Peros Banhos (n = 18), and Diego Garcia (n = 5) (Fig. 1 and Table S3).Of these, 58 were collected using a 100 µm mesh plankton net and 26 with a 200 µm mesh net.All samples were stored on ice and transported to a local vessel-or shore-based laboratory.Due to differences in equipment availability during each expedition, 66 of the samples were dried in an oven at 60 • C for 24 h and then shipped to the University of Plymouth (n = 63) or the University of Exeter (n = 3) for further processing before stable isotope analysis.The remaining 18 were filtered through pre-combusted (450 • C, 5 h) glass-fibre filters (GF/F 0.7 µm, Whatman), frozen (-20 • C) and transported to Woods Hole Oceanographic Institution, USA, for storage (Roche et al., 2022).
Due to the low number of samples obtained from Diego Garcia, plankton stable isotope data for this location were supplemented with that from fish of a similar trophic level (TL) to the planktonic prey of M. alfredi (baitfish of the family Clupeidae, TL = 2.3 -2.9, Rastgoo and Navarro, 2016).Six baitfish were regurgitated by tunas caught using baited hand-lines around Diego Garcia.Baitfish remained intact and fresh, indicative of very recent feeding and were therefore deemed uncompromised.These opportunistic samples were then measured and immediately frozen whole (− 20 • C) and shipped to the University of Windsor, where they were processed and analysed.All dried plankton and frozen baitfish samples were first homogenized before lipids were removed following the same method as described above for M. alfredi tissue, including rinsing in deionised water (DW).Once dried, 800 μg of each sample was weighed into tin capsules.The effect of lipid removal on plankton was examined (Supporting information Section 2, Fig. S2).Frozen plankton samples were thawed, rinsed with dilute HCl (5%) to remove carbonate material before being rinsed in DW, dried, and then 2000 μg was weighed into tin capsules for analysis.

Stable isotope analysis
Stable isotope analysis of all M. alfredi tissue samples (skin = 49, muscle = 54) and 76 plankton samples (66 that underwent lipid extraction and ten that did not, see Supporting information Section 2) were conducted at the University of Exeter Environmental & Sustainability Institute, stable isotope facility on the Penryn Campus, Cornwall, using a Sercon Integra 2 stable isotope analyser.Samples were run alongside in-house reference material (bovine liver; δ 13 C − 28.6, δ 15 N − 19.62 and Alanine; δ 13 C − 19.62, δ 15 N − 1.85), which was used to correct for instrument drift and to determine the δ 13 C and δ 15 N values of the samples.Analytical precision for both isotopes reported to two standard deviations for both alanine and bovine liver was ± 0.1‰.Eighteen plankton samples were analysed at the University of New Mexico Centre for Stable Isotopes (Roche et al., 2022) where they were run in triplicate using Thermo Scientific Delta V mass spectrometer with a dual inlet and Conflo IV interface connected to a Costech 4010 elemental analyser and a high-temperature conversion elemental analyser, respectively.Analytical precision was ± 0.2‰ based on an assessment of repeated analysis of internal proteinaceous reference materials, Pugel and Acetanilide (Roche et al., 2022).The six baitfish samples were analysed at the University of Windsor using an elemental analyser (Costech 4010) interfaced to a Thermo Finnigan Delta PLUS mass spectrometer.Precision was assessed by the standard deviation of replicate analyses of four standards; NIST1577c, internal lab standard (tilapia muscle), USGS 40 and Urea (n = 15 for all), measured ≤ 0.16‰ for δ 15 N and ≤ 0.15‰ for δ 13 C for all standards (Curnick et al., 2019).Both carbon and nitrogen isotope values were expressed in the delta (δ) notation as parts per million (‰), with δ 13 C and δ 15 N reported relative to Vienna Pee Dee Belemnite and atmospheric nitrogen, respectively.
Prior to data analysis, isotope δ 13 C values for plankton samples that did not undergo lipid extraction were transformed to δ 13 C normalized using the formula for aquatic organisms from Post et al., (2007).
To validate the inclusion of plankton samples with transformed δ 13 C values, the effect of lipid removal was compared to those calculated using the transform equation (Supporting information Section 2, Fig. S2).

Statistical analyses
A Welch's ANOVA followed by a Games-Howell post-hoc test using the 'oneway.test'function of the 'Rmisc' R package (Hope, 2022), and 'games_howell_test' function of the 'rstatix' R package (Kassambara, 2023), respectively, were applied to investigate the effect of sex, maturity status, size class, sampling year, and sampling location [grouped into north (Peros Bahos and Salomon) and south (Egmont)] on the values of δ 15 N and δ 13 C and the ratio of C:N in M. alfredi muscle and skin tissue separately.As samples were only collected in the north and the south of the archipelago in 2021, the test of sampling location was applied to samples collected in 2021 only.To assess the potential influence of rainfall on nutrient loading that may increase productivity via nutrient runoff, the relationship between δ 15 N and δ 13 C and the ratio of C:N of skin samples and total rainfall (30 days before sampling) was tested with a Spearman's rank correlation coefficient.The length of time for rainfall was conservatively estimated from the expected period where nutrient loading may have influenced isotopic ratios based on tissue turnover rates reported for other elasmobranch species (Carlisle et al., 2012;Logan and Lutcavage, 2010).Rainfall data were obtained from Meteoblue AG, Basel, Switzerland (www.meteoblue.com).The relationship between δ 15 N and δ 13 C and the ratio of C:N of muscle samples and rainfall were not tested as total rainfall during the estimated isotopic incorporation period [min 300 days, (Logan and Lutcavage, 2010)] did not vary for each M. alfredi due to the short period over which samples were collected.
Trophic enrichment between M. alfredi and plankton was calculated using the following equation: where X represents mean δ 13 C or δ 15 N for M. alfredi skin or muscle (consumer) and plankton (prey) (Fry, 2006).
To assess short-and long-term feeding patterns, only M. alfredi tissue samples containing both skin and muscle were included in subsequent analysis.A paired-sample t-test was used to compare the differences in the δ 15 N and δ 13 C and C:N values between skin and muscle tissue from samples that contained both tissue types, and a Spearman's rank was applied to see if there was a correlation between δ 15 N and δ 13 C of muscle and skin tissue (Li et al., 2016a).
Intraspecific differences in short-and long-term feeding patterns were investigated by first applying a clustering algorithm to identify isotopic niche (clusters) based on δ 15 N and δ 13 C of M. alfredi skin and muscle tissue separately.Ward's hierarchical clustering based on Euclidean distance was applied to scaled δ 15 N and δ 13 C values using the 'factoextra' R package (Kassambara and Mundt, 2022).The hierarchical clustering method was chosen over non-hierarchical methods (e.g., k-means) to allow for the examination of the interrelationship between short-and long-term feeding dynamics.The optimal number of clusters was identified using 30 indices via the 'NbClust' R package (Charrad et al., 2014), whereby the clustering strategy supported by the greatest number of indices was applied.For skin, the greatest number of indices identified three clusters as optimal; however, due to the low number of samples in one of the clusters, which may have implications on subsequent analysis, this was reduced to two, which was supported by the second-highest number of indices.The optimal number of clusters identified for muscle was two; therefore, Ward's hierarchical clustering algorithm was applied to a two-cluster solution for both tissue types.Clusters were investigated in terms of variation in δ 15 N, δ 13 C and C:N ratios.A Welch's ANOVA followed by a Games-Howell post-hoc test.
Generalized linear models (GLM) were used to investigate whether cluster assignment was influenced by M. alfredi sex, maturity, year of sampling or sampling location (north or south).Skin and muscle were modelled separately with a binary response variable (cluster one = 0, cluster two = 1), and all combinations of one to five explanatory variables were tested.The most parsimonious model was selected based on the corrected Akaike Information Criterion (AIC C ) test statistic, an asymptotically unbiased estimator of model quality (Burnham and Anderson, 2002).
Isotopic niche analysis of δ 13 C and δ 15 N values was used to investigate the potential difference in niche occupied by each of the identified clusters.δ 13 C and δ 15 N values for all tissue samples were plotted in biplot space.Standard ellipses corrected for small sample Table 1 Summary of mean ( ± se) values of δ 15 N, δ 13 C, and C:N ratios for plankton samples collected at Egmont, Salomon, and Peros Banhos Atolls, plankton combined with baitfish samples for Diego Garcia Atoll, and all Mobula alfredi skin and muscle tissue samples.sizes (SEA c ) with 40% confidence interval were then fitted to each cluster to estimate isotope niche via Bayesian inference using the 'SIBER' R package (Jackson et al., 2011).Pairwise overlaps in SEA c of clusters were then calculated and given as a percentage, indicating the degree of resource sharing/foraging characteristic (Letourneur et al., 2017;Thibault et al., 2021).The extent of ellipses was compared by fitting Bayesian models (SEA b ) utilising 10000 iterations (Jackson et al., 2011;Roche et al., 2022).When more than or equal to 95% of the posterior draws for one group was smaller than the other, the difference in standard ellipse area was deemed significant (Jackson et al., 2011;Roche et al., 2022).
Finally, Bayesian stable isotope mixing models (SIMM) were applied using the 'simmr' R package (Parnell 2016) to identify the most probable feeding areas from the four different study locations for each M. alfredi cluster.An isotopic baseline was established for each atoll by calculating the mean ( ± sd) of δ 13 C and δ 15 N values of bulk plankton samples pooled by location (Table 1).In the absence of controlled feeding studies for M. alfredi, a commonly applied method of incorporating trophic enrichment factor (TEFs, Δ) based on that of controlled feeding studies of similar species was used (Burgess et al., 2016;Carlisle et al., 2012;Silva et al., 2019).For muscle, ΔTEFs of 1.7 ± 0.5‰ (δ 13 C) and 3.7 ± 0.5‰ (δ 15 N) were estimated from a long-term feeding study of leopard sharks (Triakis semifasciata) (Kim et al., 2012).Due to the lack of controlled feeding studies that analyse skin tissue, the same ΔTEFs were used, but with a standard deviation of 1‰ applied for both isotopes to reflect variation among other tissue types observed in laboratory studies (Hussey et al., 2010;Kim et al., 2012;Malpica-Cruz et al., 2012).Concentration dependents given as elemental concentrations (C% and N%) of the potential food source were included in all models to account for the source mass contributions of δ 13 C and δ 15 N to the mixture.Models were fitted using 10,000 iterations, 1000 utilising 1000 burn-in iterations, with a thinning interval of ten and four Markov Chain Monte Carlo (MCMC) chains, and convergence was checked prior to interpretation (Coletto et al., 2021).

Ethics approval statement
All sampling activities were approved by the University of Plymouth Animals in Science Ethics Committee under permits ETHICS-24-2019, ETHICS-37-2020, and ETHICS-37-2022.
The which had only one sample), sampling year, or sampling location, which was tested for only 2021 samples (Welch's ANOVA, p > 0.05).Skin δ 13 C also did not differ by sex or maturity, but there was a significant difference between size classes (Welch's ANOVA, F 2,13.4 = 5.1, p < 0.05), with size class 4 being significantly more enriched in δ 13 C than size class 2 (Games-Howell post hoc, p < 0.05).There was also a significant difference between sampling years (Welch's ANOVA, F 3,17.9 = 5.8, p < 0.01), with depleted δ 13 C in 2019 compared to 2020 (Games-Howell post hoc, p < 0.05).Skin samples collected in 2021 in the north (Salomon and Peros Banhos) were significantly more enriched in δ 13 C than those at Egmont in the south (Welch's ANOVA, F 1,10.5 = 5.2, p < 0.05).There was a significant positive correlation between total rainfall 30 days prior to sampling and δ 13 C of skin tissue (r = 0.35, p < 0.05) and a negative correlation between rainfall 30 days prior to sampling and C:N ratio (r = 0.35, p < 0.05).There was no correlation between rainfall and skin δ 15 N.For muscle tissue (Fig. 2 and Fig. S4), δ 13 C and δ 15 N and C:N did not differ by sex or maturity, size class or sampling location (Welch's ANOVA, p > 0.05).δ 15 N did not differ by sampling year, but there was a significant difference between years in δ 13 C values (Welch's ANOVA, F 4,15 = 3.8, p < 0.05) and C:N ratio (Welch's ANOVA, F 4,16 = 3.5, p < 0.05).δ 13 C was significantly enriched in 2019 compared to 2018 (Games-Howell post hoc p < 0.05) and 2019 compared to 2022 (Games-Howell post hoc p < 0.05).The ratio of C:N was significantly higher in 2021 compared to 2019 (Games-Howell post hoc p < 0.05) and in 2022 compared to 2019 (Games-Howell post hoc p < 0.05).
The most parsimonious GLM model (Table S4) for cluster assignment, identified from the lowest AIC and the greatest relative weight, included location for both M. alfredi skin ( Δ AIC = 0, w AIC = 0.25) and muscle samples ( Δ AIC = 0, w AIC = 0.17); however, the effect was not significant.
The biplot of δ 13 C and δ 15 N of all paired tissue samples fitted with standard Bayesian ellipse areas (SEA c ) indicates clear segregation between the two skin clusters and the two muscle clusters with 0% overlap SEA c (Fig. 5).Both skin clusters overlapped with muscle cluster one, with skin cluster one having the highest area of overlap (total area of the skin cluster that overlapped with the muscle cluster = 75%), while skin cluster two overlapped with muscle cluster one by 40%.
A comparison of SEA b (Fig. 5) showed skin cluster one occupied a smaller isotopic niche area than skin cluster two (probability = 97%), and both skin clusters occupied a smaller isotopic niche area than both muscle clusters (probability = 100%), which both occupied a similar isotopic niche area.
The SIMM for skin (Fig. 6) indicates that short-term feeding activity of M. alfredi in cluster one was predominantly concentrated in the south of the archipelago, with > 75% of the dietary input estimated to be from the southern atolls [Diego Garcia and Egmont Atoll, 41 ± 11% and 35 ± 15%, respectively] and only a small contribution from the northern atolls [Peros Banhos and Salomon Atoll, 16 ± 10% and 9 ± 6%, respectively].In contrast, short-term feeding activity of M. alfredi in cluster two was more widespread, with similar dietary contributions from Egmont (31 ± 11%), Diego Garcia (29 ± 7%), and Peros Banhos (29 ± 11%).The contribution from  S1). Green and red trees depict M. alfredi in cluster one and cluster two that engage in local and wide-ranging foraging strategies, respectively.Lines connect each manta showing changes between short-and long-term foraging strategies.Lines are coloured according to the shortterm foraging clusters.
The SIMM for muscle (Fig. 6) indicates that long-term feeding activity in cluster one was similar to short-term feeding activity in skin cluster one, with the predominant dietary input (74%) coming from the southern atolls (Diego Garcia and Egmont, 40 ± 15% and 34 ± 20%, respectively), with only a low contribution from Perso Banhos and Salomon Atoll in the north (16 ± 11% and 11 ± 7%, respectively).Muscle cluster two also showed a similar widespread feeding pattern to skin cluster two, with a similar contribution from both the north (Perso Banhos and Salomon Atoll, 39 ± 12% and 17 ± 8%, respectively) and south atolls (Egmont and Diego Garcia, 31 ± 8% and 14 ± 5%, respectively).
To further investigate short-and long-term feeding patterns, M. alfredi in cluster one for muscle and skin is hereafter referred to as displaying a 'local foraging strategy', whereby feeding was predominantly concentrated in the south (Egmont Atoll and Diego Garcia).Mobula alfredi in cluster two for both tissue types are hereafter referred to as displaying a 'wide-ranging foraging strategy', as feeding was inferred to occur throughout the archipelago.
Of the 48 M. alfredi sampled, the majority (41.7%) displayed a wide-ranging short-and long-term foraging strategy (Fig. 7), 33% displayed a short-term local and long-term wide-ranged foraging strategy, 21% appear to engage in a wide-ranging short-term and local long-term foraging strategy, and 4.2% appear to engage in a local short-term and local long-term foraging strategy.
All demographics (Table S1 and Fig. 7) except adult females appear to engage predominantly in a short-and long-term wideranging foraging strategy (female juveniles = 40%, male adults = 50%, male juveniles = 50%).Adult females predominantly displayed a short-term local and long-term wide-ranging foraging strategy (80%).A short-and long-term local foraging strategy was only inferred for two M. alfredi, one juvenile female and one adult male.

Mobula alfredi foraging ecology
The foraging ecology of M. alfredi was investigated using δ 13 C, δ 15 N, and C:N ratio of skin and muscle samples, none of which varied significantly by demographics.Similar findings were reported in M. alfredi in Seychelles (Peel et al., 2019a) and Mobula birostris in Ecuador (Burgess et al., 2016).δ 13 C was enriched in large females (size class 4) compared to small juveniles (size class 2), and there was no a significant correlation of δ 13 C between paired M. alfredi skin and muscle samples.These findings likely reflects greater feeding efficiency with the development of filtering gill plate structures with increased size or ontogenetic differences in movement patterns as reported for other filter feeders, such as whale sharks, Rhincodon typus (Borrell et al., 2011;Prebble et al., 2018).
Skin collected in the south of the Chagos Archipelago was depleted in δ 13 C compared to the north, which may be attributed to the higher seabird densities found in the northern atolls (Carr et al., 2022;Wu et al., 2021).High seabird densities can be associated with greater nitrogen input into the marine environment from seabird guano, which can boost primary production and zooplankton biomass (Benkwitt et al., 2021(Benkwitt et al., , 2019;;Graham et al., 2018;Kim et al., 2022;Lorrain et al., 2017;McCauley et al., 2012;Mulimbwa et al., 2014).Increased rainfall may lead to higher coastal advection of guano-derived nitrogen, which is consistent with skin δ 13 C values being most depleted in 2019, when the lowest total rainfall occurred and the positive correlation between rainfall 30 days before sampling and skin δ 13 C values.
Due to higher nutrient concentrations, δ 13 C values of primary producers are generally more enriched in nearshore and benthic systems (δ 13 C >− 17‰) compared to relatively oligotrophic pelagic environments [δ 13 C = − 22‰ to ~ − 17‰ (France, 1995)].In the current study, enriched δ 13 C values for skin (− 15‰) and muscle (− 15.9‰) are consistent with nearshore feeding, suggesting M. alfredi in Chagos may have a higher affinity to nearshore environments than populations in other regions.For example, in Seychelles, δ 13 C values of M. alfredi muscle tissue were depleted (− 18.6‰), indicating a predominantly pelagic diet (Peel et al., 2019), and on the Great Barrier Reef, M. alfredi muscle δ 13 C values of − 17.4‰ are consistent with a diet derived from both pelagic and nearshore sources (Couturier et al., 2013).The apparent higher contribution of nearshore resources for M. alfredi around the Chagos Archipelago could be associated with high seabird densities; increased M. alfredi activity in nearshore habitats with higher seabird densities has previously been observed for M. alfredi in Palmyra Atoll, central Pacific (McCauley et al., 2012).It is also plausible that the seemingly higher affinity to nearshore environments could be associated with the lack of in-water disturbance within the Chagos Archipelago.For example, the absence of tourism means this population have very little interaction with snorkellers or divers, which have been shown to have a significant negative influence on M. alfredi foraging behaviour, potentially driving individuals away from productive feeding locations (Murray et al., 2019;Venables et al., 2016).Feeding in less productive environments may reduce population fitness as it could lead to reduced fecundity and offspring survivorship (Stevens, 2016).Future research into the effect of in-water disturbance on M. alfredi would benefit from the incorporation of stable isotope analysis of multiple populations from regions with varying levels of disturbance to investigate the extent to which this may affect foraging strategies and potentially population fitness.
Annual variations in δ 13 C values for skin are consistent with variability in oceanographic conditions that have been shown to influence primary production (Brewin et al., 2012).For example, the Indian Ocean Dipole an ocean-atmosphere interaction that causes variability in sea surface temperatures and anomalous wind and precipitation events (Saji et al., 1999;Huang et al., 2022).The IOD cycles through active (positive and negative) and neutral phases, active phases typically persist for less than six months (Cai et al., 2013;Vinayachandran et al., 2009).In 2019, the Indian Ocean experienced the strongest positive IOD event of the 21st Century (Shi and Wang, 2021), which led to a deepening of the thermocline, a suppression of cold-water upwelling, and reduced rainfall (Du et al., 2020;Lee et al., 2022;Ratna et al., 2021).These anomalies reduced the transport of nutrients and plankton, which was evident from significantly reduced chlorophyll-α (chl-α) concentration in the Indian Ocean of > 30% compared to that of typical years (monthly mean calculated from 2012 to 2020, Shi and Wang, 2021).Subsequently, decreased chl-α concentration likely limited the availability of zooplankton prey for M. alfredi, thus leading to the depleted δ 13 C values of skin in 2019.For muscle, which reflects δ 13 C assimilated over a much longer period, short-term changes that affect food availability are less likely to strongly influence δ 13 C values.However, the positive IOD event of 2019 was prolonged, developing in March 2019 (Du et al., 2020) with the depression of chl-α levels, continuing until May 2020 (Shi and Wang, 2021).Therefore, it may have partly contributed to the progressively depleted δ 13 C values observed between 2020 and 2022.

Short-and long-term foraging strategies
Paired skin and muscle samples showed that skin was significantly enriched in δ 13 C, and depleted in δ 15 N, which is consistent with skin and muscle differences for silky (Carcharhinus falciformis) and blue sharks (Prionace glauca) (Li et al., 2016a).It is also consistent with dorsal fin compared to muscle tissue for shortfin mako (Isurus oxyrinchus), common thresher (Alopias vulpinus), and P. glauca (Hussey et al., 2011), as well as δ 13 C values for white sharks (Carcharodon carcharias) sampled in the Pacific (Carlisle et al., 2012).It is likely that these differences are associated with the higher dietary assimilation rate of skin tissue compared to muscle or tissue-specific diet discrimination factors (Hussey et al., 2011;Nigel E Hussey et al., 2012;Li et al., 2016a).However, differences may also be associated with biochemical variations between tissue types, given that each amino acid that makes up the tissue proteins can have different δ 13 C and δ 15 N values (Carlisle et al., 2012;Li et al., 2016a).
Cluster and SIMM analysis of paired M. alfredi skin and muscle tissue identified two distinct short-and long-term foraging strategies.For both tissue types, these strategies were classified as local (foraging restricted to Diego Garcia and Egmont Atoll in the south of the archipelago) and wide-ranging (foraging in both the north and south of the archipelago) foraging.Short-and long-term foraging strategies indicate 38% and 25% of M. alfredi, respectively, predominantly engaged in local foraging.These individuals had significantly higher δ 13 C values than the remaining M. alfredi, that appear to engage in a more wide-ranging foraging strategy.While M. alfredi habitat use and the oceanographic dynamics of Diego Garcia have not been investigated, Egmont Atoll is a hotspot for aggregations where M. alfredi take advantage of tide-enhanced feeding opportunities produced by internal waves and cold-water bores that facilitate transportation of zooplankton biomass to surface waters (Harris et al., 2021).The enhanced zooplankton biomass at this location may be associated with the enriched δ 13 C values of M. alfredi inferred to have a high dietary input from this location.
Large planktivores must frequently consume vast quantities of prey (Meekan et al., 2015) and require high densities to balance the energetic cost of feeding (Armstrong et al., 2016;Meekan et al., 2015).These resources are spatially and temporally limited; therefore, M. alfredi often travel long distances to utilise optimal feeding opportunities (e.g., Couturier et al., 2011;Jaine et al., 2014;Armstrong et al., 2019;Harris, 2019;Harris and Stevens, 2021).Overall, 96% of individuals engaged in a short-and/or long-term, wide-ranged foraging strategy, with 54% of all individuals switching between local and wide-ranged strategies over short-and long-term time scales, which indicates regular habitat switching to locate prey.However, a wide-ranging foraging strategy considerably increases energetic demands, particularly for smaller individuals (Nøttestad et al., 1999;Peel et al., 2019a).Here, the majority of small juveniles (77%) displayed a long-term, wide-ranging foraging strategy, and 54% displayed both a short-and long-term wide-ranging foraging strategy, which included the single size class 1 (neonate) M. alfredi sampled.These findings are contrary to the expected size-restricted migration, whereby the extent of migration may be lower for smaller individuals (Nøttestad et al., 1999), and studies that have shown that juvenile M. alfredi tend to display high levels of residency in specific protected locations, such as lagoons and shallow water reefs (Germanov et al., 2019;Harris et al., 2021;Peel et al., 2019b;Stevens, 2016).Alternatively, due to the slow turnover rate of elasmobranch muscle tissue, the δ 15 N and δ 13 C values of juvenile M. alfredi potentially reflect that of their mother, a concept called maternal meddling, which has been suggested to occur in other elasmobranchs such as bull (Carcharhinus leucas) and Atlantic sharpnose sharks (Rhizoprionodon terraenovae) (Belicka et al., 2012;Matich et al., 2019;Olin et al., 2011).Therefore, the suggested wide-ranging foraging strategy of juvenile M. alfredi in the current study may result from maternal meddling and is consistent with all adult females in this study displaying a long-term, wide-ranged foraging strategy.However, further research into M. alfredi tissue turnover rates is required to establish how long the effects of maternal meddling may influence the isotopic signatures of juveniles.
Mixing models for muscle samples suggest a substantial proportion of dietary input from Peros Banhos Atoll, with an estimated 38% contribution from this location for 75% of M. alfredi.It should be noted that stable isotope mixing models assume that all possible dietary sources are contained within those provided to the model.Therefore, it is plausible that the estimated dietary contribution from Peros Banhos, as well as other locations, has been influenced by unknown dietary sources, for example, offshore prey.Nonetheless, the evidence presented here suggesting M. alfredi predominantly relies on nearshore resources, and previous observations of M. alfredi utilising each of the locations included in the model suggest that the species are likely to be foraging at the estimated locations.Due to the apparent low M. alfredi utilisation of Peros Banhos suggested by satellite tracking data (Andrzejaczek et al., 2020;Harris, 2019), and logistical constraints of conducting regular, systematic surveys, this location has received relatively little attention compared to other locations within the archipelago, such as Egmont Atoll.Therefore, future M. alfredi research around the Chagos Archipelago would benefit from increased survey efforts around Peros Banhos, particularly as IUU activity is predominantly concentrated in the northern regions of the archipelago (Hays et al., 2020;Jacoby et al., 2020).

Isotopic niches
Increased productivity has been suggested to increase isotope niche size as result of foraging in a more tropically complex ecosystem (Miller et al., 2019;Roche et al., 2022).Conversely, it has also been suggested that increased productivity decrease isotope niche size due to the consumer's ability to focus on preferred resources in highly productive environments (Lesser et al., 2020).While M. alfredi foraging habitat switching has been shown to be associated with increases in productivity (Harris et al., 2020), there is currently minimal evidence that they target specific prey but rather take advantage of high densities or biomass of zooplankton (Armstrong et al., 2016(Armstrong et al., , 2021)).Therefore, the larger niche width of M. alfredi that engaged in a short-term, wide-ranging foraging strategy compared to a local foraging strategy likely reflects switching between various contrasting foraging habitats and a wide range of resources (Lesser et al., 2020;Prebble et al., 2018;Roche et al., 2022), which is supported by SIMM foraging location estimates.
Low levels of productivity may increase niche overlap (Lesser et al., 2020), which has been observed to occur interspecifically with sympatric mobulid species (Stewart et al., 2017).In this study, no overlap was observed between both the skin clusters or both muscle clusters of M. alfredi, which may reflect the high productivity of the region and lack of competition for resources (Miller et al., 2019).The lack of overlap also supports the difference in dietary contributions estimated by SIMM for the identified local and wide-ranging foraging strategies.These strategies were not found to be associated with demographics, size, year, or location of sampling, but could potentially be attributed to individual or social preferences as suggested for other megafauna (Silva et al., 2019).
Broadly, there is no evidence that M. alfredi form long-term social groupings, however, studies that investigate the social dynamics of M. alfredi suggest that groups of females may maintain short-term (several weeks) social relationships, frequently utilising the same locations at the same time (Perryman et al., 2022(Perryman et al., , 2019)).There is also some evidence of cooperative or coordinated feeding behaviour within groups of M. alfredi (Armstrong et al., 2021;Stevens, 2016).Here, the range of δ 13 C and δ 15 N values of individuals that engaged in short-term local foraging strategies was lower than for those that displayed a short-term wide-ranging foraging strategy.This could indicate short-lived social relationships between individuals that engaged in the same foraging strategies, potentially cooperatively, resulting in more similar isotopic signatures.Social dynamic research also indicates that some individuals may be less likely to maintain social relationships, frequently moving between aggregation sites (Perryman et al., 2022).These individuals may play a key role in connecting habitats, potentially making them essential biological mechanisms that play a crucial role in population dynamics, such as gene flow (Lassauce et al., 2022) and ecosystem functioning via the transfer of nutrients between discrete locations (Peel et al., 2019a).Individuals with muscle δ 13 C values akin to that of a more pelagic diet (δ 13 C <− 17‰) were all inferred to engage in a long-term wide-ranging foraging strategy, which may be partially associated with feeding along pelagic migratory pathways connecting these discrete nearshore habitats.

Conservation implications
The apparent reliance on nearshore environments makes M. alfredi particularly vulnerable to activities such as IUU in the shallow reef environments surrounding the atolls.Locations particular concern include Egmont, which contributes to a high proportion to M. alfredi diet and where large aggregations regularly occur (Harris et al., 2021), and Peros Banhos, where foraging activities and IUU likely cooccur.Many M. alfredi sampled here utilise multiple discrete nearshore habitats throughout the Chagos Archipelago, indicating a high level of foraging plasticity.This plasticity may potentially offer some resilience to environmental changes that impact productivity by enabling habitat switching (Bolnick et al., 2003).It also indicates a high level of connectivity between these locations, facilitating the transport of nutrients across ecosystem boundaries and maintaining gene flow (Lassauce et al., 2022;Peel et al., 2019a).However, frequent migrations between atolls raise concerns over the species vulnerability to IUU fishing along migration corridors (Germanov and Marshall, 2014;Harris and Stevens, 2021;Stevens and Froman, 2018).While these findings can be used to enhance knowledge of the spatial impacts of IUU on M. alfredi within the Chagos Archipelago, further research is required to identify how this risk varies temporally to assist enforcement patrols.Due to the remoteness of the archipelago, it is recommended that this research employs satellite and acoustic telemetry using a dual-tagging approach to provide a detailed understanding of fine-scale habitat use, broadscale movements, including the identification of migration corridors, and behaviour of M. alfredi (Cochran et al., 2019;Knochel et al., 2022;Peel et al., 2020;Spaet et al., 2020).With the integration of biological and physical measurements, such research would also assist in further resolving the drivers of aggregations to identify how M. alfredi may be impacted by environmental change and anthropogenic disturbance in the region (Harris, 2019;Harris et al., 2021).

Conclusion
The current study investigated intraspecific variation in stable isotopes (δ 15 N and δ 13 C values and C: N ratios) of the Chagos Archipelago M. alfredi population using a multiple tissue sampling approach.Enriched δ13C values suggest that the population predominantly forage in nearshore environments, where potentially periods of heavy rainfall can boost primary production and zooplankton biomass due to the coastal advection of seabird guano.Annual variation in skin and muscle tissue δ13C values were observed, which potentially coincided with reduced nutrient transport associated with the effects of Indian Ocean Dipole oscillations, including increased thermocline depths, suppression of cold-water upwelling, and decreased rainfall.Hierarchical clustering, isotopic niche analysis, and Bayesian stable isotope mixing models applications to paired skin and muscle tissue samples identified two distinct groups of M. alfredi that engage in either local (restricted to specific locations) or wide-ranging foraging strategies (feeding throughout the archipelago).Most sampled M. alfredi (98%) switched between strategies, potentially facilitating the transport of nutrients between multiple discrete nearshore habitats throughout the archipelago, emphasising the species' role in ecosystem functioning.However, frequent migrations that may occur between atolls and the species' reliance on nearshore environments means individuals are likely to be highly vulnerable to IUU fishing activities and entanglement from lost or abandoned fishing gear throughout the MPA.Egmont and Peros Banhos Atoll are of particular concern as isotope data suggest these are key foraging locations where there may be high IUU fishing activities.Mobula alfredi are also at risk from these activities along migration pathways between all foraging locations.By proactively identifying and addressing potential threats facing M. alfredi populations at the locations identified in this study, informed conservation strategies can be developed to protect this population and their critical habitats.

Fig. 1 .
Fig. 1.The Central Indian Ocean with the Chagos Archipelago indicated within the red box (left inset).The Chagos Archipelago with sampling sites for Mobula alfredi tissue (red dots), plankton and fish (green dots).Close up of each atoll shown in the right panels; from top to bottom: Peros Banhos, Salomon, Egmont, and Diego Garcia.

Fig. 2 .
Fig. 2. Comparison of Mobula alfredi skin δ 13 C values (top three panels) and muscle δ 13 C values and C:N ratio (bottom two panels) for categories with significant differences.All other comparisons are shown in Fig. S4.Solid black circles are mean values ± se.Common letters in each category indicate no significant difference.

Fig. 3 .
Fig. 3. Comparison of Mobula alfredi muscle and skin values of δ 13 C (left), δ 15 N (middle), and C:N ratios (right).Open circles show observed values with broken lines connecting paired samples.Solid black circles are mean values ± se.

Fig. 4 .
Fig. 4. Comparison of δ 13 C and δ 15 N values and C:N ratio between Mobula alfredi skin cluster one and cluster two (top three panels) and muscle cluster one and two (bottom three panels).Solid black circles are mean values ± se.Common letters in each category indicate no significant difference.

Fig. 5 .
Fig. 5. δ 13 C and δ 15 N biplot of Mobula alfredi skin and muscle clusters with standard ellipse areas (SEA c ) estimated for each cluster shown as a solid line (top panel).Density plot of standard Bayesian ellipse areas (SEA c ) of Mobula alfredi skin and muscle clusters.Black dots show the cluster mode, and boxes of increasing size show the credibility intervals (50%, 75%, and 95%).Red crosses are equal to the mean standard ellipse areas corrected for small sample size (SEA c ). Different letters indicate a significant difference according to Bayesian inference (p < 0.05) (bottom panel).

Fig. 6 .
Fig. 6.Isotopic biplot of δ 13 C and δ 15 N values for Mobula alfredi skin and muscle clusters in relation to potential prey sources from Diego Garcia (DG), Egmont (EG), Pehros Banhos (PB), and Salomon Atolls (SA) adjusted for trophic discrimination values shown as mean ± sd (top panel) and estimates of dietary contribution from these atolls from Mobula alfredi skin (middle two panels) for cluster one (left) and cluster two (right) and muscle (bottom two panels) for cluster one (left) and cluster two (right).Coloured boxes encompass the interquartile range, the middle line denotes the median, and the whiskers show the minimum and maximum values.

Fig. 7 .
Fig. 7. Tanglegram comparison of the hierarchical clustering of Mobula alfredi skin (short-term foraging) and muscle (long-term foraging) by manta ID number (TableS1).Green and red trees depict M. alfredi in cluster one and cluster two that engage in local and wide-ranging foraging strategies, respectively.Lines connect each manta showing changes between short-and long-term foraging strategies.Lines are coloured according to the shortterm foraging clusters.