How nematode morphometric attributes integrate with taxonomy-based measures along an estuarine gradient

Nematodes are highly susceptible to environmental change and possess a wide array of morphological and functional characteristics for the assessment of the “ Good Environmental Status ” , within Marine Strategy Framework Directive. However, while the taxonomic sufficiency of nematodes in detecting spatial gradients and related ecological niche conditions is well recognized, very little is known about nematodes functional morphometric attributes in response to environmental drivers. To explore this knowledge gap, we aimed to assess the efficacy and efficiency of nematode morphometric attributes (length, width, length/width ratio, biomass) in detecting spatial patterns along a Portuguese estuarine gradient, and compare it with the taxonomic approach. We hypothesized that abundance data weighted by the morphometric attributes will have a higher explanatory power in detecting spatial patterns than using abundance of morphometric data alone. Based on the recent recommendations regarding the time and cost related efficacy of methods in biomonitoring and ecological assessments we also hypothesized that a reduced dataset based on the most common genera will suffice to capture the same distributional patterns displayed by the whole assemblage. Our results demonstrated that dataset solely based on genera abundances had consistently better explanatory power than combined datasets or morphometric datasets alone, however, combined dataset provided different spatial patterns and performed better at discriminating estuary areas. The main gradients described by the taxonomy-based dataset were related to the sediment particle size and water depth. Considering combined datasets, spatial discrimination was mainly driven by the variation in dissolved oxygen % saturation, pointing out to the importance of this variable in determining estuarine conditions substantial for nematodes morphometric distributional patterns. The same analysis repeated for the most frequent genera resulted in similar distributional patterns as for the whole assemblage dataset, clearly demonstrating that spatial estuarine gradients can be sufficiently described by using only the most frequent genera. Such information may substantially increase the efficiency of bio-assessment surveys by reducing the cost and work associated with identification and measurements of all of the individual nematode genera.


