Can the Empty Shells of Pinna nobilis Maintain the Ecological Role of the Species? A Structural and Functional Analysis of the Associated Mollusc Fauna

Abstract

The endemic Mediterranean bivalve Pinna nobilis is a typical example of an autogenic ecosystem engineer in sedimentary environments, as it modifies the substrate, creates new complex 3D habitats and concentrates various species, thus increasing local biodiversity. Despite the significant loss of P. nobilis populations due to the epidemic mass mortality that occurred in several regions of the Mediterranean, including the northern Adriatic Sea, the numerous empty shells, still firmly anchored in the sediment, may continue to play the role of basibiont of this species. The persistence of this ecological function was here investigated on dead specimens of P. nobilis through the structural and functional study of the associated epibenthic mollusc community. A large number of mollusc taxa, characteristic of both infralittoral and sublittoral hard and soft bottoms, were identified in thirteen shells collected from four different sites in the Gulf of Trieste. Multivariate analysis revealed a clear separation between sites in terms of taxonomic composition, species abundance and biomass, with differences in expressed functional traits. Overall, considerable taxonomic and functional diversity was found, suggesting that the empty shells continue to provide a suitable habitat for a variety of molluscs that are closely linked to the surrounding environment and fulfil many different ecological functions.

1. Introduction

In marine environments, ecosystem engineering and the related biogenic habitat generation is now recognised as a form of interaction that plays an important role in structuring biological communities [1,2,3]. It is a process mediated by organisms, called ecosystem engineers, that physically modify, maintain and/or create habitats, while also modulating resource availability [4]. It has been observed that the overall net positive effect of this interaction is reflected in the 25% increase in landscape-level species diversity [1,3] and, in particular, in the 53% increase in species richness observed in marine ecosystems [3], where benthic invertebrates (e.g., bivalves, corals, polychaetes, ascidians and barnacles) play an important role [3,5]. Bivalves are a typical example of autogenic ecosystem engineers in sedimentary environments, as they alter the substrate, create new complex habitats and introduce heterogeneity through their own structure [4,5]. Molluscs that form aggregations (e.g., mussel beds and oyster reefs) or live solitarily (e.g., Pinnidae) [6], create an available hard substrate that can be colonised by other benthic species and provide refuge from predators [2,7,8], leading to an increase in local diversity. The fan mussel Pinna nobilis (Linnaeus, 1758) [9] is an endemic ecosystem engineer species of the Mediterranean Sea, because its large shell, which can be over 120 cm long and buried up to one third in the sediment, provides new suitable substratum and habitats for a diverse epibenthic community [10,11,12,13]. Moreover, as a filter-feeder species, it ingests large amount of phytoplankton, zooplankton, particulate, suspended organic matter (POM and SOM, respectively) and seagrass detritus, thus playing a key role in energy transfer from the water column to the sediments and even improving local environmental water quality [14,15,16,17]. The enormous and long-lasting exposure to direct and indirect anthropogenic pressures over decades (e.g., overexploitation, habitat degradation and illegal fishing (trawling)) has led to serious declines in P. nobilis populations in the Mediterranean Sea. Therefore, it has been necessary to include this species in a full protection regime, firstly, under the Annex II of the Barcelona Convention (SPA/BD protocol 1995) and the Annex IV of the EU Habitats Directive (92/43/EEC) and secondly even in the Marine Strategy Framework Directive (2008/56/EC) as a Mediterranean species that deserves special attention [17]. However, after a consistent population recovery favoured by the protection regime, the species faced an epidemic mass mortality event (hereafter MME) that began in the early autumn of 2016 on the Spanish coast [18] and moved eastward, affecting, one after another, the densest fan mussel populations along the French coast, the Tyrrhenian coasts of Italy, the Greek and Turkish coasts and the Adriatic coast of Croatia, Bosnia and Herzegovina, finally reaching the Gulf of Trieste in the northern Adriatic in 2019 [18,19,20,21,22,23,24,25]. The numerous molecular and histological analyses then revealed that the mortality rate of almost 100% observed in most areas of the Mediterranean was mainly caused by the new protozoan species Haplosporidium pinnae [19] alone, and in combination with Mycobacterium sp. [20,25,26]. Following the MME, P. nobilis was included in the International Union for Conservation of Nature (IUCN) “Red List” of “critically endangered species”, with the aim of protecting the few unaffected populations or surviving specimens, that may be crucial for conservation and restoration efforts.

Recently, several monitoring efforts and censuses that were conducted to determine the health status of the fan mussel after the MME have revealed a discrete resistance of populations living in lagoons and transitional waters [27,28,29], while very few and rare individuals have been found alive in the sea, also in the northern Adriatic (S. Ciriaco pers. comm., 2022). However, it has been observed that the numerous shells of different sizes, now devoid of living P. nobilis animals, but still upright and perfectly anchored in the sediment, continue to act as a hard substrate suitable for the settlement and growth of a varied sessile faunal community and as a centre of aggregation and shelter for motile faunal species. Consequently, both sessile and motile species inhabiting and colonising dead P. nobilis shells may form a spatial association that can be defined as “sclerobiosis”. This term, recently proposed by Romero et al. [30], is derived from the newly coined term “sclerobiont” by Taylor and Wilson (2003) [31] which, in contrast to the term “epibiont” [32], refers to all organisms, such as encrusting, sessile, motile and boring species that inhabit any form of hard (but not necessarily living) substrate. The high abundance of sclerobionts seems to support the assumption that the role of an ecosystem engineer can persist over time despite the death of the animals, with significant long-term positive effects on local biodiversity [33]. While there is a fairly extensive literature on the role of bivalves as ecosystem engineers and the assessment of their associated macrobenthic fauna, including the fan mussel, little research has been done on the structure and composition of the epibenthic community of P. nobilis, especially in the Adriatic Sea, and in particular, no research at all has been done in the Mediterranean on how dead specimens could still contribute to structuring macrofaunal communities and enriching diversity after the MME. Biodiversity, expressed both in the taxonomic diversity of species and in the diversity of functions of these species, is an intrinsic property and value of natural communities and also contributes greatly to the most important services that ecosystems provide to humans and society [34]. It has been widely discussed that biodiversity is strongly related to ecosystem stability and ecological functioning [35,36,37,38,39]. Each species in a community is defined by a combination of morphological, behavioural, biochemical and phenological traits that reflect its interaction with the environment and its ecological role in the ecosystem [40,41,42,43]. Species with different traits are expected to perform different functions in ecosystem processes. Therefore, trait diversity in a community can be used as a measure of functional diversity (FD) [44] and is useful to describe the role of organisms in communities and ecosystems [45]. Similar to taxonomic diversity, functional diversity can be defined by specific metrics such as functional richness (FRic), functional evenness (FEve), functional divergence (FDiv) and functional dispersion (FDis) [46,47]. They attempt to characterise the extent of the trait space occupied by the community, the degree to which the abundance and/or biomass of a community is distributed within the trait space, the predominance of extreme species and functional differentiation in the trait space and the variability and position of species trait groups in the space trait, respectively. In recent years, both the FD and trait-based approaches combined with studies of taxonomic diversity have proven to be valid tools for understanding and assessing how communities and ecosystems respond to a particular environmental perturbation [40,41]. Variations in species composition and interactions can lead to changes in community structure, ecosystem functioning and associated processes. In this sense, the study of functional diversity in synergy with taxonomic diversity can shed light on how species and their interactions respond to spatial changes and environmental gradients. Therefore, this study aims to assess the remaining ecological role of empty shells of P. nobilis as ecosystem engineers by structurally and functionally investigating the associated epibenthic community. As a first preliminary investigation, this study aims to obtain a qualitative and quantitative characterisation of the mollusc community associated with dead specimens of the fan mussel in the Gulf of Trieste (northern Adriatic), with a framework also from a functional point of view, in order to confirm the hypothesis that the ecological role of these bivalves persists despite their death and that a stable microecosystem still emerges from them.

