Located on the north side of the Gargano National Park (Foggia, Italy), Lake Varano is a coastal lagoon comprising an area of ca. 65 km² with an average water depth of ca. 4 m (Spagnoli et al. 2002). The lake is connected to the Adriatic Sea via two artificial channels (Capoiale and Varano) that are located at the two ends of a coastal dune (Fig. 1). From a hydrological point of view, the lake is influenced both by exchanges with the sea through the two channels and by freshwater inputs from the catchment area through two tributaries (the Antonino and San Francesco Canals). The surface sediments are very heterogeneous, with the sandy fraction prevailing near the baymouth bar (Frontalini et al. 2013). The lagoon has been intensively used for the farming of mussels (Mytilus galloprovincialis) and other aquaculture products (Beneduce et al. 2010). Although the surrounding area is characterized by only a few small towns and no major industrial activity, the relatively high concentration of nitrogen species probably reflects the impact of agricultural activities (Specchiulli et al. 2008). Though high nutrient loads and dinoflagellate blooms might suggest eutrophic conditions, Lake Varano has been regarded as an oligo-mesotrophic ecosystem (Specchiulli et al. 2008, Belmonte et al. 2011, Frontalini et al. 2013).

Samples from 21 stations were collected using a modified model of the Van Veen grab in late March 2012. This modified model permits the insertion of a plexiglas corer from the central part of the grab, enabling portions of almost completely undisturbed sediments to be captured. A set of sub-samples of the sediments was taken from each station for the purposes of sedimentological, geochemical and nematode analyses (two and three replicates for geochemical and nematode analyses, respectively). All of the analyses were carried out on the upper surface sediment layer (0-2 cm). During the sampling process, physico-chemical parameters of the bottom water, such as temperature, pH, salinity, oxidation reduction potential (ORP) and dissolved oxygen (DO, expressed as mg/l) were measured using a conductivity, temperature and depth (CTD) probe. In particular, these environmental parameters selected as independent variables are well known to affect the structure, abundance and diversity of the nematode assemblage (Heip et al. 1985, Semprucci and Balsamo 2012).

Grain size analyses were conducted at Lille 1 University (France) following the principle of diffraction described in Trentesaux et al. (2001). Three main fractions were considered: clay (<2 µm), silt (2 to 63 µm) and sand (63 to 2000 µm).

The quantity and quality of the OM were estimated by analysing the total carbon (TC), the total organic carbon (TOC), the total sulfur (TS), the total hydrogen (TH) and the total nitrogen (TN) content in the sediment (for details, see Armynot du Châtelet et al. 2013).