Introduction
Nematodes are key components of benthic ecosystems largely contributing to organic matter mineralization and energy provision to higher trophic levels (Coull, 1999;Schmid-Araya et al., 2016;Schratzberger and Ingels, 2018). These abundant metazoans are highly susceptible to environmental variation, which turns them efficient indicators of ecological change associated to natural or anthropogenically-induced disturbance (Höss et al., 2011;Moreno et al., 2011;Semprucci and Balsamo, 2012;Semprucci et al., 2015). However as far as we started to comprehend patterns in nematodes diversity and abundance measures, very little is known about nematodes diversityfunction relationship. This information is not only crucial to develop a solid foundation for nematode-based indicators, but also to better understand the larger role that nematodes play in benthic ecosystems through the interactions with other major taxonomic groups (Cronin-O'Reilly et al., 2018). Combination of taxonomic and functional information was already recognized to provide a promising tool for the future ecological research (Schratzberger et al., 2007;Semprucci et al., 2018). However, to date only few studies investigated the diversity-function relationship of nematodes in response to environmental drivers, leaving this research field largely unexplored (Alves et al., 2014;. Functional trait approach characterizes the organism through its biological (morphological, physiological or behavioural) attributes, thereby providing a direct link between community-based processes and the underlying environmental gradients (Cadotte et al., 2011). Functional traits are also independent of species identity providing information on ecosystem functions and services (McGill et al., 2006). Body size is the most fundamental trait as it directly relates to physiological processes such as growth rate or respiration (Brown et al., 2004) and so the individuals exhibit morphometric adaptations in regard to their surrounding environment e.g., water depth and food availability (Górska et al., 2020;Losi et al., 2013). Consequently, morphometric characterization of nematodes is a useful tool for the description of their functional diversity and for comparisons of different assemblages along diverse ecological conditions. Nematodes are the principal group of meiofauna demonstrating exceptional morphometric diversity (Soetaert et al., 2002;Vanaverbeke et al., 2003Vanaverbeke et al., , 2004Warwick et al., 1998). The same species may vary remarkably in body size and proportions depending on resource availability (Grzelak et al., 2016;Vanaverbeke et al., 2003), granulometry (Tita et al., 1999) and oxygen concentration (Losi et al., 2013). For example, silt and clay favour more slender and stout nematodes, whereas short and wide ones are usually indicators of well-oxygenated environments with high availability of organic matter (Tita et al., 1999). Differences in the distribution of nematode morphotypes are also evident at different habitat types (e.g., soft and hard substrates) shaped by the above-mentioned sediment properties (Armenteros and Ruiz-Abierno, 2015). Morphometric attributes were demonstrated to be an efficient tool in depicting harboured areas with contrasting environmental conditions (Losi et al., 2013). They also proved to be a good environmental descriptor to follow up the recovery of seagrass beds (Materatski et al., 2018).
Besides these promising results of using morphometric attributes as indicators of the ecological conditions, little was done to continue this research line. Most of the meiofauna and nematode indicator-oriented research is still largely founded on traditional taxonomic based approaches (Moreno et al., 2011;Semprucci et al., 2018) or alternatively, the nematodes behavioural traits such as feeding type or life history strategy traits (Mirto et al., 2002;Semprucci et al., 2016). Taxonomic based indices are time-consuming due to genus and species-level identifications (Trigal-Domínguez et al., 2010). In the case of feeding type and life history strategy traits, a number of studies has demonstrated that they are as powerful in detecting spatial patterns as the taxonomic approach (Alves et al., 2013;Semprucci et al., 2013). Nonetheless, it was also noted that functional traits are highly complementary to abundance-based approaches, providing different type of information (Danovaro et al., 2008;Schratzberger et al., 2007). The use of morphometric attributes does not require taxonomic expertise, thereby representing a potential tool for future assessment of ecological conditions particularly related to organic matter quality and sediment characteristics (Losi et al., 2013;Soetaert et al., 2002).
Beside these advantages, there are no studies comparing the efficiency of morphometric attributes in relation to abundance approaches to detect spatial patterns in nematode assemblages. Until now, the suitability of morphometric features as ecological descriptors has been studied across sites that demonstrated contrasting differences in environmental conditions, particularly in oxygen concentration (Losi et al., 2013), organic quality (Grzelak et al., 2016) or drastic changes in granulometry associated to local seagrass extinction (Materatski et al., 2018). While these studies provide clear evidence on morphotypes as indicators of change across sites with distinct ecological conditions, they do not provide clues to incorporate morphometric attributes along environmental gradients.
An effective introduction of nematode-based morphometric attributes into biological assessment requires standardized and robust methods, with minimum cost and time-related effort. This need for efficacy in biological assessment surveys has motivated the development of a specific research line on the minimum representative sampling effort sufficient to permit the detection of changes among assemblages (e.g., Sgarbi et al., 2020). Such reasoning is based on findings that spatial diversity patterns are primarily driven by the widespread species, rather than the rare ones (Checon and Amaral, 2017;Lennon et al., 2004). For example, in aquatic communities spatial turnover along environmental and spatial gradients was demonstrated to be similar for common and rare species concluding that common species are adequate descriptors of species turnover in regard to environmental variables (Heino and Soininen, 2010). Morphometric attributes do not require a biological expertise or laboratory experiments, such as the ones necessary to assess the qualitative type traits (e.g., feeding type or a life-span). However, how far nematode morphometric-based biological assessment can be optimized in regard to representative sampling effort (common genera vs whole assemblage) has not been determined yet. Particularly, it is of great interest to understand if the most frequent species hold the same morphometric distributional patterns, along environmental gradient, as the entire nematode assemblage.
To explore these gaps of knowledge we 1) aim to assess the efficacy and efficiency of nematode morphometric attributes (length, width, length/width ratio, biomass) in detecting spatial patterns along a Portuguese estuarine gradient, and compare it with the taxonomic approach. We hypothesize that the abundance data weighted by the morphometric attributes has higher explanatory power in detecting spatial patterns than using abundance or morphometric data alone.
Our second research objective is related to recent recommendation of Sgarbi et al. (2020), regarding the time and cost related efficacy of methods in biomonitoring and ecological assessments. As a response to this recommendation, we 2) aim to replicate the above analysis for the most common genera. We hypothesize that a reduced dataset based on the most common genera will suffice to capture the same distributional patterns displayed by the whole nematode assemblage.