2. Materials and Methods

2.1. Study Area

The Gulf of Trieste is a shallow semi-enclosed basin (max depth 25 m, average depth 17 m) in the northernmost part of the Adriatic Sea (Mediterranean Sea), extending from the mouth of the Tagliamento River to Savudrija/Punta Salvore (Croatia) and covering an area of about 750 km² with a coastline of about 100 km [48,49]. The basin is characterised by the largest tidal amplitudes and the lowest winter temperature in the Mediterranean Sea, reaching values below 10 °C [50], and by great salinity (37 psu on average) and temperature seasonal variations [51]. Hydrodynamic conditions are mainly characterised by a general circulation pattern clockwise in the surface layer and counter-clockwise in the lower layer, which, in turn, is strongly influenced by prevailing winds from eastern quadrants, such as the Bora [52]. The Gulf of Trieste is characterised by quite various sediments, mainly sandy in the coastal areas and sandy pelite or with shelf sand offshore. The seabed hosts the traditionally associated biocoenoses (sensu Pérès and Picard 1964) [53] that start from SFS (Sables Fins Superficiels) and SFBC (Sables Fins Bien Calibrès) along the strictly coastal bottoms, especially in the northernmost part, to the DC (Détritique Cotier), DE (Détritique Envasé) and VTC (Vases Terrigènes Cotières) biocoenoses offshore, the latter characterising the areas facing the estuaries of the Tagliamento and Isonzo rivers. At a distance from 3 to 10 nautical miles (nm) from the coast and a depth ranging between 13 and 25 m, off the Marano and Grado lagoon, from the sandy-detritic unvaried bottoms emerge some characteristic scattered rocky outcrops of biogenic concretions [54]. These concretion outcrops, constituted mainly by encrusting coralline algae, are known as “trezze” or “grebeni” in the Gulf of Trieste and as “tegnùe” off the coast of Venice and are considered the most peculiar features of the northern Adriatic Sea [55]. These coralligenous outcrops form complex three-dimensional structures that consist of a rocky substrate suitable for the settlement and the development of rich assemblages of associated flora and fauna, thus representing a unique biodiversity hotspot [54,56]. Therefore, in order to be preserved despite any climatic and anthropogenic pressures, some outcrops have been included in the European Natura 2000 network as Sites of Community Importance (SCIs), such as the site called “Trezze San Pietro e Bardelli” (IT3330009). The Gulf of Trieste is also home to the Miramare Marine Protected Area (MPA). This small MPA is the first to be established in Italy in 1979 and is divided into two distinct zones: the inner part (30 ha), which is subject to an integral protection regime, is surrounded by a larger buffer zone (90 ha). The protected area is also part of the European Natura 2000 network and represents the only fully protected area in the Italian part of the Gulf of Trieste [57].

Prior to the MME, the Gulf of Trieste hosted one of the densest healthy populations of Pinna nobilis in the Mediterranean Sea. Relatively high average densities were recorded in the Isonzo estuary and Panzano Bay with 14.38 ind/100 m² in association with Cymodocea nodosa meadows [58], in the Miramare MPA with 17.25 ind/100 m² [59] and even in the San Pietro outcrop with 6.77 ind/100 m² [60]. The epidemic has affected all of these populations, and 100% mortality rates have been documented in the Gulf between 2019 and 2021. Currently, wide areas of dead mussels can be observed throughout the area, still firmly anchored in the sediment. Most of these now provide a stable hard substrate for epibenthic species, including molluscs species that normally are found in hard bottoms.

2.2. Sample Collection and Processing

For the study of mollusc communities, a total of thirteen dead specimens of Pinna nobilis with intact and vertically anchored shells were collected in the sediment from four sites in the Gulf of Trieste (Figure 1). The two inshore sites (less than 100 m from the coast) Barcola (BAR) and Miramare (MIR), the latter within the buffer zone of the Miramare MPA, both with a maximum depth of 5 m, were located along the inner Gulf coastline, while the other two offshore sites Panzano Bay (PAN) and San Pietro (SP) were 1 km and 9 km from the coast, respectively (Figure 1). In addition, the deepest site (13–15 m depth), San Pietro (SP), was located in the Site of Community Importance (SIC) “Trezze San Pietro and Bardelli”. All samples were collected by SCUBA diving, using an extremely accurate sampling procedure consisting of three main steps: (1) carefully approaching the erect bivalve shell against the direction of light; (2) quickly wrapping the shell with a nylon bag with a mesh size of 500 µm and finally (3) removing it from the sediment with immediate sealing of the bag to prevent the escape of motile species. With this procedure we could ensure that the collected motile organisms belonged exclusively to the “ecosystem” of P. nobilis and did not originate from other adjacent habitats. Each specimen was quickly frozen in its nylon bag while waiting for processing. In the laboratory, images were taken of the inner and outer valvar surfaces of each shell to later determine the percent coverage of sessile species. Epibenthic fauna was then collected by scraping and washing each valve and the nylon bag content over a 500 µm sieve. The sieved material was then sorted, and molluscs were separated from the other major taxonomic groups (e.g., crustaceans, polychaetes, echinoderms and other) and preserved in a 80% ethanol solution.

