Identification and Characterization of Nursery Areas of Red Mullet Mullus barbatus in the Central Tyrrhenian Sea

Abstract

Red Mullet Mullus barbatus is an important target of fishing activities in the central Tyrrhenian Sea, so it is essential to identify its critical habitats in order to manage this resource efficiently. Our research specifically focused on the identification and characterization of nursery areas. The use of spatial interpolation techniques enabled us to identify five nurseries that were highly persistent through time. Moreover, the estimate of juvenile density confirmed the strong aggregation effect of these nursery grounds, as a great portion of young individuals were concentrated in a relatively small surface of the study area. The environmental characterization of these areas showed that juveniles were mainly distributed on bottoms with a relatively high percentage of sand (>70%; P < 0.05). Shannon biodiversity index analysis indicated that the southern nurseries reached the highest values of habitat quality (P < 0.0001). Multivariate analysis showed that nursery grounds were divided into three main groups, and analysis of spatial dynamics showed that two different strategies characterized Red Mullet juveniles when density changes over time. In particular, in some areas young individuals selected habitats in a density-dependent way following the basin model scheme, while in other zones they selected habitats in a density-independent way according to the proportional density model. Results also showed that juveniles followed the proportional density model strategy into nursery areas with the highest Shannon biodiversity index values.

Received September 7, 2016; accepted January 24, 2017

Red Mullet Mullus barbatus inhabits sandy and muddy bottoms up to a 200-m depth (Machias and Labropoulou 2002; Tserpes et al. 2002). This species is distributed all around the Mediterranean basin, including the Black Sea, and also in the eastern Atlantic, from Senegal to Scandinavia (Fischer et al. 1987).

Red Mullet is an important element of the demersal assemblages living on the Mediterranean continental shelf (Demestre et al. 2000; Tserpes et al. 2002; Colloca et al. 2004). Moreover, this species has a high commercial value, so it represents one of the main targets of fishing activity operating in the Mediterranean Sea (Relini et al. 1999; Tserpes et al. 2002). Due to the commercial importance of Red Mullet, several studies have been carried out on Red Mullet biology or population dynamics (Fiorentino et al. 1998; Voliani et al. 1998; Cherif et al. 2007; Sabatés et al. 2015). Some studies have also analyzed Red Mullet spatial distribution (Abella et al. 1996; Tserpes et al. 2002; Carlucci et al. 2009), while others have focused on the effects of environmental parameters on its distribution (Lombarte et al. 2000; Machias and Labropoulou 2002; Levi et al. 2003; Maravelias et al. 2006). Few studies, however, specifically dealt with the identification of Red Mullet critical habitats as nursery grounds (Garofalo et al. 2004, 2011; Carlucci et al. 2009). Moreover, no study was previously carried out on the spatial dynamics of this species, especially in these aggregation areas.

Concerning the life cycle of a fish population, nursery areas represent “essential habitats,” as they host a high concentration of juveniles testing several advantages such as a lower predation rate or an increase in available food. These characteristics improve the survival of young specimens and support their fastest growth. The introduction of the nursery habitat concept (Beck et al. 2001) enabled the classification of different habitats in relation to their importance as nursery areas. According to Beck et al. (2001), a habitat could be defined as a nursery for a particular species if it contributes a greater than average number of individuals to the adult population on a per unit area basis compared with other habitats used by juveniles. However, Dahlgren et al. (2006) highlighted the importance of other habitats used by juveniles that fall outside the previous nursery definition but that contribute significantly to the adult population. Consequently, the concept of “effective juvenile habitat” was introduced. This is defined as a habitat for a particular species that contributes a greater proportion of individuals to the adult population than the mean level contributed by all habitats used by juveniles, regardless of the area covered. Nevertheless, according Sheaves et al. (2006), these concepts were incomplete because it is important not only to identify nursery grounds but also to evaluate their persistence through time. In this way, it is possible to establish which areas can really be considered as nurseries. Moreover, the abundance of juveniles is not always dependent on habitat suitability, but it also depends on local conditions (i.e., oceanographic patterns) which aggregate larvae in specific zones (Agostini and Bakun 2002).