Study area
The Sado estuary (38 • 27 ′ N, 08 • 43 ′ W) is located in the southwest coast of Portugal ( Fig. 1). It ranks as the second largest estuarine system in Portugal, having an area of 240 km 2 . It extends from east to northwest and is partially separated from the sea by the long (25 km) intertidal sandbank (Troia Peninsula) located on the west side. The main branch of the estuary diverges towards the northern direction forming a broad Bay, where the prevalence of low hydro-dynamism and high water residence time support the extensive aquaculture-related activities. The estuary's salinity is highly variable among seasons and tidal cycles. Salinity measured seasonally in the upstream area of the estuary might vary between 3 and 30 units, while next to the estuary mouth salinity ranges from 25 to 35 (Portela, 2016). Tidal cycles are semidiurnal with an amplitude range of 0.76-3 m (Bettencourt et al., 2004) and salinity variation along the tidal cycles can reach up to 22 units (Portela, 2016). Bottom water temperature ranges from 12 − 15 • C during the winter months, towards 18-20 • C during spring and might reach even 27 • C in the peak of summer. Daily temperature changes associated with tidal cycles vary from 2 • C in the summer to 0.5 • C in the winter (Cabeçadas et al., 1999;Portela, 2016). Major subtidal bottom types consist of both, sandy and silty clay sediments (Rodrigues and Quintino, 1993). Ecological conditions of the Sado estuary are highly heterogeneous. Northern part of the estuary is subjected to various sources of contamination: urban effluents from the city of Setúbal and its industrialized harbour area as well as runoff from upstream agricultural fields. Conversely, the southern and eastern parts of the estuary are classified as a natural reserve designated as "National Reserve of the Sado estuary".

Sampling methods
A subtidal sampling survey of the Sado estuary took place during May 2018. In order to account for the inter-tidal variations in the salinity, the sampling was conducted at three consecutive days, during the same period of the tidal cycle. Nematode samples were collected at 35 sampling stations along the estuarine gradient ( Fig. 1). Sampling included upstream locations, with a higher freshwater influence and downstream locations, closer to the estuary mouth, classified as euhaline stations.
Simultaneously with nematode samples, salinity, temperature (•C) and dissolved oxygen (DO) (mg L − 1 ) were measured near the bottom at each sampling station using a multiparameter probe (YSI Data Sonde Survey 4). Sediment samples (~100 g) were collected using a Van Veen grab and split to separately determine total organic matter (TOM) and sediment grain size. Sediment samples were oven dried for 72 h at 60 • C and subsequently combusted at 550 • C for 8 h. TOM was calculated as the difference between the total weights of dry sediment and noncombusted portion of sediment obtained through burning. TOM was expressed as a total % of organic matter. Sediment grain size was determined by sieving the collected sediments through a battery of different mesh sizes sieves. Grain sizes were assigned to five classes: gravel (>2 mm), coarse sand (0.5-2.0 mm), mean sand (0.25-0.5 mm), fine sand (0.063-0.25 mm) and silt & clay (<0.063 mm). All sediment fractions were expressed by the % of the total sediment weight (Brown and McLachlan, 2010).
Nematode samples were collected from the same Van Veen grab (0.05 m 2 ) as the sediment samples, by forcing a hand core (3.6 cm inner diameter) to a depth of 3 cm of sediment and further preserved in a 4% buffered formalin solution. In the laboratory, fixed samples were first rinsed through a 1000 µm mesh and then through a 38 µm mesh. The retained fraction was washed and centrifuged three times, using colloidal silica polymer LUDOX HS-40 (specific gravity 1.18 g cm − 3 ) (Heip et al., 1985). All the extracted nematodes were counted under Leica M205 C (100X magnification) and a set of 120 randomly picked individuals were collected from each sample for further identification (Heip et al., 1985;Vincx, 1996).
Total of 3729 nematodes were identified until the genus level under an Olympus BX50 light microscope, using identification keys Warwick, 1983, 1988;Warwick et al., 1998) and Nemys on-line database (Bezerra et al., 2019). All of the identified individuals were measured for length (L, excluding filiform tail portions) and maximum body width (W) under an Olympus BX-50 compound microscope (1000 × magnification) with Olympus Cell^D software. Individual biomass (µg dry weight [dwt]) for all specimens was determined using the Andrassy formula to estimate biomass as wet weight (Andrassy, 1956), which is based on nematodes L and W. The conversion to a dry weight was based on a dry/wet weight ratio of 0.25 (Heip et al., 1985). The length/width (L/W) ratio was calculated by dividing total individual body length (L) by its maximum body width (W) (Platt and Warwick, 1988;Vanaverbeke et al., 2004). Nematodes were further classified based on their body shape accordingly to three size classes: stout nematodes (L/W ratio < 18), slender nematodes (L/W ratio between 18 and 72) and long/ thin nematodes (L/W ratio > 72) (Schratzberger et al., 2007).

