We sampled 13 sites located in extensive seagrass beds (at least 1 km² of extension) in seven areas around the Cuban archipelago (Fig. 1 and Table 1). Six sites were exposed to the direct influence of oceanic waters (i.e. relatively close to the shelf border and no coral reef crest separating the bed from the shelf border): GG, RG, ON, OS1, OS2 and GB3. The other seven sites were separated from open waters by chains of cays and coral reef crests (i.e. sheltered): SM, AM1, AM2, AM3, AM4, GB1 and GB2. We ensured the accuracy of this classification on the basis of in situ observations of the sediment type (e.g. silty or sandy), bottom slope and biota that may indicate physical disturbance (e.g. gorgonians, sponges). The seagrass beds were composed mainly of turtle grass (Thalassia testudinum) in soft bottoms with depth range from 1 to 5 m.

Figure 1.  

Map of the Cuban archipelago indicating the 13 sampled sites in seagrass beds. Dot colours indicate different areas in the Cuban shelf.

We could not sample the whole extension of the archipelago in a single expedition; actually, sampling events were done within an interval of five years (see dates in Table 1) in seven expeditions with different goals. This is the reason for small variations in the field protocol, such as the number of collected cores or the use of preserving agent. Sediment cores (4–8 per site) of internal diameter 2.5 cm were collected by SCUBA divers in an area of ca. 36 m² within each site and at a depth inside the sediment of 6 cm. The sediment was preserved on board in either 95% ethanol or 10% formalin.

Temperature, salinity and dissolved oxygen (DO) were measured in situ (but not at all sites) at approximately 10 cm from the bottom using an oceanographic Hydrolab multiprobe 4a instrument. Samples of the uppermost 3-cm layer of sediments were taken with a 250-ml propylene container for the measurement of grain size and total organic matter (TOM). Grain size was determined with the gravimetric method using a standard sieve column and an analytical balance (Loring and Rantala 1992) and was expressed as the percentages of mud (silt + clay, < 63 µm) and sand (> 63 µm). TOM was determined by the gravimetric method, measuring the loss of weight on ignition at 450°C for eight hours (Heiri et al. 2001).

The geographical distance between sites was calculated as the shortest distance across the sea to allow for potential dispersal by currents; i.e. no land barriers were crossed. Distances were calculated in the software OpenCPN 4.2.0 with nautical charts provided by GEOCUBA Nautical Cartography Agency.

In the laboratory, sediment was sieved with filtered water (32 µm) through a nested set of sieves of 45 and 500 µm of mesh size to separate the fractions of meiofauna and macrofauna, respectively. The material retained in the sieves was preserved in ethanol 70% (Evans and Paulay 2012) and both fractions were analysed to extract the nematodes. The sediment containing the organisms was mixed gently with a solution of sugar and water (density 1.15 g cm^-3) following the procedure in Armenteros et al. (2008). This procedure allows the extraction of the organisms from the sediment by difference of density (i.e. sediment settles on the bottom and organisms float on the surface of the supernatant). The procedure was repeated three times, the supernatant being concentrated in a volume of ca. 25 ml and observed under a stereomicroscope Olympus SZX7 with magnification between 15 and 115x. All the nematodes in the samples were picked up with a sleeved needle ending in a hook and preserved; however, only a subset of nematodes were included in this study because other individuals were diverted for DNA analyses. We collected between 100 and 400 nematodes, the first to be observed in the counting chambers (i.e. no size or shape bias) and preserved them in 10% formalin. This splitting of the samples prevented a quantitative assessment of nematode abundance; so we necessarily focused on a nematode subset and not on the whole assemblage. Selected nematodes were mounted in permanent glass slides for microscopy, following the technique in Vincx (1996).

Nematodes were identified to the lowest possible taxonomic level, following the taxonomic literature, such as those by Platt and Warwick (1983), Platt and Warwick (1988), Warwick et al. (1998) and Schmidt-Rhaesa (2014). In many cases, the name of some species could not be assigned with reasonable accuracy because there were single individuals and/or they did not match with any described species within the genus (i.e. likely a new species); therefore, these species were referred as sp. in the text. If more than one morphospecies were recorded within a same genus, they were named with a number (i.e. sp. 1, sp. 2). For some complex genera with very small-sized nematodes and those, such as Desmoscolex and Tricoma, we could not discriminate beyond doubt between species and we preferred to group them as spp. (i.e. more than one species within the genus).

Nematodes were classified according to two different functional traits (feeding groups and coloniser/persister ability) (Suppl. material 1). The classification by Wieser (1953) of nematodes in four trophic categories, based on the morphology of the buccal cavity, was used: selective deposit feeders (1A), non-selective deposit feeders (1B), epigrowth feeders (2A) and predators/omnivores (2B). We also used the classification of nematodes, based on the coloniser/persister ability (Bongers 1990, Bongers et al. 1991) from c-p 1 (best colonisers) to c-p 5 (bad colonisers).

