Hard time to be parents? Sea urchin fishery shifts potential reproductive contribution of population onto the shoulders of the young adults

Study sites and sea urchin sampling

The reproductive potential of the Paracentrotus lividus populations was examined in two zones of Sardinia that differ in sea urchin fishing pressure (Fig. 1). These two zones were selected because they are at two extremes as regards the exploitation of sea urchins in Sardinia. Su Pallosu Bay, located along the Sinis peninsula (central-western Sardinia, 40°03′N; 008°25′E), is subjected to very high pressure (HP zone) which is widespread across the entire bay. Harvesting can be practised through scuba diving from November to April by 189 professional fishermen authorized by regional decree. Each one of them, when helped by an assistant, is allowed to collect up to 3,000 sea urchins (TD ≥ 50 mm) per day, according to the regional decree n. 1967/DecA/67 issued by the Regione Autonoma Sardegna (RAS, Regione Autonoma della Sardegna 2013).

On the contrary, the marine protected area of Tavolara—Punta Coda Cavallo (north-eastern Sardinia, 40°53′N; 009°40′E) is considered a low-pressure zone (LP zone). This marine reserve was legally established in 1997 and has been considered well-enforced since 2003 (Di Franco et al., 2009). This is one of the few zones in Sardinia where sea urchin fishing is strongly restricted (there are only few authorized fishermen and no more than 500 catches per day for an overall maximum of 300,000 sea urchins harvested each year) and the populations are well preserved (Guala et al., 2011). Although no data are available on annual catches, Su Pallosu Bay is among the most exploited in the whole region and the theoretical number of sea urchins that can be legally harvested in one year (calculated according to the number of licenses, allowed catches and fishing days) is over 300 times greater than that of Tavolara—Punta Coda Cavallo.

Sea urchin specimens were collected from two areas of high harvesting pressure (A and B) at Su Pallosu Bay, and low harvesting pressure (C and D) at Tavolara MPA. In each zone, the two areas were 1–3 km apart. Sampling (approved by the Regione Autonoma Sardegna through the release of the fishing license for scientific purposes n. 9727/AP SCIE/N.7 03/06/2016) involved all sea urchin size classes that have gonads and contribute to the reproductive potential of the population. Specifically, once a month, we collected 10 individuals of the commercial size (CS, TD ≥ 50 mm), five undersized individuals (US, 40 ≤ TD < 50 mm) and five smaller, undersized sea urchins (Small-US, 30 ≤ TD < 40 mm) for an entire year.

The sea urchins were collected from both zones by scuba diving over a rocky bottom covered in photophilic algal communities at a depth of five meters. Sea urchin samples were wrapped in cloth soaked in salt water, stocked in iceboxes and immediately transported to the laboratory.

Environmental features

Despite the differences in position along the Sardinian coast and in wind exposure, sampling was performed in areas with similar environmental features (i.e., depth, slope, shelter from the waves), on rocky bottoms.

Differences in food availability, potentially able to affect sea urchin gonad growth and consequently fecundity (Byrne, 1990; Minor & Scheibling, 1997; Brady & Scheibling, 2006; Scheibling & Hatcher, 2007), were assessed as temporal changes of algal cover between zones. Digital photographs were taken over three PVC quadrats of 50 × 50 cm which were randomly placed on the sea bottom in each area. This was done three times (July 2013, January and May 2014) during the surveyed year. Image analysis was carried out using Seascape^® software (Segmentation and Cover Classification Analyses of Seabed Images, Teixidó et al., 2011) to detect the percentage cover of conspicuous algal taxa or morphological groups. A permutational multivariate analysis of variance (PERMANOVA, Anderson, 2001a) was done, on the basis of a Bray–Curtis dissimilarity matrix calculated from square-root transformed data (Primer-E 6 Permanova^®), to estimate the variability of assemblages between the two zones over time. A 3-way model was used with Time (random, 3 levels) and Zone (fixed, 2 levels, Su Pallosu Bay vs. Tavolara—Punta Coda Cavallo) as crossed factors, and Area (2 levels) as a random factor nested in Zone. P-values were obtained through Monte Carlo random draws from the asymptotic permutation distribution (Anderson, 2005) and a pairwise test was used to discriminate between various levels of significant factors. A non-metric Multi-Dimensional Scaling (nMDS) ordination was used as a graphical representation of data.