Data analysis 2.3.1. Environmental data
To provide a link between environmental variables and nematode assemblage patterns, we performed a correlation analysis between each pair of abiotic variables. Since salinity is the key variable in estuaries, the classification of the estuary in salinity sections (e.g., oligo-, meso-, etc.) was used as grouping factor for better visualization of the environmental variables in Fig. 2 and morphometric attributes in Fig. 3.   Spearman correlations between pairs of abiotic variables. Abbreviation for variables are as follows: tempbottom water temperature ( • C); TOMp -Total organic matter (%); silt&clay -silt and clay fractions (%); depthwater depth (m); sandsand fraction (%); fsandfine sand fraction (%); gravelgravel fraction (%); salsalinity; oxyPabove sediment dissolved oxygen (%).
To test if abundance dataset weighted by the morphometric attributes has higher explanatory power in detecting spatial patterns than using abundance of morphometric data alone we used modified Bray-Curtis dissimilarity matrices. Each modified Bray-Curtis dissimilarity matrix was based on the "Generalized Canberra dissimilarity", proposed by Ricotta and Podani (2017): where χ Uj and χ Vj are the abundance values of genus j at sampling stations U and V, G is the total number of genera. Genera specific contributions π j were based on a given morphometric attribute that contributed to compute dissimilarity between sampling stations. The use of the modified Bray-Curtis dissimilarity index resulted in five dissimilarity matrices: Abundance weighted by length (A + L); Abundance weighted by width (A + W); Abundance weighted by the L/W ratio (A + L/W); Abundance weighted by the biomass (A + B); Abundance weighted by the LxW (A + LxW). The relationships between each dissimilarity response matrix and untransformed environmental variables were tested using Distance-based redundancy analysis (dbRDA; function "capscale" in vegan R package, Oksanen et al., 2016). This analysis allows to carry out constrained ordinations based on noneuclidean resemblance measures (Legendre and Anderson 1999). Variance inflation factors (VIF) were calculated to check for multicollinearity and to ensure that only variables with small VIFs (<10) were included.

Morphometric vs abundance efficacy using 12 most frequent nematode genera
In order to investigate the efficacy of morphometric attributes versus abundance measures we repeated the above-described analysis, with the 12 most frequent nematode genera only. The criterion for the number of genera used for the analysis was decided based on the histogram of the genera frequency (See the results section for detailed description of the selection criteria). To additionally justify that the dataset composed of only 12 genera conserves the multivariate patterns as the original dataset based on the entire 95 genera composition, we performed Procrustean test (PROTEST Jackson et al., 1995). Procrustean test measures the degree of concordance between two or more datasets having different characteristics and if statistically significant, two datasets reflect in the same way the processes that determine their association (Peres-Neto and Jackson, 2001). Procrustes analysis requires that both datasets have the same number of columns and rows. For this reason, prior to the ordination analysis, we added 83 (95 genera − 12 genera) columns of zeros to the reduced assemblage raw dataset so that the number of rows and columns will equal the original dataset. Consequently, the Protest function was applied to all of the datasets for 12 genera versus datasets for complete 95 genera composition. All the analyses were performed in R using vegan statistical package (Core R Team, 2012).

Environmental variables
All environmental variables are presented in the Appendix 1. Fig. 2 represents environmental variables that contribute to the abiotic spatial characterization of the Sado estuary. The salinity gradient of the Sado estuary increases gradually from the uppermost area (oligohaline), towards middle channels (mesohaline and polyhaline) until the estuary mouth (euhaline).
Dissolved oxygen saturation also increases towards downstream following the salinity gradient (Figs. 2 and 3). The highest water depth was recorded near the estuary mouth, but also at the three stations (S7, S8, S9) in the middle of the estuary. Water depth was demonstrated to be negatively correlated with organic matter content and bottom water temperature (Fig. 3). Consequently, sampling stations located in the proximity of the estuary mouth had low organic matter content, but high dissolved oxygen concentrations. The shallowest area of the estuary was the Bay area and few stations in the middle of the estuary. These sampling stations were rich in organic matter content and had less dissolved oxygen. Stations located at the uppermost area of the estuary presented the lowest dissolved oxygen saturation values (Fig. 2, Appendix 1). Grain size distributions were rather patchy and did not present any significant gradient along the estuary. Nevertheless, granulometry was correlated with water depth, with fine grained sediment (silt and clay fractions) being negatively correlated with the water depth, while coarser grain size (sand and fine sand) was positively correlated with water depth (Fig. 3).

