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.

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

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.

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.

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.

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.

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.

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%).

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.

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

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

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.