Combining Traditional Taxonomy and Metabarcoding: Assemblage Structure of Nematodes in the Shelf Sediments of the Eastern Antarctic Peninsula

Pantó, Gabriella; Pasotti, Francesca; Macheriotou, Lara; Vanreusel, Ann

doi:10.3389/fmars.2021.629706

ORIGINAL RESEARCH article

Front. Mar. Sci., 23 July 2021

Sec. Marine Evolutionary Biology, Biogeography and Species Diversity

Volume 8 - 2021 | https://doi.org/10.3389/fmars.2021.629706

This article is part of the Research Topic Effects of Ice Loss on Marine Biodiversity View all 11 articles

Combining Traditional Taxonomy and Metabarcoding: Assemblage Structure of Nematodes in the Shelf Sediments of the Eastern Antarctic Peninsula

$\r\nGabriella Pant*$ Gabriella Pantó^*

Francesca Pasotti

Lara Macheriotou $Ann Vanreusel\r\n$ Ann Vanreusel

Biology Department, Marine Biology Research Group, Ghent University, Ghent, Belgium

This study provides a snapshot of the largely understudied meiobenthic and nematode communities in the Prince Gustav Channel (PGC) and Duse Bay (DB). We compared five stations sampled at different water depths along the shelf and investigated their meiobenthic community structure. We approached nematode biodiversity combining traditional taxonomic identification and high throughput sequencing (HTS), with the use of Amplicon Sequence Variants (ASVs). Additionally, we characterized the environment by primary production proxies, grain size and seasonal ice conditions. Our results suggest that the availability of organic matter and its freshness are responsible for the high densities found at all depths. However, potential factors influencing the high local and regional variability of meiofauna density and biodiversity are less clear. A bathymetric transect consisting of three stations in DB (200, 500, and 1,000 m depth) showed increasing pigment concentrations in the first centimeters of the sediment vertical profile with increasing water depth, whereas the meiofauna densities showed the opposite trend. The deepest station of DB seems to function as a sink for fine material as supported by the higher silt fraction and higher organic matter concentrations. When comparing the two basins in the PGC (1,000 and 1,250 m) and the one in DB (1,000 m), differences in terms of environmental variables, meiofaunal densities, and composition were observed. The deepest basin in PGC is located further South (closer to the highly unstable Larsen area), and marked differences with the other basins suggest that it might be experiencing different conditions as a result of its presence near the summer ice margin and its more elongated topography. Both, the shallowest and the deepest stations showed the highest number of unique sequences, suggesting a more biodiverse nematode assemblage. The morphological identification did not show significant differences in the biodiversity of all stations, differently from the ASVs approach. However, the lack of reference sequences in online databases and the thickness of nematode’s cuticule are still important issues to consider as they potentially lead to underestimations of biodiversity and functional traits.

Introduction

Ecosystem dynamics of the sea-bottom are strongly influenced by the density and diversity of organisms inhabiting marine sediments. Macro- and meiofauna organisms have proven to be especially important in investigating biodiversity, as well as characterizing ecosystem dynamics both in shallow and deep waters (Schratzberger and Ingels, 2018). Benthic communities have a positive effect on ecosystem functions by modifying the environment they inhabit, for instance facilitating the biodegradation of organic matter by microbial activity stimulation and influencing the vertical redox potential of sediments by bioturbation (Zeppilli et al., 2018; Hoffmann et al., 2019; Gogina et al., 2020). Even though meiofauna (>32 μm, <1 mm) has been historically understudied because of the tediousness of sampling and identifying such small size organisms (Schratzberger and Ingels, 2018), it plays a key role in aquatic sediments. By feeding on both microalgae and bacteria and being consumed by higher trophic levels, meiofauna taxa can be both, consumers, and producers (Vinagre et al., 2012; Carpentier et al., 2013; Schmid-Araya et al., 2016; Schratzberger and Ingels, 2018). The reduced size of meiofauna makes them extremely sensitive to even the slightest changes in sediment grain size and structure (Taheri et al., 2017; Korbel et al., 2019), but also in oxygen availability and organic load within the sediment (Hoffmann et al., 2018). Because of their short generation time and direct benthic development, free-living nematodes in particular may be able to reveal environmental perturbations on a small spatial and temporal scale, as confirmed in several studies (Kennedy and Jacoby, 1999; Moreno et al., 2011; Semprucci et al., 2015). As they represent >90% of the metazoans inhabiting marine sediments, nematodes are also considered good proxies for assessing benthic biodiversity, especially in relation to meiobenthic influence on ecosystem dynamics. In light of possible climatic changes, the analysis of the biodiversity in extremely sensitive habitats might shed light on potential impacts of ongoing and future climate stressors.

Despite polar regions being particularly sensitive to such climatic alterations, Antarctic marine-shelf habitats are still amongst the least studied on Earth. Their accessibility is limited due to logistic difficulties in light of their remoteness and the regional climatic conditions of the Antarctic (Cook et al., 2005) yet the high diversity and abundance of their biological assemblages and related ecosystem services that they provide are of main interest for a number of scientific questions. The Antarctic Peninsula (AP) is known for being one of the most impacted regions by anthropogenic global warming, which results in catastrophic events such as ice-shelf collapse, glacier thinning (Martín et al., 2009; Kunz et al., 2012), acceleration (Pritchard and Vaughan, 2007) and recession (Cook et al., 2005; Davies et al., 2012). The AP has experienced rates of atmospheric warming greater than the global mean (Vaughan et al., 2003) and a loss of 28 000 km² of ice shelves since the 1960s (Cook and Vaughan, 2010). The Eastern Antarctic Peninsula (EAP) is particularly affected by variations in westerly circumpolar winds which are responsible for basal melting leading to the thinning and eventually the collapse of ice shelves (Oppenheimer, 1998; Rignot and Jacobs, 2002; Polvani et al., 2011; Pritchard et al., 2012). From a global warming perspective, it is appropriate to assume that further increase of average global temperatures will lead to dramatic and irreversible changes in the AP (Etourneau et al., 2019; Siegert et al., 2019), especially the already unstable EAP ice shelf.

The ecosystem of the AP is characterized by high primary productivity, with large phytoplanktonic blooms that have a bottom-up effect in sustaining the regional food web, strictly linked to the seasonality of sea-ice dynamics (Prezèlin et al., 2000; Nardelli et al., 2018). When the sea ice has melted and the ocean surface is exposed to solar radiation, conditions are favorable for phytoplankton to grow, and productivity reaches its maximum (Moline and Prezelin, 1996). Previously ice-covered benthic areas seasonally experience a shift to higher phytoplanktonic abundance due to the temporary system transition to more open waters (Montes-Hugo et al., 2009). Because primary productivity is of crucial importance for shaping the communities of benthic organisms living in shallow and deep sediments, ice cover and changes in ice dynamics (determining a shift to open water systems and possibly changing the amount of OM depositing on the seafloor) are essential players when studying the meiofauna community of polar ecosystems.

While morphological identification is still a common practice used in meiofauna diversity assessments, the development of molecular techniques as high throughput sequencing (HTS) has opened the possibility to analyze a larger set of samples in a much shorter timeframe, often resulting in faster and more reliable analysis (Brannock et al., 2018; Steyaert et al., 2020). Molecular approaches can overcome some of the limitations of microscope-based identification: for example, the detection of cryptic species (i.e., individuals with same morphological features but different DNA sequences) is not possible based solely on morphological identification. Traditional meiofauna identification is in fact time-consuming and requires a specific expertise, thus molecular techniques are being increasingly used as alternative diagnostic approach (Dell’ Anno et al., 2015; Fonseca et al., 2017; Macheriotou et al., 2019). Disadvantages of the molecular analysis for species identification mainly lie in the incompleteness of online databases used to assign the taxonomical label to a specific DNA sequence and the need to find the right marker for a specific community (Ahmed et al., 2015). Nevertheless, many studies of the last decade have proven this approach to be particularly useful when dealing with such abundant organisms as nematodes, particularly by means of DNA metabarcoding and the use of Amplicon Sequence Variants (ASVs) (Macheriotou et al., 2020; Steyaert et al., 2020).