Biometric measurements of the entire shell were taken, including the total length (Ht) consisting of maximum ante-posterior shell length, the maximum width (W) corresponding to the point of maximum dorso-ventral length of the shell that may coincide with the bend and the minimum width (w) coinciding with the interface shell sediment [61]. The outer surface area of the shell was determined by using ImageJ software (version 1.53T) [62] and tracing the profiles of both valves (left and right) of the completely cleaned and photographed shells after processing, to have the actual shell area. In addition, the proportion of the outer shell area covered by sclerobionts was also calculated.

2.3. Abundance and Biomass Estimation of Mollusc Community

Only molluscs alive at the time of sampling were counted and identified at the lowest possible taxonomic level, using the relevant available literature. The taxa are listed according to the updated nomenclature of the World Register of Marine Species [63]. To estimate biomass, the wet weight (WW) and shell-free wet weight (SFWW) of all individuals of each species present in each P. nobilis shell were measured. Each individual was blotted dry for a standard time of 30 s before weighing [64]. When measuring WW, the whole organism, including its shell, was considered, while the shell was removed to obtain SFWW. For small gastropods (<5 mm size) and chitons (Polyplacophora), only WW was measured as their soft tissues are firmly attached to the shell.

2.4. Structural Analysis of Mollusc Community

The diversity components of the investigated mollusc community were assessed using the following univariate indices: species richness (S), number of individuals (N), Margalef’s richness index (d), Pielou’s evenness index (J′) and Shannon–Wiener diversity index (H′(log₂)). A one-way Kruskal–Wallis ANOVA (non-parametric Analysis of Variance) was performed to assess differences in univariate diversity indices based on location (4-level factors: Miramare, Barcola, Panzano and S. Pietro) and the distance from the coast (2-level factors: inshore, offshore).

To determine the most common and the rarest epizoite molluscs of each P. nobilis community, Dajoz’s (1971) [65] constancy index was used, CA1 = nA1/n1 × 100, where nA1 was the number of shells where the species A occurred and n1 was the total number of shells sampled. According to the Dajoz index, molluscs were classified as rare (less than 12%), not common (13–24%), common (25–49%), very common (50–74%) and constant (75–100%).

To assess the differences in mollusc community structure, an agglomerative hierarchical cluster analysis (group–average linkage) was performed based on the Bray–Curtis similarity resemblance matrix from the square root transformed species abundance data. The square root transformation was performed to down-weight the importance of the highly abundant species. To better visualise the distribution patterns of the investigated faunal assemblages, Non-Metric Multi-Dimensional Scaling (nMDS) was performed based on the same Bray–Curtis similarity resemblance matrix. To test (not a priori) the consistent groups in the cluster analysis, the Similarity Profile Test (SIMPROF) was applied to both objects (Pinna community assemblage) and variables (index of association on species), while to test differences in community structure on (a priori) defined groups such as location (4-level factor: Miramare, Barcola, Panzano and S. Pietro) and distance from the coast (2-level factor: inshore, offshore), the semi-parametric one-way PERMANOVA (Permutational Multivariate Analysis of Variance) was applied with an unrestricted permutation of the raw data and the 9999 permutations test [66]. In addition, Spearman’s rank correlations (r_s) [67] were used to assess the relationship between diversity indices and shell surface area.

2.5. Functional Analysis of Mollusc Community

2.5.1. Species–Trait Analysis

Considering the interest in estimating changes in the benthic mollusc community at different sites with different ecological characteristics, a Biological Trait Analysis (BTA) was conducted. Traits related to mobility, ecological habits and biological characteristics were considered. Five traits, two continuous (maximum adult size and longevity) and three categorical (mobility, feeding and living habits), were selected to assess the role of each mollusc species within the community. Each trait comprised a number of categories, which are summarised in

Table 1. List of biological traits, type and relative categories.

Biological Trait	Trait Type	Category	Labels
Max adult size	Continuous	Very small (<5 mm)	<5 mm
		Small (5–30 mm)	5–30 mm
		Medium (30–80 mm)	30–80 mm
		Large (>80 mm)	>80 mm
Adult longevity	Continuous	Very short (≤1 year)	<1 year
		Short (1–3 years)	1–3 years
		Intermediate (3–6 years)	3–6 years
		Long (6–10 years)	6–10 years
Mobility	Categorical	Sessile	Ses
		Semi-motile	Se-mot
		Motile	Mot
Feeding Habits	Categorical	Suspension Feeder	SF
		Surface deposit feeder	SDF
		Sub surface deposit feeder	SSDF
		Grazer	G
		Scavenger	SC
		Predator	PR
		Parasite	PA
Living Habits	Categorical	Free-living	Fre-liv
		Burrow dwelling	Bur-dwe
		Crevices/hole/under stones	Cre-hol-uns
		Epi-endobiotic/Attached to a substratum	Epi-end-Att

. The categories were selected following [68]. Maximum adult size was defined as follows: under 5 mm (<5 mm), between 5 and 30 mm (5–30 mm), between 30 and 80 mm (30–80 mm) and over 80 mm (>80 mm). Adult longevity was categorised as less than 1 year (<1 year), between 1 and 3 years (1–3 years), between 3 and 6 years (3–6 years) and between 6 and 10 years (6–10 years). Mobility included the following categories: sessile (Ses), semi-motile (Se-mot) referring to infaunal species and motile (Mot). Feeding habits included the following trophic groups: suspension feeders (SF), surface deposit feeders (SDF), sub-surface deposit feeders (SSDF), grazers (G), scavengers (SC), predators (PR) and parasites (PA). Finally, living habits included free-living (Fre-liv), burrow-dwelling (Bur-dwe), those organisms that live in the crevices, holes and/or under stones (Cre-hol-uns) and epibiotic, endobiotic and/or attached to a substratum (Epi-end-Att) [68]. Each species was assigned to specific category using available information extracted from primary literature [68], books [69,70,71,72], online repositories [63,73] and expert ecological knowledge. The minimum and maximum affinity for each category was expressed with a binary code (0–1), and a species matrix of functional traits was created, and then weighted according to both abundance and biomass data.

