Diversity patterns of free-living nematode assemblages in seagrass beds from the Cuban archipelago (Caribbean Sea)

Abstract Diversity patterns of free-living marine nematodes in tropical seagrass beds are understudied. Here, we describe the species richness and assemblage composition of nematodes in 13 seagrass sites covering the whole Cuban archipelago. Nematodes were collected from Thalassia testudinum seagrass beds and identified to species level. We provide a checklist of nematode species from seagrass beds. The species richness of nematode assemblages is high, with 215 species, 138 genus, 35 families, seven orders and two classes. That γ-diversity is higher than other studies and points to seagrass beds as diversity hotspots of free-living marine nematodes. Local species richness in seagrass bed sites is about 57 ± 17 species and broadly similar across the sites, despite the environmental heterogeneity. The geographical distance plays a weak, but significant, role on the decay of similarity likely affected by limited dispersal of nematodes. The pairwise similarity values, related to poor-coloniser nematodes, were twice more affected by the distance than those related to good-colonisers, possibly due to differential success of transport and settlement.


Introduction
Nematodes constitute the fourth most diverse phylum of metazoans on Earth, after Arthropoda, Mollusca and Platyhelminthes (Zhang 2013). The number of accepted species of free-living marine nematodes is about 6219 species (Bezerra et al. 2020), but the percentage of known species is as low as 12% (Appeltans et al. 2012). This gap in knowledge is even larger in tropical regions since there are many more studies on nematode diversity in Europe and North America compared with Africa and Central America, despite tropical ecosystems harbouring some of the most diverse habitats on the Earth, such as corals reefs (Cote and Knowlton 2014) and seagrass beds (Duffy et al. 2014).
Seagrass beds, when compared with unvegetated adjacent habitats, harbour larger nematode species richness (Ruiz-Abierno and Armenteros 2017) and different species composition (Fonseca et al. 2011). However, few studies about nematode diversity in tropical/subtropical seagrass beds have been undertaken in comparison with intertidal and estuarine systems. In tropical seagrass beds, nematode assemblages have been studied by Hopper and Meyers (1967a), Hopper and Meyers (1967b), Ndaro and Olafsson (1999), Fisher (2003), Fisher and Sheaves (2003) and Liao et al. (2015). In temperate seagrass beds, other studies on nematode diversity have described relationships with food availability in Posidonia beds (Danovaro and Gambi 2002) and the assemblage response to a collapse of Zostera beds (Materatski et al. 2016, Materatski et al. 2015, Branco et al. 2018. Seagrass beds provide an array of advantageous conditions to the meiofauna (and nematodes), such as physical protection against re-suspension, food availability derived from seagrass production and diverse microhabitats (Bell et al. 1984, Decho et al. 1985. This combination of sedimentary conditions/resources provides intermediate levels of resource availability and physical disturbance which, in turn, promotes the meiofaunal richness (Armenteros et al. 2019).
The difference patterns in species composition across samples, also termed as β-diversity, is a central theme in community ecology (Anderson et al. 2011). One pattern of β-diversity, suitable to explore in our study, is the distance decay of similarity (DDS, Nekola and White 1999). For meiofauna, patterns of β-diversity have been explained by a combination of niche assembly and dispersal processes (Giere 2009). Niche assembly processes have been the subject of many studies relating the nematode assemblage structure with the environment. Some of these studies explicitly addressed the niche breadth of nematodes (e.g. Wu et al. 2019) and other related assemblage structure with environmental factors looking for ecological drivers [see Heip et al. (1985) and Moens et al. (2014) for reviews]. Nematodes have the ability to effectively disperse at small scales (< 10 m) by both active swimming and passive transport (Ullberg andÓlafsson 2003, Thomas andLana 2011).
Dispersal of marine nematodes is substantial at scales of 10-100 km as indicated by genetic data (Derycke et al. 2013). A recent review (Ptatscheck and Traunspurger 2020) concluded that nematodes are effective colonisers due to the diversity of dispersal modes, continuous immigration and rapid reproduction.
In addition to taxonomic diversity, the functional diversity of assemblages can be addressed on the basis of biological traits. These traits refer to morphological, physiological or phenological features, measurable at the individual level (Violle et al. 2007). The most used biological traits of marine nematodes have been the trophic group and the coloniser/persister ability (Schratzberger et al. 2007, Armenteros et al. 2009, Alves et al. 2014). The trophic classification was proposed by Wieser (1953), based on the structure of stoma and, later, refined by Moens et al. (2004). The c-p classification was proposed as an arbitrary scale to measure the ability to colonise/persist of nematodes, based on food preference, gonad size, metabolic rate, mobility and occurrence of larva dauer (Bongers 1990, Bongers et al. 1991.
In this study, we contribute to the knowledge of the diversity of free-living marine nematodes in seagrass beds from the Cuban archipelago, one of the hotspots of diversity in the Caribbean Sea. Therefore, the aims of our research are: (1) Quantify the species richness of free-living marine nematodes from seagrass beds at the scale of a single site (α-diversity) and of the whole archipelago (γ-diversity). We provide a checklist of species that accounts for the γ-diversity of Cuban seagrass beds; and (2) Explore the β-diversity patterns and the effects of two potential environmental drivers: the geographical distance (c.f. distance decay of similarity) between the studied sites and the differential oceanic exposure of the sites (lagoon versus open shelf). In addition, the analysis of two biological traits of nematodes (trophic group and colonising capacity) adds another dimension to the analysis of β-diversity.

Study sites
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.

Sampling
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.

Processing of samples
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 ) 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 -3 Figure 1.
Map of the Cuban archipelago indicating the 13 sampled sites in seagrass beds. Dot colours indicate different areas in the Cuban shelf.
Diversity patterns of free-living nematode assemblages in seagrass beds ...
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).

Data analyses
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 ).
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

Data resources
The data underpinning the analysis reported in this paper are deposited at GBIF, the Global Biodiversity Information Facility, http://ipt.pensoft.net/resource.do?r=xxxxxx

Abiotic factors
Salinity and dissolved oxygen in the water column had a narrow range (35-37 PSU and 6.1-8.6 mg l , 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).

Species richness
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 = 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 = 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).

Assemblage composition and β-diversity
The dominance was moderate with only 23 species (11% of total richness) accounting for 51% of the total accumulated abundance. 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).
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).

Species richness
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) Table 2.
Testing the distance decay of similarity after oceanographic regime of the sites; and after the functional traits of free-living marine nematode assemblages. Parameters of the linear regression between the Sorensen Similarity Index and the geographical distance are given. Trophic groups are based on the structure of the buccal cavity after Weiser (1953) and coloniser/persister are based on a scale from good (c-p = 2) to poor (c-p = 5) coloniser. Scales 4 and 5 were summed. Asterisk (and bold) indicates that the slope is significantly different of zero.
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 wellsheltered (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).

β-diversity
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 . 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"

Concluding remarks
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.