The aim of this research was to identify and characterize Red Mullet nursery areas in the central Tyrrhenian Sea, where no studies were previously carried out. Furthermore, in these areas we also evaluated models describing the spatial dynamics of juvenile populations under different annual densities. Density dependence factors are fundamental to studying the dynamics of fish populations, and, in this regard, intraspecific competition is able to determine a change in the way a specific habitat is exploited by individuals. Therefore, the aim of our work was also to describe the spatial relationship between Red Mullet juveniles and their habitats when the abundance of juveniles changes over time.

Since Red Mullet juveniles are mainly distributed in coastal waters, where human pressure is particularly high, an appropriate knowledge of nursery areas becomes very important to managing this species successfully.

Methods

Gruppo Nazionale Demersali (GRUND) trawl survey (Relini 2000) data (historical series 1994–2006) were employed to analyze the spatial distribution of Red Mullet nurseries in the central Mediterranean Sea. We used a GRUND data set because this survey was generally carried out in autumn, which usually represents the recruitment period of this species. Moreover, this data set was limited to 2006 because GRUND trawl survey was finished in this year. The study area of the GRUND project covered all zones subjected to trawl fishery from a depth of 10 to 800 m (Figure 1), sampled according to a random stratified sampling scheme based on five depth strata: 10–50 m, 51–100 m, 101–200 m, 201–500 m, and 501–800 m. The number of samplings was proportional to the surface of each stratum, and, from 1994 to 2006, 585 stations were sampled. However, in this research we considered only areas with a depth ranging from 10 to 100 m (42°15′N, 11°06′E; 41°07′N, 13°39′E), where young individuals were usually most abundant (Voliani et al. 1998; Machias and Labropoulou 2002). In this bathymetric range, during the above-mentioned period, 138 stations were sampled.

Figure 1. The study area. Black circles represent hauls of GRUND project.

The first step in identifying nursery areas was to define what we considered as juveniles. For this study, we considered juveniles to be all individuals that had finished their larval phase but still had not reached first sexual maturity. Sizes at first sexual maturity (L50) were calculated using the data series 1994–2006 (separately for males and females), according to the equation

which described a logistic curve. In this function, ST_L represented percentage of mature specimens for each size-class, T_L was total length of specimens, and S₁ and S₂ were two numerical parameters whose ratio (S₁/S₂) indicated L50 values. The calculation of L50 values was fundamental to distinguishing between juveniles and adults.

Identification of nurseries

The identification of nursery areas was carried out through a spatial distribution analysis of Red Mullet juveniles. Thus, georeferenced density (N/km²) indexes were employed and spatial interpolation techniques were used over the entire study period (1994–2006) by means of ArcGIS 9.3.

There are two main groupings of spatial interpolation techniques: deterministic and geostatistical. Deterministic interpolation techniques create surfaces from measured points, based on either the extent of similarity (inverse distance weighting [IDW]) or the degree of smoothing (radial basis functions). Geostatistical interpolation techniques (kriging) exploit the statistical properties of the measured points (Johnston et al. 2001). In particular, geostatistical techniques quantify the spatial autocorrelation among measured points and provide an account of the spatial configuration of the sample points around the prediction location. As geostatistics is based on statistics, these methods produce not only prediction surfaces but also error or uncertainty surfaces, giving an indication of how good the predictions are (Johnston et al. 2001).

We used spatial interpolation techniques to create a continuous surface for a parameter, the numeric density index (N/km²), starting from values measured at georeferenced sampling points. Employing these values, spatial interpolation techniques made it possible to predict values of this variable at nonsampling points. This occurred through geostatistical (kriging) and deterministic interpolation (IDW).

Kriging techniques are methods based on spatial autocorrelation shown by a variable. This property indicates the degree of relationship that exists between two or more spatial data of a variable, quantifying a basic principle of geography: things that are closer are more alike than things farther apart (Johnston et al. 2001). Thus, the values of a variable at unobserved locations depend on values observed at the nearest sampling points. Usually, it is common practice to limit the number of measured values that are used when predicting the unknown value for a location by specifying a search neighborhood. The specified shape of the neighborhood restricts how far and where to look for the measured values to be used in the prediction, and it is also influenced by the input data and the directional influences on these data (anisotropy). Moreover, kriging techniques make it possible to assess spatial autocorrelation of a variable through the empirical semivariogram, which displays the statistical correlation of nearby data points. In order to do this, it is necessary to fit a model to the empirical semivariogram (Johnston et al. 2001).

