The Devil We Don't Know: Investigating Habitat and Abundance of Endangered Giant Devil Rays in the North-Western Mediterranean Sea

The giant devil ray Mobula mobular, the only Mediterranean mobulid, is subject to mortality caused by directed and accidental captures in fisheries throughout the region. Whilst the combination of human impacts, limited range and a low reproductive potential is not inconsistent with its endangered listing, there are insufficient data to enable a quantitative assessment of trends. Without this, it is difficult to assess and prioritise threats and develop effective conservation actions. Using results from aerial surveys conducted between 2009 and 2014 over the Ligurian, Corsican, Sardinian, northern and central Tyrrhenian seas (626,228 km2), this study provides the first quantitative information on giant devil ray abundance and habitat choice in the western Mediterranean. Devil rays were observed in all seasons except winter, with their estimated abundance in the study area peaking in summer. The overall uncorrected mean density in the study area during summer was estimated at 0.0257 individuals km-2 (range: 0.017–0.044), resulting in a total abundance estimate of 6,092 (12.7%CV) individuals at the surface; once corrected for availability bias, this estimate indicates a summer presence of >12,700 devil rays in the study area. Rays were mostly observed alone even if occasionally, larger aggregations up to a maximum of 18 individuals were observed. Although observed throughout the study area, spatial modelling identified their preferred habitat to be over a broad strip connecting the Tuscan Archipelago to Eastern Sardinia, over a wide range of water depths ranging from 10 to 2000m. The observed seasonal changes in giant devil ray distribution in this study, combined with similar evidence from other areas in the Mediterranean, support the hypothesis that the species undertakes latitudinal migrations across the region, taking advantage of highly productive waters in the north during summer, and warmer southern waters during winter.