Concentrations of 37 trace elements were obtained at Activation Laboratories Ltd. (Ontario, Canada, http://www.actlabs.com), where a fraction of ~0.5g was digested with aqua regia and then analysed using Inductively Coupled Plasma Optical Emission Spectrometry. Only concentrations of As, Cd, Co, Cr, Cu, Ni, Pb, Zn and Mn were considered. The enrichment of some chemical elements whose higher concentrations may be toxic was calculated using the Pollution Load Index (PLI) following Martins et al. (2013).

Samples utilized for the nematode analysis were treated with a 7% MgCl₂ aqueous solution for narcotizing fauna, fixed in a 4% formaldehyde solution in buffered sea-water and stained with Rose Bengal (0.2 g L^–1). In laboratory, the samples were rinsed with a gentle jet of fresh water through a 0.5-mm sieve to separate the macrofauna (Danovaro et al. 2004). They were then decanted, sieved 10 times through a 42-μm mesh and centrifuged three times with Ludox HS30 (specific density 1.18) (Pfannkuche and Thiel 1988). According to Moreno et al. (2011), 100 nematodes were randomly withdrawn from each of the three replicates and subsequently mounted on permanent slides according to Seinhorst (1959). The specimens were identified under a 100× oil immersion objective using Nomarski differential interference contrast illumination (Optiphot-2 Nikon) and the NeMys online identification key (Deprez et al. 2005).

The trophic and life strategies were among the functional traits considered. Following Wieser (1953), the nematodes were divided into four trophic groups: selective (1A) and non-selective (1B) deposit feeders, epistrate feeders (2A), and predators/omnivores (2B).

Bongers (1990) and Bongers et al. (1991) proposed classifying nematodes as ‘colonizers’ (r-strategists) or ‘persisters’ (k-strategists). Therefore, the maturity index (MI) was calculated as the weighted average of the individual colonizer-persistent (c-p) values.

The biodiversity pattern of the assemblage was described by means of the Shannon (H’, using log-base 2), Pielou, Margalef and Simpson indices.

All of the statistical analyses were performed on the relative abundances of the nematode species. Following Armynot du Châtelet et al. (2004), only taxa with a relative abundance exceeding 5% in at least one sample were taken into consideration. However, though rare species often considered the most sensitive ones were down-weighted in this way, their contribution in the TE system was measured by means of MI calculation. Prior to the statistical analyses, all of the biotic and abiotic data were normalized by applying an additive logarithmic transformation log(x+1). In order to define the different assemblages, which included a group of species with a similar spatial distribution pattern, an R-mode cluster analysis (CA) was performed. Alternatively, a Q-mode CA was carried out for the ordination of the samples based on the relative abundances of species. These tests were performed by applying Ward’s linkage method and the findings were given in terms of the Euclidean distance. A principal component analysis (PCA) was used to determine the assemblage relationships to the abiotic parameters and for the ordination of sample locations.

The data on the environmental parameters are summarized in Supplementary material Table S1. The pH was slightly alkaline, ranging from 7.6 (V1, V15) to 8.0 (V29, V39). The salinity varied from 23.0‰ (V50) to 31.9‰ (V13), with a clear increase in the outer part of the lake and, in particular next to the Capoiale Channel. The DO values varied between 5.7 (V15) and 12.8 (V35), with the highest values found in the central part of the lake and the lowest near to the two channels. The ORP values ranged from 164.0 (V39 and V41) to 286.0 (V1), with the highest values in the outer part of the lagoon. The sediments were primarily composed of silt (47% on average) and sand (51.9% on average), with a low amount of clay (1.1% on average). The outer part of the lake, particularly in front of the two channels and the baymouth bar, was dominated by sand, whereas silt was predominant in the southeastern part of the lake. The CaCO₃ varied between 9.5% (V4) and 80.8% (V37), and was particularly abundant on the western side of the lake. The TN ranged between 0.1 (V4, V7, V11, V41, V48, V49 and V50) and 0.4 (V27 and V47), whereas TS was between 0.9% (V7) and 1.8% (V1). The TC values varied from 1.3 (V4) to 10.9 (V37 and V42), while TH ranged from 0.1 (V4) to 1.3 (V39 and V47). The TOC varied from 0.2% (V4) to 4.2% (V27), with higher values being found in the central and southwestern parts of the lake. A clear spatial distribution of most of the trace elements was visible, with the highest concentrations being found in both the central and southeastern parts of the lake. The PLI values ranged from 1.0 (V4) to 11.4 (V41), and revealed a general enrichment of the concentrations of trace elements in both the central and southeastern parts of the lake.

Fifty-five nematode species, belonging to 37 genera and 17 families, were recorded at Lake Varano (Supplementary material Table S2). The richest families were Chromadoridae (10 species), Linhomoeidae (9) and Xyalidae (8). Almost half of the selected specimens belonged to either the Desmodoridae (35%) or Linhomoeidae (16%) families, followed by Leptolaimidae (12%), Microlaimidae (9%) and Xyalidae (9%). Three species had very high abundances and a wide distribution: Molgolaimus allgeni (35%), Terschellingia longicaudata (13%) and Leptolaimus luridus (11%). They were followed by Daptonema normandicum (6%), Sabatieria pulchra (5%), Aponema torosa (4%) and Calomicrolaimus honestus (4%). Aponema torosa was more abundant in the outer (particularly in front of the two channels) and central parts of the lake, whereas D. normandicum and S. pulchra were more abundant in the stations close to the San Francesco output and in the central area.

The richness of the species ranged from 4 (V29) to 27 (V7). The H’ ranged from 0.9 (V29) to 4.2 (V7), the Margalef index from 0.7 (V29) to 5.6 (V7), the Pielou evenness from 0.4 (V29) to 0.9 (V7), and the Simpson evenness from 0.33 (V29) to 0.94 (V7 and V50) (Supplementary material Table S3). These indices revealed that the nematode assemblage in the central part of the lake was less diversified.

The trophic structure of the assemblage was mainly characterized by 1A and 1B (60.1% and 21.8%, respectively, on average), followed by 2A (16.2%) and 2B (2%) (Supplementary material Table S3). The group 1A was widely distributed, whereas the group 1B was particularly associated with freshwater inputs and the Varano Channel. The group 2B and, in particular, 2A were more abundant in the outer part of the lake.

The MI ranged from 2.1 (V4) to 2.9 (V27), with higher values mainly in the central part of the lake. All of the classes of the colonizers-persisters (c-p) were revealed, with the only exception being class c-p5 (Supplementary material Table S3). On average, the c-p3 species were the dominant component of the nematode assemblage (58.5%) followed by c-p2 (38.4%), c-p1 (2.8 %) and c-p4 (0.4%).

The results of the R-mode CA enabled two main (A and B) clusters and two different subclusters (A1 and A2) to be recognized (Fig. 2). The most abundant taxa in subcluster A1 were A. torosa, D. normandicum and S. pulchra, whereas M. allgeni, L. luridus and T. longicaudata were the only three representatives of the subcluster A2. Cluster B grouped together several taxa with very low relative abundances, including Daptonema fistulatus, O. wieseri and P. longicaudata as the main representatives. The Q-mode CA separated two different clusters (A and B) and two subclusters (B1 and B2) (Fig. 3). Cluster A was represented by station 4, which was in front of the baymouth bar, where the main taxa characterizing the nematode assemblages were O. wieseri, P. longicaudatus, P. pellucida and A. caudostriatus. This cluster was characterized by the lowest richness and MI, and was dominated by c-p2 taxa and 1B feeders. This subcluster was dominated by sand and characterized by the lowest values of CaCO₃, TOC, TN, TC, TH and TS. The subcluster B1 included the stations of the southeastern part of the lake: it was mainly dominated by M. allgeni, L. luridus, T. longicaudata, D. normandicum and S. pulchra. The same subcluster was characterized by higher levels of richness and the highest H’ value. It was dominated by c-p2 and c-p3 taxa and by 1A and 1B feeders. The subcluster B1 also showed the lowest salinity and ORP values and was mainly characterized by silty and, secondarily, sandy substrates. It also had higher values of TC and TOC and the highest PLI value.

Subcluster B2 included all the other stations, and can be further subdivided into two subclusters, B2a and B2b. The subcluster B2a grouped a few stations located in front of the Varano Channel, in the southern margin and in the central part of the lake. It was dominated by M. allgeni, C. honestus and D. normandicum, followed by S. pulchra. This subcluster showed the lowest values of Pielou and richness, if cluster A was excluded. It was also dominated by 1A and secondarily by 1B and 2A feeders and by c-p2 and c-p3 taxa. It was also characterized by sandy substrates with a high abundance of silt. It had higher values of OM and the highest value of TS. B2b included plenty of stations located in both the outer margins and the central parts of the lake, which are mainly represented by M. allgeni, along with T. longicaudata, L. luridus and A. torosa. This subcluster showed the highest values of sand and MI, and was also dominated by 1A and, to a lesser extent, by 1B and 2A feeders. It was also dominated by c-p3 taxa. This subcluster was characterized by the highest salinity values and sand-silt dominated sediments with the highest percentages of CaCO₃. It also had the highest values of TOC, TN, TC, and TH.

The eigenvalues of the PCA revealed that the first two factors explained ca. 49.3% of the total variance, whereas the factor variable plans showed how the species (those used for the CA and representing the main variables) and secondary variables (abiotic) were related to the different factors (Fig. 4). Factor 1, which explained 33.3% of the data variability, appeared to be related to the sediment characteristics in terms of grain size. Accordingly, silt and sand were placed on the opposite sides of the factor plan.

Strong relationships among the fine sediment fractions of both silt and clay were evidenced with T. longicaudata, L. luridus and the PLI. On the other hand, an opposite behaviour was highlighted among the finer fractions and the relative abundances of some taxa, including A. caudostriatus, P. pellucida, P. longicaudatus and O. wieseri. Factor 2, which justified 16.0% of the data variability, was mainly related to some of the physicochemical parameters of the bottom water, namely salinity, and the OM quality and quantity within the sediments (TOC, TC, TN, TS and TH). In particular, this factor seemed to be positively related to Microlaimus sp. 2, C. honestus, M. allgeni, and the MI, and negatively with D. normandicum, Thalassomonhystera parva, S. pulchra, N. papillosa and the 1B feeders.

The description of the distribution patterns of species is one of the fundamental starting blocks in the ecology of biological assemblages. However, there is very little information available on nematode taxonomic composition in Italian TEs and in the entire southern part of the Adriatic Sea (Semprucci 2013).

The nematode assemblage of the study area appeared to be mainly comprised of representative species of the fine sediments (Heip et al. 1985, Travizi and Vidakovic 1997) or of TEs (Villano and Warwick 1995, Pallo et al. 1998, Fabbrocini et al. 2005, Barnes et al. 2008). The overall richness revealed values that were highly comparable with those reported for other Italian TEs (Guerrini et al. 1998, Fabbrocini et al. 2005), but were generally lower than those reported for European brackish-water systems (e.g. Hendelberg and Jensen 1993, Barnes et al. 2008, Ferrero et al. 2008). This is probably due to the low salinity gradient of Lake Varano (Frontalini et al. 2013), which does not determine a transition from a marine to typical freshwater fauna, leading to a lower level of richness than in the European TEs.

Studies carried out in both marine and coastal lagoon ecosystems have revealed that the quantity and quality of the OM and salinity are among the most important factors influencing meiofaunal and nematode assemblages (Gambi et al. 2003, Barnes et al. 2008). However, the distribution of the assemblages is often more complex, reflecting multiple interactions between the distance to the sea, water depth, grain size, hydrodynamic turnover time, sediment oxygen and nutrient availability (Frontalini et al. 2013).

As shown in the PCA, the sediment grain size appeared to be one of the most important environmental variables affecting the nematode species distribution of Lake Varano. This parameter may, in fact, exert significant influence on the structure of the nematode assemblage, which appears to be very sensitive, even to slight sedimentological variations (Steyaert et al. 2003, Vanaverbeke et al. 2002, Semprucci et al. 2010a). In particular, the increasing silt content appeared to have the greatest effect on the lake’s nematodes. The presence of silt, even in small amounts, can actually reduce the permeability of sediment. Furthermore, a greater retention of OM is often observed in silty sediments (Albertelli et al. 1999, Frontalini et al. 2011, Semprucci et al. 2013), leading to OM accumulation and oxygen deficiency (Graf 1992). Accordingly, total OM and, in particular, its refractory component, increased in the study area in parallel with the silt fraction, while the highest ORP values detected were only associated with sand (Frontalini et al. 2013). As suggested by Vanaverbeke et al. (2011), this may prove that grain size only indirectly shapes nematode assemblages, and the biogeochemical environment resulting from the interplay between hydrodynamic features and grain size is a key factor structuring the assemblages. As shown by the PCA, the positive correlation between silt and the PLI is noteworthy. Spagnoli et al. (2002) and Fabbrocini et al. (2010) have already evidenced a correlation between trace elements and the fine sediment in Lake Varano. This greatly depends on the adsorptive properties of clay and silt particles, and supports the notion that grain size is also important in influencing pollutant accumulation in sediment (Bernardello et al. 2006).

As suggested by the CA and the PCA, the species closely associated with the silt were mainly M. allgeni, T. longicaudata and L. luridus. The high abundance of Desmodoridae in the study area was mainly due to M. allgeni. The Molgolaimus genus is present in all of the marine environments, ranging from shallow waters to the deep sea. It is generally associated with muddy sediment (e.g. Lambshead et al. 2000, Muthumbi et al. 2004) and its opportunistic life style has been documented (Schratzberger et al. 2003). In the Adriatic Sea, the presence of M. allgeni has been recorded in the Venice Lagoon, as well as in the offshore area of the Rovinj-Po River mouth and the Pored-Venezia transect (Travizi and Vidakovic 1997).

Terschellingia longicaudata is a typical inhabitant of estuarine, intertidal and subtidal areas, with its highest abundances being in muddy sediment (Warwick and Gee 1984, Barnes et al. 2008, Semprucci et al. 2010b). It is also well known as tolerant to a variety of natural and anthropogenic stressors (e.g. Armenteros et al. 2009, Gollner et al. 2010, Guilini et al. 2012). The physiological and behavioural adaptations of T. longicaudata to poorly oxygenated environments include a low respiratory rate and slow movement (Warwick and Gee 1984, Guerrini et al. 1998). The deposition of insoluble metal sulphides in intracellular inclusions in T. longicaudata has been suggested as a mechanism of detoxification of sulphide (Nicholas et al. 1987). Leptolaimus is another genus that is frequently and widely distributed in TE systems (Heip et al. 1985, Pallo et al. 1998, Ferrero et al. 2008). In particular, L. luridus has been found in muds characterized by the presence of organic detritus (Timm 1963).

Sabatieria pulchra inhabits all types of sediment, but its highest densities are in mud (Hendelberg and Jensen 1993, Steyaert et al. 2007). It is generally present in very shallow down to deeper waters. In anoxic and strongly reduced sediment, it is often the only remaining species. It dwells deep in the sediment, and is often found in high densities near the redox potential discontinuity layer (Wetzel et al. 2002, Steyaert et al. 2003, 2005). Data from a heavily polluted, tropical semi-closed bay (Havana Bay) revealed an enhancement of S. pulchra over T. longicaudata, probably due to the increase in grain size (Armenteros et al. 2009). This could explain the lack of a close correlation between S. pulchra and silt. However, the same authors also referred to possible interspecific competition for food or space between S. pulchra and T. longicaudata.

As highlighted by the statistical analyses, the enhancement of ORP and the sand percentage values led to a higher abundance of A. caudostriatus, O. wieseri, P. pellucida and P. longicaudatus, especially in the area in front of the baymouth bar. With a progressive increase in grain size (and, above all, a decrease in the silt-clay content), Axonolaimidae and Cyatholaimidae, and some representative species of Xyalidae, become increasingly abundant (Heip et al. 1985). Indeed, the presence of these taxa has previously been documented, particularly in intertidal or sublittoral sandy substrata (Warwick 1971, Semprucci 2013). Diversity indices, unlike the trends observed in the literature (Vanaverbeke et al. 2002, Semprucci et al. 2010a,b), were not significantly affected by grain size, maintaining levels that were comparable to those recorded for other TEs. On the other hand, the relative proportions of the feeding guilds appeared to be more related to the sediment features. This depends on the nature of the available food, which in turn depends on the nature of the sediment. In particular, the 1A and 2A feeders were significantly affected by silt and sand increases, respectively. The increase of 1A, mainly bacteriovorus, in the silty sediment may be related to the higher abundances of microorganisms in the finest sediments (Albertelli et al. 1999, Semprucci et al. 2010b). In sand, a high benthic primary production, a great abundance of diatoms, and wide surfaces suitable for scraping off the algal and bacterial biofilms can be found (Semprucci et al. 2010a,b), explaining the association between the 2A feeders and well-sorted fine sand. The only c-p class correlated with sediment type was c-p3 (intermediate colonizers), whose numbers increased with the silt amount, along with the PLI and OM enrichment (Bongers et al. 1991).

In the second factor of the PCA, TN, TOC, TH, TS and salinity were the most important environmental variables. Salinity was not strongly correlated with the nematode assemblages. As stated above, this could be due to the small variations in salinity detected in the area, and also confirms the macrofaunal spatial distribution observed by Specchiulli et al. (2010). However, Barnes et al. (2008) suggested that to document the controlling action of salinity on nematode assemblages, longer-term salinity data are fundamental. As highlighted by the CA and the PCA, C. honestus, which was reported by Ferrero et al. (2008) as being abundant in mixoeuhaline zones, was one of the species more positively affected by salinity. In contrast, T. parva, D. normandicum and S. pulchra appeared to be more negatively affected by this parameter. This is in line with the results of Barnes et al. (2008), who reported the first two species in an oligohaline salinity regime and the third in a polyhaline regime. Daptonema normandicum was also reported to be widely distributed in the Thames estuary by Ferrero et al. (2008) and in the ‘Valli di Comacchio’ complex (Northern Adriatic Sea) by Guerrini et al. (1998).

Pielou and Simpson indices seemed to be more affected by the second factor plan of the PCA, increasing with the reduction in the OM supply and salinity. This could be related to the discharge from the San Antonino and Francesco Canals, as documented by Specchiulli et al. (2008). As also suggested by the CA, the 1B feeders were widely distributed in the lake. However, the PCA showed that they were, as an opportunist trophic group, more associated with the freshwater inputs and secondly to sand (Semprucci and Balsamo 2012). Moreover, the MI, along with the other c-p classes, was more affected by the second factor plan, exhibiting a notable correlation with the quality and quantity of OM. Indeed, it appeared to be positively correlated with OM content and salinity enhancements. The lower values of the MI and the increases in abundance of the extreme and general colonizers (c-p1 and 2) in relation to the lower salinity values could corroborate the negative effects of freshwater inputs on the Varano system.

The nematode assemblage in Lake Varano appeared to be mainly represented by colonizer species. This group has been reported to increase in abundance under stress conditions (Bongers et al. 1991). However, we must highlight that the study area did not show higher values of trace elements when compared with the Effect Range-Low and Effect Range-Median values (Long et al. 1995), as well as high OM enrichment.

Although the MI may generally give good results in ecological assessments (Balsamo et al. 2012, Semprucci et al. 2013), Armenteros et al. (2009) reported unclear results in Cienfuegos Bay (Caribbean Sea), especially for the interaction of anthropogenic and natural factors. A clear relationship between anthropogenic activity (PLI) and the MI was not found in Lake Varano. This could be explained by the environmental complexity of transitional habitats as well as the necessity to update the life style of some species.

In order to classify the EQS of Lake Varano in accordance with the Water Framework Directive (WFD, Directive 2000/60/EC), we used the thresholds for nematode descriptors proposed by Moreno et al. (2011). This is the only tool available to classify the EQS by means of marine nematodes and it was integrated considering the percentage of all the c-p classes (Semprucci et al. 2014) (Fig. 5). On the basis of these thresholds, the EQS of most of the stations can be mainly classified as moderate, followed by good and poor, which is in agreement with the level of organic enrichment and the trace element data (Fig. 5). In detail, moderate conditions of the assemblage were generally revealed in the central area of the lake, while in front of the San Antonino and Francesco Canals, where Specchiulli et al. (2010) also documented the negative influence of the water discharges of the canals, a poor EQS was detected. The notable difference revealed between the western and eastern part of the lake could be related to the stronger influence of the Capoiale Channel and the poor mixing action of the Varano Lake, and this fits perfectly with the findings of Specchiulli et al. (2008) and Frontalini et al. (2013). Although temporal replications are lacking, this research enables there to be a preliminary evaluation of the current ecological conditions of the lake. It may also provide a baseline for the future monitoring of this TE and for assessing long-term changes therein.