We used indicator kriging to assess where juveniles were annually aggregated. This geostatistical technique predicts, for a study area, the probability that values of a variable are higher than a fixed threshold (Isaaks and Srivastava 1989). In particular, if values are above the threshold, they become 1, and if they are below the threshold, they become 0 (Johnston et al. 2001). As these values are 0 or 1, the spatial interpolation of sampling points produces values between 0 and 1. Specifically, indicator kriging was employed to predict the probability that the values of the density numeric index (N/km²) were higher than a specified threshold value. Since we are not able to define a specific biological threshold to identify aggregation areas of Red Mullet juveniles, we set this threshold at the 75th percentile of logarithmic distribution of density numeric indexes (log N/km²) to include only areas with an important aggregation of young individuals (Figure 2). Once we specified this threshold, abundance numeric indexes in the sampling points were converted in indicator values. The specified search neighborhood was an ellipse, with the major axis approximately parallel to the shoreline because Red Mullet juveniles are mainly aggregated in the shallow waters. Moreover, the spherical model was fitted to the semivariogram to analyze the spatial autocorrelation of the indicator values. It showed a progressive decrease of spatial autocorrelation until a distance of approximately 10 km, beyond which autocorrelation was zero. Finally, the continuous surfaces created by indicator kriging were probability maps which illustrated expressly the probabilities (included between 0 and 1) that the considered variable (N/km²) was higher than the fixed threshold. This enabled us to identify where juveniles were annually aggregated. To validate these results, we also evaluated the prediction error statistics.

Figure 2. Estimation of threshold value set at 75th percentile (black arrow) of logarithmic distribution of density indexes (log N/km²) to include only areas with important aggregations of young individuals.

In order to define the main location of Red Mullet nursery areas, we analyzed persistence through time of annual aggregation areas, as temporal stability represents a fundamental requirement in designating them as nursery areas. As mentioned above, probability maps showed values between 0 and 1, exceeding the defined abundance index threshold. The specification of five probability thresholds (0.3, 0.4, 0.5, 0.6, and 0.8) on these probability surfaces was used to obtain presence–absence maps according to these different thresholds tested. In particular, in these presence–absence maps value 1 identified juvenile annual hotspots, i.e., areas where probabilities were higher than the investigated threshold probability, and value 0 indicated the other zones of the study area. Finally, the overlap of these presence–absence maps, by means of a spatial analysis technique, enabled us to identify nursery areas as a function of their persistence through time and according to the five different scenarios tested.

Inverse distance weighting was used to estimate density indexes values (N/km²) in the study area divided in cells of 1 km². Inverse distance weighting has been employed in this study because it is an exact interpolator, where the maximum and minimum values in the interpolated surface occurred only at sampling points. The overlay of density index values with probability map values made it possible to estimate how many Red Mullet juveniles (shown as a percentage) were aggregated in the nurseries area, according to the different scenarios analyzed in this study.

The inverse distance weighting interpolation technique determines cell values using a linearly weighted combination of a set of sample points and weight as a function of inverse distance. In order to predict a value for any unmeasured location, IDW uses the measured values surrounding the prediction location, assuming that each measured point has a local influence that decreases with distance. This technique weights the points closer to the prediction location greater than those farther away, hence the name IDW (Johnston et al. 2001). In particular, weights are proportional to the inverse distance raised to the power value, p. Therefore, as the distance increases between the measured sample locations and the prediction location, the weight that the measured point will have on the prediction will decrease exponentially (Johnston et al. 2001). How fast the weights decrease is dependent on the value of p. As for kriging interpolation, it was possible to specify a search neighborhood. We used an ellipse, with the major axis approximately parallel with the shoreline, also in IDW interpolation.

Characterization of nurseries