By analyzing the meiofauna and specifically nematodes along a depth gradient in Duse Bay (DB) and by comparing two adjacent deep basins in the Prince Gustav Channel (PGC), the following hypothesis were tested: (i) the community composition and density of meiofauna and nematodes along a bathymetric transect in DB is dependent on changes with depth and the related environmental variables; (ii) meiobenthic community structure in the three deep basins is dependent on local small-scale variability drivers and seasonal sea-ice cover; (iii) the use of ASVs efficiently provides accurate information on the connectivity of nematode communities along the shelf and thus possible gene flow between stations.

Materials and Methods

Study Sites and Sampling Methodology

The focus of this study are the meiobenthic communities of the PGC and DB, both situated in the EAP. The PGC is a 6–7 km narrow seaway enclosed between the Trinity Peninsula and the James Ross and Vega Islands, northeast of the AP. The bathymetry of the channel is characterized by a “U-shaped” 1,000–1,250 m deep trough with steep walls interrupted by shallow ledges, which separate the deep basins (Camerlenghi et al., 2002). As part of the AP, the PGC has been characterized by the collapse of its ice-shelf since the late 1980s. The glacier dynamics in the region are particularly sensitive to ice-shelf back pressure (which slows down the ice-shelf movement toward the ocean); in fact, the rapid retreat of glaciers after a collapse event seems to be common in the PGC and Larsen A and B districts (Rott et al., 2002). This is confirmed by the disintegration of the PGC ice shelf and Larsen A embayment in 1995 and the northern and central sections of the Larsen B embayment in March 2002, which resulted in considerable ice losses due to intensified ice discharge (Shuman et al., 2011). Mass reduction has been recorded until recent times, and frontal retreat and intensification of calving fluxes is expected to occur in correlation to glaciers which showed significant temporal variations of ice flow during recent years (Rott et al., 2018). DB is a sheltered inlet situated on the south side of the Trinity Peninsula, between Tabarin Peninsula and View Point (Figure 1). It is characterized by extremely rapid glacier loss compared to the rest of the region, together with Larsen Inlet, Sjögren Inlet, and Eyrie Bay (Davies et al., 2012). The topography of the area is characterized by a steep slope. Understanding that limited information is available on the benthic ecology of DB, this study represents an initiative for further studies in the area. The designated sites for this study were two deep basins along the PGC at 1,250 and 1,000 m, and a bathymetric transect in DB (200, 500, and 1,000 m depth; Table 1). The study area and sampling stations are shown in Figure 1.

FIGURE 1

Figure 1. Sampling stations in the Prince Gustav Channel (PGC, 1,000 and 1,250 m) and Duse Bay (DB, 1,000, 500, and 200 m). Bathymetric map courtesy of: “British Antarctic Survey” and “Bathymetry Group of the Alfred-Wegener-Institute.”

TABLE 1

Table 1. Overview of stations sampled during the RRS James Clark Ross cruise JR17003a (February–March 2018).

Sampling was carried out during the RRS James Clark Ross cruise JR17003a in February-March 2018¹. Sediment samples were collected using a 12-core Multicorer (MUC) type Oktopus, each core having an inner diameter of 10 cm. At each station three to four cores were collected from the same or different MUC deployments, as shown in Table 1. Limitations on the amount of deployments were dictated during the cruise by the ice conditions and time constraints. Overlying water from each intact core was filtered over a 32-μm sieve and preserved with the corresponding sediment sample. For each core, the sediment was divided in equal horizontal slices of 2 cm, down to 8 cm depth. In turn, each slice was equally divided in two subsamples and preserved in different conditions according to the prospected analysis. Specifically, the first half was stored in a 250 mL container and fixed with 4% formaldehyde (buffered with pre-filtered seawater) for morphological meiofauna community analysis (density, biodiversity). The second half was placed in plastic bags and frozen at −80°C for molecular analysis of meiofauna (i.e., metabarcoding) and for the examination of environmental variables (pigments, organic matter, sediment grain size).

Environmental Parameters: Pigments and Granulometry Data

From the frozen sliced cores, 250 mg of sediment was subsampled from each layer for phytopigment analysis. Environmental parameters included phaeopigments, organic matter and sediment grain size (Table 2). Pigments were extracted with acetone (90%, 8 mL), ultrasonicated for 30 s and further extracted overnight at 4°C. Subsequently, samples were centrifuged for 10 min at 4,000 rpm (4°C) and the supernatant was filtered through 0.2 μm Polytetrafluoroethylene (PTFE) membrane filters (VWR International). As described by Van Heukelem and Thomas (2001), pigments were then separated on a HPLC equipped with a cooled auto-sampler, column heater, a photodiode array detector and a fluorescence detector as described. Authentic pigment standards together with compounds isolated from reference cultures were used to calibrate the system following the Scientific Committee on Oceanic Research (SCOR) protocols (DHI, Denmark). Finally, the concentration of individual pigments was calculated using the response factor of standard pigments (Van Heukelem and Thomas, 2001) and it is reported as μg/g. Sample preparation for total organic carbon (TOC; expressed as% C_org w/w – weight of carbon to total weight of sediment) and total nitrogen (TN) (expressed as % N) was performed by first drying and homogenizing the sediment. Subsequently, 37% HCl (purity grade) was added to TOC samples for decalcification. TOC and TN contents were measured using Flash 2000 NC Sediment Analyser from Interscience. Grain size was estimated using the Mastersizer Hydro 2000G (MALVERN) particle size analyser: the measurable particle range size is between 0.01 and 2,000 μm in volume percent. The grain size classes used were the following: clay (<4 μm), silt (4–63 μm), and sand (63–500 μm). Grain size distributions were evaluated with the Mastersizer 2000 5.4 program.

TABLE 2

Table 2. Concentration of all the environmental variables analyzed.

Meiofauna Extraction and Characterization

The volume of each sample was recorded after sedimentation in a 250 mL cylinder before extracting the meiofauna fraction. Samples fixed in 4% formaldehyde were then sieved over a 1,000-μm sieve on top of a 32 μm sieve. Organisms larger than 1,000-μm (macrofauna) were dyed with Rose Bengal (0.5 g L^–1) and stored in 4% formaldehyde for further analysis. The meiofauna fraction (between 1,000 and 32 μm) was extracted by density-gradient centrifugation (3 × 12 min, at 3,000 rpm) with the colloidal silica polymer LUDOX HS40, with a specific density of 1.18 g/cm³ (protocol described in Heip et al., 1982; Vanhove et al., 1999). The resulting extracted sample was finally dyed with Rose Bengal (0.5 g L^–1) and stored in 4% formaldehyde. The meiofauna was counted at the stereomicroscope. The counting of the most abundant taxon (Nematoda) was performed with the use of an eight-chambered meiofauna sample-splitter (Jensen, 1982). Three out of the eight chambers were randomly chosen as sub-replicates to count nematodes, copepods and copepoda nauplii. Organisms were initially identified to higher taxon level, and results were standardized to individuals per 10 cm^–2. We conducted a preliminary inspection on the entire sediment vertical profile at each station to establish which layer hosted the majority of the organisms. The first 2 cm of the sediment are expected to harbor the majority of meiofaunal taxa, which has been reported in previous studies (Ingels et al., 2012) and confirmed by the results of the vertical abundance analysis in this study (Supplementary Annex A). Quantification of meiofauna was performed on three replicates (A–C) per station: the abundance analysis of the first replicate (A) was done along the complete depth gradient whereas for the other two replicates only the first 0–2 cm layer was considered. Given the high patchiness and variability of nematode communities within the sediments, it might be necessary to assume a high standard deviation also in the deeper layers of the sediment. For these reasons, the results of the vertical gradient analysis (Supplementary Annex A) must be interpreted with care and such considerations taken into account.