Morphometric attributes
Long nematodes almost exclusively inhabited euhaline and polyhaline sampling stations and accounted for only 4.92% of the total individuals. The most abundant group was slender nematodes accounting for 86.64% of the total abundance along the entire estuary. However, even within the range of this body type a clear separation was observed between bigger nematodes, with higher L/W ratio variability, with greater extent at the polyhaline and euhaline stations and smaller ones, less variable in L/W ratio, occurred at the mesohaline and oligohaline areas. Stout nematodes characterized by short and robust body were again predominant at higher salinity stations restricted to the middle estuary and estuary mouth, accounting for 8.44% of the total abundance (Fig. 4).

Morphometric vs abundance efficacy for entire nematode genera assemblages
The results from the multidimensional scaling and dbRDA analysis for datasets for the complete nematode genera assemblages are summarized in Table 1A (See the Appendix 2A for the full dbRDA results for all of the datasets). The ordination accounting for the highest amount of variance in the first two ordination axes (33%) was reported for the dataset based solely on genera abundance composition. A similar result was observed for the dataset based on L/W ratio (30% of the variance was accounted by the first two ordination axes). Datasets based on biomass as well as modified Bray-Curtis datasets weighted by the morphometric attributes captured less multivariate variation in the first two axes (Table 1A, Appendix 2A).
Significant environmental vectors from dbRDA analysis demonstrated that silt and clay fractions, depth and temperature were the best explanatory variables of nematode genera abundances, with the highest adjusted R 2 reported of 0.15 (Table 1A, Fig. 5A). The second dbRDA component was clearly discriminating between shallow areas with the dominance of fine sediments of silt and clay fractions (such as in the Bay area) and the deeper areas confined to estuary mouth, but also to the two deep sampling stations located in the upper estuary (Fig. 5A).
For L/W ratios dataset salinity was a major explanatory variable together with % of organic matter and temperature. Whereas depth, temperature and pH were major contributors to explain the biomassbased assemblage distributional patterns (Table 1A).
The % of dissolved oxygen was indicated as the most influential variable for the combined datasets (A + L/W and A + B), where abundance was weighted by morphometric attributes (Table 1A, Fig. 5B, C, D, E, F). It is worth noting that oxygen was not indicated as a significant factor for any other dataset analysed. Although oxygen was positively correlated with salinity, this last variable was not significantly associated with this distribution. Therefore, for combined datasets, the oxygen content determined the sampling station's distribution patterns where high oxygen content corresponded to stations in the mid and low estuary while stations with low oxygen concentrations were located mainly in the upper estuary (Fig. 5B, C, D, E, F). Ordinations based on abundance-based dataset discriminated shallow and rich in fine sediment areas from deeper stations, with prevalence of coarse substrate located next to the estuary mouth (Fig. 5A). For the combination of A + B, depth was not significant and only oxygen % and sand were significant variables in explaining the distribution of the assemblages (Fig. 5E).  Schratzberger et al. (2007). Each symbol represents different salinity section.