Accumulation curves of richness versus individuals were built in the software EstimateS 9.0 (Colwell 2013). Richness was assessed using the observed number of species and the non-parametric Chao 2 estimator that accounts for the unseen taxa. The observed number of species and the rarefied richness for 100 specimens were calculated for each site. The associated 0.95 confidence intervals (CI) were calculated, based on 100 permutation in EstimateS 9.0. Overlapping of CIs was used as a relaxed indication of non-significant differences between sites. The association between richness and abundance across all the sites was calculated with the Spearman rank correlation coefficient (R_S).

We used the complement of the Sorensen Similarity Index as a measure of β-diversity (i.e. 100 – Sorensen) between pairs of sites. Based on the triangular matrix of similarities, we did a numerical ordination of the sites by non-metric multidimensional scaling (NMDS) to visualise potential patterns in β-diversity. We used the Sorensen Index, which relies on presence/absence data, because only a subset of nematodes in the samples were analysed. Statistical significance between the oceanographic regimes (two states: exposed versus sheltered) was undertaken using an Analysis of Similarity (ANOSIM) in the software PRIMER 6.1.15 (Clarke and Gorley 2006).

An analysis of linear regression was carried out to test the dependence of the Sorensen pairwise similarity with the geographical distance. Here, we used the Sorensen similarity measure (instead of dissimilarity) because previous studies have tested the distance decay model using similarity measures. A Mantel-type test of significance, using a non-parametric approach, was made using the routine RELATE in the software PRIMER 6.1.15 to test if both triangular matrices were significantly related.

Salinity and dissolved oxygen in the water column had a narrow range (35–37 PSU and 6.1–8.6 mg l^-1, respectively), typical of marine and well-oxygenated shallow waters (Table 1). The fraction of fine sediment (grain size < 63 µm) in sediment had broad differences between exposed (mean: 6.5%) and sheltered sites (mean: 55%). The content of organic matter was also broadly different with exposed sites having poorer content of carbon in sediments (mean: 6%) compared to sheltered sites (mean: 17%) (Table 1).

We identified 2678 nematodes belonging to 215 species, 138 genera, 34 families, seven orders and two classes. The observed species richness at local scale (α-diversity) had a median (± interquartile range, n = 13) of 57 ± 17 species (range: 31–88 species). Species richness did not show significant differences between most of the sites as indicated by the broad overlapping of the confidence intervals; but GB3 and OS2 had the lowest and highest values of α-diversity, respectively (Fig. 2a). Since species richness was strongly influenced by the number of identified specimens (Spearman correlation, R_S = 0.82, p < 0.001, n = 13), we calculated the expected number of species (ES) for a sample size of 100 specimens, using rarefaction. ES₍₁₀₀₎ was not correlated with the abundance (R_S = 0.05, p = 0.87, n = 13). The median ES₍₁₀₀₎ was 42 ± 6 species (range: 27–49 species). The broad overlapping of the 0.95 confidence intervals for ES₍₁₀₀₎ suggests no significant differences in the species richness between sites; albeit the significant differences remained between GB3 and OS2 (Fig. 2b).

Figure 2.  

Species richness of free-living marine nematode assemblages in seagrass beds. (a) Observed species richness per site (α-diversity) with 0.95 confidence intervals. (b) Expected number of species (ES) rarefied to a sample of 100 nematodes per site. (c) Accumulation curves of observed species richness and the non-parametric Chao 2 estimator with all sites combined (γ-diversity).

The species richness at regional scale (γ-diversity) of free-living marine nematodes in seagrass beds was estimated from the combination of the 13 studied sites. The curve of accumulation of observed species richness did not approach to an asymptote (Fig. 2c) and the total richness (± SD) was 215 ± 6 species (0.95 CI: 204–226 species). However, the curve of Chao 2 non-parametric estimator approached to an asymptote and the estimated species richness was 253 ± 13 species (0.95 CI: 235–287 species). The checklist of species contributing to the regional inventory and their occurrence across the 13 sites is given as supplementary material (Suppl. material 2).

The more abundant trophic groups (median ± interquartile range) were predator/omnivores (2B, 35% ± 15%) and epigrowth feeders (2A, 34% ± 7%). The third most abundant trophic group was non-selective deposit feeders (1B, 19% ± 11%) and the less abundant was selective deposit feeders (1A, 14% ± 5%). According to the coloniser/persister scale, the nematodes with intermediate colonising ability were the most abundant (c-p 3, 47% ± 13%), followed by nematodes with high colonising ability (c-p 2, 40% ± 13%). Nematodes with low colonising abilities had lower abundance (c-p 4, 12% ± 4%) and nematodes, with the lowest colonising ability, were the least abundant (c-p 5, 0.5% ± 0.9%).