The second step of our research was to characterize nurseries by evaluating what models could describe the spatial dynamics of fish under different population densities. First, we assessed how some environmental parameters, such as benthic community and sediment composition (considered as percentage of sand), affected juvenile distribution. These two parameters were selected given the biological features of Red Mullet, as this species is a benthic carnivore (Lombarte et al. 2000). In particular, for each nursery ground we estimated the mean number of juveniles associated with a specific benthic community and sediment type (expressed as percentage of sandy component). As mentioned above, the estimation of the number of Red Mullet juveniles for a unit area (1 km²) was made by means of IDW. As a result, these values were related to a specific benthic community and sediment composition, and, consequently, the average number of juvenile specimens was so calculated for each of these parameters. We then estimated the Shannon biodiversity index to define the habitat quality of these areas. This species richness index was calculated for each sampling of the entire data set (1994–2006) by using the species composition of catches, according to the equation

where p is the proportion (n/N) of individuals of one particular species found (n) divided by the total number of individuals found (N), ln is the natural log, Σ is the sum of the calculations, and s is the number of species. The spatial interpolation of these values was made in order to estimate the biodiversity index values for the whole study area divided, as mentioned above, in cells of 1 km². In this way, we obtained Shannon index values for each single cell, and, consequently, we calculated the mean value of this biodiversity index for the nursery grounds.

In order to verify how nursery areas were related to each other, we used data multivariate analysis. Specifically, we employed the Bray–Curtis similarity index (Clifford and Stevenson 1975), so a similarity matrix between the nursery areas was calculated using square root transformation, and cluster linkage was determined by group average sorting. Then, hierarchical agglomerative classification (cluster) and ordination by nonmetric multidimensional scaling (MDS) based on the Bray–Curtis similarity matrix were performed. These analyses were carried out through the PRIMER statistical package (Clarke and Warwick 1994).

Successively, we analyzed how juvenile populations responded, from a spatial dynamic point of view, to different annual densities in these nursery areas. Usually, when a fish population increases, some individuals settle in previously unoccupied marginal habitats to elude intraspecific competition. This produces an extension in the spatial distribution of a fish population. The presence of this relationship in the marine fish population has been considered as evidence of density-dependent habitat selection (Marshall and Frank 1995; Shepherd and Litvak 2004; Hiddink et al. 2005), and it has been attributed to a specific ecological model of spatial dynamics—the basin model (see below). However, this is not the only response to an increase in total abundance that we could observe. The nature of the relationship between local density, total abundance, and the area over which a fish population is distributed can assume different forms (Shepherd and Litvak 2004). Nevertheless, these are fundamentally explained by three ecological models (Table 1) of spatial dynamics of marine species (Hilborn and Walters 1992; Petitgas 1998). The first, the proportional density model, is appropriate when the area occupied by a fish population stays constant over time and the local density varies proportionally to the abundance (Petitgas 1998). In this model, there is no relationship between abundance and area occupied by individuals. Moreover, habitats show different carrying capacities in response to the abundance variation. This model has been referred to as density-independent habitat selection (Marshall and Frank 1995).

Table 1. Properties of models describing the spatial dynamics of marine species (adapted from Shepherd and Litvak 2004).

Relationship	Proportional density model	Constant density model	Basin model
Abundance–area	Independent	Positive	Positive
Abundance–maximum density	Positive	Independent	Positive
Abundance–average density	Positive	Independent	Positive

The second model, the constant density model, predicts that a fish population increases its spatial range in relation to total abundance, but maximum density appears constant throughout the spatial distribution of the fish population. In this model, the relationship between abundance and area occupied by individuals is positive. Under the constant density model, habitats show common carrying capacities. Furthermore, this model has been attributed to density-dependent habitat selection (Shepherd and Litvak 2004).