2.5.2. Functional Diversity Analysis

Metrics of functional diversity such as Functional Richness (FRic), Functional Evenness (FEve), Functional Divergence (FDiv) and Functional Dispersion (FDis) were calculated using FDiversity Software [74] on combined species-scaled matrices by traits, weighted by abundance and biomass. Each index provides information on a single facet of functional diversity and is useful in its entirety for understanding the overall framework from a functional perspective [46,75] of the mollusc communities associated with P. nobilis shells. A one-way Kruskal–Wallis ANOVA was performed to assess differences in functional diversity metrics based on location (4-level factor: Miramare, Barcola, Panzano and S. Pietro) and the distance from the coast (2-level factor: inshore, offshore). In particular, FRic allows the study of the functional range occupied by the community; FEve defines how the species of a community are functionally distributed in the trait space; FDiv indicates how the species abundance is distributed within the volume of the functional trait space occupied by the species [46] and finally, FDis measures the difference between species in the community by estimating the distribution of species in a functional trait space [40,47].

Functional redundancy was estimated by testing linear and logarithmic regression models between the species richness (S) and the FRic index, as the extent of functional redundancy depends on the shape of the regression line [40,76]. In addition, functional identity was calculated as the community-weighted means (CWM) of trait category expression. The CWM values largely indicate the trait strategies given by the species pool and environmental conditions of a site [77]. The CWM values were weighted by species abundance (CWM-abu) and biomass (CWM-bio). The functional patterns of communities sampled at four sites were displayed performing two Canonical Analyses of the Principal coordinates (CAP) on the Euclidean distance matrix from CWM-abu and CWM-bio values [40,78]. As with the structural community data, a semi-parametric one-way (4 levels location and 2 levels Inshore-offshore factors) PERMANOVA was used for the functional data with an unrestricted permutation of the raw data and a 9999 permutations test.

3. Results

3.1. Community Structure and Composition

A total of 99 mollusc taxa belonging to 3 different classes (57 Gastropoda, 38 Bivalvia and 4 Polyplacophora) and 42 families were found in thirteen sampled shells of P. nobilis, collected from the four sites investigated in the Gulf of Trieste (Table S1). Of these, 96 were identified to species level (97% of the total). A total of 1172 individuals were sampled, with bivalves being the dominant group with 712 individuals, followed by gastropods with 425 individuals and polyplacophorans with 35 individuals.

The majority of the associated molluscs (68%) proved to be “rare” and “not common”, while the remaining 32% comprised “constant”, “very common” and “common” species. The constant species included three bivalves (Rocellaria dubia, Hiatella arctica, and Mimachlamys varia) and one gastropod (Hexaplex trunculus), which are typical of hard bottoms. The associated mollusc species of the families Mytilidae (Modiolula phaseolina, Musculus subpictus), Cerithiidae (Bittium latreilli, Bittium scabrum) and Pyramidellidae (genera Odostomia, Spiralina, Parthenina, Megastomia and Auristomia) proved to be very common and frequent species (Table S1).

All univariate indices (except evenness) were highest for P. nobilis shells collected from BAR (

Table 2. Univariate diversity indices: species richness (S), number of individuals (N), Margalef richness index (d), Pielou index of equitability (J′), Shannon–Wiener index (H′) and the biomass (g SFWW) of mollusc community found in P. nobilis shells with the mean and standard deviation (SD) values per site (Barcola, Miramare, Panzano, S. Pietro).

Site	Station	S	Mean	SD	N	Mean	SD	d	Mean	SD	J′	Mean	SD	H′	Mean	SD	Biomass	Mean	SD
Barcola (BAR)	B1	39	34.00	5.57	99	151.67	55.15	8.27	6.68	1.46	0.93	0.85	0.06	4.89	4.33	0.52	2.71	4.84	4.50
	B2	28			147			5.41			0.80			3.86			1.79
	B3	35			209			6.36			0.83			4.24			10.01
Miramare (MIR)	M1	18	23.67	6.03	103	101.00	34.04	3.67	4.95	1.16	0.67	0.78	0.10	2.79	3.54	0.65	10.14	12.31	4.92
	M2	23			66			5.25			0.87			3.95			8.85
	M3	30			134			5.92			0.79			3.87			17.95
Panzano (PAN)	P1	16	18.33	2.08	38	76.67	36.69	4.12	4.08	0.05	0.93	0.82	0.10	3.70	3.41	0.26	5.48	5.44	0.66
	P2	19			81			4.10			0.79			3.35			6.08
	P3	20			111			4.03			0.74			3.19			4.77
S. Pietro (SP)	S1	15	20.50	6.56	27	45.25	23.88	4.25	5.15	1.07	0.91	0.89	0.03	3.57	3.82	0.31	6.00	9.77	5.08
	S2	15			28			4.20			0.91			3.54			17.24
	S3	28			78			6.20			0.85			4.08			7.34
	S4	24			48			5.94			0.90			4.11			8.49

, Figure 2). The number of individuals ranged from a minimum of 45.25 ± 23.88 in shells from SP to a maximum of 151.77 ± 55.15 in BAR, while the number of species ranged from a minimum of 18.33 ± 2.08 in the shells from PAN to a maximum of 34.00 ± 5.57 in BAR. The associated Margalef richness index (d) also showed the same trend as the number of species, varying from a minimum of 4.08 ± 0.05 in PAN and a maximum of 6.68 ± 1.46 in BAR. The evenness Pielou index (J′) shows relatively high values at all sites and was lowest at MIR with 0.78 ± 0.10 and highest at SP with 0.89 ± 0.03. Based on the richness and evenness observed in the sampled shells, the Shannon–Wiener diversity index (H′) shows values clearly related to the richness component of diversity, with the highest values recorded at BAR with 4.33 ± 0.52 and the lowest at PAN with 3.41 ± 0.26. The associated mollusc biomass expressed as SFWW ranged from a minimum of 4.84 ± 4.54 g in shells from BAR to a maximum of 12.31 ± 4.92 g in specimens from MIR. The shell surface areas and biometric measurements were similar in the specimens collected from the four sites, as shown in