Taxonomical Identification of Nematoda

At each station, 80 nematodes specimens from three replicates were picked out from 0 to 2 cm layer, transferred to anhydrous glycerol (Seinhorst, 1959) and mounted on glass slides for identification to genus level. The relatively low number of nematodes identified per replicate (in other studies the number ranges between 100 and 200 specimens), and the choice to analyze only the first 0–2 cm layer is related to the time-consuming identification process and the weight of the morphological analysis in this study. Further, the most abundant genera and their species would be represented in correct proportions in the small subsample used for morphological identification, and their relative importance within the sample would be still the highest. This would constitute meaningful support information to the molecular dataset (ASVs). The molecular analysis allows the identification of a higher number of individuals in a shorter timeframe and (hypothetically) with more accuracy, hence metabarcoding analysis (ASVs) was prioritized for biodiversity assessment over the use of the morphological dataset, which was instead used to integrate any missing information. The Taxonomical classification was carried out following the online key for free living marine nematodes (NeMysKey©, Bezerra et al., 2021) and the pictorial key from Warwick et al. (1998), based on relevant morphological traits (buccal cavity, tail, amphid, and reproductive organs). Nematode were then divided into trophic guilds following the subdivision suggested by Wieser (1953), considering the structure of their buccal cavities. Specifically, the guilds were: selective deposit feeders (1A) with small buccal cavities mainly feeding on bacteria; non-selective-deposit feeders (1B), characterized by a larger unarmed buccal cavity expected to feed on detritus and particulate organic matter; epistrate feeders (2A), with one or more small teeth to perforate the outer shell of unicellular organisms such as microalgae; predators/omnivores (2B), usually larger nematodes with cuticularized teeth (and possibly strong mandible structures) that predate on other meiofauna organisms and/or ingest detritus particles.

Metabarcoding Analysis

DNA Extraction

The molecular analysis of meiofauna was performed on the entire vertical profile (0–8 cm) of all frozen samples (three to four replicates per station). Samples frozen at −80°C were initially defrosted and meiofauna was extracted following the density-gradient centrifugation procedure using LUDOX HS40. The sieved samples were stored in CTAB in 1.5 mL Eppendorf tubes and frozen at −20°C. Genomic DNA extraction was performed by adding 6 μL proteinase-K (10 mg/mL) to the defrosted samples and centrifuging for 5 min at 14,000 rpm at room temperature. The formed pellet was grounded and bead-beaten for 2 min at 30 cycles/s in order to break nematodes’ membranes and tubes were then incubated for 1 h at 60°C before adding ammonium acetate (250 μL, 7.5 M). The samples were then centrifuged for 10 min at 14,000 rpm and 750 μL of the supernatant were transferred into a sterile 2.0 mL Eppendorf tube; subsequently 750 μL of cold 80% isopropanol was added. Tubes were incubated for 30 min at room temperature and centrifuged for 10 min at 14,000 rpm at 4°C. After the removal of the supernatant, samples were incubated on ice for 30 min with 1 mL of washing buffer (solution of 10 mM ammonium acetate in 76% ethanol) to remove salts. Tubes were centrifuged for 5 min at 14,000 rpm at 4°C, the supernatant was removed and the pellet was stored in 20 μL of sterile milliQ water. The samples were then stored frozen at −20°C.

Library Preparation and Sequencing

Using the “16S Metagenomic Sequencing Library Preparation” protocol (Amplicon et al., 2013) SSU_F_04-SSU/22_R primers (GCTTGTCTCAAAGATTAAGCC and TCCAAGGAAGGCAGCAGGC, respectively) were assembled with Illumina overhang adapters. The above mentioned primers were used to amplify the 18S (V1–V2 region) ribosomal locus of each sample in triplicate, using the PCR conditions shown in Annex B (Supplementary Material). The PCR mix comprised 8.4 μL of PCR-grade H2O, 4 μL Phusion Buffer, 4 μL Dye, 0.4 μL dNTP (10 mM), 1 μL forward and reverse primer (10 μM), 0.2 μL Phusion Hot Start II High Fidelity Polymerase (New England BioLabs, United States) and 1 μL DNA template, for a total sample volume of 20 μL. All DNA templates were initially diluted 1/10,000: failed amplifications were readjusted with lower ratios (1/1,000, 1/100, and 1/10 dilutions depending on the sample). To confirm amplification of the DNA, PCR products were run on a 1% agarose electrophoresis gel and triplicates pooled. Samples were then purified using Agencourt AMPure XP beads and run on Bioanalyzer 2100 to confirm length and size distribution of the PCR fragments. The FC131-1002 NexteraXT Index Kit (Illumina, United States) and Kapa High Fidelity PCR kit (Kapa Biosystems, United States) were used to implement library indexing. The mix consisted of 11.25 μL PCR-grade water, 5 μL Buffer, 0.75 μL dNTPs (10 mM), 2.5 μL Index1 and Index2, 0.5 μL Kapa Hot Start High Fidelity Polymerase and 2.5 μL PCR product. Index PCR conditions are shown in Annex B (Supplementary Material). Indexed PCR products were then purified a second time using Agencourt AMPure XP beads; they were run on a Bioanalyzer 2100 to confirm successful indexing and DNA concentration of each sample was measured using Qubit^® dsDNA High Sensitivity Assay Kit. Finally, biological replicates were pooled and sequenced at Edinburgh Genomics on an Illumina MiSeq-v3 2 × 300 bp paired-end read run².

Sequence Data Analysis

The removal of gene-specific adapter sequences was performed using Cutadapt (v1.12) where sequences were truncated from the 5′ and 3′ read ends. By truncating the forward and reverse reads (respectively at 225 and 250 bp), ASVs were generated, following the Dada2 tutorial³. Subsequently, taxonomic assignment of ASVs was performed in Qiime1 (assign_taxonomy.py) using the Naive Bayesian RDP28 classifier, with an estimate confidence of 0.80 (Macheriotou et al., 2019). The process was completed in two steps: the first step included the employment of a large eukaryotic reference training set (Silva release 123 for Qiime1, 99% OTUs and UGent nematode Sanger sequences; n = 20,201); the second step consisted in the extraction and isolation of all the ASVs which were identified as “Nematoda,” in order to have a comprehensive list of nematode sequences separated from the rest of the meiofauna species. The taxonomic assignment (until genus level) was then completed using another, smaller training set (2,178 sequences) containing nematodes exclusively. Finally, the samples were rarefied to the lowest count (n = 25,258) (Qiime1, alpha_rarefaction.py), for further statistical analyses.

Statistical Analysis

Environmental Variables and Morphological Dataset