The dominance was moderate with only 23 species (11% of total richness) accounting for 51% of the total accumulated abundance. The most abundant species were Paradesmodora immersa Wieser, 1954 (4%); Desmodora pontica Filipjev, 1922 (4%); Viscosia abyssorum (Allgén, 1933) (3%); Dorylaimopsis punctata Ditlevsen, 1918 (3%); Daptonema sp. (3%); Marylynnia sp. (3%), Euchromadora vulgaris Bastian, 1865 (3%); Halichoanolaimus chordiurus Gerlach, 1955 (3%); and Zalonema ditlevseni (Micoletzky, 1922) (3%). These species belonged to five orders and seven families widely distributed in marine habitats, namely: Chromadoridae, Comesomatidae, Cyatholaimidae, Desmodoridae, Oncholaimidae, Selachinematidae and Xyalidae.

The values of pairwise dissimilarity (β-diversity) between sites ranged from 40 to 67% with an average of 55%. The ordination of the sites, based on the presence/absence of species, did not indicate any clustering of sites, based on the geographical area or oceanographic regime (i.e. exposed vs. sheltered) (Fig. 3a). The multivariate assemblage structure between regime conditions was not significantly different (ANOSIM, R = 0.08, p = 0.25, 999 permutations).

Figure 3.  

Species composition of free-living marine nematode assemblages in seagrass beds. (a) Ordination of 13 sites in seagrass beds of the Cuban archipelago, based on presence/absence of free-living nematode species. (b) Relationships between pairwise Sorensen similarity and geographical distance. (c) Relationships between Sorensen similarity and geographical distance with sites separated by the oceanographic regime.

Geographical distance may be a driver of β-diversity, namely, the model distance decay of similarity (DDS). To explore the adjustment of DDS model of our data, we computed 78 shortest geographical distances amongst the 13 sites. The geographical distance amongst sites had a median of 599 km (range: 5–1269 km). The Sorensen pairwise similarity between sites was significantly, but weakly, related with the geographical distance (linear regression, slope = -0.01 ± 0.003, p < 0.001, R² = 0.22, n = 78) (Fig. 3b). A non-parametric approach using a Mantel Test indicated that both triangular matrices (Sorensen similarity and geographical distance) were significantly associated (RELATE, Rho = 0.46, p = 0.004, 999 permutations).

We also explored the relationships between Sorensen similarity with geographical distance independently for exposed and sheltered sites. Exposed sites lacked a DDS as indicated by the non-significant relationship between Sorensen similarity and geographical distance. However, sheltered sites had a DDS relationship with a significant relationship between similarity and distance (Fig. 3c and Table 2).

The DDS occurred for the four tested trophic groups: deposit feeders (1A and 1B), epigrowth feeders (2A) and predator/omnivore (2B) (Fig. 4, left column and Table 2). The slope of the relationship similarity-distance for these groups was statistically significant and similar in magnitude with a decay of ca. 1% of similarity for each 100 km. The DDS was also significant for all the categories of nematodes classified after the coloniser/persister scale (Fig. 4, right column and Table 2). However, the slope of poor colonisers (i.e. c-p 4 & 5) was two times greater (ca. 2% change of similarity each 100 km) than good colonisers (ca. 1% change each 100 km).

Figure 4.  

The distance decay of similarity for pairwise similarity values (Sorensen Index) of free-living marine nematode assemblages and shortest geographical distance between sites. Sorensen similarities are based on matrices of functional traits. Left column: Trophic groups (selective deposit feeders = 1A, non-selective deposit feeders = 1B, epigrowth feeders = 2A and predator/omnivores = 2B). Right column: Coloniser/persister scale (c-p 2 good coloniser/poor persister to c-p 5 bad coloniser/good persister).

The γ-diversity of nematode assemblages was high as expected in seagrass beds. We have reported higher species richness than other studies in tropical seagrass beds [e.g. 100 species in Hopper and Meyers (1967b) from four sites, 152 species in Fisher and Sheaves (2003) from four sites]. Most of the studies of nematode diversity in seagrass beds have used genus as the lowest taxonomic level. The diversity at genus level reported in our study (138 genera) is higher than other reports in temperate seagrass beds [e.g. 88 genera in Danovaro and Gambi (2002) within a single Posidonia bed, 58 genera in Materatski et al. (2015) from two sites]; and also higher than other estimates in tropical seagrass beds [e.g. 100 genera in Ndaro and Olafsson (1999)from three sites, 63 genera in Liao et al. (2015) within a single Thalassia bed]. The most plausible explanation is the large spatial coverage of our sampling scheme that included 13 sites in seagrass beds from all the sub-regions in the Cuban archipelago (i.e. separated in median 599 km). Another factor that could promote the high species richness is the variety of oceanographic regimes in the sampling sites, ranging from widely exposed (OS1 and OS2) to well-sheltered (AM1 and AM2). The non-asymptotic shape of the curves of accumulated number of species suggests that there is still undiscovered diversity. The non-parametric Chao 2 estimator indicates an as yet unseen diversity as high as 253 species that can be interpreted as the lower bound of diversity (Gotelli and Colwell 2011) in the studied seagrass beds.