Introduction
The giant devil ray Mobula mobular is one of the largest elasmobranchs, with a maximum disc width greater than 4.5m and a weight that can exceed 1.5 tonnes [1]. While the low reproductive potential common to all mobulids makes these rays vulnerable to fishing [2,3], this may be particularly true for M. mobular with its limited range (the Mediterranean, possibly extending into the Northeast Atlantic [4,5]). Until recently, no major fisheries specifically targeting giant devil rays were known to occur in the Mediterranean. Only small (3-4 animals annually) catches were known to be taken in a few specific locations (e.g. Algeria [6] and Sicily [7]). However, this picture substantially changed after the discovery of a directed seasonal fishery in the Levantine Sea off Gaza. Between January and April from 2005 to 2014, an annual mean catch of 83.2 rays was taken using purse seines, with peaks of 363 and 370 rays, respectively, in 2006 and 2013 [8]. Fishery bycatch further adds to giant devil ray mortality in the Mediterranean, in particular in the now illegal (since 2002) large-scale pelagic driftnets [8,9,10,11,12,13], and to a smaller extent in bottom-set nets [14,15,16], purse seines for small pelagics [2], bottom trawls [15], pelagic trawls [17], longlines [16,18,19,20], fixed traditional tuna traps [21,22,23]. In addition, there is evidence of sport fishing most often followed by release (e.g., youtube. com). Today, however, bycatch appears to be a lesser concern than directed catches off Gaza, particularly as a consequence of the progressive phasing out of pelagic driftnets in the Mediterranean, triggered by the 2002 ban. All sources of mortality considered, the conservation status of the giant devil ray was assessed as "Endangered, A4d" in IUCN's Red List [24]. The species is also explicitly protected by Croatia and Malta, as well as being listed in Annex II ("List of endangered and threatened species") to the Barcelona Convention Protocol on Specially Protected Areas and Biological Diversity in the Mediterranean (SPA/BD), in Appendix II ("Strictly protected fauna species") to the Bern Convention on the Conservation of European Wildlife and Natural Habitats, and in both Appendices I and II to the Convention on the Conservation Migratory Species of Wild Animals.
Notwithstanding their status and the existing international legal framework for their protection, there are no enforced management measures in the Mediterranean or dedicated monitoring programmes [25]. Two scientific issues may be hampering proper management and conservation actions for this endangered elasmobranch: (1) uncertainty over its range caused by taxonomic difficulties; and (2) the intrinsic complexities of undertaking quantitative investigations of the population ecology and status of free-ranging elasmobranchs.
With respect to taxonomy, M. mobular can be confused with the spinetail devil ray M. japanica, a circumtropical species which is found in the tropical and warm-temperate North Atlantic, although not the Mediterranean [4]. Since M. mobular and M. japanica differ mainly in size (the latter being smaller) and in a number of subtle morphological and morphometric characters that cannot be detected during sightings at sea, it is possible that past reports of giant devil rays from the Atlantic-e.g., in coastal African waters from Morocco to Senegal, the Canary Islands, Madeira, the Azores, Portugal, and southern Ireland [1]-might have been in fact spinetail devil rays [4,26]. Given that the only recorded mobulid in the Mediterranean is M. mobular, this uncertainty really only concerns its presence (or not) in Atlantic waters. Furthermore, the possibility that M. mobular is endemic to the Mediterranean should add to the need to develop effective conservation actions and monitoring.
Investigating the abundance and distribution of pelagic elasmobranch populations is challenging, due to their range, somewhat inaccessible habitat and the fact that, unlike mammals and reptiles, they do not have to come to the surface to breathe [27,28]. Traditionally therefore, most of the existing knowledge of their biology and demographics has been obtained from fisheries data (e.g., [26]), which is problematic for species such as giant devil rays that are only subject to direct fishing on a small scale. Recent technical advances, such as the use of telemetry, have provided important new insights into the lives of larger elasmobranchs such as basking sharks, Cetorhinus maximus (e.g., [29]) and whale sharks, Rhincodon typus (e.g., [30]). Telemetry data have revealed previously unsuspected details about the extent and timing of horizontal and vertical movements of mobulids such as M. japanica [31], M. tarapacana [32], and Manta alfredi [33]. Of most relevance to this paper are the results of a telemetry experiment describing the movements and diving habits of M. mobular tagged in the Strait of Messina [34]. Telemetry data are invaluable in interpreting the results of aerial survey data as discussed later in this paper.
A separate development has been the increased use of aerial surveys to investigate the distribution, density and abundance of marine species such as cetaceans over large and remote areas, using the well-established approach of 'Distance sampling' e.g. [35]. They have been used in studies of large elasmobranchs such as manta rays [27], whale sharks [36,37,38], basking sharks [39], and great white sharks, Carcharodon carcharias [40]. Most recently, aerial surveys of the Adriatic Sea focussing primarily on cetaceans have generated valuable data on M. mobular in that Mediterranean sub-region, including the provision of novel information on giant devil ray distribution, density, and absolute abundance [25].
The present paper presents results on giant devil rays obtained from a series of aerial surveys, as in the case of [25], primarily aimed at investigating cetacean ecology, providing new insights into giant devil ray distribution, densities, and seasonality. Surveys were conducted between 2009 and 2013 in the North-Western Mediterranean Sea (Fig 1) [41,42], including the waters of the Pelagos Sanctuary for Mediterranean Marine Mammals (hereafter`Pelagos Sanctuary`), a marine protected area established in 1999 by a treaty amongst France, Italy and the Principality of Monaco [43]. The establishment of this sanctuary to protect cetaceans, by helping the phasing out of the use of driftnets [44], has extended umbrella protection to giant devil rays.