Table 3. Biometric measures: total length (Ht), maximum width (W), minimum width (w) and coverage by sclerobionts (%) of P. nobilis shells with the mean and standard deviation (SD) values per site (Barcola, Miramare, Panzano and S. Pietro).

Site	Station	Suface Shell Area (cm2)	Mean	SD	Ht (cm)	Mean	SD	W (cm)	Mean	SD	w (cm)	Mean	SD	% Coverage	Mean	SD
Barcola (BAR)	B1	952.48	962.49	51.02	48.60	50.93	2.30	16.90	17.30	0.61	11.00	10.63	0.55	96.50	95.53	3.92
	B2	917.23			53.20			17.00			10.90			98.87
	B3	1017.78			51.00			18.00			10.00			91.21
Miramare (MIR)	M1	1216.53	1120.56	185.79	55.50	54.30	1.66	17.50	17.53	1.85	13.50	12.23	1.70	99.13	98.95	0.52
	M2	906.41			52.40			15.70			10.30			98.37
	M3	1238.74			55.00			19.40			12.90			99.35
Panzano (PAN)	P1	980.94	932.76	172.95	55.20	54.67	2.44	15.00	15.00	1.00	9.80	9.97	0.96	74.58	79.92	9.21
	P2	1076.50			56.80			16.00			11.00			74.62
	P3	740.83			52.00			14.00			9.10			90.55
S. Pietro (SP)	S1	1394.84	1273.03	180.42	59.00	59.38	3.35	20.50	18.78	1.37	12.00	10.68	0.91	99.85	96.31	5.11
	S2	1004.63			55.00			17.40			10.00			96.54
	S3	1342.02			60.50			18.00			10.20			89.01
	S4	1350.62			63.00			19.20			10.50			99.84

. The surface area occupied by sclerobionts was very large at BAR, MIR and SP with more than 95%, while a proportion of 80% was recorded for the shells from PAN (Table 3). According to the results of the Spearman rank correlation test, no correlations were found between all diversity indices and shell surfaces, nor between biomass and shell surfaces.

The multivariate analysis performed on abundance with the Bray–Curtis dissimilarity matrix showed that the mollusc communities were clearly separated according to the sampling sites. The n-MDS plot shows the sampled P. nobilis shells grouped by site, which is highlighted by the SIMPROF analysis (p < 0.05) in which groups are superimposed on the ordination plot (Figure 3). The PERMANOVA analysis performed on the mollusc community also confirms the differences between sites (Pseudo-F = 4.45; p < 0.01).

To have an indication of which molluscs species are most involved in the observed differences between sites, we performed an agglomerative k-R clustering shade plot (Figure 4). In the shade plot the four significantly different groups of sites (I–IV) resulting from the SIMPROF test are represented by vertical lines, and the sites’ grouping can be seen to be driven by seven groups (A–G) of species (index of association among standardised square root transformed abundance data and SIMPROF test) (Figure 4). The shade plot clearly identifies the key groups of species that typify the four clusters (sites) of P. nobilis shells in the mollusc community (Figure 4). The P. nobilis shells collected at Barcola (BAR1, BAR2 and BAR3) were characterised by groups B, C, and D, where gastropods predominate over other molluscs. The samples from Miramare (MIR1, MIR2 and MIR3) were mainly characterised by group F, in which bivalves dominated and the species Rocellaria dubia was present in greater abundance than others. In the pen shell samples from Panzano (PAN1, PAN2 and PAN3) and Barcola, there was slight overlap in community composition, as group E was mainly represented by Mytilidae and Modiolula phaseolina was evenly distributed at both sites. In addition, shells from Panzano were also characterised by low abundances by group G, which included a total of five species, of which three gastropods (Odostomia turrita, Cerithidium submammillatum and Pseudofusus rostratus) and two bivalves (Chama gryphoides and Dosinia lupinus). Finally, S. Pietro (SP1, SP2, SP3 and SP4) was mainly characterised by group A, which included gastropods found only at this most offshore and deepest site (e.g., Bolma rugosa, Caecum trachea, and Calliostoma zizyphinum).

3.2. Functional Diversity Indices

The multi-trait indices of functional diversity calculated on a Gower distance matrix of species traits weighted by abundance and biomass are summarised in

Table 4. Comparison of functional diversity indices weighted by abundance (A) and biomass (B) of mollusc species found in P. nobilis shells with mean and standard deviation values per site. Functional Richness (FRic), Functional Divergence (FDiv), Functional Evenness (FEve) and Functional Dispersion (FDis).

Site	Station	FRic (A,B)	Mean	SD	FDiv (A)	Mean	SD	FDiv (B)	Mean	SD	FEve (A)	Mean	SD	FEve (B)	Mean	SD	FDis (A)	Mean	SD	FDis (B)	Mean	SD
Barcola (BAR)	B1	8.14	6.50	1.52	0.92	0.89	0.03	0.83	0.92	0.08	0.62	0.53	0.08	0.27	0.31	0.05	4.71	4.47	0.22	3.79	3.06	1.65
	B2	5.15			0.88			0.93			0.50			0.37			4.43			4.21
	B3	6.20			0.86			0.99			0.48			0.29			4.28			1.17
Miramare (MIR)	M1	5.67	4.20	1.43	0.94	0.93	0.02	0.92	0.93	0.04	0.68	0.61	0.07	0.44	0.30	0.12	4.45	4.54	0.14	3.77	3.53	0.27
	M2	2.81			0.94			0.97			0.60			0.25			4.70			3.24
	M3	4.12			0.90			0.89			0.55			0.22			4.46			3.57
Panzano (PAN)	P1	5.12	4.61	0.66	0.85	0.86	0.03	0.76	0.86	0.10	0.81	0.66	0.16	0.31	0.36	0.06	4.21	4.28	0.32	0.91	2.63	1.49
	P2	3.86			0.89			0.86			0.68			0.43			4.62			3.43
	P3	4.84			0.83			0.96			0.50			0.35			4.00			3.56
S. Pietro (SP)	S1	4.51	4.46	1.01	0.92	0.90	0.02	0.84	0.87	0.08	0.76	0.70	0.07	0.49	0.37	0.12	4.88	4.68	0.26	3.53	3.23	0.52
	S2	3.56			0.90			0.81			0.72			0.21			4.89			3.73
	S3	3.92			0.89			0.83			0.60			0.36			4.34			2.55
	S4	5.85			0.87			0.98			0.70			0.42			4.61			3.11