The third model, the basin model, predicts that both density and area are positively related to abundance. In the basin model, density changes in all areas with abundance as a result of a positive relationship between habitat suitability and local density. Therefore, an extension of fish population spatial distribution and an increase in local density is observed. As stated previously, this model has been referred to as density-dependent habitat selection. In our current research, we tried to identify what models could better describe the spatial dynamics in Red Mullet nursery areas in response to annual changes in juvenile total abundance. The identification of the spatial dynamic model followed by Red Mullet in these areas was carried out by applying and analyzing linear regressions to the three relationships reported in Table 1. The methodological approach could be described as follows. First, for each nursery we located, we related the annual total abundance of Red Mullet juveniles to the area occupied year by year by each nursery. Specifically, the area was calculated by counting how many cells of 1 km² were annually included into each nursery. As a consequence of the number of cells constituting a nursery, juvenile total abundance in each nursery was estimated by adding together the values of density indexes (N/km²) calculated for each single cell by means of IDW. As shown in Table 1, if these two variables were independent of each other, the spatial dynamics model followed by juveniles in a nursery area would be the proportional density model. However, when the linear regression between total abundance of Red Mullet juveniles and the area occupied by a nursery was positive, it became very important to analyze the other relationships as well (see Table 1). Total abundance of juveniles was thus related to the maximum density reached by these specimens in a nursery and to the average density observed within the same nursery ground. If both these interactions were positive, then the spatial dynamics model would correspond to the basin model. Conversely, if these variables were independent of each other, then we would observe the constant density model as the spatial dynamics model in a nursery area. To validate these results, ANOVA tests were performed. Lastly, we checked if the observed spatial dynamics were affected by habitat quality by analyzing the values of the Shannon biodiversity index calculated for the different nursery areas.

Results

Since the length at first maturity (L50) was 10-cm TL for males and 12-cm TL for females, we used these sizes to distinguish juveniles from the adult population.

We set the threshold value at the 75th percentile of logarithmic distribution of density numeric indexes (Figure 2). The value observed corresponded to 3,300/km², so sampling stations with the number of individuals equal to or higher than 3,300/km² were considered important aggregation areas for Red Mullet juveniles. Then, the specification of different probability thresholds on the probability maps, created by the indicator kriging technique, enabled us to evaluate different scenarios. Thus, we obtained presence–absence maps which were overlaid by means of the spatial analysis technique, making it possible to identify Red Mullet nursery areas according to the different scenarios proposed. In particular, persistence analysis showed five nursery grounds along the study area sited mainly between 10- and 50-m depths (Figure 3, a–e). The first (named C1) is located nearby the Cape Linaro promontory; the second and the third are positioned, respectively, north (C2) south (C3) of the Tiber mouth; the forth is located off Cape Circeo (S1); and the last is positioned off the Gaeta Coast (S2). When we tested 0.3 as the probability threshold, the surface occupied by Red Mullet nurseries was higher than when we investigated the other probability thresholds (from 0.4 to 0.8). Specifically, when testing 0.8 as the probability threshold, we noted that nursery S2 became very small. However, changing probability thresholds, the total surface occupied by the nurseries did not vary in a consistent way (Figure 4), as it ranged from approximately 1% to 8% of the whole study area. At the same time, the estimation of juvenile abundance showed that a high percentage of juveniles were concentrated in this relatively small surface (Figure 4; Table 2). For example, using 0.5 as the probability threshold, we observed that about 60% of juveniles were aggregated in about 4% of the whole study area (Table 2). Moreover, we observed a low increase in juveniles into these areas for probability thresholds smaller than 0.5 (Figure 4). Last, we used IDW to estimate the mean density of juveniles per unit area (number of individuals within a cell of 1 km²) for each nursery ground. In particular, S1 and S2 were the areas with highest mean density of juveniles (over 10,000 individuals/km²), as shown in Figure 5.

Table 2. Estimate of young individuals (expressed as percentage) in nursery areas according to five different probability thresholds. The surfaces occupied by nurseries were reported as percentage of the study area.

Probability thresholds	Area (percentage of the study area)	Percentage of juveniles
0.8	0.89	23.25
0.6	2.87	48.02
0.5	4.34	58.34
0.4	6.09	65.42
0.3	8.35	71.17

Figure 3. (a–e) Nursery areas of Red Mullet in the central Tyrrhenian Sea, according to the different scenarios investigated. Persistence analysis substantially showed five nursery zones along the study area.

Figure 4. Relationship between the percentage of study area occupied by nurseries and the percentage of juveniles aggregated in these zones, according to the five different probability thresholds.

Figure 5. Mean density abundance (N/km²) in the study area.