The meiofauna community dataset and the environmental parameters were initially subjected to univariate data exploration to assess the dispersion of the data. Afterward, a multivariate analysis on the 0–2 cm layer of the meiofauna community data and all environmental variables was performed. This was done with a non-parametric permutational analysis of variance (PERMANOVA) in both cases. Subsequently, a PermDisp analysis (permutational analysis of multivariate dispersion) was used to test for differences in meiofauna densities between groups (Stations). A pairwise comparison was then performed on significant results. For the analysis of morphologically identified meiofauna, a two-factor design was used (the fixed factor “Station” and the factor “Replicate” nested into “Station”), whereas a three-factor design was adopted in the analysis of the environmental variables (the fixed factor “Station,” the factor “Replicate” nested into “Station” and the fixed “Layer” factor). For the environmental parameters’ analysis, a resemblance matrix was calculated based on Euclidean distances, whereas a Bray–Curtis similarity matrix was used for the meiofauna community. Meiofauna abundance data were fourth root transformed prior to the analysis in order to decrease the influence of most abundant taxa (i.e., Nematoda). Environmental parameters – Chlorophyll C2, Pheophorbide a, Fucoxanthin, Diadinoxanthin, Diatoxanthin, Chlorophyll a, Pheophytin, Beta Carotene, CPE, TN, TOC, total organic matter (TOM) – with different units were analyzed together (Table 2) hence normalization was carried out prior to the analyses. Multicollinearity of environmental parameters was tested by means of correlation matrix (Pearson’s method). A distance-based redundancy analysis (dbRDA) was used to investigate the relationship between meiofauna and the environmental variables. This was based on Bray–Curtis distance and included automatic data standardization using the “metaMDSdist” argument. Environmental parameters were considered as predictor variables. To graphically visualize the different variation of the environmental variables between stations, a principal component analysis (PCA) was performed both on the complete dataset and subsets considering the significant factors of the dbRDA (Table 2). For the analysis of biodiversity, individual numbers of nematode species identified by morphological identification were standardized to ind/10 cm² prior to further statistical analysis. Biodiversity indexes were used to test the degree of diversity at each station, particularly the Shannon–Wiener test and the Simpson’s test were performed on the dataset. All analysis were performed with the “vegan” and “ggplot2” packages on Rstudio (v 3.6.1) versions 2.5–5 and 3.2.1.

Amplicon Sequence Variants Dataset

A multivariate analysis on the entire depth gradient (0–8 cm) of the ASVs dataset was performed (all replicates). The non-parametric PERMANOVA was used on the nematodes’ subset, followed by a PermDisp analysis (permutational analysis of multivariate dispersion). A three-factor design was adopted in the analysis of nematodes’ ASVs (the fixed factor “Station,” the factor “Replicate” nested into “Station” and the fixed “Layer” factor) using the Bray–Curtis similarity matrix. Data were fourth root transformed prior to the analysis in order to decrease the influence of the most abundant genera. A non-metric Multi-dimensional Scaling (nMDS) was subsequently performed on the relative abundance of nematode genera. The dbRDA was calculated to study the relationship between the nematodes community and the environmental variables at each station. The distance matrix was based on Bray–Curtis dissimilarity matrix and data were standardized using the “metaMDSdist” argument. An “adonis” analysis confirmed the significance of the model used. Environmental parameters were considered as independent variables while nematode community was the response variable. All above mentioned analysis were performed in Rstudio (version 3.6.1) using the “vegan,” and “ggplot2” packages (versions 2.5–5 and 3.2.1, respectively). To graphically visualize the intersection of the five ASVs datasets (where each station is a set), UpSet plots were generated using the “UpsetR” package (version 1.4.0) in Rstudio. The Upset plot shows the commonality between overlapping datasets (intersection between groups): it gives information both on the unique sequences present at a specific site and on the number of sequences shared between selected datasets (i.e., stations). Further biodiversity tests were performed by means of diversity indexes, particularly the Shannon–Wiener (H), Simpson (D), and inverse Simpson’s (1/D) indexes. Although both indexes provide similar information on communities’ diversity, they are differently influenced by species abundance, particularly the Shannon index is more affected by richness and rare species, while Simpson’s index gives more weight to common species. The use of both indexes can provide a more complete set of information on diversity. A PERMANOVA analysis was performed on the biodiversity indexes to test for significant differences between stations.

Results

Grain Size

Clay content of the 0–2 cm layer was rather constant among stations, ranging between 7 and 11%. In contrast, silt percentage was higher in DB 1,000 m (73%) and the two deep basins in PGC where it represented about 73 and 50% of the sediment, respectively. The three stations in DB clearly show a decreasing silt gradient with decreasing depth, from 73% (1,000 m) to 64% (500 m) to 46% (200 m). Interestingly, the two deep basins in PGC had a sand concentration more similar to DB 200 m than DB 1,000 m, although the silt fraction was high at all stations (Supplementary Annex C). In general, the grain size composition did not differ much along the sediment depth gradient, with all four layers showing similar concentrations (Table 2).

Pigments and Organic Material

The highest values of CPE in the 0–2 cm layer were observed in PGC 1,000 m (367.5 ± 172.63 μg/g), greatly exceeding that of the other stations (values ranging between 32.7 and 148.5 μg/g, see Supplementary Annex C). A clear increase along the bathymetric transect in DB was visible: a mean value of 32.7 ± 9.97 μg/g was observed at 200 m depth, 122.4 ± 7.59 μg/g at 500 m and 148.5 ± 40.05 μg/g at 1,000 m. The most southern and deepest station at PGC (1,250 m) showed the lowest concentration of CPE (27.3 ± 3.8 μg/g) and Chlorophyll a (9.45 ± 1.45 μg/g). The results of the Chl-a: CPE ratio were similar across stations, ranging between 0.34 and 0.43 (Supplementary Annex C). TOM showed high values in all the stations except PGC 1,250 m, where the lowest value was recorded (3.28 ± 0.44%), and followed no specific vertical gradient distribution within the sediment layers (Supplementary Annex D). The Chl-a vertical gradient did not decrease with sediment depth in DB in contrast to the PGC stations (Figure 2). In fact, the Chl-a was most abundant in the 2–4 cm layer in all three stations of DB. The vertical gradient of CPE followed the same trend as Chl-a, except for the DB 200 m station, where the amount of Chl-a decreased with depth (Supplementary Annex E and Figure 2, respectively).

FIGURE 2

Figure 2. Comparison between Duse Bay stations (from the top 200 m, 500 m and 1,000 m depth) and the Prince Gustav Channel (from the top 1,000 m and 1,250 m depth). The barplot shows Chl-a concentrations (μg/g, bottom axis, plus Sd) within the sediment layers (0–2, 2–4, 4–6, and 6–8 cm). Meiofauna distribution within the same layers expressed in ind 10 cm^–2 (upper axis) was plotted over the barplot. One replicate per layer was used to calculate meiofauna densities. Note the difference in scale between DB and PGC stations.

Environmental Characterization

The multivariate PERMANOVA analysis on environmental parameters resulted in a significant p-value (p < 0.05) for the “Station” factor, indicating statistically significant differences between stations. However, the permutation test for the homogeneity of dispersion (PermDisp) was significant, indicating that the differences in the dispersion of the group variance may influence the main test result and thus require caution in interpreting its significance. The PCA (Figure 3) shows separation of the group ellipses (confidence interval 0.95), pointing at differences in environmental variables between the different stations and depths. Some degree of overlap was only observed for the shallowest DB station and the nearest deep PGC station (1,000 m). The sum of the two main principal components (PC1 and PC2) explains together 73.22% of the variance. The main variables contributing along the PC1 axis are the ratio between Chlorophyll-a and CPE (Chloroplastic Pigment Equivalent; Chl-a:CPE), TOC and TN ratio (TOC:TN) and TOM concentration. The dbRDA resulted in a significant value (p < 0.05) for some of the pigments, phaeopigments, and TOM (Table 2) suggesting a correlation between the amount of fresh and/or degraded organic material and the meiofauna structure at each station.

FIGURE 3

Figure 3. Principal component analysis (PCA) of a subset of environmental characteristics of the sediment of all layers combined (i.e., 0–2, 2–4, 4–6, and 6–8 cm). Particularly, Chlorophyll a, CPE, Chl:CPE ratio, TOM, TOC:TN ratio, Clay, Silt and Sand parameters were considered. Confidence ellipses were plotted around group mean points, represented as larger symbols on the plot.