and shown in Figure 5. The Functional Richness index (FRic) values, which for intrinsic characteristics show the same values when both mollusc abundance and biomass are considered, ranged from 2.81 to 8.14, with the highest values observed in the samples from Barcola (mean value 6.50 ± 1.52). Panzano showed lower values than Barcola (mean value 4.61 ± 0.66) but similar values to the Miramare and S. Pietro samples (4.20 ± 1.43 and 4.46 ± 1.01, respectively) (Table 4, Figure 5A). Functional Divergence and Functional Evenness (FDiv and FEve) are usually expressed as values between 0 and 1. Comparing FDiv values between abundance and biomass, they were close to 1 in both cases at all sites studied, with the same values in Miramare (0.93 ± 0.02), where the highest values were observed (Figure 5B,E; Table 4). The samples from Panzano had the same lowest value for both abundance and biomass (0.86 ± 0.03), while an opposite trend was observed in the samples from Barcola and S. Pietro. In the former sample, there was a slight increase in abundance data compared to biomass data, while the latter slightly decreased (Figure 5B,E). Clear differences were found in the FEve values between abundance and biomass, but overall no differences were found among the sites (Table 4, Figure 5C,F). In particular, FEve values were clearly above 0.5 for abundance data and below 0.5 for biomass data. Furthermore, different trends were observed when comparing sites by abundance and biomass: Barcola samples had the lowest values weighted by mollusc abundance (0.53 ± 0.08), while Miramare had the lowest values weighted by mollusc biomass (0.30 ± 0.12). The highest values for FEve were observed in S. Pietro (Table 4, Figure 5C,F). Finally, Functional Dispersion (FDis) showed overall higher values considering the abundance of mollusc species than mollusc biomass (Figure 5D,G). The mollusc species found in the P. nobilis shells from S. Pietro had the highest FDis weighted by abundance (4.68 ± 0.26), while the biomass of mollusc species from Miramare had the highest values (3.53 ± 0.27). Both mollusc abundances and biomasses had the lowest FDis values in Panzano (4.28 ± 0.32 and 2.63 ± 1.49, respectively) (Table 4, Figure 5D,G).

Regarding functional redundancy, the results showed a non-significant correlation between species richness and functional richness according to both the linear and logarithmic models (Table S4, Figure S1).

3.3. Functional Analysis of Mollusc Community

An ordination plot using CAP on the Euclidean distance matrix from normalised CWM data weighted on mollusc abundance and biomass is shown in Figure 6A,B. A clear variation in the functional pattern of epifaunal mollusc abundance (CWM-abu) was observed among the four sampled sites. The functional expression of the mollusc community of SP was mainly separated from the shells of PAN by CAP1, while the samples of BAR, but also of PAN and SP, were mainly separated from the shells of MIR by CAP2 (Figure 6A). PERMANOVA tests revealed a significant difference between sites (Pseudo-F: 4.91, p < 0.01) but not between inshore and offshore sites (MIR and BAR vs. PAN and SP) (Pseudo-F = 1.08; p < 0.36). With a Pearson correlation coefficient of >0.4 with the axes of CAP, twenty vectors of the trait categories were identified and plotted on the CAP (Figure 6A). Of these, nine categories (e.g., Ses, Epi-end-Att, Cre-hol-uns, SF, SDF, SSDF, 3–6 years, 6–10 years and >80 mm) generally pointed towards samples from Miramare (MIR), six categories (e.g., 30–80 mm, G, SC, PR, Se-mot, and Bur-dwe) were pointed towards samples from the offshore site S. Pietro (SP),while two categories (e.g., Fre-liv and Mot) and the remaining three (e.g., PA, <5 mm and 1–3 years) were pointed towards Barcola (BAR) and Panzano (PAN), respectively. Regarding the mollusc biomass (CWM-biom) on the shells of P. nobilis, CAP revealed a different functional pattern between the sites (Figure 6B). Samples from Barcola (BAR) were separated from those from Panzano (PAN) and S. Pietro (SP) mainly by CAP1 and from those from Miramare (MIR) mainly by CAP2; in addition, samples from S. Pietro (SP) and Panzano (PAN) were grouped together. PERMANOVA tests revealed significant differences in the functional structure (biomass weighted) of mollusc community between sites (Pseudo-F = 2.82; p < 0.05) and also between the distance from shore factor (Pseudo-F = 6.48; p < 0.01). With a correlation coefficient of >0.4 with the CAP axes, eighteen trait categories were identified (Figure 6B). Of these, seven categories (e.g., Bur-dwe, Cre-hol-uns, G, SDF, SSDF, <5 mm, 5–30 mm) pointed to samples from Barcola (BAR), while six categories (e.g., Epi-end-Att, Ses, SF, 3–6 years, 6–10 years) and the remaining five (e.g., Fre-liv, Mot, SC, PR and 30–80 mm) pointed to samples from S. Pietro (SP) and Panzano (PAN) (Figure 6B).

4. Discussion

An epidemic that spread throughout the Mediterranean Sea triggered a Mass Mortality Event (MME) of Pinna nobilis in several areas, including the northern Adriatic Sea, affecting one of the densest populations of this fan mussel. The aim of this study was to assess the remaining role of dead specimens of P. nobilis as ecosystem engineers by describing the associated mollusc sclerobiontic fauna through structural and functional analyses. As far as we know, this is the first attempt to evaluate the remaining ecological role of P. nobilis shells, still anchored in the sediment, after the occurrence of the MME in the Mediterranean Sea. This study was conducted on shells collected from four different sites in the Gulf of Trieste. The results show that a large number of mollusc species inhabit the shells, regardless of the sampling location. The overall diversity, expressed as univariate diversity metrics (e.g., Shannon diversity, species richness and evenness) was even higher compared to the findings of the studies conducted by Giacobbe (2002), Cosentino and Giacobbe (2007; 2008) [12,79,80] on mollusc assemblages in live specimens of P. nobilis in the Strait of Messina and Rabaoui et al. [13] on the macrofaunal community associated with P. nobilis in live specimens collected at three sites along the Tunisian coast.

