Data were collected during a survey for which the primary target species were cetaceans and sea turtles. Given the ecological characteristics of the giant devil ray, two strata were chosen: (1) the entire Adriatic Sea (133000 km²); and (2) the “central and southern Adriatic”, i.e. the pelagic portion of the Adriatic Sea (about 88400 km²). Both strata were covered by 20 km-spaced parallel transects with a random start point (Fig. 1) providing equal coverage probability using the Distance 6.0 software (Buckland et al. 2004, Thomas et al. 2010; http://www.ruwpa.st-and.ac.uk/distance/). Survey design and data collection followed the standard line-transect distance sampling methodology (Buckland et al. 2004).

The survey took place between 31 July and 16 August 2010 when weather conditions were expected to most likely allow coverage of the predetermined track-lines under acceptable conditions for the target species (i.e. sea states of ≤3 and good visibility). A twin engine high-wing aircraft (Partenavia P-68) equipped with two bubble windows was used as the research platform. The observer team consisted of three researchers: two searching for animals through opposite lateral bubble windows and a third recording data onto a semi-automated logging system. Based on experience from surveys elsewhere, the plane flew at an altitude of approximately 200 m (650 ft) with a ground speed of about 185 km/h (100 knots)—precise height and speed was obtained from the GPS records. Only ‘on-effort’ encounters (searching along the predetermined track-lines under acceptable conditions) were used for the final abundance estimation. The distance between transects (20 km), survey speed and the average distance covered daily by giant devil rays (2.5-3.7 km/day, Canese et al. 2011) minimize the risk of double counting animals/groups.

Data recorded, in addition to the GPS data, were (1) time of sighting when abeam; (2) species identification (discussed below); (3) school size; (4) declination angles measured using Suunto PM-5/360PC inclinometers on the observed animals when abeam (to allow calculation of perpendicular distance of the sighting from the track-line); (5) environmental information (including sea state, turbidity, cloud cover, etc.); (5) effort (e.g. whether plane is 'on-effort', i.e. circling to confirm data, or 'off-effort', i.e. transiting or in unacceptable conditions) and positional/time data; and (6) approximate estimated disc width (DW) category ("0.5-1 m", "2-3 m" and ">3 m") using average lengths of sea turtles, striped dolphins and bottlenose dolphins as a reference.

Conventional distance sampling (CDS) analysis was applied and data were analysed using the Distance 6.0 software. In summary, density was estimated from:

$$\widehat{D}=\frac{n\overline{s}}{2\widehat{w}L}$$

where n is the number encounters on-effort; s, the school size; L, the length of transect; and w, the effective strip half width.

The effective strip width was estimated by fitting a detection function to the distribution of perpendicular distances of sightings from the track-line. Since the reliability of the detection function analysis is dependent upon sample size, we chose to use all available sighting data for Myliobatiformes. Myliobatiformes other than M. mobular are somewhat smaller and more difficult to detect and their inclusion will not introduce any major bias into the detection function and hence w. The best model fit was chosen using Akaike’s Information Criterion (AIC), as proposed in Buckland et al. (2004). Total uncorrected abundance ($\widehat{N}$ ) was then calculated from the density estimate and the surface of surveyed area.

Total abundance corrected ($\widehat{N}$ ) for availability bias (â) and its coefficient of variation

$$CV({\widehat{N}}_c=\frac{\widehat{N}}{\widehat{a}})$$

were calculated as follows: $${\widehat{N}}_c=\frac{\widehat{N}}{\widehat{a}}$$ and

$$CV({\widehat{N}}_c)=\sqrt{CV{({\widehat{N}}_c)}^2+CV{(\widehat{a})}^2}$$ .

Whenever possible, Myliobatiformes were identified to species. In the Adriatic, adult giant devil rays (DW=2-4 m) are easily identifiable from the air at 200 m in sea state ≤3 (Fig. 2). In addition to the dark stripe on the head and the two obvious cephalic fins extending beyond its head, this species is characterized by a white lateral ventral colouration apparent when it is filter feeding. When observed dorsally, it appears to be a bluish black to brown colour, with a distinctive darker and lighter pattern between the eyes.

Smaller animals (DW≤2 m) may include young giant devil rays, bull rays (Pteromylaeus bovinus), common eagle rays (Myliobatis aquila) or the dasyatid pelagic sting ray (Pteroplatytrygon violacea), which are also present in the Adriatic Sea (Fortuna et al. 2010). The bull and common eagle rays, although generally similar in shape to the giant devil ray, are smaller and lack cephalic fins (adult DW 100-170 cm and 40-80 cm, respectively; Jardas et al. 2004, Dulčić et al. 2008). The sting ray is also smaller (80 cm maximum) and has a more discoid shape (Bauchot 1987, Mollet 2002).

Abundance estimates are presented for (1) adult giant devil rays in central-southern Adriatic; and (2) all Myliobatiformes combined, in the entire Adriatic Sea.

The CDS assumes that all animals on the track-line are seen (commonly stated as g(0)=1). This is rarely the case in any survey and certainly not the case for animals that do not spend all their time at the surface. There are two primary categories of bias affecting g(0): (1) availability bias, where the diving behaviour of the animal renders it unavailable to be seen by the observers; and (2) perception bias caused by observers missing animals that in principle are available to be seen (i.e. they are at the surface as the plane passes over). However, if appropriate data are available these biases can be accounted for and thus an unbiased estimate of absolute abundance can be obtained. It should be noted that if a time series of relative abundance estimates is required (e.g. to estimate trends), then the g(0) issue can be ignored if it can be assumed that it does not change over time.

Availability bias was corrected using time depth recorder data from Canese et al. (2011). who tagged three animals within the Tyrrhenian Sea and Sicily Channel between June and October 2007. These animals spent approximately 49% (SD=0.25) of their time between 0 and 10 m depth. In the Adriatic Sea in summer, the limpidity of the northern offshore Italian waters normally ranges between 2 and 10 m (Casellato et al. 2005, Anonymous 2006). On the Croatian side it normally ranges from 11.7 to 17.5 m in summer 2010 (IOR 2011). These water characteristics allow for the potential sighting of animals below but near the surface. As an initial correction for the availability bias, 0.49 (SD=0.25) was applied. While the availability bias correction appreciably increases the point estimate, it also substantially increases the CV.

Perception bias can be corrected by using independent observer approaches (e.g. Buckland et al. 2004). However, the limitations of the sighting platform prevented collection of independent observer data. Not correcting for this bias will result in an underestimate of true abundance, to an unknown degree.

The sparse information on the biology of the giant devil ray is summarized in Table 1. Disc width maximum size and average DW length at birth were 520 cm and 152 cm, respectively. Length at maturity is unknown as is frequency of reproduction (a single pup is believed to be the norm). As for many fish species, natural mortality (M) is unknown.

As an initial step towards developing a more quantitative assessment approach, growth was assumed to follow the von Bertalanffy (1938) growth curve:

$$L_t=L_{\infty }\left(1-{\exp }^{\left(-k\left(t-t_0\right)\right)}\right)$$

where k is the growth rate, L_∞ is the asymptotic length and t₀ is a parameter of the von Bertalanffy growth formula, which represents the time at which the size is equal to 0. We therefore used L₀, which can be derived from t₀, i.e.

$$t_0=\log \left(\left(1-\frac{L_0}{L_{\infty }}\right)/k\right)$$

In the following analyses all lengths and length-based parameters (e.g. L_∞) were based on disc width (DW). Empirical studies have shown that life history parameters, such as L_∞, k, maturity and natural mortality, are correlated and can be used to obtain estimates for data-poor stocks (e.g. Roff, 1984, Jensen 1996, Reynolds et al. 2001). Therefore, due to the lack of data on life history parameters for M. mobular, length-based parameters were derived by fitting linear regressions to data from related species and families. L_∞ was derived from maximum size (Fig. 3A) and k from L_∞ (Fig. 3C) using data on the families Myliobatidae and Rajidae in FishBase (Froese and Pauly 2012) and length-at-maturity from L_∞ (Fig. 3B) using data on the Mobula genus from Couturier et al. (2012). These relationships (see linear regressions fitted in Figure 3A, B and C) provide estimates for L_∞, length at maturity and k of 880 cm, 320 cm and k=0.09 per year, respectively. Owing to the large size of M. mobula, in some cases this required extrapolation, but the extrapolations did not result in excessively large confidence limits and so were not rejected (Gislason 2010). Life history parameters also allow an approximate rate of natural mortality to be derived from life history ‘invariants’ (e.g. Beverton and Holt 1959, Jensen 1996), i.e.

M=1.5k

where k is the growth coefficient from the von Bertalanffy growth equation. Methods for estimating approximate size-dependent natural mortality vectors are also available (e.g. Peterson and Wroblewski 1984, Lorenzen 2000). However, an increase in M₀ and a reduction in fecundity have the same effect, i.e. they reduce recruitment, so to simplify this already exploratory analysis we considered recruitment, rather than M₀ and fecundity separately.

The population growth rate (r) was estimated in the standard manner from the Leslie Matrix (A), a transition matrix with columns representing age classes (Caswell 1989), i.e. N_t+1=AN_t, where N_t is a vector describing the

$$A=\left(\begin{array}{ccccc}{s}_0{f}_1& s_0{f}_2& \cdots & s_0{f}_{m-1}& s_0{f}_m\\ s_1& 0& \cdots & 0& 0\\ 0& s_2& \cdots & 0& 0\\ ⋮& \cdots & \ddots & 0& 0\\ 0& 0& \cdots & s_{m-1}& 0\end{array}\right)$$

r can be derived from λ, the dominant eigenvalue of A:

r = log (λ)

For a given value of r the population doubling time is:

T_d = log(2)/log (1+r)

Since nothing is known about the breeding frequency of M. mobular, 1.5 years, the mean value provided for the closely-related Manta birostris (Dulvy et al. 2008) was used in this exploratory approach.

Assuming the relationships described above, we conducted a sensitivity analysis in which parameters were varied in turn, to examine the importance of the levels of uncertainty for the various parameters.

When designing monitoring programmes, it is important to determine what magnitude of population changes, should they occur, is detectable by proposed survey approaches over given time scales. Power analyses can be undertaken to evaluate the survey CV and frequency required to be reasonably certain to detect a given change in population size (Gerrodette 1987). Here we chose to examine a significance level of 5% and a detection power of 60% based on conventional practice. Survey CVs considered were 20% and 30%.

Figure 1 shows the transect lines covered on-effort (Beaufort ≤3).

A total of 68 sightings (56 on-effort) of all Myliobatiformes were obtained. The mean group size of Myliobatiformes was 1.1 animals (SD=0.49; range=1-4). Within this, 42 (39 on-effort) were confirmed as adult giant devil rays. The mean group size was 1.2 animals (SD=0.6; range=1-4).

Adult giant devil rays were mostly (88 %) observed in the central and southern Adriatic Sea, whereas sightings of other “Myliobatiformes” were more frequent (65%) in the northern portion of the basin (Fig. 1). The distribution of sightings was more clumped offshore in the central stratum than in the southern stratum.

Of the 23 (17 on-effort) ‘other’ Myliobatiformes encounters, only 7 could be identified to the level of family (6 myliobatids and 1 dasyatids).

Figure 4 shows the distribution of perpendicular distances for all Myliobatiform sightings combined. There were fewer sightings close to the track-line than would be expected but the maximum sighting distance was relatively small, with a severe drop in sightings after about 236 m (at 40° declination angle). The best fit was to a uniform model with a truncation distance of 225 m; this effectively equates to a strip census. The variance in density was affected more by group encounter rate (90%) than by group size (10%), as apparent from the clustered distribution of sightings (Fig. 1). The uncorrected CDS abundance estimates for confirmed adult giant devil rays for the central and southern Adriatic are shown in Table 2, along with the uncorrected estimate for all rays for the whole Adriatic Sea. When we applied the availability correction factor (0.49; SD=0.25) to the uncorrected estimate (N=1595; CV=0.23) for confirmed adult giant devil rays, the estimated abundance increased to 3255 with a substantially larger CV(CV=0.56).

Table 3 summarizes the results of the sensitivity analysis for the exploratory population modelling exercise. The ‘base case’ gave a population growth rate of 0.036 and a population doubling time of 19.5 years. Increasing M or k has a similar effect since M is a function of k, i.e. population doubling time increases to 135 years. The analysis shows that the population growth rate and doubling time are very sensitive to the assumed parameters, i.e. a 25% increase in recruitment leads to a 7 fold increase in r.

Figure 5 shows the results of the power analysis for two levels of survey CV (20% and 30%) and two survey frequencies (1, 2 and 3 years). CVs of this level are plausible if uncorrected estimates are considered as relative abundance estimates, and if it can be assumed that availability and perception bias remains constant over the survey periods investigated. The curved lines represent a power of 0.6 to detect a trend in the population for each combination of survey CV and frequency. For a given r (y-axis), the time taken to detect a significant change in the population can be found from the intercept of the corresponding horizontal line with the power curves. For example, for a CV of 20% and an r of 0.1 with surveys every 1, 2 or 3 years it would take 6, 7 and 10 years, respectively, to detect a change in the population. In contrast, detection of a growth rate of –0.1 would take longer, i.e. 7, 10 and 14 years, respectively. With growth rates of between 0.05 and –0.05 it appears difficult to detect any change in the population within a reasonable timeframe.

On 14 August 2010, a group of four giant devil rays was observed northeast of the Otranto strait. Three animals had an estimated DW of about 3 m and one appeared to be smaller (Fig. 2). Animals were observed for about 6 minutes swimming in anticlockwise circles. During this observation the small individual and one of the larger animals engaged in belly-to-belly swimming and flipping over while circling. Such behaviour in Myliobatiformes is usually connected with courtship and mating (Tricas 1980, Yano et al. 1999).

The present study is the first to provide a more quantitative and systematic description of the summer distribution of the giant devil ray (and information on some other Myliobatiformes) in the Adriatic. The results can be used to consider its distribution within certain biogeographic partitions of a Mediterranean sub-region: i.e. coastal, deep pelagic and intermediate ecosystems (Bianchi 2007). We recognize that this survey represents a seasonal snapshot for a single year. However, this is the first step towards having sufficient data to undertake a spatial modelling approach to better understand the key oceanographic and ecological factors that influence its distribution or an ERA to assess the impacts of fisheries.

Data from the present survey (Fig. 1) generally agreed with the information from earlier opportunistic data with more animals in the central and southern parts of the Adriatic (Notarbartolo di Sciara et al. 2006, Holcer et al. 2013).

Sighting records for other Myliobatiformes [i.e. the benthopelagic bull and common eagle rays (Bauchot 1987)] were more common in the shallower northern Adriatic, where local studies often detected them at high densities (Fortuna et al. 2010, Manfredi et al. 2010). It should be noted that these may include some sub-adult Mobula as referred to earlier under species identification.

The northwestern Adriatic generally has more turbid waters (e.g. Justić 1988, Bernardi Aubry et al. 2004). It is possible that this turbidity could have diminished the ability of observers to spot the key morphological features used for identification of the giant devil ray, especially for younger specimens. While in principle this may render the observed spatial segregation between giant devil rays and other Myliobatiformes (Fig. 1) an artefact, we do not believe this to be the case. Given the number of sightings identified to species and the clear differences in average size, we believe that the observed difference in distribution pattern is likely to be real and reflect the ecology of this species.

Both eagle and bull rays feed mainly on bivalves, gastropods and demersal crustaceans (Jardas et al. 2004, Dulčić et al. 2008). Despite the limited direct information on giant devil ray diet, the genus Mobula is believed to have a highly specialized diet, feeding on macro-zooplankton and predominantly euphausiids (Orsi-Relini and Cappello 1992, Sampson et al. 2010). This fits well with the distribution data available from the survey (Šipos 1977).

As concerns the quantitative estimate of summer abundance of giant devil rays (and indeed Myliobatiformes) in the Adriatic Sea, the limitations of the study have been described and explain the high CV for the availability-corrected estimate. The estimate is only partially corrected and underestimates the actual abundance to an unknown extent for two main reasons: (1) the lack of data to correct for “perception bias” (which may be high given the need for good condition to see and identify them); and (2) the fact that possible unidentified sub-adult giant devil rays are among the 26 smaller ray sightings. The partially corrected estimate of 3255 (CV=0.56) may be considered to challenge the existing general perception of the ‘rarity’ of the giant devil ray. In absolute terms, a middle estimate of over 3000 animals would not be considered ‘rare’ as suggested in Holcer et al. (2013). Note that while the species may not be considered ‘rare’, this does not mean that the population has no conservation problems.

In addition, it should be noted that given the total surface area of the central-southern Adriatic (about 88400 km²), our estimate represents a maximum density of about 4 adults per 100 km². It is difficult to evaluate how this density relates to the IUCN qualitative evaluation that this species is present in ‘very low densities’ throughout its range (Notarbartolo di Sciara et al. 2006). Finally, any full consideration of status would need some knowledge of ‘pristine’ densities.

The giant devil ray is ‘protected’ in most Mediterranean countries and thought to be of limited commercial value (Bauchot 1987). That being said, as discussed above, there is information suggesting a recent directed fishery (e.g. in Israeli and Palestinian waters) and, possibly considerable ongoing incidental catch with various types of fishing gear for other target species.

There is considerable uncertainty over the status and present trends of the giant devil ray population, despite its presence in the appendices of various conventions and its ‘Endangered’ IUCN Red List classification. While the present study suggest that the species may not be not as rare as previously thought (the ban on pelagic driftnets since the early 1990s may have been successful in promoting population recovery), it is not known whether the population(s) is/are increasing, decreasing or static under the present conditions. As yet there is no research or management regime in place to evaluate this. If commitments to international agreements are to be kept, progress must be made, especially in the light of newly identified source of mortalities (e.g. Couturier et al. 2013). In addition, there is no information on what the original size of the population was in the 1970s before the “driftnet era” and no agreement on an appropriate conservation objective e.g. a biomass reference point corresponding to ‘good conservation status’ or an agreed desired (and reliably measured) trend in abundance.

A fuller assessment of the impact of fishery-related mortality should be carried out by conducting an ERA in which the susceptibility analysis is stratified by management sub-regions and fishing métiers. To conduct the ERA and evaluate the likely and actual success of any mitigation management measures requires: (1) improved information on recruitment; (2) a time series of relatively precise indices of abundance; (3) information on stock structure (if any) and the geographic and temporal distribution of the species; (4) robust estimates of removal rates from potentially impacting fisheries; and (5) appropriate reference points.

An interesting approach to examining recruitment (i.e. taking into account ‘natural’ biological parameters and anthropogenic removals) is to monitor population indices directly through, for example, aerial surveys. The present study has shown that this is possible in principle for the Adriatic Sea. While effective monitoring of population trends in devil rays is possible, this will require a long-term commitment to a monitoring programme and aerial surveys are probably the most cost-efficient approach. Given the observed densities of the giant devil ray, an increase in the research effort (e.g. by increasing the number of transects) will decrease the CVs of the estimates and thus shorten the time needed to detect trends should they occur.