Meiofauna Densities and Composition

Results for meiofauna densities are displayed in Table 3. The highest densities were recorded in DB, with similar values (in the order of 6,000 individuals (ind) 10 cm^–2) at 500 and 200 m. Slightly lower values were detected in the two 1,000 m deep basins (one in DB and one in PGC), with abundances ranging in the order of 4,000 ind 10 cm^–2). The lowest values were observed in PGC 1,250 m (2,339.524 ± 785.2732 ind 10 cm^–2 ± SD), but it is noteworthy to mention that the SD value is high. The univariate PERMANOVA analysis on densities showed a significant p-value (p < 0.05) for the “Station” factor and the PermDisp was non-significant for the same. The vertical distribution of total meiofauna within the sediment followed a depth gradient analogous across all stations, where the 0–2 cm layer showed the highest densities. At each sampling site, nematodes were the major component of the meiofauna, constituting between 69 and 91% of the total abundance (Supplementary Annex F). Besides nematodes, the main groups present in the samples were nauplii (3–12%) and Copepoda (4–11%) in the first 0–2 cm layer, with similar relative densities at each of the stations of the bathymetric transect in DB. A different structure could be noticed in PGC, where both groups showed slightly lower densities at 1,000 m and particularly lower at 1,250 m (Supplementary Annex G).

TABLE 3

Table 3. Number of individuals per 10 cm^–2 of the 0–2 cm layer showed as mean, followed by standard deviation (SD) values.

Nematode Biodiversity

Amplicon Sequence Variants Analysis

The metabarcoding analysis resulted in a total of 7,175 ASVs belonging to the Eukarya domain, of which 2,564 were assigned the label “Nematoda.” It was possible to assign a genus to 40.6% of the sequences, whereas the 59.4% of them remained “Unassigned.” Nevertheless, all ASVs (including the “Unassigned”) were included in the biodiversity analysis (biodiversity tests and UpSet plot), whereas the NMDS and trophic analysis was only based on the genus-assigned ASVs. A total of 49 genera were identified across all stations, of which the most abundant were Molgolaimus, Daptonema, Halalaimus, Chromadorita, Sabatieria, and Bathylaimus (>11%). The PERMANOVA analysis performed on the ASVs dataset resulted in a significant p-value for the fixed factor “Station,” as well as the nested factor “Layer” (Table 4), suggesting that differences were also present in the vertical profile. The PermDisp test performed after the PERMANOVA analysis was non-significant for both factors. A pairwise comparison was performed subsequently, however, no significant results were detected. Hence, the interpretation of data must consider this issue. The intersection between datasets (all layers pooled together) is visualized in Figure 4. Here we observe that the majority of sequences (27.4%) were unique to PGC 1,250 m when all layers were considered followed by the shallowest one in DB 200 m (17.4%). As shown in the NMDS performed on the HTS dataset (Figure 5), the PGC 1,250 m and DB 200 m ellipses are markedly separated from each other and from the other stations, suggesting high variation in the assemblages of these two sites compared with the other stations. The DB 200 m ellipse was observed at a similar distance from PGC 1,250 m and DB 1,000 m, although some of DB 200 m sample points are closer to the ones in PGC 1,250 m than DB 1,000 m. This vicinity implies some degree of similarity between the two stations (DB 200 m and PGC 1,250 m) despite being located in different areas. The middle stations, regardless of area, were less isolated and closer to each other in the plot, with PGC 1,000 m and DB 500 m overlapping. This proximity implies a certain amount of similarity between the community composition of the two stations. The dbRDA analysis to investigate relationships between environmental variables and nematode diversity resulted in a significant value (p < 0.05) for the main pigments and TOM (Table 2), indicating some degree of association between the organic matter in the sediments and nematodes community structure. The results of the Shannon–Wiener (H), Simpson (D) and inverse Simpson’s (1/D) tests are shown in Table 5. Results show the mean value of each station (all replicates, all layers considered). The Shannon diversity index showed a high degree of biodiversity at all stations, with DB 500 m having the highest values (4.172) and PGC 1,250 the lowest (3.450). These are confirmed by both the Simpson’s test and its inverse, which is positively correlated to the Shannon index. Inverse Simpson’s values are in fact particularly high at all station, suggesting high biodiversity at all sites, and especially at DB 500 (35.150); the lowest value was recorded at PGC 1,250 (15.188). However, the PERMANOVA analysis performed on biodiversity indexes did not result in significant values for any of the ones considered.

TABLE 4

Table 4. Permutational analysis of variance results for the metabarcoding dataset, considering the three factors “Station,” “Layer,” and “Replicate” and the interaction between significant values (Station*Layer).

FIGURE 4

Figure 4. UpSet plot of all ASVs assigned a “Nematoda” label (all layers considered, three to four replicates). The left barplot (Set Size) represents the total number of ASVs per station; the bottom plot shows every possible intersection between sets and the top barplot displays their occurrence.

FIGURE 5

Figure 5. Non-metric Multi-dimensional Scaling analysis performed on the ASVs dataset (three to four replicates per station, all layers considered). The lighter ellipses represent the Duse Bay area, while the darker ones belong to the PGC stations. The further segregated ellipses (lower part of the graph) belong to DB 200 m and PGC 1,250 m, whereas an overlap between DB 500 m and PGC 1,000 m is clearly visible, with DB 1,000 m located closer to the “middle stations.”

TABLE 5

Table 5. Significant values for biodiversity indexes: Shannon–Wiener (H′), Simpson (D), inverse Simpson (1/D) on the metabarcoding dataset (ASVs).

Morphological Identification of Nematodes

A total of 107 nematode genera and 121 morphospecies were identified across the set of subsamples of 80 individuals per sample (0–2 cm, 3 samples per station). The greatest biodiversity was detected at 1,000 PGC and PGC 1,250 m (48 and 47 genera respectively in two replicates per station combined), although a substantial number of genera were observed at each sampling site (43 genera in DB 200 m; 39 in DB 500 m; and 45 in DB 1,000 m). The NMDS performed on this dataset (Supplementary Annex H) showed a large distance between the assemblages of the stations in DB (with DB 200 m and DB 500 m particularly isolated), while PGC 1,000 m was located at a similar distance from DB 1,000 m and PGC 1,250 m, when considering all data points. At each station, a few genera (e.g., Desmodora in DB 200 and Chromadorita in PGC 1,250 m) seem to dominate the assemblage (18.7 and 13.2%, respectively), while the majority of taxa represented less than 2% of the community. Figure 6 shows the comparison between the most abundant genera detected via morphological identification and ASVs analysis. Although the number of genera detected per station was almost identical, the community composition at each station differed. When ASVs were used, the genus Daptonema was dominant, constituting the 32.6% of all nematode genera at DB 500 m, whilst it represented only 11.4% in the morphological dataset at the same station. Similarly, in PGC 1,250 m Chromadorita was the most abundant genus using both identification methods: however, it represented 29.3% of all nematode genera when ASVs were used and only the 13.2% in the morphological dataset. It is worth noting that the genus Desmodora was absent in the ASVs at DB 200 m, while it represented the most abundant genus (18.7%) when the sample was identified morphologically.

FIGURE 6

Figure 6. Nematodes community composition for the most abundant genera (>3%). Comparison between ASVs (A) and morphologically identified dataset (B); less abundant genera (<3%) were pooled in the “Other” item. From left to right, image sequence of ice melting in the Eastern Antarctic Peninsula during the 4 months preceding the JR17003a research campaign. Satellite images source: https://worldview.earthdata.nasa.gov.

Trophic Guilds Analysis