In contrast to other observations, no significant correlation was found between diversity indices and shell surface area in our study. This result seems to contradict the ecological theory that states the larger the available surface area, more species it can host [81,82] and also the results of Giacobbe (2002) [79] who found a positive correlation between the shell size of P. nobilis and the Shannon-Wiener diversity index of molluscs. The lack of correlation between our shell size and the observed diversity metrics is probably related to the small differences in the dimensions of the P. nobilis shells sampled in this study. Indeed, when selecting the fan shells, we tried to sample similar sizes to avoid attributing differences in mollusc community structure between sites to differences in shell size. In our P. nobilis samples, the small differences in size do not seem to lead to remarkable differences in species number and diversity. The considerable number of mollusc species found in our samples was quite comparable, if not higher, than in the few previous studies carried out on live specimens of P. nobilis in the Mediterranean Sea [12,13,79,80]. This suggests that even the dead specimens that remain fixed in the sediment continue to serve as available substrate for a variety of mollusc species, just like the living P. nobilis. To assess how long this role may last requires further investigation and analysis which are beyond the scope of this paper. However, it can be speculated that the hydrodynamic conditions of the surrounding environments in which the P. nobilis samples were collected may have positively influenced colonisation by preventing fan mussel shells from dislodging and scattering onto the surrounding soft bottoms, and instead favouring their structural integrity. As has been observed in other dead marine bivalves [83,84] and other macro- and megabenthic species [85,86,87,88,89,90,91,92,93,94,95], the shells of P. nobilis should represent an available space for a typical hard bottom epibenthic community [13,96,97,98,99]. However, the particular nature of the vertical shell surface and the surrounding soft sediment habitat may influence the associated species assemblage and thus represent potential selection factors for the epifauna [13,80], which also seems to be confirmed by our results.

Indeed, when looking at the overall taxonomic composition, clear differences were found between sites and distance from the shore as clearly shown by the n-MDS ordination (see Figure 3) and confirmed by both SIMPROF and PERMANOVA analyses. Species composition varied between the sites studied and a monotonic difference between sites emerged in some species groups. The community of S. Pietro, the most offshore site, was characterised by different species compared to those found in the other sites, such as the gastropod molluscs Bolma rugosa and Calliostoma zizyphinum, typical of hard bottoms and common in the adjacent coralligenous outcrops habitats [60]. The other sites, in turn, showed further differences among themselves, with a community in Miramare characterised mainly by suspension-feeding bivalves such as Ostrea edulis and Mytilus galloprovincialis, the community in Panzano by small mussels such as Modiolula phaseolina and Musculus subpictus and the community in Barcola by a variety of species, mainly gastropods. Some of the species were found to be common epibionts of Pinna nobilis, such as the sessile bivalves species Anomia ephippium, Arca noae, Mimachlamys varia and Hiatella rugosa, as well as some gastropods of the genus Bittium [12,79,80]. Conversely, some species such as the endolithic boring bivalves Rocellaria dubia and the endolithic Hiatella arctica, both present in high densities mainly in Barcola, never occurred in the mollusc assemblages found in live specimens of P. nobilis in the Strait of Messina. In particular, R. dubia was found in a calcareous abiotic structure nearby the living northern Adriatic populations of the fan mussel [100]. The presence of these sessile, endolithic and hard bottom species at the inshore sites of Miramare and Barcola could be attributed to their proximity to artificial structures and reef-like breakwaters, which are known to be easily colonised by a large number of epibionts [101] and may facilitate a mechanism for larval export and active spillover towards the surrounding environments [102]. Furthermore, the high frequency found in all samples regardless of the location may indicate that the now inert calcareous shell may represent an excellent substrate even compared to the calcareous shell of living organisms.

The influence of the particular characteristics of the surrounding bottoms and the overall potentially different conditions at the individual sites on the structuring of the various mollusc communities observed at the sampled sites is also confirmed by the large number of rare or not common species, which account for exactly 69% of the total (see Table S1). This large number of rare or not common species underlines a clear turnover component of biodiversity between sites and could reflect species sorting caused by environmental differences (location factor) or dispersal limitation leading to the selective differentiation of species pools between communities (inshore-offshore factor) [103]. A notable difference between the dead shells of P. nobilis compared to those sampled alive before the MME is the considerable amount of infaunal bivalve species (always in small numbers) found in our samples, many of them even typical of muddy bottoms. Indeed, some of characteristics infaunal bivalves such as Abra alba, Chamelea gallina, Dosinia lupinus, Gastrana fragilis, Gouldia minima, Kurtiella bidentata, Loripes orbiculatus, Lucinella divaricata, Nucula nucleus, Pitar rudis, Timoclea ovata, Varicorbula gibba, Venus casina and Venus verrucosa were found at all sites. This apparent contradiction with the characteristics and the nature of hard bottom of the shell of P. nobilis is probably related to the amount of sediment filling the buried part of the shell, which has now been freed from the living organism. In fact, each P. nobilis shell was found to be filled with very fine sediments in its part anchored in the soft bottom. These sediments inside the shells, which are clearly different from the detritic sand (SP) and the pelitic sand (BAR, MIR and PAN) of the surrounding bottoms, probably act as a fine soft bottom enclave, leading to a peculiar biotope associated with P. nobilis dead shells. This could also be related to the observed significant presence of smaller gastropods below 5 mm in size, considered micromolluscs, and belonging to the Pyramidellidae family, which are normally an important infaunal component of soft bottoms [104]. These observations are also confirmed by the results of Gutierrez and Iribarne (1999) [105], who found that the dead but life-positioned shells of the stout razor clam Tagelus plebeius not only provided hard substrate for epifaunal species in muddy bottoms, but also were filled with sediment and hosted a relevant number of infaunal species. Similar to the fan mussels’ empty shells, this contributed to a positive enhancement of local biodiversity in muddy bottoms.