Materials and Methods
The study area (626,228 km 2 ) covered by the aerial surveys is shown in Fig 1. Overall, between January 2009 and January 2014, 88 days were spent surveying 5 different sub-areas, some twice (Table 1).
Surveys covered areas including broad bathyal plains at depths between 2000 and 2500 m associated with narrow continental shelves (sub-areas B, C, D, E, and the western portion of A), as well as shallower areas with extended shelves (eastern portion of A).
The aircraft used for all surveys was a two-engine high-wing Partenavia P-68 equipped with bubble windows to allow direct observation of the trackline under the plane, flying at the altitude of 750 feet (229 m) at a ground speed of 100 knots (185 km -1 ). Observations were performed by three experienced observers, two seated in the rear searching for animals through the opposite lateral bubble windows, the third in the co-pilot seat recording the data on a semiautomatic logging system. The plane flew along equally-spaced (15 km apart except in winter and summer 2009 when it was 10km) parallel transects, designed to provide equal coverage probability, along a total length of 60,895 km. Data collected for each sighting included: GPSderived position, date and time of sighting, declination angle when the sighted animal(s) was seen abeam, aggregation size, and observer name. Primary search effort (searching carried out in acceptable conditions, i.e., with Beaufort 3 and good visibility) and altitude were recorded directly onto a laptop through data-logging software connected to the GPS. Additional information that may affect sightability (e.g., sea state, glare, cloud cover) was entered at the beginning of each transect and/or whenever changes occurred along the trackline. The perpendicular distance of each sighting from the trackline was derived trigonometrically from the altitude of the aircraft at the time of sighting and the declination angle to the sighting, measured with a Suunto clinometer.
Based upon the hypothesis that the peak presence of rays in the surface waters occurred throughout the study area during summer (as inferred from the high summer density of giant devil rays at the surface in sub-area A, starkly contrasting with their absence during winter; see Table 2), density estimates and spatial analyses have been carried out only with data coming from surveys conducted during those months (sub-areas A, B, and C; shaded in Fig 1).
Two methods to estimate giant devil ray abundance were used: (1) a design-based method [35], based on a survey design that ensures equal coverage probability across the study area; and (2) a model-based method [45,46], in which line-transect sampling is combined with spatial analysis. The former has been the standard method used to estimate density and abundance from systematically collected survey data. Such estimates therefore enable comparison with many previous surveys and other areas. The latter approach potentially increases precision and also allows the production of density maps which can be especially valuable for developing conservation-and management-oriented actions. A further advantage is that abundance can be estimated for any subarea within the study area [45].

Design-based method
A standard line transect distance sampling approach was implemented to obtain abundance and density estimates [35,47]; survey design and sampling strategy were obtained using the latest available version of the dedicated software Distance [48]. Specific details and a full description of the methods and data collection protocols can be found in [41,42].
Estimates were obtained using both the Conventional Distance Sampling, CDS [35] and the Multi Covariate Distance Sampling, MCDS [49] analyses engines in the Distance Software. In CDS, the detection probability is modelled solely as a function of the perpendicular distance to sightings.
The abundance in each area is estimated as: Table 1. Subareas (letters indicate corresponding sectors in Fig 1), listed in chronological order of survey timing, with subarea surfaces and transect lengths (on average over 90% of the planned transect length was covered). Rows with bold lettering correspond to summer surveys, which were used for density and spatial analyses. where A is the surface (in km 2 ), L is the total effort (in km), n is the number of sightings,m is the estimate of the effective strip width (esw), andÊ½s is the estimate of the average group size. In MCDS, the probability to detect objects depends not only on the distance from the trackline, but also on additional covariates that are directly incorporated in the detection probability estimation (e.g., [49,50,51]). Provided appropriate covariate data are collected, MCDS can potentially provide more precise estimates than CDS and reduce heterogeneity e.g., when regional differences in relevant environmental variables may differentially affect detection probability and lead to biased estimates if ignored when modelling the detection probability [51]. In MCDS the abundance is obtained using the Horvitz-Thompson estimator as it follows:

Sub-areas
where the symbols are as for (Eq 1) and z i represents the covariates. A global detection function g(y), the probability of detecting an object, given that it is at distance y from the random line [35] was obtained by pooling all the summer data (2009,2010,2013), stratified by area. The use of a pooled detection function is justified by the consistency among the surveys in terms of observers, aircraft and data collection protocols. Two additional detection functions were obtained: Pelagos summer stratified by year; and Central Tyrrhenian summer stratified by year, to compare the density estimates across years for the different areas. The Sardinian Sea was only surveyed in 2010, so its detection function was obtained from the global detection function.
Perpendicular distance data were right truncated prior to the analysis after visual exploration of the frequency histograms of the observations plotted against the perpendicular distances themselves [35]. The truncation distance for all detection functions was set at 400m ( Fig  2); this resulted in elimination of only one observation (at 564 m). A total of 266 observations remained: 21 for the Ionian Sea; 118 for the Central Tyrrhenian Sea; and 130 for Pelagos (Fig  3a, 3b and 3c). The additional explanatory covariates considered to model the detection probability in MCDS were: Beaufort, an assessment of overall sighting conditions and observer, all of them treated as factor variables. For both CDS and MCDS, the final selection of the model was based on the minimisation of the Akaike Information Criterion, AIC [35,52], and examination of the performance of the qq-plot and the goodness of fit tests (chi-square, Kolmogorov-Smirnov and Cramer-von Mises).

Model-based method
With this method, the perpendicular distance data are used to estimate a detection function, which allows abundance to be modelled as a function of physical and environmental data associated with the surveyed transects. Abundance then can be estimated for the entire study area through extrapolation and maps of created densities, using features of the environment to predict abundance, which can increase precision.
For the spatial analysis, a grid of cells was created covering sub-areas A, B and C, characterised according to potential spatial and environmental variables [53,54]. A total of 26,944 cells were generated with a resolution of 17 km 2 . The covariates examined were: 1. spatial-latitude and longitude; 2. fixed-depth, distance to coast, distance to 200, 1000 and 2000m isobaths, slope, and aspect; and 3. dynamic: seasonal and annual means of sea surface temperature (SST), and chlorophyll concentration (Chl).
Depth was extracted from ETOPO [55], and its derivates were obtained using ArcGis 9. The SST and Chl-a values covered the period between 2009 and January 2014 and were obtained from SeaWiFS (between 2000 and 2010) and MODIS-Aqua sensors (from July 2002), and the SST of MODIS-Terra (from 2000) and MODIS-Aqua according to the operational lifespan of these three satellite platforms. Chl-a data from MODIS-Terra were not used due to the current low quality of the data. SST and Chl-a have wide coverage and both are available synchronously on a daily base (at the scale of the processes involved, i.e. within 12 h) at a medium resolution (geo-projected data at 4.6 km for MODIS and 9.2 km for SeaWiFS) [56].
All on-effort transects (Beaufort 3) were classified and divided into segments (mean = 3.1 km; max = 5.9 km) with homogeneous effort types, and under the assumption that little variability in the physical and environmental features occurred within each segment. These segments were subsequently associated to the attributes of the specific cell in which they fell, based on their spatial relationship. Using the count of animals in each segment as the response variable, the abundance of animals was modelled using a Generalized Additive Model (GAM) with a logarithmic link function and a Poisson error distribution. The general structure of the model was: where the offset a i is the effective search area for the i th segment (calculated as the length of the segment multiplied by twice the effective strip width, esw), θ 0 is the intercept, f k are smoothed functions of the explanatory covariates and z ik is the value of the k th explanatory covariate in the i th segment. The esw was obtained from each of the detection functions, according to the spatial model to be created in terms of areas/years. Abundance of animals was modelled directly, rather than using a two-step approach (modelling abundance of groups and modelling group sizes) given the low mean group size (see Table 2); almost 82% of the encounters (237) were of single animals, 9% of 2 animals (25), 8% from 3 to 8 animals (24) and 1% of 16 to 18 animals (2)). Therefore, the few encounters with more than one animal were treated as individual encounters with the same perpendicular distance.
Models were fitted using the R package 'mgcv' version 1.7-22 [57], and manually selected using three diagnostic indicators: (a) the Generalised Cross Validation score, GCV [58]; (b) the percentage of deviance explained; and (c) the probability that each variable was included in the model by chance. The decision to include/drop a term from the model was adopted following the criteria proposed in [59]. The point estimate of total abundance was obtained by summing the abundance estimates in all grid cells over the study area. Finally, to obtain the coefficient of variation and percentile-based 95% confidence intervals, using day as the resampling unit, 400 non-parametric bootstrap re-samples were applied to the whole modelling process. In each bootstrap replicate, the degree of smoothing of each model term was selected by the statistical package, thus incorporating some model selection uncertainty in the variance.