The second step of our research was to describe the nursery grounds from an environmental point of view. According to the coastal seabed maps (SIBM 2004), we identified four benthic communities: (1) biocoenosis of well-sorted fine sands (SFBC), (2) biocoenosis of coastal terrigenous muds (VTC), (3) biocoenosis with intermediate features between these two (SFBC–VTC), and (4) biocoenosis of the muddy detritic bottoms (DE). These biocoenoses were not equally distributed in Red Mullet nursery areas (Figure 6), and it was only in nursery S2 that we observed all the biocoenoses identified above. Thus, we noticed a heterogeneous pattern of juvenile distribution in relation to the benthic community. Moreover, depending on the estimation of juvenile density, we did not observe a pattern of juvenile distribution strictly related to a particular benthic community except for nursery C3, where juveniles were significantly distributed (ANOVA: F_2.267 = 14.238, P < 0.05) on SFBC.

Figure 6. Distribution of juveniles into nursery areas in relation to the benthic communities (top panel) and the percentage of sandy component of the sediment (bottom panel).

With regard to sediment composition, we considered the percentage of sand contained on superficial sediments and identified five categories (Figure 6), where the first category included sediments with a very small percentage of sand (0–5%) and the fifth, sediments with a very high percentage of sand (>95%). Except for C1 and S1, all five categories were observed in Red Mullet nurseries (Figure 6). We also usually noticed that juveniles were more greatly distributed where bottoms had a high percentage of sand (>70%). While nursery C3 juveniles were significantly distributed on these bottoms (ANOVA: F_4.265 = 12.884, P < 0.05), only in nursery C1 was sediment composition characterized by a small percentage of sand (<30%).

A Shannon biodiversity index was estimated for the entire study area. With regard to Red Mullet nursery grounds, we found that this index was absolutely highest in the southern areas S1 and S2 (ANOVA: F_4.605 = 576.62, P < 0.0001). Consequently, we considered these zones as having the best habitat quality among all nurseries for M. barbatus identified along the study area. We did observe, however, that nursery C1 also had Shannon index values significantly higher than that for nurseries C2 and C3 (ANOVA: F_2.496 = 342.73, P < 0.0001).

The Bray–Curtis index showed similarity between nursery areas (Figure 7), and subsequently, cluster analysis and MDS identified three separate groups: one composed of nursery C1 (group 1), one of nurseries S1 and S2 (group 2), and one of nurseries C2 and C3 (group 3). Group 1 showed a low similarity percentage with groups 2 and 3 (47.9%), while groups 2 and 3 were more similar to each other, reaching a similarity value of about 64%. Moreover, regarding group 2, we observed nurseries S1 and S2 to be very similar, with a similarity percentage of 80%; for group 3, we also observed that nurseries C2 and C3 reached a high level of similarity (77%).

Figure 7. Clustering of nursery areas by means of multidimensional scaling.

We sought to understand how juveniles responded to different densities in the nursery grounds, and we substantially detected two different spatial dynamics models (Table 3). In particular, in the zones located near the mouth of the Tiber River (C2 and C3), we observed that an annual increase in young individuals’ abundance caused an expansion of spatial distribution of Red Mullet juveniles and at same time an increase in local density. In these zones, juveniles followed spatial dynamics typical of the basin model (Table 3). In the nursery areas south of the study area (S1 and S2), however, we observed that the surface occupied by juveniles did not change according to annual abundance variation, while local density varied proportionally to juvenile total abundance. These spatial dynamics are typical of the proportional density model (Table 3). In response to an increase in annual total abundance, we did not observe an expansion of spatial distribution, but only an increase of local density, and consequently we assumed an increase of habitat carrying capacity. This model also described the spatial dynamics observed in the nursery C1.

Table 3. Spatial dynamics models followed by juveniles in nursery areas in response to different densities. Models were described analyzing the relationships shown below.

	Relationship			Model
Nursery areas	Abundance–area	Abundance–maximum density	Abundance–average density	Spatial dynamics
C1	Independent (p > 0.05)	Positive (p < 0.05)	Positive (p < 0.01)	Proportional density
C2	Positive (p < 0.01)	Positive (p < 0.05)	Positive (p < 0.01)	Basin
C3	Positive (p < 0.05)	Positive (p < 0.05)	Positive (p < 0.01)	Basin
S1	Independent (p > 0.05)	Positive (p < 0.01)	Positive (p < 0.01)	Proportional density
S2	Independent (p > 0.05)	Positive (p < 0.05)	Positive (p < 0.01)	Proportional density