The analysis of the trophic structure of nematode communities is shown in Supplementary Annex J. The first barplot (Supplementary Annex JA) shows the trophic composition of the ASVs dataset based on the assigned genera, where non-selective deposit feeders (Supplementary Annex JB) comprised between 46 and 61% of all the trophic groups both in DB and the PGC. Selective deposit feeders were the second most abundant trophic group in DB 200 m and DB 1,000 m (26% in both the stations), while epistrate feeders (2A) had similar percentages in DB 500 m, PGC 1,000 m and PGC 1,250 m (22, 22, and 30%, respectively). On the right (Supplementary Annex JB), the trophic analysis performed on the morphological dataset is shown, where epistrate feeders were dominant at all stations (between 31 and 50%).

Discussion

This study investigates the meiobenthos of DB and the PGC for the first time and sets the foundation to understand the benthic ecological dynamics in these areas by applying a unique combination of taxonomical identification and DNA analysis of organisms, with a focus on nematodes. In general, high levels of biodiversity were detected across stations in both areas (Table 5 and Figure 4), which are most likely associated to the abundance and freshness of food.

The DB Bathymetric Transect

Environmental Characterization

The local warming of the AP is responsible for increased algal blooms in the region, related to ice-shelf melting and calving events, which affect the abundance, structure, and dynamics of benthic organisms to a great extent (Veit-Köhler et al., 2018). Because of their tight connection with the sediment organic pool, meiofauna reflects the local dynamics of benthic-pelagic coupling processes. The structure of meiofauna communities in terms of densities and biodiversity are strictly linked to the amount of available organic matter. This is in turn dependent on deposition dynamics occurring in the upper water column. In general, the Chl-a and CPE values in DB were higher than those observed at similar depths in the AP (Vanhove et al., 2004; Hauquier et al., 2015), but relatively low in the DB 200 m station, where they showed an average concentration four times lower than that of the other DB stations (Table 2). The lower Chl-a values at the shallowest DB station, coincides with lower TOM concentrations and a higher sand fraction compared to the deeper DB stations. Sediment grain size also plays a crucial role in determining the infaunal assemblage structure and diversity since it relates to the organic matter accumulation/remineralization and oxygen penetration through the sediment. Sediment oxygenation has been observed to have a significant effect on meiofauna abundance and composition considering that different taxa have specific tolerances for oxygen availability (Schratzberger et al., 2006; Vanreusel et al., 2010). The sediment composition is highly affected by the lithology of parent rocks, the local mixing with recycled detritus (Dinis et al., 2017) together with the deposition and stratification of glacial and subglacial meltwater finer particles (Meslard et al., 2018). The latter are regularly resuspended in the water column because of local storms in shallow waters which is less the case in the deeper parts. The finer sediment composition in the deeper DB stations (DB 500 and 1,000 m, Table 2) was possibly a consequence of subglacial ablation and the associated disaggregation of the overlying sediment and its mixing with debris. Coupled with local hydrodynamics and geomorphology, this likely determined the transport of finer particles, i.e., silt, to deeper parts of the bay (Lawson, 1982; Nokes et al., 2019; Perolo et al., 2019), hence resulting in a muddier, soft sediment habitat compared to the shallowest station.

Meiofauna Community Structure

Differences in the sediment food availability of the 200 m station were not reflected in the meiofauna densities of the first 0–2 cm of the sediment, since they were equally high in the 500 m station and even higher at 1,000 m despite the lower chlorophyll-a values. In general, the meiofauna densities in the shallower stations (200 and 500 m) in DB were within the range of those found at other Antarctic shelf locations (Vanhove et al., 1995; Ingels et al., 2006; B_West station, Hauquier et al., 2015). High meiofauna densities are often related to the accumulation of phytoplanktonic organisms in the form of detritus on the surface of shelf sediments. Due to the cold temperatures of Antarctic waters, the process of degradation of the phytodetritus is slower compared to other areas, thus organic matter may accumulate in the form of “foodbanks” (Mincks et al., 2008; Smith et al., 2008). These organically rich sediments which are able to sustain the benthic densities throughout the year, preventing a reduction in biomass during austral winter months. High Chl-a and TOM concentration might be explained by such dynamics. The high densities found at DB 200 m corresponding to lower Chl-a and TOM concentration, might be due to the grazing action of benthic organisms on the food source prior to sampling, which would explain how such high abundances were sustained despite the low values of organic matter recorded. This is confirmed by the dbRDA performed on densities, which highlighted a significant correlation between organic matter and meiofauna abundances (Table 2). Copepods were the second most abundant taxon after nematodes (∼9.03% of the total meiofauna density in DB 200 m and ∼6.16% in DB 500 m; Supplementary Annex G), which is not surprising, considering that they are able to feed on both microalgae and bacteria (De Troch et al., 2005). The meiofauna densities seemed to decrease with increasing sediment depth, thus hosting the majority of the total meiofauna in the first 0–2 cm layer (Supplementary Annex A). The TOM followed the same decreasing vertical gradient, with higher concentrations in the first 0–2 cm layer. However, the Chl-a was distributed in a contrasting pattern in all DB stations, with the second sediment layer (2–4 cm) being the one with higher concentrations (Figure 3). The subsurface maximum of Chl-a in the DB stations might be in indication of the presence of larger benthic organisms (macrofauna) which are often responsible for burial of fresh phytodetritus into the deeper sediment layers (Furukawa, 2001; Lohrer et al., 2004; Meysman et al., 2006; Barsanti et al., 2011; van de Velde and Meysman, 2016). Biological interactions, such as the occurrence of macrobenthic organisms, are known to have a strong effect on the small-scale variability of nematodes, due to predator-prey dynamics and competition for food. The diversity of larger consumers has in fact been shown to be a crucial aspect in regulating the access of meiofauna to shared food sources (Nascimento et al., 2011). Although no macrofauna data are available in this study, it is still central to highlight the potential impact of these organisms. The values of CPE and Chl-a were highest at DB 1,000 m and corresponded to lower average meiofauna densities (although with a higher standard deviation). Macroalgal decay and accumulation of organic matter can lead to hypoxic conditions of the sediment surface (Pasotti et al., 2014; Hoffmann et al., 2018), influencing the abundance and composition of meiofauna assemblages. The possible hypoxic conditions of the sediment at this station was also suggested by the presence of the genus Anoplostoma, which was recorded in high numbers (15% of all nematode genera at this station) by morphological identification but not in the ASVs dataset (probably linked to the limited sequences deposited in online databases, see section “Morphological Identification Versus Metabarcoding: Drawbacks and Interpretation of Data”). This nematode genus is categorized as a selective deposit feeder often recorded in organically enriched sediments of mangrove ecosystems and muddy estuarine sediment types (Surey-Gent, 1981; Li and Guo, 2016), which are usually associated with hypoxic conditions, hence lower meiofauna abundances (Zeppilli et al., 2018). Moreover, one of the most abundant genera according to the ASVs dataset (Supplementary Annex H) was Sabatieria, which is also known for being tolerant of hypoxic/anoxic conditions in several environments (Wetzel et al., 2002; Steyaert et al., 2007; Franzo et al., 2019).

Nematode Biodiversity in the DB Area

