Front. Mar. Sci. Frontiers in Marine Science Front. Mar. Sci. 2296-7745 Frontiers Media S.A. 10.3389/fmars.2021.629706 Marine Science Original Research 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 Biology Department, Marine Biology Research Group, Ghent University, Ghent, Belgium

Edited by: Katrin Linse, British Antarctic Survey, United Kingdom

Reviewed by: Gritta Veit-Köhler, Senckenberg am Meer Wilhelmshaven, Germany; Federica Semprucci, University of Urbino Carlo Bo, Italy

*Correspondence: Gabriella Pantó, Gabriella.Panto@UGent.be

This article was submitted to Marine Evolutionary Biology, Biogeography and Species Diversity, a section of the journal Frontiers in Marine Science

23 07 2021 2021 8 629706 15 11 2020 18 06 2021 Copyright © 2021 Pantó, Pasotti, Macheriotou and Vanreusel. 2021 Pantó, Pasotti, Macheriotou and Vanreusel

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

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.

Eastern Antarctic Peninsula benthos nematodes Prince Gustav Channel Duse Bay metabarcoding ASVs biodiversity

香京julia种子在线播放

    1. <form id=HxFbUHhlv><nobr id=HxFbUHhlv></nobr></form>
      <address id=HxFbUHhlv><nobr id=HxFbUHhlv><nobr id=HxFbUHhlv></nobr></nobr></address>

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

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

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

      Date Area Depth (m) Latitude (S) Longitude (W) # cores
      2018–03-07 Duse Bay 200 63.62481 57.48192 3
      2018–03-07 Duse Bay 500 63.61546 57.49913 4
      2018–03-01 Duse Bay 1000 63.57557 57.29855 2
      2018–03-01 Duse Bay 1000 63.56889 57.29919 2
      2018–03-06 PGC mid 1000 63.76125 57.96736 3
      2018–03-05 PGC South 1250 63.97651 58.42942 2
      2018–03-05 PGC South 1250 63.397658 58.842954 2
      The table shows the date of each sampling event (each row represents a deployment), the depth of the station, coordinated and number of cores deployed.

      Sampling was carried out during the RRS James Clark Ross cruise JR17003a in February-March 20181. 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% Corg 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.

      Concentration of all the environmental variables analyzed.

      Chlorophyll Pheophorbide Fucoxanthin Diadinoxanthin Diatoxanthin Chlorophyll Pheophytin Beta carotene CPE TN TOC TOM Clay Silt Sand
      C2 (mg/g)
      a (μg/g)
      (μg/g)
      (μg/g)
      (μg/g)
      a (μg/g)
      a (μg/g)
      (μg/g)
      (μg/g)
      (%)
      (%)
      (%)
      (%)
      (%)
      (%)
      Station Layer Mean SD Mean SD Mean SD Mean SD Mean SD Mean SD Mean SD Mean SD Mean SD Mean SD Mean SD Mean SD Mean SD Mean SD Mean SD
      Duse Bay 200 0—2 0.98 0.54 9.53 2.22 3.85 1.46 0.28 0.1 1.66 0.57 13.24 5.96 2.36 0.22 0.77 0.24 32.67 9.97 0.2 0.09 0.68 0.25 5.66 0.16 7.71 0.75 49.46 2.05 42.83 2.74
      2—4 1.21 0.3 8.66 0.92 2.66 0.4 0.17 0.04 2.15 0.28 15.93 1.79 1.9 0.58 0.97 0.21 33.65 3.18 0.13 0.02 0.53 0.07 5.6 0.06 7.21 3.05 47.5 8.78 45.29 11.83
      4—6 1.25 0.5 7.39 1.35 2.5 1.13 0.15 0.13 1.99 0.14 10.11 3.35 1.63 0.32 0.99 0.07 26.02 6.7 0.14 0.09 0.52 0.24 5.51 0.16 7.15 1.77 46.83 5.44 46.02 6.46
      6—8 1.28 0.67 6.64 2.64 1.94 0.78 0.1 0.03 1.89 0.8 8.28 3.28 1.38 0.3 0.93 0.23 22.43 8.57 0.12 0.09 0.41 0.27 5.41 0.29 7.13 0.11 43.08 4.1 49.79 4.17
      Duse Bay 500 0—2 4.93 2.02 31.02 2.09 16.96 0.51 1.58 0.51 5.47 2.13 52.64 2.29 7.51 0.55 2.33 0.13 122.44 7.59 0.29 0.03 1.32 0.09 6.52 0.27 10.26 0.9 66.71 6.05 23.03 6.92
      2—4 6.33 3.8 21.8 4.14 12.84 7.36 1.26 0.38 6.97 1.03 59.51 18.66 5.62 1.33 2.37 0.51 116.7 36.68 0.29 0.02 1.2 0.05 5.91 0.2 10.58 3.35 62.08 15.49 27.34 18.76
      4—6 3.89 1.51 16.97 0.37 9.09 1.34 0.77 0.12 5.86 0.54 41.43 5.95 5.07 0.9 2.04 0.17 85.12 9.65 0.29 0.02 1.2 0.05 5.66 0.13 10.03 4.03 59.43 16.05 30.54 20.01
      6—8 4 0.85 15.11 1.42 7.11 1.01 0.63 0.12 5.7 0.45 33.54 5.66 4.97 0.09 1.92 0.14 72.97 9.39 0.29 0.05 1.1 0.1 5.6 0.13 11.15 2.61 68.25 7.81 20.6 10.42
      Duse Bay 1,000 0—2 6.94 1.99 48.98 7.23 19.61 6.76 2.34 1.09 5.06 1.19 56.57 19.33 6.7 3.64 2.33 0.62 148.53 40.05 0.26 0.03 2.17 0.35 7.99 1.08 10.49 0.58 70.08 2.73 19.43 2.73
      2—4 5.46 1.8 43.56 17.41 17.85 5.01 1.89 0.82 6.39 0.45 66.26 17.69 5.86 3.93 2.5 0.44 149.77 44.45 0.31 0.1 2.15 0.49 6.77 0.36 11.73 0.9 74.35 0.43 13.92 1.12
      4—6 6.67 2.03 26.99 1.62 13.94 2.31 1.35 0.18 6.4 0.64 57.53 5.16 4.91 1.28 2.06 0.3 119.85 10.3 0.23 0.05 1.85 0.34 6.34 0.24 12.55 0.25 74.67 1.84 12.78 1.97
      6—8 5.09 1.67 23.95 2 10.91 1.48 0.91 0.26 5.44 0.37 39.68 8.54 4.46 0.4 1.98 0.07 92.42 11.72 0.25 0.04 1.72 0.11 5.85 0.18 12.68 0.34 73.99 1.71 13.33 1.94
      PGC 1,000 0—2 12.4 6.44 116.46 56.85 52.39 29.26 8.07 3.72 15.52 15.06 131.42 48.47 26.71 14.71 4.56 1.78 367.53 172.63 0.33 0.09 1.29 0.26 6.28 0.88 8.22 0.33 60.86 5.24 30.92 5.14
      2—4 5.23 2.94 21.86 3.25 18.66 7.42 1.41 0.76 4.45 1.42 60.76 21.69 5.54 2.45 2.12 0.62 120.02 32.37 0.43 0.36 0.91 0.18 4.47 0.44 6.81 1.56 47.4 6.06 45.79 7.58
      4—6 1.89 1.65 10.01 3.61 5.73 1.86 0.3 0.2 3.07 0.5 20.05 6.67 2.21 0.4 1.19 0.15 44.46 14.39 0.48 0.5 0.69 0.05 3.19 0.25 9.15 0.46 51.24 2.15 39.61 2.61
      6—8 1.01 0.84 5.6 4 2.71 1.12 0.14 0.09 2.01 0.54 9.84 3.71 1.89 0.49 0.92 0.22 24.11 10.57 0.15 0.06 0.53 0.08 3.31 0.23 8.33 0.28 47.26 2.74 44.41 3.01
      PGC 1,250 0—2 1.32 0.44 9.76 1.04 3.35 0.61 0.27 0.09 1.36 0.25 9.45 1.45 1.22 0.22 0.6 0.1 27.33 3.8 0.17 0.05 0.83 0.14 3.28 0.44 11.46 0.45 61.06 3.65 27.48 3.46
      2—4 0.98 0.44 5.86 1.88 2.42 0.4 0.1 0.01 1.21 0.08 7.48 0.97 0.72 0.28 0.51 0.07 19.29 3.1 0.16 0.04 0.68 0.08 3.02 0.61 12.82 0.87 60.23 2.19 26.95 3.01
      4—6 0.24 0.22 3.11 1.4 0.66 0.48 0.01 0.02 0.55 0.22 2.24 0.9 0.71 0.18 0.35 0.07 7.88 3.28 0.13 0.05 0.62 0.12 3.14 0.37 9.75 2.95 57.93 4.68 32.32 7.23
      6—8 0.08 0.07 1.78 1.25 0.48 0.34 0 0 0.42 0.19 1.52 0.57 0.76 0.12 0.32 0.05 5.36 2.48 0.13 0.02 0.6 0.07 3.32 0.42 8.63 3.43 53.67 2.95 37.7 5.95
      dbRDA ASVs (p-value) 0.008** 0.035* 0.017* 0.062 0.085 0.156 0.069 0.178 - 0.707 0.08 0.010** 0.195 0.693 0.152
      dbRDA meiofauna (p-value) 0.049* 0.008** 0.131 0.082 0.081 0.031* 0.004** 0.072 - 0.142 0.033 0.007** 0.225 0.06 0.059
      The table reports the mean and standard deviation (SD) of three to four replicates per layer of each station. At the bottom, significant values for dbRDA analysis are reported for each variable: the first row reports values of the analysis performed on the metabarcoding dataset – dbRDA ASVs (p-value); the second row shows values for the analysis performed on morphological dataset-dbRDA meiofauna (p-value). The number of * indicates the degree of significance.
      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/cm3 (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 run2.

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

      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.

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

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

      Area Depth Replicates Mean SD
      Duse Bay 200 m 3 5,943.33 1,367.21
      Duse Bay 500 m 3 6,008.89 814.82
      Duse Bay 1,000 m 3 4,551.43 2,577.60
      PGC 1,000 m 3 3,993.97 1,045.78
      PGC 1,250 m 3 2,339.52 785.27
      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.

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

      Factors Df Sums of sqrs Mean sqs F. model R2 Pr (>F)
      Station 4 5.33 1.33 4.94 0.20 0.001***
      Layer 3 2.46 0.82 3.04 0.09 0.001***
      Replicate 67 20.28 0 0 0.79 1
      Station*Layer 12 3.81 0.32 1.17 0.15 0.023*
      The significant values are indicated in bold. The number of * indicates the degree of significance.

      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.

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

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

      Station H D 1/D
      Duse Bay 200 m 3.630 0.933 16.549
      Duse Bay 500 m 4.172 0.969 35.150
      Duse Bay 1,000 m 3.795 0.943 22.683
      PGC 1,000 m 3.845 0.947 20.350
      PGC 1,250 m 3.450 0.933 15.188
      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.

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

      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.

      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.

      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.

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

      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: /articles/10.3389/fmars.2021.629706/full#supplementary-material

      References Ahmed M. Sapp M. Prior T. Karssen G. Back M. (2015). Nematode taxonomy: from morphology to metabarcoding. Soil Discuss. 12 11751220. 10.5194/soild-2-1175-2015 Amplicon P. C. R. Clean-Up P. C. R. Index P. C. R. (2013). 16s Metagenomic Sequencing Library Preparation. Available online at: https://support.illumina.com/documents/documentation/chemistry_documentation/16s/16s-metagenomic-library-prep-guide-15044223-b.pdf (accessed March 15, 2019). Barsanti M. Delbono I. Schirone A. Langone L. Miserocchi S. Salvi S. (2011). Sediment reworking rates in deep sediments of the Mediterranean Sea. Sci. Total Environ. 409 29592970. 10.1016/j.scitotenv.2011.04.025 21561644 Bezerra T. N. Eisendle U. Hodda M. Holovachov O. Leduc D. Mokievsky V. (2021). Nemys: World Database of Nematodes. Available online at: http://www.marinespecies.org/imis.php?module=dataset&dasid=66 (accessed February 10, 2019). Boeckner M. J. Sharma J. Proctor H. C. (2009). Revisiting the meiofauna paradox: dispersal and colonization of nematodes and other meiofaunal organisms in low- and high-energy environments. Hydrobiologia 624 91106. 10.1007/s10750-008-9669-5 Brannock P. M. Learman D. R. Mahon A. R. Santos S. R. Halanych K. M. (2018). Meiobenthic community composition and biodiversity along a 5500 km transect of Western Antarctica: a metabarcoding analysis. Mar. Ecol. Prog. Ser. 603 4760. 10.3354/meps12717 Buffen A. Leventer A. Rubin A. Hutchins T. (2007). Diatom assemblages in surface sediments of the northwestern Weddell Sea, Antarctic Peninsula. Mar. Micropaleontol. 62 730. 10.1016/j.marmicro.2006.07.002 Bunker A. J. Hirst A. G. (2004). Fecundity of marine Planktonic copepods: global rates and patterns in relation to chlorophyll a, temperature and body weight. Mar. Ecol. Prog. Ser. 279 161181. 10.3354/meps279161 Camerlenghi A. Domack E. Rebesco M. Gilbert R. Ishman S. Leventer A. (2002). Glacial morphology and post-glacial contourites in northern Prince Gustav Channel (NW Weddell Sea, Antarctica). Mar. Geophys. Res. 22 417443. 10.1023/A:1016399616365 Carpentier A. Como S. Dupuy C. Lefrançois C. Feunteun E. (2013). Feeding ecology of Liza spp. in a tidal fl at: evidence of the importance of primary production (bio fi lm) and associated Meiofauna. J. Sea Res. 92 8691. 10.1016/j.seares.2013.10.007 Cook A. J. Fox A. J. Vaughan D. G. Ferrigno J. G. (2005). Retreating glacier fronts on the Antarctic Peninsula over the past half-century. Science 308 541544. 10.1126/science.1104235 15845851 Cook A. J. Vaughan D. G. (2010). Overview of areal changes of the ice shelves on the Antarctic Peninsula over the past 50 years. Cryosphere 4 7798. 10.5194/tc-4-77-2010 Coomans A. (2002). Present status and future of nematode systematics. Nematology 4 573582. Danovaro R. Carugati L. Corinaldesi C. Gambi C. Guilini K. Pusceddu A. (2013). Multiple spatial scale analyses provide new clues on patterns and drivers of deep-sea nematode diversity. Deep Sea Res. Part II Top. Stud. Oceanogr. 92 97106. Darnis G. Wold A. Falk-Petersen S. Graeve M. Fortier L. (2019). Could offspring predation offset the successful reproduction of the arctic copepod Calanus hyperboreus under reduced sea-ice cover conditions? Prog. Oceanogr. 170 107118. 10.1016/j.pocean.2018.11.004 Davies B. J. Carrivick J. L. Glasser N. F. Hambrey M. J. Smellie J. L. (2012). Variable glacier response to atmospheric warming, Northern Antarctic Peninsula, 1988-2009. Cryosphere 6 10311048. 10.5194/tc-6-1031-2012 De Troch M. Steinarsdóttir M. Chepurnov V. Ólafsson E. (2005). Grazing on diatoms by harpacticoid copepods: species-speciWc density-dependent uptake and microbial gardening. Aquat. Microb. Ecol. 39 135144. Deagle B. E. Clarke L. J. Kitchener J. A. Polanowski A. M. Davidson A. T. (2018). Genetic monitoring of open ocean biodiversity: an evaluation of DNA metabarcoding for processing continuous plankton recorder samples. Mol. Ecol. Resources 18 391406. 10.1111/1755-0998.12740 29171158 Dell’ Anno A. Carugati L. Corinaldesi C. Riccioni G. (2015). Unveiling the biodiversity of deep-sea nematodes through metabarcoding: are we ready to bypass the classical taxonomy? PLoS One 10:e0144928. 10.1371/journal.pone.0144928 26701112 Derycke S. Backeljau T. Moens T. (2013). Dispersal and gene flow in free-living marine nematodes. Front. Zool. 10:1. 10.1186/1742-9994-10-1 23356547 Dinis P. Garzanti E. Vermeesch P. Huvi J. (2017). Climatic zonation and weathering control on sediment composition (Angola). Chem. Geol. 467 110121. 10.1016/j.chemgeo.2017.07.030 Duprat L. P. A. M. Bigg G. R. Wilton D. J. (2016). Enhanced Southern Ocean marine productivity due to fertilization by giant icebergs. Nat. Geosci. 9 219221. 10.1038/ngeo2633 Eskin R. A. Palmer M. A. (1985). Suspension of marine nematodes in a turbulent tidal creek: species patterns. Biol. Bull. 169 615623. Etourneau J. Sgubin G. Crosta X. Swingedouw D. Willmott V. Barbara L. (2019). Ocean temperature impact on ice shelf extent in the eastern Antarctic Peninsula. Nat. Commun. 10 815. 10.1038/s41467-018-08195-6 30659177 Fabiano M. Danovaro R. (1999). Meiofauna distribution and mesoscale variability in two sites of the Ross Sea (Antarctica) with contrasting food supply. Polar Biol. 22 115123. 10.1007/s003000050398 Fonseca V. G. Sinniger F. Gaspar J. M. Quince C. Creer S. Power D. M. (2017). Revealing higher than expected meiofaunal diversity in Antarctic sediments: a metabarcoding approach. Sci. Rep. 7 111. 10.1038/s41598-017-06687-x 28733608 Franzo A. Asioli A. Roscioli C. Patrolecco L. Bazzaro M. Del Negro P. (2019). Influence of natural and anthropogenic disturbances on foraminifera and free-living nematodes in four lagoons of the Po delta system. Estuar. Coast. Shelf Sci. 220 99110. 10.1016/j.ecss.2019.02.039 Furukawa Y. (2001). Biogeochemical consequences of macrofauna burrow ventilation. Geochem. Trans. 2:83. Gambi C. Pusceddu A. Benedetti-Cecchi L. Danovaro R. (2014). Species richness, species turnover and functional diversity in nematodes of the deep Mediterranean S ea: searching for drivers at different spatial scales. Glob. Ecol. Biogeogr. 23 2439. Gogina M. Zettler M. L. Vanaverbeke J. Dannheim J. Van Hoey G. Desroy N. (2020). Interregional comparison of benthic ecosystem functioning: community bioturbation potential in four regions along the NE Atlantic shelf. Ecol. Indic. 110:105945. 10.1016/j.ecolind.2019.105945 Gonçalves A. M. M. Pardal M. A. Marques S. C. De Troch M. Azeiteiro U. M. (2010). Distribution and composition of small-size zooplankton fraction in a temperate shallow estuary (Western Portugal). Fresenius Environ. Bull. 19 31603176. Hauquier F. Leliaert F. Rigaux A. Derycke S. Vanreusel A. (2017). Distinct genetic differentiation and species diversification within two marine nematodes with different habitat preference in Antarctic sediments. BMC Evol. Biol. 17:120. 10.1186/s12862-017-0968-1 28558672 Hauquier F. Suja L. D. Gutt J. Veit-Köhler G. Vanreusel A. (2015). Different oceanographic regimes in the vicinity of the Antarctic Peninsula reflected in benthic nematode communities. PLoS One 10:e0137527. 10.1371/journal.pone.0137527 26355457 Heip C. Vincx M. Smol N. Vranken G. (1982). The systematics and ecology of marine free-living nematodes. Helminthological Abstracts Series B, Plant Nematology 51 131. Hoffmann R. Al-Handal A. Y. Wulff A. Deregibus D. Zacher K. Quartino M. L. (2019). Implications of glacial melt-related processes on the potential primary production of a microphytobenthic community in potter cove (Antarctica). Front. Mar. Sci. 6:655. 10.3389/fmars.2019.00655 Hoffmann R. Pasotti F. Vázquez S. Lefaible N. Torstensson A. MacCormack W. (2018). Spatial variability of biogeochemistry in shallow coastal benthic communities of Potter Cove (Antarctica) and the impact of a melting glacier. PLoS One 13:e0207917. 10.1371/journal.pone.0207917 30566444 Holland P. R. Kwok R. (2012). Wind-driven trends in Antarctic sea-ice drift. Nat. Geosci. 5 872875. 10.1038/ngeo1627 Ingels J. Tchesunov A. V. Vanreusel A. (2011). Meiofauna in the gollum channels and the whittard canyon, celtic margin-how local environmental conditions shape nematode structure and function. PLoS One 6:e020094. 10.1371/journal.pone.0020094 21629829 Ingels J. Vanhove S. De Mesel I. Vanreusel A. (2006). The biodiversity and biogeography of the free-living nematode genera Desmodora and Desmodorella (family Desmodoridae) at both sides of the Scotia Arc. Polar Biol. 29 936949. 10.1007/s00300-006-0135-4 Ingels J. Vanreusel A. Brandt A. Catarino A. I. David B. De Ridder C. (2012). Possible effects of global environmental changes on Antarctic benthos: a synthesis across five major taxa. Ecol. Evol. 2 453485. 10.1002/ece3.96 22423336 Isla E. Gerdes D. Palanques A. Gili J. M. Arntz W. (2006). Particle fluxes and tides near the continental ice edge on the eastern Weddell Sea shelf. Deep. Res. II 53:866. Jensen P. (1982). A new meiofauna sample splitter. Ann. Zool. Fenn. 9 233236. Kang S.-H. Kang J.-S. Lee S. Chung K. H. Kim D. Park M. G. (2001). Antarctic phytoplankton assemblages in the marginal ice zone of the northwestern Weddell Sea. J. Plankton Res. 23 333352. Kennedy A. D. Jacoby C. A. (1999). Biological indicators of marine environmental health: meiofauna - a neglected benthic component? Environ. Monit. Assess. 54 4768. Korbel K. L. Stephenson S. Hose G. C. (2019). Sediment size influences habitat selection and use by groundwater macrofauna and meiofauna. Aquat. Sci. 81 110. 10.1007/s00027-019-0636-1 Kunz M. King M. A. Mills J. P. Miller P. E. Fox A. J. Vaughan D. G. (2012). Multi-decadal glacier surface lowering in the Antarctic Peninsula. Geophys. Res. Lett. 39 15. 10.1029/2012GL052823 Lawson D. (1982). Mobilization, movement and deposition of active subaerial sediment flows, matanuska glacier, Alaska. J. Geol. 90 279300. 10.1086/628680 Lee H. Vanhove S. Peck L. S. Vincx M. (2001). Recolonisation of meiofauna after catastrophic iceberg scouring in shallow Antarctic sediments. Polar Biol. 24 918925. 10.1007/s003000100300 Li Y. Guo Y. (2016). Two new free-living marine nematode species of the genus Anoplostoma (Anoplostomatidae) from the mangrove habitats of Xiamen Bay, East China Sea. J. Ocean Univ. China 15 1118. 10.1007/s11802-016-2896-x Lohrer A. M. Thrush S. F. Gibbs M. M. (2004). Bioturbators enhance ecosystem function through complex biogeochemical interactions. Nature 431 10921095. 10.1038/nature03042 15470385 Macheriotou L. Guilini K. Bezerra T. N. Tytgat B. Nguyen D. T. Phuong Nguyen T. X. (2019). Metabarcoding free-living marine nematodes using curated 18S and CO1 reference sequence databases for species-level taxonomic assignments. Ecol. Evol. 9 12111226. 10.1002/ece3.4814 30805154 Macheriotou L. Rigaux A. Derycke S. Vanreusel A. (2020). Phylogenetic clustering and rarity imply risk of local species extinction in prospective deep-sea mining areas of the clarion-clipperton fracture zone. Proc. R. Soc. B 287:20192666. 10.1098/rspb.2019.2666 32228410 Martín C. Gudmundsson G. H. Pritchard H. D. Gagliardini O. (2009). On the effects of anisotropic rheology on ice flow, internal structure, and the age-depth relationship at ice divides. J. Geophys. Res. Earth Surf. 114 118. 10.1029/2008JF001204 Meslard F. Bourrin F. Many G. Kerhervé P. (2018). Suspended particle dynamics and fluxes in an Arctic fjord (Kongsfjorden, Svalbard). Estuar. Coast. Shelf Sci. 204 212224. 10.1016/j.ecss.2018.02.020 Meysman F. J. R. Galaktionov O. S. Gribsholt B. Middelburg J. J. (2006). Bioirrigation in permeable sediments: advective pore-water transport induced by burrow ventilation. Limnol. Oceanogr. 51 142156. 10.4319/lo.2006.51.1.0142 Mincks S. L. Smith C. R. Jeffreys R. M. Sumida P. Y. G. (2008). Trophic structure on the West Antarctic Peninsula shelf, detritivory and benthic inertia revealed by d13C and d15N analysis. Deep. Res. 55:2502. Moline M. A. Prezelin B. B. (1996). Long-term monitoring and analyses of physical factors regulating variability in coastal Antarctic phytoplankton biomass, i n situ productivity and taxonomic composition over subseasonal, seasonal and interannual time scales. Mar. Ecol. Prog. Ser. 145 143160. Montes-Hugo M. Doney S. C. Ducklow H. W. Fraser W. Martinson D. Stammerjohn S. E. (2009). Recent Changes in Phytoplankton communities associated with rapid regional climate change along the Western Antarctic Peninsula. Science 323 14701473. Moreno M. Semprucci F. Vezzulli L. Balsamo M. Fabiano M. Albertelli G. (2011). The use of nematodes in assessing ecological quality status in the Mediterranean coastal ecosystems. Ecol. Indic. 11 328336. 10.1016/j.ecolind.2010.05.011 Nardelli S. Kohut J. Schofield O. Waite N. Brown M. Saba G. (2018). Changes in the upper ocean mixed layer and phytoplankton productivity along the West Antarctic Peninsula. Philos. Trans. R. Soc. A Math. Eng. Sci. 376 20170173. 10.1098/rsta.2017.0173 29760116 Nascimento F. J. A. Karlson A. M. L. Elmgren R. (2011). Diversity of larger consumers enhances interference competition effects on smaller competitors. Oecologia 166 337347. 10.1007/s00442-010-1865-0 21161548 Nokes C. R. Bostock H. C. Strachan L. J. Hadfield M. G. (2019). Hydrodynamics and sediment transport on the North Canterbury Shelf, New Zealand. New Zeal. J. Mar. Freshw. Res. 55 120. 10.1080/00288330.2019.1699584 Oppenheimer M. (1998). Global warming and the stability of the West Antarctic Ice Sheet. Nature 393 325332. Pasotti F. Convey P. Vanreusel A. (2014). Potter Cove, west Antarctic Peninsula, shallow water meiofauna: a seasonal snapshot. Antarct. Sci. 26 554562. 10.1017/S0954102014000169 Perolo P. Bakker M. Gabbud C. Moradi G. Rennie C. Lane S. N. (2019). Subglacial sediment production and snout marginal ice uplift during the late ablation season of a temperate valley glacier. Earth Surf. Process. Landf. 44 11171136. 10.1002/esp.4562 Polvani L. M. Waugh D. W. Correa G. J. Son S. W. (2011). Stratospheric ozone depletion: the main driver of twentieth-century atmospheric circulation changes in the Southern Hemisphere. J. Clim. 24 795812. 10.1175/2010JCLI3772.1 Prezèlin B. B. Hofmann E. E. Mengelt C. Klinck J. M. (2000). The linkage between Upper Circumpolar Deep Water (UCDW) and phytoplankton assemblages on the west Antarctic Peninsula continental shelf. J. Mar. Res. 58 165202. Pritchard H. D. Ligtenberg S. R. M. Fricker H. A. Vaughan D. G. Van Den Broeke M. R. Padman L. (2012). Antarctic ice-sheet loss driven by basal melting of ice shelves. Nature 484 502505. 10.1038/nature10968 22538614 Pritchard H. D. Vaughan D. G. (2007). Widespread acceleration of tidewater glaciers on the Antarctic Peninsula. J. Geophys. Res. Earth Surf. 112 110. 10.1029/2006JF000597 Raes M. De Troch M. Ndaro S. G. M. Muthumbi A. Guilini K. Vanreusel A. (2007). The structuring role of microhabitat type in coral degradation zones: a case study with marine nematodes from Kenya and Zanzibar. Coral Reefs 26 113126. 10.1007/s00338-006-0184-8 Raes M. Rose A. Vanreusel A. (2010). Response of nematode communities after large-scale ice-shelf collapse events in the Antarctic Larsen area. Glob. Chang. Biol. 16 16181631. 10.1111/j.1365-2486.2009.02137.x Rignot E. Jacobs S. S. (2002). Rapid Bottom melting Widespread near Antarctic Ice sheet grounding lines. Science 296 20202023. 10.1126/science.1070942 12065835 Rott H. Abdel Jaber W. Wuite J. Scheiblauer S. Floricioiu D. Van Wessem J. M. (2018). Changing pattern of ice flow and mass balance for glaciers discharging into the Larsen A and B embayments, Antarctic Peninsula, 2011 to 2016. Cryosphere 12 12731291. 10.5194/tc-12-1273-2018 Rott H. Rack W. Skvarca P. De Angelis H. (2002). Northern Larsen Ice Shelf, Antarctica: further retreat after collapse. Ann. Glaciol. 34 277282. 10.3189/172756402781817716 Scambos T. Fricker H. A. Liu C.-C. Bohlander J. Fastook J. Sargent A. (2009). Ice shelf disintegration by plate bending and hydro-fracture: satellite observations and model results of the 2008 Wilkins ice shelf break-ups. Earth Planet. Sci. Lett. 280 5160. 10.1016/j.epsl.2008.12.027 Schmid-Araya J. M. Schmid P. E. Tod S. P. Esteban G. F. (2016). Trophic positioning of meiofauna revealed by stable isotopes and food web analyses. Ecology 97 30993109. Schratzberger M. Ingels J. (2018). Meiofauna matters: the roles of meiofauna in benthic ecosystems. J. Exp. Mar. Bio. Ecol. 502 1225. 10.1016/j.jembe.2017.01.007 Schratzberger M. Warr K. Rogers S. I. (2006). Patterns of nematode populations in the southwestern North Sea and their link to other components of the benthic fauna. J. Sea Res. 55 113127. 10.1016/j.seares.2005.07.002 Seinhorst J. W. (1959). A rapid method for the transfer of nematodes from fixative to anhydrous glycerin. Nematologica 4 6769. Semprucci F. Colantoni P. Baldelli G. Rocchi M. Balsamo M. (2010). The distribution of meiofauna on back-reef sandy platforms in the Maldives (Indian Ocean). Mar. Ecol. 31 592607. 10.1111/j.1439-0485.2010.00383.x Semprucci F. Sbrocca C. Rocchi M. Balsamo M. (2015). Temporal changes of the meiofaunal assemblage as a tool for the assessment of the ecological quality status. J. Mar. Biol. Assoc. U.K. 95 247254. 10.1017/S0025315414001271 Shuman C. A. Berthier E. Scambos T. A. (2011). 2001-2009 elevation and mass losses in the Larsen A and B embayments, Antarctic Peninsula. J. Glaciol. 57 737754. 10.3189/002214311797409811 Siegert M. Atkinson A. Banwell A. Brandon M. Convey P. Davies B. (2019). The Antarctic Peninsula under a 1.5°C global warming scenario. Front. Environ. Sci. 7:102. 10.3389/fenvs.2019.00102 Smith C. R. Mincks S. DeMaster D. J. (2008). The FOODBANCS project, Introduction and sinking fluxes of organic carbon, chlorophyll-a and phytodetritus on the western Antarctic Peninsula continental shelf. Deep. Res. II 55:2404. Stark J. S. Mohammad M. McMinn A. Ingels J. (2020). Diversity, abundance, spatial variation, and human impacts in marine Meiobenthic nematode and copepod communities at casey station, East Antarctica. Front. Mar. Sci. 7:480. 10.3389/fmars.2020.00480 Steyaert M. Moodley L. Nadong T. Moens T. Soetaert K. Vincx M. (2007). Responses of intertidal nematodes to short-term anoxic events. J. Exp. Mar. Bio. Ecol. 345 175184. 10.1016/j.jembe.2007.03.001 Steyaert M. Priestley V. Osborne O. Herraiz A. Arnold R. Savolainen V. (2020). Advances in metabarcoding techniques bring us closer to reliable monitoring of the marine benthos. J. Appl. Ecol. 57 22342245. Stuecker M. F. Bitz C. M. Armour K. C. (2017). Conditions leading to the unprecedented low Antarctic sea ice extent during the 2016 austral spring season. Geophys. Res. Lett. 44 90089019. 10.1002/2017GL074691 Surey-Gent S. C. (1981). Distribution of Anoplostoma viviparum (Nematoda) in Southampton Water sediments. Mar. Biol. 62 157160. Taheri M. Bastami K. D. Yazdani M. (2017). “The role of the sediment conditions in shaping meiofauna spatial distribution in the shallow water of the south Caspian Sea,” in Proceedings of the First International Conference on Oceanography for West Asia (Tehran: Civilica). Tarjan A. C. Keppner E. J. (2008). Illustrated Key to the Genera of Free-Living Marine Nematodes in the Superfamily Chromadoroides Exclusive of the Chromadoridae. Available online at: https://entnemdept.ufl.edu/creatures/nematode/marine_nematodes.htm (accessed December 8, 2019). van de Velde S. Meysman F. J. R. (2016). The influence of bioturbation on iron and sulphur cycling in marine sediments: a model analysis. Aquat. Geochem. 22 469504. Van Heukelem L. Thomas C. S. (2001). Computer-assisted high-performance liquid chromatography method development with applications to the isolation and analysis of phytoplankton pigments. J. Chromatogr. A 910 3149. 10.1016/S0378-4347(00)00603-4 Van Leeuwe M. A. Tedesco L. Arrigo K. R. Assmy P. Meiners K. M. Rintala J. (2018). Microalgal community structure and primary production in Arctic and Antarctic sea ice: a synthesis. Elem. Sci. Anthropocene 6:4. Vanhove S. Arntz W. Vincx M. (1999). Comparative study of the nematode communities on the southeastern Weddell Sea shelf and slope (Antarctica). Mar. Ecol. Prog. Ser. 181 237256. 10.3354/meps181237 Vanhove S. Vermeeren H. Vanreusel A. (2004). Meiofauna towards the South sandwich trench (750-6300 m), focus on nematodes. Deep. Res. Part II Top. Stud. Oceanogr. 51 16651687. 10.1016/j.dsr2.2004.06.029 Vanhove S. Wittoeck J. Desmet G. Van den Berghe B. Herman R. L. Bak R. P. M. (1995). Deep-sea meiofauna communities in Antarctica: structural analysis and relation with the environment. Mar. Ecol. Prog. Ser. 127 6576. Vanreusel A. Fonseca G. Danovaro R. Cristina M. Ferrero T. Gad G. (2010). The contribution of deep-sea macrohabitat heterogeneity to global nematode diversity. Mar. Ecol. 31 620. 10.1111/j.1439-0485.2009.00352.x Vaughan D. G. Marshall G. J. Connolley W. M. Parkinson C. Mulvaney R. Hodgson D. A. (2003). Recent rapid regional climate warming on the Antarctic Peninsula. Clim. Chang. 60 243274. Veit-Köhler G. Durst S. Schuckenbrock J. Hauquier F. Durán Suja L. Dorschel B. (2018). Oceanographic and topographic conditions structure benthic meiofauna communities in the Weddell Sea, Bransfield Strait and Drake Passage (Antarctic). Prog. Oceanogr. 162 240256. 10.1016/j.pocean.2018.03.005 Vermeeren H. (2002). Biogeografie van Antarctische Diepzeenematoden: Species Turn-Over in Dominante Genera van de Familie Chromadoridae. MSc thesis. Gent: Universiteit Gent. Vinagre C. Salgado J. P. Mendonça V. Cabral H. Costa M. J. (2012). Isotopes reveal fl uctuation in trophic levels of estuarine organisms, in space and time. J. Sea Res. 72 4954. 10.1016/j.seares.2012.05.010 Warwick R. M. Platt H. M. Somerfield P. J. (1998). “Freeliving marine nematodes part III Monhysterids: pictorial key to world genera and notes for the identification of British species,” in Synopses of the British Fauna (New Series) 53, eds Barnes R. S. K. Crothers J. H. (Shrewsbury: Field Studies Council), 296. Wetzel M. A. Weber A. Giere O. (2002). Re-colonization of anoxic/sulfidic sediments by marine nematodes after experimental removal of macroalgal cover. Mar. Biol. 141 679689. 10.1007/s00227-002-0863-0 Wieser W. (1953). Die Beziehung zwischen Mundhohlengestalt, ernahrungsweise und vorkommen bei freilebenden marinen Nematoden. Ark. Zool. 4 439484. Zeppilli D. Fuchs S. Sanchez N. Sarrazin J. Leduc D. Fontanier C. (2018). Characteristics of meiofauna in extreme marine ecosystems: a review. Mar. Biodivers. 48 3571. 10.1007/s12526-017-0815-z

      https://www.bodc.ac.uk/resources/inventories/cruise_inventory/report/16954/

      https://genomics.ed.ac.uk/

      https://benjjneb.github.io/dada2/index.html

      https://worldview.earthdata.nasa.gov

      ‘Oh, my dear Thomas, you haven’t heard the terrible news then?’ she said. ‘I thought you would be sure to have seen it placarded somewhere. Alice went straight to her room, and I haven’t seen her since, though I repeatedly knocked at the door, which she has locked on the inside, and I’m sure it’s most unnatural of her not to let her own mother comfort her. It all happened in a moment: I have always said those great motor-cars shouldn’t be allowed to career about the streets, especially when they are all paved with cobbles as they are at Easton Haven, which are{331} so slippery when it’s wet. He slipped, and it went over him in a moment.’ My thanks were few and awkward, for there still hung to the missive a basting thread, and it was as warm as a nestling bird. I bent low--everybody was emotional in those days--kissed the fragrant thing, thrust it into my bosom, and blushed worse than Camille. "What, the Corner House victim? Is that really a fact?" "My dear child, I don't look upon it in that light at all. The child gave our picturesque friend a certain distinction--'My husband is dead, and this is my only child,' and all that sort of thing. It pays in society." leave them on the steps of a foundling asylum in order to insure [See larger version] Interoffice guff says you're planning definite moves on your own, J. O., and against some opposition. Is the Colonel so poor or so grasping—or what? Albert could not speak, for he felt as if his brains and teeth were rattling about inside his head. The rest of[Pg 188] the family hunched together by the door, the boys gaping idiotically, the girls in tears. "Now you're married." The host was called in, and unlocked a drawer in which they were deposited. The galleyman, with visible reluctance, arrayed himself in the garments, and he was observed to shudder more than once during the investiture of the dead man's apparel. HoME香京julia种子在线播放 ENTER NUMBET 0016fg500.com.cn
      hknbsoft.com.cn
      www.formycon.com.cn
      dyxdhs.com.cn
      lcchain.com.cn
      sxyhhk.org.cn
      www.simibaby.com.cn
      www.thrjn.com.cn
      pjrbtd.com.cn
      woobuy.com.cn
      处女被大鸡巴操 强奸乱伦小说图片 俄罗斯美女爱爱图 调教强奸学生 亚洲女的穴 夜来香图片大全 美女性强奸电影 手机版色中阁 男性人体艺术素描图 16p成人 欧美性爱360 电影区 亚洲电影 欧美电影 经典三级 偷拍自拍 动漫电影 乱伦电影 变态另类 全部电 类似狠狠鲁的网站 黑吊操白逼图片 韩国黄片种子下载 操逼逼逼逼逼 人妻 小说 p 偷拍10幼女自慰 极品淫水很多 黄色做i爱 日本女人人体电影快播看 大福国小 我爱肏屄美女 mmcrwcom 欧美多人性交图片 肥臀乱伦老头舔阴帝 d09a4343000019c5 西欧人体艺术b xxoo激情短片 未成年人的 插泰国人夭图片 第770弾み1 24p 日本美女性 交动态 eee色播 yantasythunder 操无毛少女屄 亚洲图片你懂的女人 鸡巴插姨娘 特级黄 色大片播 左耳影音先锋 冢本友希全集 日本人体艺术绿色 我爱被舔逼 内射 幼 美阴图 喷水妹子高潮迭起 和后妈 操逼 美女吞鸡巴 鸭个自慰 中国女裸名单 操逼肥臀出水换妻 色站裸体义术 中国行上的漏毛美女叫什么 亚洲妹性交图 欧美美女人裸体人艺照 成人色妹妹直播 WWW_JXCT_COM r日本女人性淫乱 大胆人艺体艺图片 女同接吻av 碰碰哥免费自拍打炮 艳舞写真duppid1 88电影街拍视频 日本自拍做爱qvod 实拍美女性爱组图 少女高清av 浙江真实乱伦迅雷 台湾luanlunxiaoshuo 洛克王国宠物排行榜 皇瑟电影yy频道大全 红孩儿连连看 阴毛摄影 大胆美女写真人体艺术摄影 和风骚三个媳妇在家做爱 性爱办公室高清 18p2p木耳 大波撸影音 大鸡巴插嫩穴小说 一剧不超两个黑人 阿姨诱惑我快播 幼香阁千叶县小学生 少女妇女被狗强奸 曰人体妹妹 十二岁性感幼女 超级乱伦qvod 97爱蜜桃ccc336 日本淫妇阴液 av海量资源999 凤凰影视成仁 辰溪四中艳照门照片 先锋模特裸体展示影片 成人片免费看 自拍百度云 肥白老妇女 女爱人体图片 妈妈一女穴 星野美夏 日本少女dachidu 妹子私处人体图片 yinmindahuitang 舔无毛逼影片快播 田莹疑的裸体照片 三级电影影音先锋02222 妻子被外国老头操 观月雏乃泥鳅 韩国成人偷拍自拍图片 强奸5一9岁幼女小说 汤姆影院av图片 妹妹人艺体图 美女大驱 和女友做爱图片自拍p 绫川まどか在线先锋 那么嫩的逼很少见了 小女孩做爱 处女好逼连连看图图 性感美女在家做爱 近距离抽插骚逼逼 黑屌肏金毛屄 日韩av美少女 看喝尿尿小姐日逼色色色网图片 欧美肛交新视频 美女吃逼逼 av30线上免费 伊人在线三级经典 新视觉影院t6090影院 最新淫色电影网址 天龙影院远古手机版 搞老太影院 插进美女的大屁股里 私人影院加盟费用 www258dd 求一部电影里面有一个二猛哥 深肛交 日本萌妹子人体艺术写真图片 插入屄眼 美女的木奶 中文字幕黄色网址影视先锋 九号女神裸 和骚人妻偷情 和潘晓婷做爱 国模大尺度蜜桃 欧美大逼50p 西西人体成人 李宗瑞继母做爱原图物处理 nianhuawang 男鸡巴的视屏 � 97免费色伦电影 好色网成人 大姨子先锋 淫荡巨乳美女教师妈妈 性nuexiaoshuo WWW36YYYCOM 长春继续给力进屋就操小女儿套干破内射对白淫荡 农夫激情社区 日韩无码bt 欧美美女手掰嫩穴图片 日本援交偷拍自拍 入侵者日本在线播放 亚洲白虎偷拍自拍 常州高见泽日屄 寂寞少妇自卫视频 人体露逼图片 多毛外国老太 变态乱轮手机在线 淫荡妈妈和儿子操逼 伦理片大奶少女 看片神器最新登入地址sqvheqi345com账号群 麻美学姐无头 圣诞老人射小妞和强奸小妞动话片 亚洲AV女老师 先锋影音欧美成人资源 33344iucoom zV天堂电影网 宾馆美女打炮视频 色五月丁香五月magnet 嫂子淫乱小说 张歆艺的老公 吃奶男人视频在线播放 欧美色图男女乱伦 avtt2014ccvom 性插色欲香影院 青青草撸死你青青草 99热久久第一时间 激情套图卡通动漫 幼女裸聊做爱口交 日本女人被强奸乱伦 草榴社区快播 2kkk正在播放兽骑 啊不要人家小穴都湿了 www猎奇影视 A片www245vvcomwwwchnrwhmhzcn 搜索宜春院av wwwsee78co 逼奶鸡巴插 好吊日AV在线视频19gancom 熟女伦乱图片小说 日本免费av无码片在线开苞 鲁大妈撸到爆 裸聊官网 德国熟女xxx 新不夜城论坛首页手机 女虐男网址 男女做爱视频华为网盘 激情午夜天亚洲色图 内裤哥mangent 吉沢明歩制服丝袜WWWHHH710COM 屌逼在线试看 人体艺体阿娇艳照 推荐一个可以免费看片的网站如果被QQ拦截请复制链接在其它浏览器打开xxxyyy5comintr2a2cb551573a2b2e 欧美360精品粉红鲍鱼 教师调教第一页 聚美屋精品图 中韩淫乱群交 俄罗斯撸撸片 把鸡巴插进小姨子的阴道 干干AV成人网 aolasoohpnbcn www84ytom 高清大量潮喷www27dyycom 宝贝开心成人 freefronvideos人母 嫩穴成人网gggg29com 逼着舅妈给我口交肛交彩漫画 欧美色色aV88wwwgangguanscom 老太太操逼自拍视频 777亚洲手机在线播放 有没有夫妻3p小说 色列漫画淫女 午间色站导航 欧美成人处女色大图 童颜巨乳亚洲综合 桃色性欲草 色眯眯射逼 无码中文字幕塞外青楼这是一个 狂日美女老师人妻 爱碰网官网 亚洲图片雅蠛蝶 快播35怎么搜片 2000XXXX电影 新谷露性家庭影院 深深候dvd播放 幼齿用英语怎么说 不雅伦理无需播放器 国外淫荡图片 国外网站幼幼嫩网址 成年人就去色色视频快播 我鲁日日鲁老老老我爱 caoshaonvbi 人体艺术avav 性感性色导航 韩国黄色哥来嫖网站 成人网站美逼 淫荡熟妇自拍 欧美色惰图片 北京空姐透明照 狼堡免费av视频 www776eom 亚洲无码av欧美天堂网男人天堂 欧美激情爆操 a片kk266co 色尼姑成人极速在线视频 国语家庭系列 蒋雯雯 越南伦理 色CC伦理影院手机版 99jbbcom 大鸡巴舅妈 国产偷拍自拍淫荡对话视频 少妇春梦射精 开心激动网 自拍偷牌成人 色桃隐 撸狗网性交视频 淫荡的三位老师 伦理电影wwwqiuxia6commqiuxia6com 怡春院分站 丝袜超短裙露脸迅雷下载 色制服电影院 97超碰好吊色男人 yy6080理论在线宅男日韩福利大全 大嫂丝袜 500人群交手机在线 5sav 偷拍熟女吧 口述我和妹妹的欲望 50p电脑版 wwwavtttcon 3p3com 伦理无码片在线看 欧美成人电影图片岛国性爱伦理电影 先锋影音AV成人欧美 我爱好色 淫电影网 WWW19MMCOM 玛丽罗斯3d同人动画h在线看 动漫女孩裸体 超级丝袜美腿乱伦 1919gogo欣赏 大色逼淫色 www就是撸 激情文学网好骚 A级黄片免费 xedd5com 国内的b是黑的 快播美国成年人片黄 av高跟丝袜视频 上原保奈美巨乳女教师在线观看 校园春色都市激情fefegancom 偷窥自拍XXOO 搜索看马操美女 人本女优视频 日日吧淫淫 人妻巨乳影院 美国女子性爱学校 大肥屁股重口味 啪啪啪啊啊啊不要 操碰 japanfreevideoshome国产 亚州淫荡老熟女人体 伦奸毛片免费在线看 天天影视se 樱桃做爱视频 亚卅av在线视频 x奸小说下载 亚洲色图图片在线 217av天堂网 东方在线撸撸-百度 幼幼丝袜集 灰姑娘的姐姐 青青草在线视频观看对华 86papa路con 亚洲1AV 综合图片2区亚洲 美国美女大逼电影 010插插av成人网站 www色comwww821kxwcom 播乐子成人网免费视频在线观看 大炮撸在线影院 ,www4KkKcom 野花鲁最近30部 wwwCC213wapwww2233ww2download 三客优最新地址 母亲让儿子爽的无码视频 全国黄色片子 欧美色图美国十次 超碰在线直播 性感妖娆操 亚洲肉感熟女色图 a片A毛片管看视频 8vaa褋芯屑 333kk 川岛和津实视频 在线母子乱伦对白 妹妹肥逼五月 亚洲美女自拍 老婆在我面前小说 韩国空姐堪比情趣内衣 干小姐综合 淫妻色五月 添骚穴 WM62COM 23456影视播放器 成人午夜剧场 尼姑福利网 AV区亚洲AV欧美AV512qucomwwwc5508com 经典欧美骚妇 震动棒露出 日韩丝袜美臀巨乳在线 av无限吧看 就去干少妇 色艺无间正面是哪集 校园春色我和老师做爱 漫画夜色 天海丽白色吊带 黄色淫荡性虐小说 午夜高清播放器 文20岁女性荫道口图片 热国产热无码热有码 2015小明发布看看算你色 百度云播影视 美女肏屄屄乱轮小说 家族舔阴AV影片 邪恶在线av有码 父女之交 关于处女破处的三级片 极品护士91在线 欧美虐待女人视频的网站 享受老太太的丝袜 aaazhibuo 8dfvodcom成人 真实自拍足交 群交男女猛插逼 妓女爱爱动态 lin35com是什么网站 abp159 亚洲色图偷拍自拍乱伦熟女抠逼自慰 朝国三级篇 淫三国幻想 免费的av小电影网站 日本阿v视频免费按摩师 av750c0m 黄色片操一下 巨乳少女车震在线观看 操逼 免费 囗述情感一乱伦岳母和女婿 WWW_FAMITSU_COM 偷拍中国少妇在公车被操视频 花也真衣论理电影 大鸡鸡插p洞 新片欧美十八岁美少 进击的巨人神thunderftp 西方美女15p 深圳哪里易找到老女人玩视频 在线成人有声小说 365rrr 女尿图片 我和淫荡的小姨做爱 � 做爱技术体照 淫妇性爱 大学生私拍b 第四射狠狠射小说 色中色成人av社区 和小姨子乱伦肛交 wwwppp62com 俄罗斯巨乳人体艺术 骚逼阿娇 汤芳人体图片大胆 大胆人体艺术bb私处 性感大胸骚货 哪个网站幼女的片多 日本美女本子把 色 五月天 婷婷 快播 美女 美穴艺术 色百合电影导航 大鸡巴用力 孙悟空操美少女战士 狠狠撸美女手掰穴图片 古代女子与兽类交 沙耶香套图 激情成人网区 暴风影音av播放 动漫女孩怎么插第3个 mmmpp44 黑木麻衣无码ed2k 淫荡学姐少妇 乱伦操少女屄 高中性爱故事 骚妹妹爱爱图网 韩国模特剪长发 大鸡巴把我逼日了 中国张柏芝做爱片中国张柏芝做爱片中国张柏芝做爱片中国张柏芝做爱片中国张柏芝做爱片 大胆女人下体艺术图片 789sss 影音先锋在线国内情侣野外性事自拍普通话对白 群撸图库 闪现君打阿乐 ady 小说 插入表妹嫩穴小说 推荐成人资源 网络播放器 成人台 149大胆人体艺术 大屌图片 骚美女成人av 春暖花开春色性吧 女亭婷五月 我上了同桌的姐姐 恋夜秀场主播自慰视频 yzppp 屄茎 操屄女图 美女鲍鱼大特写 淫乱的日本人妻山口玲子 偷拍射精图 性感美女人体艺木图片 种马小说完本 免费电影院 骑士福利导航导航网站 骚老婆足交 国产性爱一级电影 欧美免费成人花花性都 欧美大肥妞性爱视频 家庭乱伦网站快播 偷拍自拍国产毛片 金发美女也用大吊来开包 缔D杏那 yentiyishu人体艺术ytys WWWUUKKMCOM 女人露奶 � 苍井空露逼 老荡妇高跟丝袜足交 偷偷和女友的朋友做爱迅雷 做爱七十二尺 朱丹人体合成 麻腾由纪妃 帅哥撸播种子图 鸡巴插逼动态图片 羙国十次啦中文 WWW137AVCOM 神斗片欧美版华语 有气质女人人休艺术 由美老师放屁电影 欧美女人肉肏图片 白虎种子快播 国产自拍90后女孩 美女在床上疯狂嫩b 饭岛爱最后之作 幼幼强奸摸奶 色97成人动漫 两性性爱打鸡巴插逼 新视觉影院4080青苹果影院 嗯好爽插死我了 阴口艺术照 李宗瑞电影qvod38 爆操舅母 亚洲色图七七影院 被大鸡巴操菊花 怡红院肿么了 成人极品影院删除 欧美性爱大图色图强奸乱 欧美女子与狗随便性交 苍井空的bt种子无码 熟女乱伦长篇小说 大色虫 兽交幼女影音先锋播放 44aad be0ca93900121f9b 先锋天耗ばさ无码 欧毛毛女三级黄色片图 干女人黑木耳照 日本美女少妇嫩逼人体艺术 sesechangchang 色屄屄网 久久撸app下载 色图色噜 美女鸡巴大奶 好吊日在线视频在线观看 透明丝袜脚偷拍自拍 中山怡红院菜单 wcwwwcom下载 骑嫂子 亚洲大色妣 成人故事365ahnet 丝袜家庭教mp4 幼交肛交 妹妹撸撸大妈 日本毛爽 caoprom超碰在email 关于中国古代偷窥的黄片 第一会所老熟女下载 wwwhuangsecome 狼人干综合新地址HD播放 变态儿子强奸乱伦图 强奸电影名字 2wwwer37com 日本毛片基地一亚洲AVmzddcxcn 暗黑圣经仙桃影院 37tpcocn 持月真由xfplay 好吊日在线视频三级网 我爱背入李丽珍 电影师傅床戏在线观看 96插妹妹sexsex88com 豪放家庭在线播放 桃花宝典极夜著豆瓜网 安卓系统播放神器 美美网丝袜诱惑 人人干全免费视频xulawyercn av无插件一本道 全国色五月 操逼电影小说网 good在线wwwyuyuelvcom www18avmmd 撸波波影视无插件 伊人幼女成人电影 会看射的图片 小明插看看 全裸美女扒开粉嫩b 国人自拍性交网站 萝莉白丝足交本子 七草ちとせ巨乳视频 摇摇晃晃的成人电影 兰桂坊成社人区小说www68kqcom 舔阴论坛 久撸客一撸客色国内外成人激情在线 明星门 欧美大胆嫩肉穴爽大片 www牛逼插 性吧星云 少妇性奴的屁眼 人体艺术大胆mscbaidu1imgcn 最新久久色色成人版 l女同在线 小泽玛利亚高潮图片搜索 女性裸b图 肛交bt种子 最热门有声小说 人间添春色 春色猜谜字 樱井莉亚钢管舞视频 小泽玛利亚直美6p 能用的h网 还能看的h网 bl动漫h网 开心五月激 东京热401 男色女色第四色酒色网 怎么下载黄色小说 黄色小说小栽 和谐图城 乐乐影院 色哥导航 特色导航 依依社区 爱窝窝在线 色狼谷成人 91porn 包要你射电影 色色3A丝袜 丝袜妹妹淫网 爱色导航(荐) 好男人激情影院 坏哥哥 第七色 色久久 人格分裂 急先锋 撸撸射中文网 第一会所综合社区 91影院老师机 东方成人激情 怼莪影院吹潮 老鸭窝伊人无码不卡无码一本道 av女柳晶电影 91天生爱风流作品 深爱激情小说私房婷婷网 擼奶av 567pao 里番3d一家人野外 上原在线电影 水岛津实透明丝袜 1314酒色 网旧网俺也去 0855影院 在线无码私人影院 搜索 国产自拍 神马dy888午夜伦理达达兔 农民工黄晓婷 日韩裸体黑丝御姐 屈臣氏的燕窝面膜怎么样つぼみ晶エリーの早漏チ○ポ强化合宿 老熟女人性视频 影音先锋 三上悠亚ol 妹妹影院福利片 hhhhhhhhsxo 午夜天堂热的国产 强奸剧场 全裸香蕉视频无码 亚欧伦理视频 秋霞为什么给封了 日本在线视频空天使 日韩成人aⅴ在线 日本日屌日屄导航视频 在线福利视频 日本推油无码av magnet 在线免费视频 樱井梨吮东 日本一本道在线无码DVD 日本性感诱惑美女做爱阴道流水视频 日本一级av 汤姆avtom在线视频 台湾佬中文娱乐线20 阿v播播下载 橙色影院 奴隶少女护士cg视频 汤姆在线影院无码 偷拍宾馆 业面紧急生级访问 色和尚有线 厕所偷拍一族 av女l 公交色狼优酷视频 裸体视频AV 人与兽肉肉网 董美香ol 花井美纱链接 magnet 西瓜影音 亚洲 自拍 日韩女优欧美激情偷拍自拍 亚洲成年人免费视频 荷兰免费成人电影 深喉呕吐XXⅩX 操石榴在线视频 天天色成人免费视频 314hu四虎 涩久免费视频在线观看 成人电影迅雷下载 能看见整个奶子的香蕉影院 水菜丽百度影音 gwaz079百度云 噜死你们资源站 主播走光视频合集迅雷下载 thumbzilla jappen 精品Av 古川伊织star598在线 假面女皇vip在线视频播放 国产自拍迷情校园 啪啪啪公寓漫画 日本阿AV 黄色手机电影 欧美在线Av影院 华裔电击女神91在线 亚洲欧美专区 1日本1000部免费视频 开放90后 波多野结衣 东方 影院av 页面升级紧急访问每天正常更新 4438Xchengeren 老炮色 a k福利电影 色欲影视色天天视频 高老庄aV 259LUXU-683 magnet 手机在线电影 国产区 欧美激情人人操网 国产 偷拍 直播 日韩 国内外激情在线视频网给 站长统计一本道人妻 光棍影院被封 紫竹铃取汁 ftp 狂插空姐嫩 xfplay 丈夫面前 穿靴子伪街 XXOO视频在线免费 大香蕉道久在线播放 电棒漏电嗨过头 充气娃能看下毛和洞吗 夫妻牲交 福利云点墦 yukun瑟妃 疯狂交换女友 国产自拍26页 腐女资源 百度云 日本DVD高清无码视频 偷拍,自拍AV伦理电影 A片小视频福利站。 大奶肥婆自拍偷拍图片 交配伊甸园 超碰在线视频自拍偷拍国产 小热巴91大神 rctd 045 类似于A片 超美大奶大学生美女直播被男友操 男友问 你的衣服怎么脱掉的 亚洲女与黑人群交视频一 在线黄涩 木内美保步兵番号 鸡巴插入欧美美女的b舒服 激情在线国产自拍日韩欧美 国语福利小视频在线观看 作爱小视颍 潮喷合集丝袜无码mp4 做爱的无码高清视频 牛牛精品 伊aⅤ在线观看 savk12 哥哥搞在线播放 在线电一本道影 一级谍片 250pp亚洲情艺中心,88 欧美一本道九色在线一 wwwseavbacom色av吧 cos美女在线 欧美17,18ⅹⅹⅹ视频 自拍嫩逼 小电影在线观看网站 筱田优 贼 水电工 5358x视频 日本69式视频有码 b雪福利导航 韩国女主播19tvclub在线 操逼清晰视频 丝袜美女国产视频网址导航 水菜丽颜射房间 台湾妹中文娱乐网 风吟岛视频 口交 伦理 日本熟妇色五十路免费视频 A级片互舔 川村真矢Av在线观看 亚洲日韩av 色和尚国产自拍 sea8 mp4 aV天堂2018手机在线 免费版国产偷拍a在线播放 狠狠 婷婷 丁香 小视频福利在线观看平台 思妍白衣小仙女被邻居强上 萝莉自拍有水 4484新视觉 永久发布页 977成人影视在线观看 小清新影院在线观 小鸟酱后丝后入百度云 旋风魅影四级 香蕉影院小黄片免费看 性爱直播磁力链接 小骚逼第一色影院 性交流的视频 小雪小视频bd 小视频TV禁看视频 迷奸AV在线看 nba直播 任你在干线 汤姆影院在线视频国产 624u在线播放 成人 一级a做爰片就在线看狐狸视频 小香蕉AV视频 www182、com 腿模简小育 学生做爱视频 秘密搜查官 快播 成人福利网午夜 一级黄色夫妻录像片 直接看的gav久久播放器 国产自拍400首页 sm老爹影院 谁知道隔壁老王网址在线 综合网 123西瓜影音 米奇丁香 人人澡人人漠大学生 色久悠 夜色视频你今天寂寞了吗? 菲菲影视城美国 被抄的影院 变态另类 欧美 成人 国产偷拍自拍在线小说 不用下载安装就能看的吃男人鸡巴视频 插屄视频 大贯杏里播放 wwwhhh50 233若菜奈央 伦理片天海翼秘密搜查官 大香蕉在线万色屋视频 那种漫画小说你懂的 祥仔电影合集一区 那里可以看澳门皇冠酒店a片 色自啪 亚洲aV电影天堂 谷露影院ar toupaizaixian sexbj。com 毕业生 zaixian mianfei 朝桐光视频 成人短视频在线直接观看 陈美霖 沈阳音乐学院 导航女 www26yjjcom 1大尺度视频 开平虐女视频 菅野雪松协和影视在线视频 华人play在线视频bbb 鸡吧操屄视频 多啪啪免费视频 悠草影院 金兰策划网 (969) 橘佑金短视频 国内一极刺激自拍片 日本制服番号大全magnet 成人动漫母系 电脑怎么清理内存 黄色福利1000 dy88午夜 偷拍中学生洗澡磁力链接 花椒相机福利美女视频 站长推荐磁力下载 mp4 三洞轮流插视频 玉兔miki热舞视频 夜生活小视频 爆乳人妖小视频 国内网红主播自拍福利迅雷下载 不用app的裸裸体美女操逼视频 变态SM影片在线观看 草溜影院元气吧 - 百度 - 百度 波推全套视频 国产双飞集合ftp 日本在线AV网 笔国毛片 神马影院女主播是我的邻居 影音资源 激情乱伦电影 799pao 亚洲第一色第一影院 av视频大香蕉 老梁故事汇希斯莱杰 水中人体磁力链接 下载 大香蕉黄片免费看 济南谭崔 避开屏蔽的岛a片 草破福利 要看大鸡巴操小骚逼的人的视频 黑丝少妇影音先锋 欧美巨乳熟女磁力链接 美国黄网站色大全 伦蕉在线久播 极品女厕沟 激情五月bd韩国电影 混血美女自摸和男友激情啪啪自拍诱人呻吟福利视频 人人摸人人妻做人人看 44kknn 娸娸原网 伊人欧美 恋夜影院视频列表安卓青青 57k影院 如果电话亭 avi 插爆骚女精品自拍 青青草在线免费视频1769TV 令人惹火的邻家美眉 影音先锋 真人妹子被捅动态图 男人女人做完爱视频15 表姐合租两人共处一室晚上她竟爬上了我的床 性爱教学视频 北条麻妃bd在线播放版 国产老师和师生 magnet wwwcctv1024 女神自慰 ftp 女同性恋做激情视频 欧美大胆露阴视频 欧美无码影视 好女色在线观看 后入肥臀18p 百度影视屏福利 厕所超碰视频 强奸mp magnet 欧美妹aⅴ免费线上看 2016年妞干网视频 5手机在线福利 超在线最视频 800av:cOm magnet 欧美性爱免播放器在线播放 91大款肥汤的性感美乳90后邻家美眉趴着窗台后入啪啪 秋霞日本毛片网站 cheng ren 在线视频 上原亚衣肛门无码解禁影音先锋 美脚家庭教师在线播放 尤酷伦理片 熟女性生活视频在线观看 欧美av在线播放喷潮 194avav 凤凰AV成人 - 百度 kbb9999 AV片AV在线AV无码 爱爱视频高清免费观看 黄色男女操b视频 观看 18AV清纯视频在线播放平台 成人性爱视频久久操 女性真人生殖系统双性人视频 下身插入b射精视频 明星潜规测视频 mp4 免賛a片直播绪 国内 自己 偷拍 在线 国内真实偷拍 手机在线 国产主播户外勾在线 三桥杏奈高清无码迅雷下载 2五福电影院凸凹频频 男主拿鱼打女主,高宝宝 色哥午夜影院 川村まや痴汉 草溜影院费全过程免费 淫小弟影院在线视频 laohantuiche 啪啪啪喷潮XXOO视频 青娱乐成人国产 蓝沢润 一本道 亚洲青涩中文欧美 神马影院线理论 米娅卡莉法的av 在线福利65535 欧美粉色在线 欧美性受群交视频1在线播放 极品喷奶熟妇在线播放 变态另类无码福利影院92 天津小姐被偷拍 磁力下载 台湾三级电髟全部 丝袜美腿偷拍自拍 偷拍女生性行为图 妻子的乱伦 白虎少妇 肏婶骚屄 外国大妈会阴照片 美少女操屄图片 妹妹自慰11p 操老熟女的b 361美女人体 360电影院樱桃 爱色妹妹亚洲色图 性交卖淫姿势高清图片一级 欧美一黑对二白 大色网无毛一线天 射小妹网站 寂寞穴 西西人体模特苍井空 操的大白逼吧 骚穴让我操 拉好友干女朋友3p