In order to investigate any potential thermic anomalies, seasonal variations of the coastal water temperature were monitored with the “Mediterranean Sea—High Resolution and Ultra High Resolution L3S Sea Surface Temperature” product (http://marine.copernicus.eu/web/69-myocean-interactive-catalogue.php). Daily sea surface temperatures (SST) were extrapolated from the catalogue after choosing an intermediate point between the two sampling areas within the zones, then mean monthly temperatures were obtained and used to represent the annual trend.

Gonadosomatic Index and fertility

The gonadosomatic index (GSI) was examined every month from June 2013 to May 2014 (with the exception of November due to adverse weather conditions) for all of the three sampled size classes. Sea urchins were allowed to drip dry for some minutes and then weighed. The test diameter (TD) without spines was measured and the gonads were successively extracted and weighed as well.

GSI was calculated by the formula: GSI = [gonad wet weight/total wet weight] × 100 as reported by Lawrence, Lawrence & Holland (1965).

The fertility was tested according to Brundu et al. (2016) during the months of maximum gonadal development (i.e., from December to April). Gonads from females were extracted and gently shaken in filtered seawater to allow the mature ova to come out. Sperm was then added and finally, the fertility of the individuals was assessed as a percentage of the effectively fertilized eggs after the appearance of the fertilization membrane (fertilization was considered successful if it took place in at least 80% of the eggs according to Falugi & Angelini, 2002). The percentage of achievement of the first larval stage (development of four-arm echinopluteus) was also measured. The fertility test was carried out on the US individuals but not on Small-US specimens because of the paucity of their gonads (see “Results”). Nor was it done on the commercial-sized individuals (CS) because they were assumed to be fertile (Ouréns, Fernández & Freire, 2011).

Monthly mean values of GSI were calculated for the two zones and areas of sampling. Since sex ratio was nearly 1:1 for both sampling zones (Data S8), females and males were pooled together to obtain a single mean GSI value per month. A four-way ANOVA (using Statistica 6.0, Statsoft Inc.) was performed to highlight the differences in GSI values in different months (fixed and orthogonal factor, 11 levels), zones (fixed and orthogonal factor, 2 levels), areas (random and nested in zones, 2 levels), and fertile size classes (fixed and orthogonal factor, 2 levels, CS and US). Eight replicates for each size class were haphazardly selected from among those available in order to get a balanced design. Cochran’s C test was used to check for the assumption of homogeneity of variances and a posteriori SNK tests were performed to find alternative hypotheses (Underwood, 1997).

Population structure and potential reproductive contribution

Abundance and size-frequency distribution of the populations were estimated for both zones by counting all sea urchins found in the PVC plots and measuring them with calipers. The plots were 50  × 50 cm and were placed randomly as many times as necessary to cover two replicates of 25 m² in each area (100 m² for each zone).

Sea urchin abundance was estimated as the total density (individuals m⁻²) and the density of each 10 mm size class: TD < 10 mm, 10 ≤ TD < 20 mm, 20 ≤ TD  < 30 mm, 30 ≤ TD < 40 mm, 40 ≤ TD < 50 mm and 50 ≤ TD < 60 mm, TD ≥ 60 mm. Size-class densities were then translated into frequency percentages and used to compose the population structure for each zone. The gamete output and the spawning magnitude of each spawning event were calculated for the fertile CS and US classes according to Brewin et al. (2000). The highest and the lowest mean monthly GSI recorded during the year of sampling corresponded to the period just before the beginning (pre-spawning) and after the end (post-spawning) of the spawning events. We defined the mean individual gamete output (IGO), in units of gamete wet weight per urchin per spawning event (g g⁻¹ se⁻¹), as the difference of the mean monthly pre-spawning GSI and the mean monthly post-spawning GSI. The spawning magnitude was defined as the percentage ratio of the mean individual gamete output and the mean monthly pre-spawning GSI.

For each spawning event, we calculated the gamete output per m² (GO, g g⁻¹ m⁻² se⁻¹), released by fertile size classes in both zones, as their IGO multiplied for the respective natural density. The total gamete output (TGO) and the mean gamete output (MGO), used to estimate the reproductive contribution of each size class per m² per year, were defined respectively as the sum and the average of GO (g g⁻¹ m⁻² yr⁻¹). Finally, the potential reproductive contribution of the whole population in both zones (popTGO and popMGO, respectively) were calculated for the investigated year as the sum of the contributions of both fertile size classes.

Environmental features

Multivariate analysis on algal cover has detected a significant interaction between Time and Zone. More specifically, pairwise tests highlighted that algae where sea urchins live changed significantly over time but not between the two zones (Table S1). The nMDS showed that plots of the HP zone were interspersed on the graph with those of the LP zone (Fig. S2).

Moreover, the water temperature changed over time but no difference was observed in the annual trend between the two zones studied (Fig. S3).

Gonadosomatic Index and fertility

At each area, 220 specimens of the commercial-size class (CS) and 110 undersized and smaller undersized individuals (US and Small-US) were randomly collected over a year to compare reproductive potential between populations. Regarding CS individuals, sampled sea urchins with TD ≥ 60 mm were 24 of 440 (5.5%) in the HP zone and 219 of 417 (53%) in the LP zone.

Fertility tests showed that 100% of the undersized individuals (US) checked were fertile and contributed to the reproductive potential of the populations, with high percentages of fertilized eggs (ranging from 87 to 96%) and developing larvae (ranging from 95 to 100%) (Table S4). Conversely, the Small-US individuals had reduced gonads (Fig. 2) and, even during the period of maximum development, they never released gametes, therefore their contribution to the reproductive potential of the population can be considered negligible.

The GSI trend over the year was generally higher for CS individuals than US ones for both zones (Fig. 2). At the HP zone, we recorded a single large spawning period from March to May for both fertile classes (see Fig. 2A). In pre-spawning time, GSI values reached 6.6 ± 0.3% (mean ± standard error) and 4.4 ± 0.4%, while in post-spawning time values were 1.3 ± 0.2% and 1.6 ± 0.1% for the CS and US classes respectively (Fig. 2A).

At the LP zone, a spawning event was observed twice (Fig. 2B). The first one was recorded from June to December with a pre-spawning GSI of 6.7  ± 0.3% and 5.2 ± 0.5%, and a post-spawning GSI of 2.5  ± 0.3% and 1 ± 0.2% for the CS and US classes respectively. The second spawning event occurred from February to April with a pre-spawning GSI of 5.4 ± 0.4% and a post-spawning GSI of 2.5 ± 0.2% for the CS class. From February to May we observed a pre-spawning GSI of 4 ±  0.6% and a post-spawning GSI of 1.2 ± 0.2% for the US class (Fig. 2B).

ANOVA highlighted significant differences during the year and between the different size classes, while there were no major statistical differences found between sampling areas and zones. The CS class had higher GSI values in the LP zone while no differences were found for the US individuals (Table 1). A significant interaction between Zone and Size Class was detected, and the GSI of the CS class was significantly higher than that of the US individuals for both zones. A significant interaction was also found between Month and Size Class with SNK pointing out significantly higher GSI values for the CS individuals during the whole sampling year excluding April (Table 1).

Population structure and potential reproductive contribution

Size class distribution was consistently different between the two zones (Fig. 3). In the HP zone, sea urchin density was almost two-fold higher than in the LP zone: 10 ± 1.4 and 5.4 ± 0.5 individuals per m² respectively. Sea urchins with TD ranging from 0 to 20 mm were 1.7 ± 0.1 per m² and 1.1 ± 0.6 per m² in the HP and LP zones respectively, and they represent 17% and 21% of their populations. The most abundant size classes were those ranging from 20 to 50 mm diameter (77%) with a density of 3 ± 0.6 and 3 ± 0.3 individuals per m² for the Small-US and US classes respectively (30 ≤ TD < 40 mm and 40 ≤ TD < 50 mm). The proportion of individuals of the CS class with respect to the entire population was 6% (0.6 ± 0.2 individuals per m²) and all the individuals were included in the range of 50 ≤ TD < 60 mm (Fig. 3). In the LP zone, the individuals between 20 and 50 mm represented 28% of the population while the CS class was 52% with 1 ± 0.4 and 2 ± 0.2 individuals per m² for 50 ≤ TD < 60 mm and TD ≥ 60 mm respectively (Fig. 3).

In relation to the population structure, the reproductive contribution was compared between the two zones, but no comparisons were made between the sampling areas because no differences were found (see Table 1). The potential reproductive contribution was calculated according to the number of spawning events during the surveyed year and the natural density of the fertile size classes (US and CS). In the HP zone, a single spawning event occurred and the spawning magnitude for the year was 73% on average (Table 2) with an individual gamete output (IGO) of 0.03 and 0.05 g g⁻¹ se⁻¹ for the US and CS individuals, respectively (Table 2). The gamete output (GO), calculated in relation to the natural density of this zone, was 0.08 g g⁻¹ m⁻² se⁻¹ for the US sea urchins (2.7 ± 0.3 individuals per m²) and 0.03 g g⁻¹ m⁻² se⁻¹ for the CS class (0.6 ± 0.2 individuals per m²). Because of the single spawning event, the total gamete output per m² (TGO) overlaps the mean gamete output per m ² (MGO) (Table 2). Accordingly, the total gamete output of the whole population (i.e., sum of TGO of the two fertile size classes) corresponded to the mean gamete output for m² (popMGO) with a value of 0.11 g g⁻¹ m⁻² yr⁻¹ (Table 2).

Conversely, in the LP zone, two spawning events were observed (Fig. 2B). Spawning magnitude varied from 54 to 81% with higher values for US individuals (Table 2). IGO was similar for both size classes with values ranging from 0.03 to 0.04 g g⁻¹ se⁻¹ according to the spawning period. The GO of the US individuals was 0.02 g g⁻¹ m⁻² se⁻¹ during the first spawning event and 0.01 g g⁻¹ m⁻² se⁻¹ during the second one, while it was 0.11 and 0.08 g g⁻¹ m⁻² se⁻¹ for the CS class. Total gamete output (TGO) and mean gamete output (MGO) of the US individuals, whose density was 0.4 ± 0.1 individuals per m², were 0.03 and 0.01 g g⁻¹ m⁻² yr⁻¹ respectively. Meanwhile they were 0.19 and 0.10 g g⁻¹ m⁻² yr⁻¹for the CS class, whose natural density was 2.7 ± 0.3 individuals per m². Consequently, the total gamete output of the whole population was estimated to be 0.22 g g⁻¹ m⁻² yr⁻¹ and the total mean gamete output was 0.11 g g⁻¹ m⁻² yr⁻¹ (Table 2).