Morphometric vs abundance efficacy for the most frequent nematode genera
The criterion for the number of genera used in the analysis was decided based on the histogram of the genera frequency. Genera observed at more than 19 sampling stations (12 genera out of total 95) were used for the analysis. This cut-off line was clearly represented in the histogram frequency of the genera distribution (Appendix 3). While genera with more restricted distribution were present at one to five sampling stations (66 genera) and five to 15 sampling stations (17 genera), more widespread genera were present in 19 sampling stations with some genera occurring at 33 sampling stations.
The results from the multivariate analysis for the 12 most frequent genera demonstrated a higher explanatory power of the first two components, particularly regarding morphometric attributes, in comparison to the dataset based on the entire nematode assemblage (Table 1B). First two principal components explained 37% of the variance in the abundance dataset while the first two ordination axes accounted for 35% of the variability when using the L/W matrix (Table 1B). Similarly, other datasets based on morphometric attributes allowed to increase the amount of variability captured by the first two ordination axes by 5% in comparison to datasets based on the full nematode assemblages (Table 1B, Appendix 2B). An interesting aspect became apparent when reduced genera datasets were subjected to dbRDA. While nothing changed for the abundance dataset, morphometry-based datasets  Fig. 6. dbRDA for the most frequent genera for abundance dataset and morphometric attributes datasets. Only the significant variables (after forward selection) are displayed as vectors. Abbreviation for variables are as follows: sil&cla -silt and clay fractions (%); gravelgravel fractions (%); sandsand fractions (%); depthwater depth (m); tempbottom water temperature ( • C); salsalinity; TOMp -Total organic matter (%); TOMg -Total organic matter (g); oxyPabove sediment dissolved oxygen (%).
showed a significant increase in the variation explained by the fitting of environmental variables, indicated by the high adjusted R 2 in comparison to the full genera dataset. Dataset based on L/W ratios demonstrated that mainly four environmental variables played an important role in explaining the distribution of the assemblages (Fig. 6B). The first axis differentiated the sampling stations rich in organic matter content from those located next to the estuary mouth, poor in organic matter, but rich in oxygen and respective high pH. By turn, the stations with prevalence of sand were distributed mainly along the second ordination axis, associated also to the upper estuary and discriminating between deeper stations with gravel prevalence. Similar patterns were observed for the "length" and "width" datasets ( Fig. 6 C, D). In the case of the "length × width" dataset, salinity turned to be an important environmental explanatory variable (Fig. 6F). For the "biomass" dataset, sand and pH were again the main explanatory variables with a high number of stations being tightly distributed along the pH vector (Fig. 6E). For combined datasets, the same patterns as for the entire assemblage datasets were observed, with dissolved oxygen % and depth as the most important explanatory variables. Nevertheless, both explanatory power and adjusted R 2 were the highest for the reduced datasets, only based on the 12 most frequent genera. The Procrustean test showed a significant similarity between all of the ordinations (Appendix 4), satisfactorily indicating congruence in multivariate ordination among analyses based on the full assemblage dataset and on the most frequent genera dataset.