Correcting for biases
In line transect surveys, two main biases-perception bias and availability bias-represent often substantive violations of a primary assumption of the method, i.e., that all animals on the trackline are detected [35], and both can cause an underestimation of abundance. Perception bias, when an observer misses an animal that is available to be seen, is customarily evaluated using an independent observer approach, i.e., with >1 observer independently searching the same area and recording data separately. This was not possible during the present surveys but the limitation of the surveys to good sighting conditions and the use only of experienced observers suggests that perception bias is likely to be small. By contrast, availability bias may be significant, especially for species such as the giant devil rays that unlike cetaceans do not need to come to the surface to breathe [27]. To account for the availability bias, we follow the approach of [25] for the giant devil rays in the Adriatic Sea. The authors of [25] derived a correction factor of 0.49 based upon telemetry data for three animals in the central Mediterranean between June and October 2007, which showed that devil rays spent some 49% of their daylight time (SD = 0.25) at a depth 10m, where they can be sighted under good conditions from the air [34]. Accordingly, corrected abundance estimates and CVs were obtained applying the following formula [60]: Whereâ is 0.49 (from [25].

Species identification
Due to their large size, distinctive body shape (including a very broad head with two prominent cephalic fins), a long, whip-like tail, and greyish dorsal pigmentation with a marked dark "collar" [4], giant devil rays can be unmistakeably identified from the air (Fig 3). Furthermore, Mobula mobular is the only member of the family Mobulidae found the Mediterranean Sea [4], which eliminates the problem of species misidentification.

Ethical statement
No permit or approval was sought for this study, given that no takes were involved and the presence of observers is not considered to cause disturbance considering the fleeting and distant occurrence of the aircraft in the individual animals' space.

Results
The sightings data (n = 298) are summarised in Table 2. As noted earlier, almost 82% of the sightings were of single animals, with only a few sightings comprising larger numbers (up to 18). Rays were sighted in all sub-areas (Fig 1) and reported sightings were more frequent in the summer. The complete absence of sightings of devil rays during the 2009 winter survey over the Pelagos Sanctuary compares remarkably with the frequent sightings made in the same area using the same transects the following summer.

Design-based approach
In examining the fits for the CDS and MCDS, no improvement was found for the MCDS and thus only CDS estimates were taken further. The best model for pooled data and for the Pelagos Sanctuary data was the hazard-rate function with no adjustment terms. For the Central Tyrrhenian Sea the best model was the half-normal with no adjustment terms. Fig 2 shows the perpendicular distance distribution (histograms), with fitted detection functions (lines) for the three models. The hazard rate function allows for a "shoulder" in the closer distances to the track line, indicating that the probability of detecting animals only decreases at certain distance from the track line (after around 200 m in the Pelagos Sanctuary and 80 m in the Tyrrhenian Sea). The abundance estimates are shown in Table 3.

Model-based approach
The best models for all years pooled together (all areas), for 2009 (Pelagos) and for 2010 (all areas) incorporated an interaction between the geographical covariates (Latitude and Longitude) and depth, while the model for 2013 (Central Tyrrhenian Sea) only incorporated Longitude and depth. In all cases covariates were highly significant. The total deviance explained was 14.8% for the global model, 30.5% for 2009, 24% for 2010, and 14.8% for 2013. The abundance estimates are shown in Table 3.
The results from the design-based and the model-based estimates correspond well, with an overall uncorrected abundance estimate for the study area (sum of sub-areas A, B and C) of just over 6, 000 (CV 13.4% for the design-based method and 12.7% for the model based method). Once corrected for availability bias, the total giant devil ray abundance estimate in the study area resulted to be around 12,500, with a CV of around 53% (12, 396, CV = 52.75% design-based; 12,722 CV = 52.57% model-based).
The predicted summer abundance of giant devil rays for the area occupied by the Pelagos Sanctuary, the Sardinian Sea and the Central Tyrrhenian Sea (all years pooled) is shown in Fig  4. The model predicts the presence of areas of high abundance along a wide band connecting the coastline of central Italy with the northeastern portion of Sardinia, spilling slightly through the strait between Corsica and Sardinia into the Gulf of Asinara (northwest Sardinia). Smaller areas of high abundance are predicted off northwest Sardinia, extending to the west with decreasing density, as well as along the coast of southern Latium and Campania. Based on the model, the Pelagos Sanctuary is interested by areas of highest predicted devil ray abundance mostly in its southernmost part, i.e., along its border with the Central Tyrrhenian Sea and off northern Sardinia.

Distribution and movements
In the present study, although sightings were scattered throughout the study area, a summer concentration of rays was apparent along a wide strip connecting the Tuscan Archipelago to Eastern Sardinia, over a wide range of water depths (Fig 1). Giant devil rays are believed to range widely throughout the Mediterranean Sea, both in neritic [2,15,16,17] and offshore waters [34], over depths from few tens of metres to several thousand metres (giant devil rays are epipelagic batoids, which can dive up to at least 600-700 m [34]). However, in addition to the range given in the formal scientific literature, a search of other sources (e.g. images posted on the Internet by marine users), reveals these rays' occurrence in many additional areas (e.g., Crete, Malta, Sicily Channel, Eastern Ionian Sea). This suggests that the species might be distributed more widely than previously thought, which is important when moving towards a full assessment of the species and evaluation of threats.
Strong seasonal changes in the presence of giant devil rays from specific Mediterranean locations suggest that they may undertake regular migrations. Migrations have been reported from several mobulid species [27,31,32]. Our data from the northern Pelagos Sanctuary show regular sightings in summer whereas despite extensive effort in good conditions, no animals were seen in winter ( Table 2). The species is reported to aggregate in the Strait of Messina between late spring and summer [7], but is no longer seen there from around the middle of Table 3. Density and uncorrected abundance estimates for the three sub-areas, obtained with design-based and model-based methods. Esw: effective strip width.

Sub-areas
Year Sightings on track autumn [34]. In late autumn and winter, our data revealed low densities in the southern Tyrrhenian Sea (Table 2). Seasonal occurrence has also been reported for the Adriatic, where reports begin in April-May [16]. Observations of large aggregations of giant devil rays in winter in the southeastern corner of the Mediterranean, off the coast of Gaza [8], also suggests seasonal migrations. One hypothesis is that they migrate between the northern Mediterranean in summer and the southernmost shores in winter. Migratory behaviour of giant devil rays is consistent with these rays' propensity for long-distance movements. The available tagging data for three animals showed movements from the Strait of Messina of 298 km (in 120 days), 337 km (in 120 days), and 278 km (in 60 days) [34]. Further tagging studies throughout the Mediterranean are required to address the question of giant devil ray migration, and how it might relate to seasonal changes in water circulation and primary productivity in the concerned areas. Seasonal movements may be related to possible energetic advantages deriving from warmer waters if this species thermoregulates like other mobulids [61], or to the localized availability of high densities of prey [34] (e.g., mesopelagic and epipelagic fish in the Strait of Messina [7] and epipelagic fish in the central and southern Adriatic [16]).
During summer, the rays' presence inside the Pelagos Sanctuary was higher along its southern boundaries, particularly in the Tyrrhenian Sea. This might be considered surprising given that the deep (down to 1,000m) and relatively productive [62] pelagic waters of the north- western Pelagos Sanctuary are a favoured habitat of one of giant devil rays' prey species, Meganyctiphanes norvegica [63,64], an euphausiid crustacean known to undergo diurnal vertical migrations in excess of 1 km [65]. One possible explanation might be that in the more northern Sanctuary waters, the giant devil rays spend more time at the greater depths frequented by euphausiids during the day, making them less available to sighting from the air.
Another driver of migratory behaviour in giant devil rays could be the species' still unknown breeding habits [66]. Some evidence of sexual segregation is apparent from the catches (n = 30) of adults caught from 21 March to 10 April 2014 within the large aggregations observed off Gaza, 90% of which were males, some with sperm on the claspers [8]. By contrast, all the few (n = 5) specimens caught in the Strait of Messina between 1990 and 2003 were females [7].

Density and abundance
The corrected absolute abundance estimate of >12,000 giant devil rays in the study area of the north western Mediterranean that covers approximately 10% of the Mediterranean is the first abundance estimate of the species in this region, and, despite the discussed uncertainties surrounding availability bias (see [25]), provides a possible basis for a future assessment of the conservation status and trends of M. mobular in a significant portion of its range.
Densities varied both in time (e.g., between successive years in the same area: 0.016 individuals km -2 in 2009 and 0.026 in 2010 in the Pelagos Sanctuary) and in space (consistently highest in the central Tyrrhenian Sea) ( Table 3). The variation in time may reflect a real slight variation in density, or more likely, simply normal non-significant interannual variations due to environmental differences and/or the random effect of the sample size. Variation in space is reflected in the covariates retained in the model, depth being the only one with a clear ecological meaning. The geographic covariates (latitude and longitude) would be proxies for some unknown ecological factors. Both depth and other unknown factors most probably affect the distribution of prey, a major element affecting the spatial distribution of animals. To be able to detect a real trend in densities, a much lower CV would be needed during a longer term density estimates series. A power analysis is recommended to assess whether a real trend can be detected and under which circumstances. Mean density values calculated for the Adriatic Sea (0.018 indiv. km -2 , 95% C.I.: 0.011-0.028, 22.6% CV) were slightly lower than in this study (mean around 0.25) but of the same order of magnitude [25].
Presudoreplication is unlikely to be an issue when pooling estimates from different regions and years. On the one hand, and stated above, movements of M mobular in this region appear to be very limited (around 3.3 km/day on average [34]), which reduces considerably the probability of duplicate sightings. On the other hand, as long as a duplicate sightings occur in different tracks (considered the independent sampling units), distance sampling theory shows that it does not result in bias; var(n) is estimated from variation in encounter rate between independent replicate lines [35]. Given the speed of the aircraft, it is not possible to have potential duplicates on the same sampling unit. Similarly, the fact that same animals may be present for example in Pelagos one survey year and in the Tyrrhenian in another sampling year again is not a problem as they are different sampling units.

Conclusion
The present corrected abundance estimate of >12,000 is comparable (similar methods and time of year) with that for the Adriatic [25] of over 3,200 with similar large CVs (around 57%) which infers a total population size for the two areas of over 15,200. Despite the difficulties with respect to availability bias ( [25]), this information provides for the possibility for the first time of undertaking a more quantitative assessment of the status of the giant devil rays in the Mediterranean. However, considerable uncertainties remain with respect to the biological parameters necessary for such an assessment (see discussion in [25]), as well as information on population structure, abundance throughout its range and levels of past and present removals (directed and bycaught). Information on the original population size e.g. from before the "driftnet era"-is unavailable [25].
Aerial surveys such as these, if carried out regularly, will provide the possibility for developing a comparable relative abundance series and thus an indication of trend and of the need for mitigation measures. However, to do this, further efforts (e.g. in telemetry and genetics, both for better estimates of time at the surface and an understanding of possible migrations and population structure) are required to reduce CV values and determine the proportion of the overall population being covered.
Without a better understanding of the levels of threats such as directed catches and bycatch, it is not possible at this stage to try to evaluate their likely effect on the status of the giant devil ray. Despite its formal protected status by a number of international and national regulations, there has been little systematic research into potential threats (including chemical pollutants and habitat degradation) and mitigation measures, or into monitoring trends in the population. If this protected status is to be meaningful, such information is essential.