Looking at the trophic guilds, sessile and suspension feeders’ species were most abundant in all P. nobilis shells examined in our study, especially at the inshore sites (BAR and MIR). This result differs markedly from the findings of Giacobbe (2002) [79], who observed a clear dominance of omnivores (including non-selective grazers) in the molluscs’ community associated with P. nobilis in the Strait of Messina. In the P. nobilis shells from Sicily, the omnivore molluscs reach more than 60% on average, and suspension feeders were also very low compared to herbivore-grazers and parasites [12,80]. It is well known that space, i.e., available substrata, is therefore considered a limiting resource in sessile suspension feeders’ communities [106], and these organisms are considered to be limited primarily by space, while competition for other resources is considered to be minor or non-existent [107,108,109]. The remarkable differences observed in the number of sessile suspension feeders between the molluscs communities of P. nobilis in the Strait of Messina and our samples collected in the northern Adriatic Sea seem to indicate that competition for space may not be the main limiting resource in these sessile communities, but probably if there is not an adequate food supply for suspension feeders the competition for space should be ineffective, at least for the sessile suspension feeders. Indeed, the Strait of Messina, like the Tyrrhenian and Ionian Seas it divides, has average values of primary productivity that are much lower than in the northern Adriatic [110], and the availability of food for suspension feeders is limited. Since resources are considered limiting when an insufficient supply limits the growth, abundance or distribution of an organism or a population of organisms in an ecosystem [111,112], in an oligotrophic environment such as the Straits of Messina, food availability could be the limiting factor for sessile suspension feeders, while in mesotrophic seas such as the northern Adriatic, available substrata could be one of the most important limiting factors. This seems to be strongly supported by the high percentage of shell surfaces occupied by sclerobionts (ranging from about 80% to over 95%) (see Table 3), mainly sessile suspension feeders, observed in all P. nobilis specimens examined (see Table S2).

Functional diversity indices (e.g., FRic, FDiv, FEve and FDis) showed no significant differences across location and distance from the shore. Functional richness (FRic) was slightly higher at Barcola than at other sites because its collinearity with species richness [113] depends on the addition or removal of species with unique combinations of trait categories [46]. The high values of functional divergence (FDiv) and functional dispersion (FDis) observed at all sites indicate the occurrence of assemblages with a high relative abundance of species with unique categories of functional traits [114]. Different trends were observed when comparing abundance-based FEve and biomass-based FEve. Medium-high values of abundance-based FEve were observed at all sites, indicating good uniformity in the expression of community traits, whereas a general decline was observed in biomass-based FEve. While traits expressed by abundance may be numerically comparable (e.g., SF vs. SDF, SSDF, G, etc.), biomass distribution and expression could be biased only for certain dominant traits (e.g., SF vs. all others) (see Tables S2 and S3). Considering the discrepancy between the values obtained, a combination of the functional abundance analysis with the functional biomass analysis would be extremely useful. Overall, all functional metrics indicate a high degree of functional diversity, which is observed in a similar way in all molluscan communities, regardless of the location factor and the distance-from-shore factor (inshore-offshore). Regarding functional redundancy, the correlation between functional richness and species richness was not significant in the linear and logarithmic models tested. Since a linear correlation suggests that as species richness increases, functional richness also increases, implying a low level of functional redundancy as each new species introduces new features [76], this pattern is not confirmed in our results, indicating some functional redundancy in the molluscan community inhabiting the dead fan shells. Functional redundancy refers to a situation where the same role or function is shared in an ecosystem, which increases stability and enables positive responses to specific stressors [76].

The general functional patterns described by the abundance-based analysis are partly different from the biomass-based patterns in our study. The mollusc communities at the four sites described by the functional aspect were clearly separated according to the expression of specific traits when the analysis was based on abundance-based values, but the distribution of many traits changed between the sites when the analysis was performed on biomass-based values, showing that there are no differences in the expression of traits at the offshore sites of PAN and SP. These differences between abundance-based and biomass-based functional analysis are probably related to the numerical dominance of small-sized species (e.g., Pyramidellidae, Musculus subpictus and Modiolula phaseolina), which do not represent an expressive part of the biomass as observed in other studies [115]. However, other authors consider that abundance-based analysis in functional assessment may be more suitable for revealing functional differences, especially when used to highlight differences in polluted situations due to the importance of opportunistic species, which are usually small-sized and can be very abundant in stressed/polluted environments [113]. In this context, biomass-based analysis should also be very useful in describing environmental gradients. Indeed, ecosystem functioning could be strongly influenced by changes in size structure and biomass rather than by changes in abundance patterns of taxonomic communities [116]. Changes in the size structure of multispecies assemblages can have significant effects on interactions between organisms and the energy flow in the environment. On the other hand, it has been observed in other benthic studies that significant differences in benthic size spectra can occur in response to various environmental changes. For example, lower biomass values due to the reduced presence of large organisms along bathymetric gradients across continental margins have been attributed to decreasing food availability and food predictability. Similarly, an increasing proportion of smaller size classes in benthic size spectra has been found with a marked reduction in biomass due to a high sedimentation rate and sediment instability in gradients related to rivers’ influence [117]. Body mass is a fundamental organism trait that affects metabolic rate, energy demand and uptake rate [118] and is an important feature of overall population and community structure through density-mass relationships [119], especially when considering trophic guilds. For these reasons, we believe that biomass-based analysis should also be considered, even though it is often neglected, when arguing about different environmental gradients at least at the same level as abundance-based analysis.

From the results of this study, it can be concluded that the empty shells of Pinna nobilis retain an important part of the ecological role of the living shells, despite the death of the animal. A first investigation conducted in this study confirmed that a rich community of molluscs lives on the shells of the fan mussel, which not only continue to provide additional hard substrate for colonisation but also accumulate a certain amount of very fine sediment inside, which can become a receptacle for infaunal species. The ability to still supply hard substrate and provide a sheltered and protected muddy microhabitat inside leads to a significant aggregation of different species that fulfil different ecological roles, as revealed by the functional trait analysis. In addition, a considerable number of mollusc species tend to share common traits, resulting in some functional redundancy within the community. Moreover, this lends a certain stability to the rich “micro-ecosystem” formed by the empty shells of Pinna nobilis.

5. Conclusions

In this study, the mollusc community on dead Pinna nobilis specimens was examined to assess their remaining ecological role as ecosystem engineers and biodiversity enhancers. Ecosystem engineers are organisms that modify and create new habitats, providing substrate and protection for many other species. The results of this study suggest that the ecological role of this fan mussel persists after the animal’s death. Community structure analysis revealed that the empty shells of the fan mussel can harbour and aggregate a variety of molluscs, providing valuable hard substrate while also serving as a soft-bottom enclave. Furthermore, functional analysis revealed a high degree of functional diversity in the mollusc community inhabiting the P. nobilis shells collected in the Gulf of Trieste. Although BTA has proven to be a valid tool for understanding ecological dynamics, the discrepancies between abundance-based results and biomass-based results suggest that a combination of the two approaches, and especially knowledge of where most of the biomass is allocated in the trophic guilds, could be useful to provide a more comprehensive framework for trophic and ecological patterns within a community.