Ability of nematode morphometric attributes in detecting spatial patterns in relation to the taxonomic approach
Ecological data are often complex and multi-dimensional and thus require multifaceted approaches to explore the underlying mechanisms that drive community assembly (Dray et al., 2012). Differences between communities can be either reflected through morphological differences, such as different species composition or through differences in their functional attributes, or both (Webb et al., 2002).
Large majority of ecological studies deal with species-based datasets separately from trait-based datasets (Bremner et al., 2006). And while the use of functional traits has demonstrated to be highly complementary in relation to taxonomy-based approaches (Alves et al., 2014;Schratzberger et al., 2007), its combined usage to explain spatial patterns has never been addressed for meiofauna. Our initial hypothesis was that genera-based abundance dataset weighted by the morphometric traits would better reflect the assemblage distribution patterns along an estuarine gradient, when compared to the traditional taxonomic approach based solely on genera abundances. This hypothesis was based on the prediction that functional traits are more sensitive to different sets of environmental variables (Schratzberger et al., 2007) than taxonomic attributes and thereby will have ability to detect additional spatial gradients contributing to an overall higher explanatory power of the multivariate ordination. Our results demonstrated that the dataset solely based on genera abundances had consistently a better explanatory power than combined datasets or morphometric datasets alone. Main gradients described by the taxonomy-based dataset were related to the sediment particle size and water depth, where shallow areas with mud predominance had different nematode communities than the deeper areas of coarser grain size. Both salinity and granulometry play the major role in structuring nematode communities (Adão et al., 2009;Wu et al., 2019). Salinity did not emerge as an important grouping factor for this dataset and communities were demonstrated to be more responsive to the overruling differences in granulometry. For example, the Bay section (stations: S31, S32, S33, S34) was dominated by the deposit-feeding generalist type of taxa (Terschellingia, Sabatieria) that are known for their preferences for fine sediments, rich in organic matter, where they can outcompete other genera (Austen et al., 1989;Soetaert et al., 1995). However, their distribution patterns do not seem to be related with the salinity, as they can be found within high salinity ranges (Soetaert et al., 1995), contrary to other genera that usually live within their optimum salinity regime (Hourston et al., 2009). Consequently, high affinity of these two most abundant genera to silt and fine sediments might have obscured other important gradients that influence nematodes distributional patterns. Exceptionally high abundance of these two genera might be also related to seasonal differences in organic matter and microphytobenthos availability (Hourston et al., 2009). Sampling surveys were performed in spring corresponding to a peak of food availability (Hourston et al., 2005(Hourston et al., , 2009, especially in fine sediments, thereby promoting higher abundances of generalist taxa. A temperature related gradient was also detected, though its ecological importance remains more difficult to explain.
Spatial discrimination, based on combined datasets, allowed delineating ecologically distinct gradients in the estuary, than those captured by the abundance dataset. The variable that consistently emerged as significant for this spatial discrimination was dissolved oxygen % saturation, where sections with higher oxygen saturation corresponded to those with a higher marine influence in opposition to the stations located in the uppermost areas of the estuary, with freshwater predominance.
Consequently dissolved oxygen % saturation seemed to be an efficient proxy of estuarine conditions important for nematode distribution patterns. In a similar study using macrofauna communities, oxygen % saturation was highly correlated with pore water salinity and the presence of seagrass (Beard et al., 2019), thereby indicating an important hotspot for the macrofauna diversity.
Oxygen levels are diminishing worldwide both in open ocean and coastal areas (Breitburg et al., 2018) having direct negative impacts on all different aspects of aquatic life (Diaz and Rosenberg, 2008). However, unlike macrofauna, which is absent in anoxic sediments, meiofauna can have much wider adaptations to reduced oxygen concentrations (Norkko et al., 2019;Wetzel et al., 2001), being a potential important indicator of low oxygen associated to different type of disturbances (e.g., organic pollution, wastewater discharge etc.). Thereby, functional response of nematodes to oxygen saturation levels provide important information to advance the development of a nematode index that would be indicative of low oxygen levels (and thus potential organic pollution). The use of morphometric attributes as indicators of oxygen requirements are actually not new as Wieser and Kanwisher (1961) observed that the size of buccal cavity of different nematode genera can be a proxy for the amount of oxygen consumption to support their metabolic functions and thereby withstand the surrounding environment (Wieser and Kanwisher, 1961). This topic deserves to be further developed, by investigating which genera could be used as indicators of "low oxygen conditions" and what are their characteristics.
While the oxygen concentration along the sediment vertical profile is crucial for nematode distribution patterns and their morphometric adaptations , not much is known concerning the relationships between the % dissolved oxygen saturation above the sediment and oxygen availability inside the sediments. While it is well beyond the scope of this paper to interpret biological significance of the oxygen % saturation for nematode assemblages distribution patterns, it is important to draw attention to the role of this variable in the definition of estuarine habitats, and the fact that it only emerges as significant when functional community components are considered. Oxygen turned to be a significant factor also for other nematode functional traits (lifehistory strategy and functional feeding groups composition) as well as for taxonomic and functional diversity indices distribution patterns (Alves et al., 2014;Sroczyńska et al., 2020) and undoubtedly plays a significant role in spatial determination of habitat conditions in Sado estuary.
Taxonomy-based dataset captures not only phylogenetic differences, but also ecological signals resulting from the niche-driven processes and dispersal (Padial et al., 2014). Our data shows that the most explanatory dataset was the one based on the genera abundances, followed by the dataset of genera L/W ratios. Morphometric traits represent a physical adaptation of the organisms to live in a certain type of conditions and L/ W ratios of assemblages were previously demonstrated to significantly differ between adjacent tropical habitat types (Armenteros and Ruiz-Abierno, 2015). The combined dataset had less explanatory power than above mentioned abundance and L/W ratio datasets, nevertheless, when displayed on dbRDA plots, combined datasets clearly discriminated ecologically different estuarine conditions. The variance captured by ecological models may be low due to a variety of ecological processes that operate at spatial and time scales (Conde et al., 2018). The latter may happen when using multivariate matrices based on counts or alternative measurable response variables such as the organism's biomass, length or width. In this study, the lower fit of morphometric traits in the model may rely on at least two different reasons. Firstly, the values in morphometric-based dataset were mean of measurements of different individuals, which increased the overall variability of morphometric values in comparison to the absolute count values of abundance-based dataset. This lower variability of counts data in the abundance-based dataset caused higher signal to noise ratio, making its response to environmental variables more predictable than when using morphometric attributes. Secondly, while the value of count data is undisputed, the measurement of morphometric traits has an associated error depending on both accuracy and precision. Thus, these two sources of error imply a lower signal to noise ratio when using morphometric traits in multivariate models, likely lowering the explained variance (including combined datasets). Despite the lower explanatory variance, combined datasets have potential to indicate different ecological gradients of the nematodes distribution patterns within an estuary in comparison to the gradients captured by each dataset separately. Therefore, we argue that combined datasets constitute a robust alternative for nematode assemblage data analysis, in comparison to solely abundance datasets.

Distributional patterns of the whole vs reduced nematode genera dataset
Common and rare species possess different sets of functional traits and are known to respond differently to environmental variables (Cornwell and Ackerly, 2010;Magurran and Henderson, 2003). Therefore, it would be expected that they show different ordination patterns in relation to environmental variables. However, we observed that patterns displayed by the reduced versus complete dataset were the same for abundance and abundance weighted by the morphometric attribute's datasets. This finding demonstrates that rare genera do not display higher levels of specialization than common genera. On the other hand, for the dataset based solely on morphometric attributes, new environmental variables contributed to a greater spatial discrimination of the estuary habitats. For example, the dataset based solely on L/W ratio displayed clear differentiation along both dbRDA axes distinguishing four estuarine areas (1. Of sand prevalence; 2. Of high pH; 3. Dominated by organic matter deposits and 4. Deep areas with gravel prevalence). One could argue that the analysis based on functional traits instead of taxonomic entities would only retain sufficient information when analysed across strong environmental gradients or disturbance sites (De Castro et al., 2018). Other authors also demonstrated that the removal of rare species negatively affected the functional diversity of fish, birds and trees (Leitão et al., 2016). Considering that our study area did not present any sharp gradient in environmental conditions, we argued that the morphometric traits of the most common genera retained sufficient ecological information on the estuary spatial patterns being ecologically more informative than dataset based on the complete set of genera. Furthermore, the repeated analysis for the most frequent genera resulted in higher explanatory power and concomitant higher adjusted R 2 . Rare species might influence the analysis by introducing a higher noise relatively to a signal ratio producing confounding results (Cunningham and Lindenmayer, 2005;Gauch and Gauch, 1982). On the other hand, removing rare species from the analysis can be responsible for underestimation of the differences between assemblages (Cao et al., 1998) or the reduction of the significant correlations between the ordination space and environmental variables (Faith and Norris, 1989). In the present study, the removal of 83 genera from the analysis improved the ability for detection of spatial patterns for all the datasets, but particularly for datasets solely based on morphometric attributes. This improvement in the signal ratio was achieved without compromising changes in the original (whole assemblage) ordination patterns, as indicated by the results of the procrustean test. Similar results, but based on the presence-absence and abundances of aquatic macroinvertebrates, were found by Draper et al. (2019) and Sgarbi et al. (2020) after removing a substantial amount (over 80%) of rare species.
Such information is of particular importance for the bio-assessment surveys as it may substantially increase the efficiency of such programs by reducing the cost and work associated with identification and measurements of all of the individual nematode genera. It demonstrates that selective adaptation to specific estuarine conditions expressed by different morphotypes of 12 most abundant genera yield sufficient ecological information for the estuarine spatial gradient delineation. In the context of the bio-assessment, our finding also implies that at poor biodiversity sites, morphometric attributes can constitute an important and promising surrogate for delineating gradients associated to different ecological conditions. Nevertheless, we tested this relation only for the morphometric attributes (quantitative type of traits) and for only one phylum (nematodes) and thus it might not be applicable to other taxonomic groups/trait types or ecosystems.

Remarks on the calculation of functional diversity: What can we learn from nematodes
The functional trait approach has been increasingly used to answer a variety of ecological questions, particularly how changes in environmental conditions might impact biodiversity and ecosystem functioning (Gagic et al., 2015). Functional trait approach relies on the computation of functional diversity (FD) derived from different trait types (Villéger et al., 2008). Most of the FD metrics are calculated based on genera or species-specific trait measure and thus represent the weighted average of a measured trait type for all of the individuals belonging to the same taxonomic entity. Therefore, the assumption behind FD computation is that the interspecific trait variability is larger than intraspecific variability (Auger and Shipley, 2013). Our results emphasize that morphometric attributes of nematodes display changes in relation to the local environmental variables, implying that the same genus might differ in their morphometric features depending on the environment (estuary areas) it lives in. Therefore, morphometric attributes in this case are not necessarily genera specific, but also highly dependent on the surrounding environment. Furthermore, if these morphometric attributes are included into computation of functional diversity metrics, where each genus is assigned to its average trait value, it will obscure important differences derived by the individual local-specific morphometric characteristics and will not assists in the discrimination of the estuary areas (Pakeman and Quested, 2007;Pakeman, 2014). Therefore, our study based on nematode assemblages highlights the importance of considering intraspecific trait variability to obtain a better insight into a species response to the environmental conditions.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.