Based on the molecular analysis of the ASVs, the degree of biodiversity in DB was observed to be very high at all stations, as confirmed by the Shannon (H), Simpson (D), and inverse Simpson’s (1/D) indexes (Table 5). The most diverse station was found to be DB 500 m in contrast to the highest number of unique ASVs (n = 446) recorded at 200 m when considering all sediment layers. The visual results of the UpSet plot (Figure 4) show a decrease of unique sequences in relation to depth, which suggests less biodiversity at the deepest site of the area (DB 1,000 m, 320 ASVs). This is in contrast with the H and 1/D values, which show little difference between stations, with slightly lower values for DB 200 m. Since the UpSet plot visualization is based on a presence/absence matrix, the reason behind such contrasts is likely due to the influence of richness (for H) and samples number (for D and 1/D), which rely on a distance-based similarity matrix. Moreover, the high values of the Shannon and inverse Simpson indexes suggest a uniform distribution of the species at all sites, confirmed by the non-significant PERMANOVA performed to test differences in biodiversity between stations. Previous studies show that nematode diversity is more often dependent on local environmental characteristics and hydrodynamics than bathymetry (Danovaro et al., 2013; Gambi et al., 2014). Also in this study, there was no evidence of increasing biodiversity in relation to depth, whereas differences between stations were more associated to the availability and freshness of food source (Table 2). The non-selective deposit feeder Daptonema was detected at all stations in DB, with highest concentrations at DB 200 m (Supplementary Annex J). This cosmopolitan genus has a long history of records in Antarctic waters, ranging from shallow to deeper sediments (Vanhove et al., 1999; Hauquier et al., 2015; Stark et al., 2020). Surprisingly only 11 ASVs were uniquely shared among all the DB stations (Figure 4), whereas the majority were unique to each station, thus it seems likely that the local biodiversity might have been influenced by gene flow external to DB, or by local hydrodynamics and seasonal sea-ice cover. A crucial element in determining the diversity of nematode communities are in fact local currents and the potential colonization from relatively close areas (Boeckner et al., 2009).

The Deep Basins

Ice Cover and Influence on Sediment Characteristics

The changes in sea-ice cover, sea-ice melting, and ice-shelf collapse have a major impact on the intensity and frequency of phytoplankton blooms, and as a consequence, on the primary productivity of the ocean surface in Antarctica (Buffen et al., 2007; Duprat et al., 2016; Van Leeuwe et al., 2018). Satellite images⁴ show that in the period preceding the JR17003a campaign, the sea ice started melting in proximity of the DB area at the end of December 2017, continuing the retreat southwards along the narrow seaway between James Ross Island and the Trinity Peninsula (Figures 7, 8). In fact, this melting allowed DB to be ice-free by the end of January 2018, almost contemporarily in all the stations of the area, together with the 1,000 m station in PGC. However, it looks like the southernmost station (PGC 1,250 m), was free of ice just by mid-February, almost 1 month later. This melting pattern might partly explain the over 10-fold higher Chl-a content in the first 0–2 cm layer in the PGC 1,000 m station, compared to the PGC 1,250 m station which showed the lowest Chl-a. Studies on the Larsen A and B Area, south of the PGC, have shown that newly ice-free regions stimulate phytoplanktonic blooms, affecting the total meiofauna density (Hauquier et al., 2015) and often increasing the abundance and biodiversity of nematodes (Raes et al., 2010). Although the two PGC stations are geographically close and of similar depth, the much lower quantities of Chl-a in PGC 1,250 m resulted in only 1.7 times lower meiofauna densities (Supplementary Annex G). Given the longer ice cover at this station it is possible that at the time of sampling a potential phytoplanktonic bloom had not yet taken place, with the meiofauna still thriving on food reserves from the previous year. To confirm these assumptions, a more detailed study on the timing of sea-ice retreat and the related spring phytoplankton bloom would be of significant importance. Other variables to consider are sea-ice thickness, water turbidity and the amount of freshening, which determines the stability of the water column and has a strong influence on phytoplanktonic blooms (Kang et al., 2001; Buffen et al., 2007). The distribution of sea ice depends on a number of factors which fluctuate seasonally and annually, namely wind intensity, wind direction and coastal currents (Scambos et al., 2009; Holland and Kwok, 2012; Stuecker et al., 2017). In light of the large variability of sea-ice extent and duration, dynamics relating to sea-ice cover are indeed particularly difficult to assess.

FIGURE 7

Figure 7. Satellite image of the sampling Area in January 2018, with PGC 1250 still covered by sea-ice while the rest of the stations are ice-free. The different colors indicate each station in the PGC and DB, as shown in the legend. Satellite images source: https://worldview.earthdata.nasa.gov.

FIGURE 8

Figure 8. From left to right, image sequence of ice melting in the Eastern Antarctic Peninsula during the four months preceding the JR17003a research campaign. Satellite images source: https://worldview.earthdata.nasa.gov.

Meiofauna Densities and Structure in the Deep Basins

The meiofauna abundance in the deepest PGC station (1,250 m) were lower compared to the other two basins (DB 1,000 m and PGC 1,000 m), but still high or similar in comparison to other findings at similar depths (Vanhove et al., 1995; Vermeeren, 2002; Ingels et al., 2011). Even though the lowest values of CPE and Chl-a were recorded at this deep site, the amount of food and its freshness, seemed in fact sufficient to sustain relatively high meiofauna abundances, with pigment concentrations within the range of those found in previous studies from adjacent areas (Hauquier et al., 2015). The two 1,000 m basins in DB and PGC seem to be acting as a sink for organic matter as they display the highest abundance of TOM in the whole region, and particularly PGC 1,000 m which exhibited almost double the Chl-a concentration compared to DB 1,000 m. The relative abundance of nauplii and copepods was also high in some stations (particularly copepods at DB 1,000 m and PGC 1,000 m, and nauplii at DB 1,000 m, Annex G). Nauplii are early larval stages of mainly copepods, and their assemblage can be a proxy for environmental changes. They are highly sensitive to environmental stressors, and indicate the recruitment of copepods, hence their reproductive success (Bunker and Hirst, 2004; Gonçalves et al., 2010; Darnis et al., 2019). The variability of larval stages is typically linked to the availability and quantity of food, thus potentially representing a proxy of surface primary productivity in DB. The difference in density of nauplii between the three deep basins (DB 1,000 m, PGC 1,000 m, and PGC 1,250 m) with four times lower densities in the PGC compared to DB might be due to a number of variables, such as differences in bottom currents, resuspension potential and vertical mixing of the two areas, since they have been shown to alter oxygen and food availability in the deep sediments (Isla et al., 2006). Based on the UpSetR plot analysis (Figure 4) and the NMDS performed on the first 0–2 cm layer of the three basins (DB 1,000 m, PGC 1,000 m, and PGC 1,250 m, Figure 5), marked differences in nematode community structure and abundance were visible. The dominant genera were shared among stations (i.e., Sabatieria, Molgolaimus, and Daptonema), albeit in different proportions. When all the sediment layers were considered, the deep PGC 1,250 m showed the highest number of unique ASVs (n = 703, Figure 4) and only 63 ASVs were shared among all stations. Shannon and Simpson’s indexes showed higher values for PGC 1,000 m, suggesting increased biodiversity at this site. However, biodiversity indexes were relatively high at all stations, and showed no significant differences between stations. As frequently reported in other areas worldwide (Raes et al., 2007; Semprucci et al., 2010), also in this study it seems that local small-scale variability has more influence on community structure than depth or even larger spatial distances. The differences in biodiversity through the sediment layers might be of great importance in order to explain the gene flow dynamics in the EAP (Derycke et al., 2013). In fact, endobenthos and especially nematodes mainly depend on passive transportation of individuals (i.e., via bottom currents) for their dispersal, given that they do not have a pelagic larval stage (Hauquier et al., 2017). Hence, their resuspension in the water column is expected to be influenced by their vertical segregation in the sediment (Eskin and Palmer, 1985). This might explain the relatively low abundance of shared ASVs and the distance of ellipses in the NMDS analysis even between geographically close stations (i.e., PGC 1,000 m and PGC 1,250 m), which are expected to be influenced by similar hydrodynamics and sea-ice melting patterns. The deepest basin, more isolated from the rest of the stations but closer to the open waters of the Weddell Sea, shows lower amounts of TOM, CPE and Chl-a and lower densities but the highest number of unique nematode ASVs. As the dbRDA analysis was significant for organic matter and chlorophyll a, suggesting one more time some degree of correlation between food source and nematodes community structure, phytoplanktonic dynamics regulated by local ice-cover are to be considered. It is appropriate to assume that differences between the 1,000 m stations and PGC 1,250 m might be linked to the different ice cover influence on the basins (from sedimentation rates to phytoplanktonic bloom effects), probably together with local physical variables such as bottom currents. Raes et al. (2010) found that higher nematode biodiversity in the B_South station (located south of the Larsen B area) was linked to a longer ice-free period, and open to potential phytoplanktonic blooms. According to their study, quick colonizers might have invaded the relatively new space, free of ice after the retreat of the ice shelf in 1995. The deep basin in PGC 1,250 m might have experienced similar colonization dynamics, also because of its proximity to more open waters of the Weddell Sea, which might have acted as a source for colonizers and food.