Shannon biodiversity index analysis, as mentioned previously, showed that the highest values were registered in nursery areas S1 and S2, where Red Mullet juveniles followed spatial dynamics typical of the proportional density model.

Discussion

We analyzed the spatial distribution of Red Mullet juveniles in order to accurately identify nursery grounds and to characterize them according to biological parameters. Moreover, in these areas we described spatial dynamics as a response to the annual fluctuation of juvenile total abundance.

Nursery area identification was carried out through the geostatistical technique called indicator kriging, a geostatistical model that has also been used by others to identify nursery grounds (Carlucci et al. 2009). In particular, nursery areas of Red Mullet, European Hake Merluccius merluccius, and deepwater rose shrimp Parapenaeus longirostris were identified in the southern Adriatic Sea by means of indicator kriging. On a larger scale, this specific geostatistical approach was also used to pinpoint European Hake nursery grounds in Italian waters (Murenu et al. 2010), and other researchers have employed different kriging methods (Lembo et al. 2000) or techniques such as geostatistical aggregation curves (Colloca et al. 2009). However, the common point of all these studies was the employment of abundance thresholds to identify aggregation areas of the investigated species.

The first consideration regarding the identification of Red Mullet nursery grounds was the high persistence of aggregation areas, which indicated environmental conditions quite constant through time. As observed by Colloca et al. (2009), nurseries that are highly persistent through time should provide a fish population with a larger input than areas which sporadically function as nurseries. Thus, temporal stability plays an essential role and should be an indispensable reference point for designing protected areas, as emphasized by Early et al. (2008).

The estimation of juvenile abundance in the nursery areas showed that a large fraction of juveniles were concentrated in relatively small surfaces of the study area, confirming the tendency to strong aggregation displayed by young individuals; this has also been described in other research mentioned previously (Lembo et al. 2000; Carlucci et al. 2009; Colloca et al. 2009). Moreover, regarding the five scenarios investigated in this study, we observed that the probability threshold of 0.5 was a good compromise to approach an appropriate management of Red Mullet juveniles in the study area. The exclusion of these zones from fishing activities through the implementation of marine reserves could therefore protect a notable fraction of juveniles without excessively reducing fishing areas, not to mention the fact that trawl fishing activity is already forbidden at depths <50 m, where Red Mullet nurseries are mainly located. Nevertheless, in our opinion, the complete closure of nursery areas to trawl fishing activity and their monitoring are also necessary to prevent the illegal trawling which is quite common in these zones. Such measures, as also mentioned in other studies (Fogarty 1999; Gell and Roberts 2003; Roberts et al. 2005), could address the exploitation of marine resources and direct efforts towards sustainability. A positive effect on the population of M. barbatus, deriving from the exclusion of a marine zone from fishing activities, was noted in a study conducted in a zone of Sicily (Castellamare Gulf), where fishing activities have been prohibited for 14 years (Fiorentino et al. 2008). While this research specifically showed a general increase of juvenile abundance and spawner biomass, another study (Machias and Labropoulou 2002) revealed that the exclusion of fishing activities in coastal zones to protect juveniles was not a sufficient measure to preserve the entire fish population from overfishing. However, the identification of nursery areas represents, without a doubt, an important element in the management of marine resources from a spatial point of view. The need to protect these habitats depends on the necessity to reduce the impact of fishing activities on commercial species, especially in the critical phases of the life cycle such as recruitment, where individuals are aggregated and more vulnerable.

As reported by some (Vinagre et al. 2006; Le Pape et al. 2007), the spatial distribution of juvenile fish is generally explained by environmental parameters. In this regard, the analysis of some factors, such as sediment composition and benthic communities, showed how Red Mullet juveniles generally preferred areas with a high percentage of sand, as was typical of the shallowest waters in this study area, but a clear trend was not observed regarding benthic community preference, except for nursery C3. As suggested by Machias and Labropoulou (2002), it is possible in our study that since juveniles appeared to have a relatively broad trophic niche, it was difficult to relate them only to a specific benthic community. Conversely, for other fish species, some studies have emphasized the importance of benthic fauna assemblage in the spatial distribution of juveniles (Nagelkerken et al. 2000; Letourneur et al. 2003; Nicolas et al. 2007).