The significant differences in α-diversity between sites could be explained by differences in the number of identified specimens and/or environmental conditions at local scale. However, when we standardised to an equal level of abundance, the species richness was the same across most of the sites. It suggests that the seagrass beds around the Cuban archipelago harbour a rather similar species richness, despite the environmental heterogeneity. This is consistent with the lack of clustering of the samples, based on assemblage species composition, by geographical area or exposure to oceanic influence.

The recorded broad phylogenetic scope of the nine dominant nematode species in our study (belonging to five orders and seven families), high species richness and even occurrence of all trophic groups, reflect the combination of: (i) heterogeneous composition of food sources (Danovaro and Gambi 2002), (ii) the provision of sheltered environment by the seagrass blades (Leduc and Probert 2011) and (iii) the diversity of microhabitats in seagrass beds (Hall and Bell 1993). This combination of enhancing factors supports the view that seagrass beds are diversity hotspots for free-living marine nematodes with higher richness, compared to other habitats (Armenteros et al. 2019).

The geographical distance influenced the pairwise similarity in the studied assemblages, albeit with weak effect. The DDS can be caused by either a decrease in environmental similarity with distance or by limits to dispersal and niche breadth differences amongst taxa (Nekola and White 1999). The separation between seagrass beds in our study (median 599 km) likely constitutes limits to dispersal. Distances in the order of more than 100 km constitute limits to the dispersal of marine nematodes as evidenced from genetic (Derycke et al. 2013) and ecological data (Pérez-García et al. 2019). However, aquatic nematodes may disperse to long distances by the combination of diverse modes of dispersal (e.g. passive dispersal by currents, zoochory) with rapid reproduction (Ptatscheck and Traunspurger 2020). The differences in the DDS pattern between exposed and sheltered sites suggest that dispersal by hydrodynamics plays a significant role. Stronger hydrodynamics in exposed sites putatively resulted in a weaker DDS when compared with sheltered sites where weaker hydrodynamics may pose limits to the passive dispersal. However, alternative explanations, such as the difference in the pool of organic carbon in sediments (lower in exposed sites) and/or grain size (higher in exposed sites), cannot be ruled out and deserve further study using abundance data.

The DDS patterns occurred in similar rates for the four trophic groups (i.e. diminution of similarity with the increase in geographical distance). According to Danovaro and Gambi (2002), the nematode trophic group composition in seagrass beds was mainly determined by the quality and quantity of available food in sediments. In our study, it is unclear if DDS patterns, based on trophic groups, were affected by the distance itself or was mirroring the decay of similarity of species composition. Two further issues, related to the use of trophic group classification, could affect our results. First, Wieser’s classification does not characterise effectively the high feeding selectivity and flexibility of nematodes (Moens et al. 2014). Second, swimming behaviour of the species, even within the same Wieser’s feeding group, affects the susceptibility to passive re-suspension as indicated by experimental evidence (Thomas and Lana 2011).

The differential colonising ability of nematodes significantly affected the distance decay of similarity. The similarity pattern of poor coloniser nematodes (c-p 4 and 5) was twice as much affected by the geographical distance than the good colonisers (c-p 2 and 3). This suggests that the organismal features used for c-p scale, namely, generation time, reproduction rate and body size (Bongers 1990, Bongers et al. 1991) may affect the distance decay patterns of nematodes. Likely, large body-sized nematodes with low generation times, few offspring and low metabolic rate (e.g. Cylicolaimus magnus, Desmoscolex spp., Enoplus sp., Leptosomatum sp. and Polygastrophora maior) were less effectively dispersed and/or settled with lower success to adjacent sites. On the other hand, good colonisers with small body size, producing many offspring and with higher metabolic rate, show a weaker distance decay likely due to more efficient dispersal by currents and/or successful settlement. Thus, according to Poulin (2003): “distance matters, but it is not always the key determinant of similarity”

The regional species richness of free living nematode assemblages accounted for 215 species and 34 families. This γ-diversity was higher than other estimates from tropical and temperate regions and points to seagrass beds as diversity hotspots of free-living marine nematodes. Local species richness in seagrass sites was about 57 ± 17 species. The geographical distance played a weak, but significant, role in the decay of similarity and was likely affected by the limited dispersal of nematodes. Pairwise similarity of bad-coloniser nematodes was twice as much affected by distance than good-colonisers possibly due to differential success of transport and settlement.