Morphological Identification Versus Metabarcoding: Drawbacks and Interpretation of Data

The documentation of nematodes has traditionally been performed via morphological identification (Vanhove et al., 1995; Fabiano and Danovaro, 1999; Lee et al., 2001; Tarjan and Keppner, 2008; Vanreusel et al., 2010). However, such work is known for being particularly time consuming and tedious, and it requires specific skills and deep knowledge of the taxa at hand. With the development of new molecular techniques, it has been proven that DNA metabarcoding is able to increase detection sensitivity, which in turn improves the ability to document dispersion of organisms, their connectivity and biodiversity (Fonseca et al., 2017). However, even new techniques such as metabarcoding are subject to drawbacks, and need to be perfectioned (Deagle et al., 2018). In this study we combined both approaches, using the morphological identification to integrate the data collected from the bio-molecular analysis. In general, a large number of genera per station was observed with both methods: the 0–2 cm layer seemed to host the majority of nematode diversity, with slightly more genera observed in the morphological identification (Supplementary Annex I). However, the composition of each station differed notably between the two methods, with many genera undetected in the ASVs dataset. In fact, a total of 107 genera was identified by morphological identification, and only 49 via metabarcoding. It is important to emphasize that differing results of nematode diversity may lead to different interpretation of the data, especially in regard to the assignment of nematode trophic groups and their relative abundance throughout the assemblage (Figure 6 and Supplementary Annex J). For example, the dominance of the morphologically identified genus Desmodora in DB 200 m (18.7%), indicated that the main trophic group at this station was represented by epistrate feeders (2A), suggesting that the much lower Chl-a concentration recorded at that station might be the result of grazing pressure by these organisms (probably coupled with macrobenthic grazing). However, the absence of Desmodora ASVs, leads to an opposite result where non-selective deposit feeders were the majority, probably influenced by the marked presence of Molgolaimus and Daptonema (Figure 6 and Annex J). Macheriotou et al., 2020 (under review) also highlighted a similar problem: although Desmodora was identified as dominant by morphological identification in a Mozambique channel site, it was absent from the ASVs analysis. This was most likely explained by failed PCR amplification, which might have been hindered by the thick cuticle typical for this genus. Given the inclusion of eight Desmodora sequences in the database, it is likely that these specimens would have been identified in the event of successful PCR amplification. While this explanation is in general satisfactory, it is worth considering also another matter: the possible lack of this particular Desmodora sequence in online databases. The majority of known species are indeed based on descriptions from few specimens originating from a small number of localities, which poses the problem of natural variation (Coomans, 2002). Although we would expect Desmodora to be identified at least to genus level, it is plausible to assume that one or more Desmodora sequences from Antarctica had not been included in ASVs databases yet, thus resulting in a lack of information on this organism. The lack of nematode sequences in online databases results in an incomplete coverage of nematode diversity, which is confirmed in this study also by the amount of non-identified sequences in the dataset (59.4% “Unassigned” ASVs). Finally, the mismatch between morphologically identified organisms and ASVs might also be linked to subjective bias, based on personal skills, eyesight, knowledge of the taxonomist, and the differentiation between related genera is not always easy (Vanreusel et al., 2010).

Conclusion

The meiobenthic community of the three stations in DB showed a distinct differentiation following the bathymetric gradient. The meiofauna densities decreased with depth along the transect, in contrast with the concentrations of Chl-a and TOM, which were instead increasing. The lower meiofauna density, despite the high food availability at the deepest DB station, was possibly due to hypoxic conditions of the sediments, as suggested by the presence of the genus Anoplostoma, and the higher concentration of TOM. Nematode diversity was higher at the shallowest stations, probably due to heterogeneity of the habitat and food source availability. The three deep basins in DB (1,000 m) and in the PGC (one at 1,000 m and one at 1,250 m depth) were different both in terms of meiofauna densities, pigments, TOM and nematode biodiversity, despite similar depths. Although meiofauna densities in the two 1,000 m basins were similarly high, the deeper PGC 1,250 m showed the lowest densities compared to the other stations. High biodiversity was detected at all stations. The dominant taxa were shared among the stations, however, PGC 1,250 m showed the highest number of unique nematode sequences. The PGC 1,250 m station is indeed located in a more elongated, less isolated basin and likely more impacted by the sea-ice system and the Weddell Sea, hence subject to external colonization. Our data suggest that meiofauna community structure was influenced by the simultaneous action of different variables, including food availability, habitat heterogeneity and seasonal sea-ice dynamics. In general, local physical processes might have had a strong influence both on the densities and on the connectivity of the stations considered. The application of morphological identification to assess nematode biodiversity proved to be a valuable tool, especially to identify less abundant or rare specimens, but nevertheless time consuming and influenced by subjective bias. HTS offers a faster, more sensitive, less laborious option: however, there is a necessity to improve the coverage of reference sequence databases and the potential limitation concerning thicker nematode’s cuticules must be taken into consideration. These issues suggest it is premature to completely abandon the traditional taxonomic identification of meiofauna organisms. On the contrary, future studies might benefit from an integrated approach of the two methods. Moreover, given the morphological natural variation of meiofauna in Antarctica and the presence of rare species, it would be of interest to create a genetic sequence database solely for Antarctic meiobenthos, especially nematodes. The data presented in this study constitute an initial investigation of meiofauna communities in these areas and their link to abiotic variables. In this regard it provides precious baseline information which could feed into future investigations for the creation of time series of temporal variability of the EAP.

Data Availability Statement

The data presented in this study are deposited in the Dryad repository, access number 10.5061, doi: 10.5061/dryad.n8pk0p2tr. The data is publicly available.

Author Contributions

AV performed the field-sampling, designed, and coordinated the study. GP processed the samples, performed the statistical analysis, and drafted the manuscript. FP participated in the statistical design and writing. LM participated in the molecular analysis and bioinformatics. All authors contributed to the reviewing process.

Funding

This work contributes to project BR/154/A1-Recto of the Belgian Science Policy (BELSPO/BRAIN). The environmental and molecular research was carried out with infrastructure funded by EMBRC Belgium (FWO project GOH3817N).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

The authors wish to thank the Principal Scientist Katrin Linse from the British Antarctic Survey for the coordination of the scientific program and the captain and crew of RRS James Clark Ross for the logistic support during the JR17003a expedition. Bruno Vlaeminck, Bart Beuselinck, Annelien Rigaux, Annick Van Kenhove, and Guy De Smet are thanked for their help with analysis of abiotic variables and preparation of biological samples.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2021.629706/full#supplementary-material

Footnotes

References

Ahmed, M., Sapp, M., Prior, T., Karssen, G., and Back, M. (2015). Nematode taxonomy: from morphology to metabarcoding. Soil Discuss. 12, 1175–1220. doi: 10.5194/soild-2-1175-2015