Multivariate analysis essentially showed that nurseries could be divided into three different groups. It is probable that dissimilarity between group 1 (nursery C1) and the other two groups could be explained in terms of sediment composition and benthic community. Nursery C1 was characterized by the absence of a well-sorted fine sands community and by a low percentage of sand into sediments composition. In contrast, differences detected between group 2 (nurseries S1 and S2) and group 3 (nurseries C2 and C3) could be due mainly to what we observed about species richness values (Shannon index) and mean density of juveniles. Notably, in nurseries S1 and S2 the values of these parameters were significantly higher than those calculated in nurseries C2 and C3.

Analysis of spatial dynamics in nursery areas showed that Red Mullet juveniles followed two different strategies when density changes over time. In some areas (C2 and C3), they selected habitats in a density-dependent way following the basin model scheme. This strategy is probably linked to the need to avoid intraspecific competition inside areas where a high number of juveniles are present and resources, such as food availability, are not sufficient to sustain a local increase of their abundance. So, in response to this occurrence, some individuals settled in previously unoccupied marginal habitats of inferior quality, thereby expanding, in this way, their spatial distribution. The same methodological approach of this study enabled previous descriptions of the basin model for other species distributed all around the world, for example, the anchovy (family Engraulidae) and the sardine (family Clupeidae) in the southern Benguela and off Japan, California, and Perù (Barange et al. 2009), or the flatfish (family Pleuronectidae) in the eastern Bering Sea (Spencer 2008). The basin model scheme has also been reported for Haddock Melanogrammus aeglefinus (Hiddink et al. 2005). In nursery areas C1, S1, and S2, Red Mullet juveniles selected habitats in a density-independent way following a different spatial dynamics model (the proportional density model). In these areas, spatial distribution remained unvaried despite an increase in juvenile total abundance, so we assumed an increase in habitat carrying capacities. This strategy suggested that advantages offered by some local factors, such as lower predation rate or higher food availability, were greater than possible disadvantages caused by intraspecific competition for some of these resources. Also, as indicated by Shannon index values, the proportional density model was linked to high habitat quality, which could better support juvenile survival and growth. Consequently, we concluded that in nurseries S1, S2, and C1, quality played an essential role in habitat selection, so juveniles remained in these highly suitable habitats despite an increase in total abundance. The proportional density model has also been observed in some gadoids and Pleuronectiformes (Myers and Stokes 1989; Petitgas 1994, 1997). In conclusion, we observed that in nurseries C2 and C3, the spatial dynamics of Red Mullet juveniles were substantially controlled by density-dependent factors, while in nurseries S1, S2, and C1, spatial dynamics were mainly affected by density-independent factors. A similar effect was also observed in Atlantic Cod Gadus morhua (Tamdrari et al. 2010), where spatial dynamics were controlled by density-dependent or density-independent factors as a function of different zones.

Our study showed that an analysis of spatial dynamics in aggregation areas (i.e., nurseries) could be very important in evaluating what protection measures should be applied in these zones. This is a key element to effectively managing marine resources such as Red Mullet, which are subject to a high fishing effort, but further exploration of the particular efforts needed to efficiently protect a specific resource is necessary. In this regard, research development could concern the analysis of connectivity between juveniles and spawners according to a spatial population model. It is also important to understand how juveniles move toward adult and spawning areas in order to manage, in an efficient way, all of the critical habitats of a fish population, not only a part such as nursery grounds. Nevertheless, the study of nursery areas is important for the conservation of marine resources because it suggests that habitat management should be included in fishery management plans, as also recommended by the ecosystem approach to fisheries management.

Acknowledgments

We thank the Italian Ministry of Agriculture, Food, and Forestry (Fisheries Department) for the authorization to use the bottom trawling survey data (GRUND), without which this paper would not have been possible. We also thank Cristina Fontana for improving the English language quality of this manuscript.