Species and Functional Diversity of Deep-Sea Nematodes in a High Energy Submarine Canyon

Introduction

Gaoping Submarine Canyon (GPSC) is the largest underwater canyon off SW Taiwan (Yu et al., 2009). The head of GPSC connects to a small mountain river, the Gaoping River, which originates from the Mt. Jade at more than 3,000 m above sea level. The thalweg of the canyon meanders approximately 260 km to the south and eventually connects to the Manila Trench, which reaches depths of over 5,000 m. The high bedrock erodibility, heavy precipitations, steep topography and high tectonic activities in the Gaoping River catchment make the GPSC prone to receive extremely high fluvial sediment loads (∼10 kt km^–2 y^–1) and frequent turbidity currents (Liu et al., 2013, 2016; Gavey et al., 2017). The rugged, steep, and tectonically active margin off SW Taiwan also subjects the GPSC to high risk of slope failure, producing debris flows which cascade down the canyon axis, occasionally resulting in breakage of telecommunication cables (Hsu et al., 2008; Su et al., 2012; Gavey et al., 2017). Besides these episodic events, the GPSC is close to the source of the world’s largest internal solitary waves in the Luzon Strait (Jan et al., 2008; Alford et al., 2015). The energy of internal tides propagates along the axis of the GPSC following tidal cycles to drive bottom and headward intensified currents with observed current velocity exceeding 1 m/s at the canyon head (Wang et al., 2008; Chiou et al., 2011). These unique physical settings make the GPSC an example of source-to-sink sediment transports from the high mountains to the deep sea and also an ideal location to study submarine geohazards (Liu et al., 2013, 2016). However, despite a wealth of literature and interests in the GPSC’s unique geology, information is lacking on how benthic communities, especially the very small but numerically dominant meiofauna, respond to such extreme environmental conditions. This study thus attempts to conduct the first meiofauna diversity investigation in the GPSC and to contribute to the growing interests in understanding the ecological processes in the benthic communities in highly dynamic conditions.

Many studies have documented enriched benthic communities in submarine canyons and suggested that the organic accumulation in these prominent topographic depressions might be the primary driver of enhanced benthic production (Baguley et al., 2006; Ingels et al., 2009; De Leo et al., 2010; Kiriakoulakis et al., 2011; Leduc et al., 2014; Amaro et al., 2016; Román et al., 2016; Gambi et al., 2019). Due to the higher food availability within canyons, the abundance and biomass of meiofauna are generally higher than on the adjacent slope (Danovaro et al., 2009; Ingels et al., 2009, 2011; De Leo et al., 2010; Leduc et al., 2014; Gambi et al., 2019), although the opposite has also been found in active canyons with substantial sediment disturbances (Garcia et al., 2007; Van Gaever et al., 2009). The elevated density and biomass of benthos in the submarine canyon are thus considered as the typical “canyon effect” (Vetter and Dayton, 1998; Fernandez-Arcaya et al., 2017; Román et al., 2019). Also, opportunistic species may become dominant in the canyons so the local diversity could be comparable to or lower than in slope communities (Van Gaever et al., 2009; Ingels et al., 2011; Leduc et al., 2014; Gambi et al., 2019; Wei and Rowe, 2019; Wei et al., 2020). In the upper GPSC off SW Taiwan, however, Liao et al. (2017) found that strong, continuous, and recurrent bottom currents driven by the internal tides may cause sediment erosion which may negatively impact the densities and alter the taxonomic composition of meiofaunal taxa. The faunal patterns observed in the GPSC probably resemble those observed under high physical disturbance by strong near-bed flows and sediment erosion, as in the upper Nazaré Canyon and Congo Channel (Garcia et al., 2007; Van Gaever et al., 2009). The strong bottom shear not only may erode sediment and remove the local benthic communities but also likely re-suspend or prevent the organic-rich particles from settling on the seafloor and thus reduce the sedimentary organic carbon contents and the food supply to the meiofauna (Liao et al., 2017). While meiofauna abundance and the number of taxa were found to be depressed in the canyon (Liao et al., 2017), it is still unknown how species diversity, species composition, and functional traits may be affected. These biodiversity attributes are known to predict the ecosystem functioning of the benthic community, reflecting the ability of meiobenthos in moderating carbon and nutrient cycling and converting the available energy into biomass and transferring this energy up the food chain (Covich et al., 2004; Ieno et al., 2006; Radwell and Brown, 2008; Karlson et al., 2010).

Moreover, Liao et al. (2017) hypothesized that the extreme conditions in the GPSC might result in the environmental filtering of the less-fitted species or even local extinction if the disturbance is too strong. In such a case, meiofauna should exhibit varying degrees of spatial turnover (i.e., replaced by better-fitted species) or nestedness pattern (i.e., species-poor as subsets of species-rich assemblages) (Baselga, 2010, 2012). From the relative proportion of the turnover and nestedness inferences can be made regarding the process of community assembling under strong environmental perturbations. Therefore, the high-resolution, species-level data are required for quantitatively assessing these ecological mechanisms.

In this study, we used free-living marine nematodes as a model organism to examine the species and functional diversity of meiofauna in the GPSC and adjacent slope habitat. Nematodes are the most abundant metazoan in sediments, accounting for more than 90% of the meiofaunal abundance in deep-sea environments (Giere, 2009; Danovaro, 2012). They are also one of the most diverse organisms in the ocean. A conservative estimation suggests that there may be 6,900 described species to up to 50,000 species if the undescribed species are also included (Appeltans et al., 2012). Not only are nematodes abundant and diverse, they are also ubiquitous, not highly mobile, and sensitive to environmental changes and pollution; therefore, nematodes have been widely used as an indicator for environmental monitoring (Coull and Chandler, 1992; Bongers and Ferris, 1999; Schratzberger et al., 2000; Giere, 2009; Moreno et al., 2011). Due to their small sizes, short generation times, and high metabolic rates, nematodes process a considerable amount of energy in the benthic ecosystem (Giere, 2009). The grazing, excretion, bioturbation, and mechanical breakdown of detrital particles by nematodes can stimulate bacterial activities (Giere, 2009). These interactions may affect the organic carbon remineralization, oxygen dynamics, and nutrient cycling in the sediment (Aller and Aller, 1992; Rysgaard et al., 2000; Moens et al., 2005; Nascimento et al., 2012; Bonaglia et al., 2014; Schratzberger and Ingels, 2018). More importantly, the feeding guilds and morphological traits affecting the locomotion and life-history strategies of nematodes are well-developed (Wieser, 1953; Jensen, 1987; Bongers et al., 1991; Thistle et al., 1995; Soetaert et al., 2002; Schratzberger et al., 2007). These studies allow us to quickly assemble the species functional traits from the known species list and compute nematode functional diversity.

Our objective was to examine how species and functional diversity and community composition of deep-sea nematodes may respond to a high energy environment in the GPSC. We compared the GPSC and the adjacent slope by surveying parallel transects at equivalent depths across the continental margin off SW Taiwan. We hypothesized that: (1) the extreme physical conditions in the GPSC may depress the nematode species, functional and trophic diversity, and maturity, as well as alter their life strategies (e.g., feeding guilds, morphology, and lifestyles) and species composition; (2) nestedness patterns of species loss may be the dominant mechanism shaping the GPSC community if the adverse environment caused local extinction; otherwise, spatial turnover of species replacement may play a more dominant role indicating the effect of environmental filtering. Ultimately, with the species-level information on nematodes, we hope to contribute to the understanding of the responses of meiofauna to the strong internal tides and thus fast bottom currents and consequent physical disturbances in a high-energy submarine canyon.

Materials and Methods

Environmental Setting

GPSC is a major conduit on the active margin offshore SW Taiwan. The collision between the Luzon Volcanic Arc and the South China Sea margin created the Taiwan central mountain belt, foreland shelf and slope, and Manila Trench, which significantly affects the geomorphology of GPSC (Liu et al., 2016). The head of GPSC is deeply incised into the continental shelf (i.e., vertical relief >400 meters, Figure 1) with a clear bathymetric connection to a small mountain river, the Gaoping River (Chiang et al., 2020). The canyon head is believed to have formed by subaerial erosion; the head later submerged during the Holocene transgression following the last glacial maximum (∼18,000 years BP) and then was re-excavated by sediment movement processes (Yu et al., 1991). The meandering, V-shaped, and entrenched thalweg with deep-cutting outer bends are characteristic features of the GPSC (Figure 1). These erosive features are believed to be maintained by turbidity currents triggered by (1) flood events (Liu et al., 2013, 2016), (2) turbulence mixing related to internal and surface waves (Wang et al., 2008; Lee et al., 2009; Chiou et al., 2011), (3) submarine groundwater discharges (Su et al., 2012), and (4) sediment collapses from over- steepening canyon walls following earthquakes (Hsu et al., 2008; Su et al., 2012; Gavey et al., 2017). The GPSC is filled with the effluent of Gaoping River, Kuroshio Current Water, and South China Sea Water, in which the mixing of these water masses is controlled by internal tides (Liu et al., 2002, 2006). The internal tides in the GPSC follow semidiurnal cycles (M₂), in which the intensity increases toward the canyon head, or in other words, decreases with depth and with distance from the canyon head (Wang et al., 2008). The energy flux of internal tide in the GPSC is estimated to be 3–7 times greater than that in the Monterey Canyon (Lee et al., 2009); therefore, the isothermal displacement by the internal tides can be easily over 200 meters in the GPSC (Liu et al., 2016). Another consequence of the strong internal waves is the year-round existence of a benthic nepheloid layer as thick as 100 m with the suspended sediment concentration reaching 30 mg/l (Liu et al., 2010). The turbulence mixing associated with internal tides and high suspended-sediment concentrations near the seabed are believed to result in significant sediment transport (Liu et al., 2016).

Figure 1. Map of sampling stations visited in August and November 2015. (A) The sampling area off the continental margin of SW Taiwan. (B) Sampling transects along the upper Gaoping Submarine Canyon (GC) and Gaoping Slope (GS).

Field Sampling

Two cruises (cruises 1114 and 1126) were conducted off SW Taiwan by R/V Ocean Researcher 1 in August and November 2015 (i.e., summer and fall, respectively). A total of eight stations were repeatedly sampled from 200 to 1,100 m in the upper Gaoping Submarine Canyon (GPSC) and the adjacent slope (Figure 1 and Table 1). The sampling was designed to collect biological and environmental information on four equally spaced depth strata (200–400 m, 400–600 m, 600–800 m, and 800–1,100 m) along the canyon thalweg and on the adjacent slope. At each station, a CTD-Rosette cast was made and surface sediment was sampled by means of a multicorer (OSIL megacorer). One multicore tube (i.d., =105 mm) per station was selected for meiofaunal analysis and three subsamples were taken with a cut-off syringe (i.d., 28 mm; area 616 mm²) from the top 5 cm of the sediment. The subsamples were immediately fixed on board in formalin solution with Rose Bengal. During Cruise 1126, only one subsample for meiofauna analysis was retrieved from the deepest canyon station (GC4). Altogether, 46 meiofauna subsamples were retrieved during the two cruises.

Table 1. Dates, coordinates, and water depths of sample sites during August and November 2015.

The fixed sediment samples were washed with tap water through a 1,000-μm sieve and a 40-μm sieve. The meiofauna samples retained on a 40-μm sieve were extracted by Ludox HS40 (gravity = 1.18 g/cm³; centrifuged for 10 min at 8,000 rpm and repeated three times) (Danovaro, 2010), after which meiofauna were picked and counted under a stereomicroscope (Olympus SZ61). From each subsample, 100 nematode individuals (or all the individuals if fewer than 100 were present) were randomly picked out and transferred to a solution of 5% glycerol and 5% ethanol in water. The mixture was allowed to evaporate gradually on a warm hotplate, leaving the nematodes in pure glycerol. The nematodes were mounted onto permanent slides. Putative morphospecies were identified using pictorial keys (Platt and Warwick, 1983, 1988; Warwick et al., 1998; Schmidt-Rhaesa, 2014). Based on the morphology of the buccal cavities, the nematodes were classified to four feeding traits, namely selective deposit feeders (1A), non-selective deposit feeders (1B), epigrowth feeders (2A), and omnivores/predators (2B) (Wieser, 1953). Of the gutless chemotrophic nematode genus Astomonema, allocated to a fifth feeding trait by Ingels et al. (2011), only one specimen was identified, which was not included in the analysis. Furthermore, nematodes were classified to four classes of tail shape, namely short/round, elongated/filiform, conical, and clavate, because these traits might relate to the movement and vertical distribution of nematodes, which depends on the sediment type (Riemann, 1974; Thistle et al., 1995). Also, following Bongers et al. (1991), nematodes were assigned to four life-history strategy classes, using c-p values ranging from 1 for colonizers (r-selected species with short generation time, high colonization ability, and adaptation of unstable environments) to 4 for persisters (K-selected species with strong competition ability in crowded niches).

Temperature, salinity, dissolved oxygen concentration, fluorescence, and light transmission of the bottom water were measured by a conductivity-temperature-depth (CTD) recorder (Sea-Bird SBE 911) and other attached sensors. Surface sediment grain sizes (non-carbonate fraction) were measured by a laser diffraction particle size analyzer (Beckman Coulter LS13 320). Total organic carbon (TOC) and total nitrogen (TN) were analyzed with a Flash EA 1112 elemental analyzer. The detailed descriptions of hydrological and sediment sampling can be found in Liao et al. (2017). Hourly bottom current velocity at each site was derived from a 3-D, hydrodynamic internal tide model from Chiou et al. (2011). Based on the internal tide model, we calculated the hourly mean velocity of bottom currents and the duration for which bottom current velocity exceeded 20 cm/s for 1 month preceding the sampling campaign to evaluate the possible disturbances of near-bottom currents on nematodes. The environmental data were derived from Liao et al. (2017, Supplementary Figure S1), in which the canyon head and upper canyon had the highest modeled bottom current velocity and the most prolonged duration of sediment erosion. The lowest bottom water light transmission and surface sediment TOC coincided with the strongest hydrodynamic energy near the canyon head. While the temperature and dissolved oxygen concentration both declined with depth, surface sediment TOC increased with depths in both the canyon and slope transects (Liao et al., 2017). Nevertheless, the dissolved oxygen concentration never dropped below the threshold of hypoxia (2 mg/L), presumably, due to turbulence mixing by internal tides. To avoid autocorrelations, redundant environmental variables (correlations > 0.9) were removed. For example, the bottom water temperature, density, and dissolved oxygen measured in this study were highly correlated (Supplementary Figure S1). Therefore, the dissolved oxygen was not retained because the bottom water was not hypoxic in all stations. The density was not retained because it’s not an important factor known to affect nematode assemblages. Only the ecologically relevant variables, including bottom water temperature, salinity, light transmission, percent sand, silt and clay, sediment TOC and C/N ratio, mean bottom current velocity, and duration of bottom current velocity exceeding 20 cm/s were retained for the analysis. These variables were logarithm (base of 10) transformed, centered (subtracted from the mean), and normalized (divided by the standard deviation) prior to use in statistical analyses involving environmental variables.

Data Analysis

Species diversity was determined using Hill numbers (Hill, 1973), or the effective numbers of equally abundant species in a hypothetical assemblage (Chao and Jost, 2010). The Hill number, ^qD, is expressed as a function of the q-th power sum of the relative species abundance. Thus, the order q controls the sensitivity of Hill numbers to the relative species abundance. For example, when q = 0, ⁰D weights all species equally and reduces to the number of species. For q = 1, ¹D weights species in proportion to their abundance and can be interpreted as the effective number of abundant species. For q = 2, ²D then measures the effective number of highly abundant species. The Hill numbers of order q = 0, 1, 2 are mathematically equivalent to the species richness, the exponential of Shannon entropy, and the inverse of Simpson concentration index (Chao et al., 2014). We, therefore, adopted these three most commonly used diversity indices and referred to them as the species richness, Shannon diversity, and Simpson diversity following Chao et al. (2014). To make a fair comparison across samples, the Hill numbers of order q = 0, 1, 2 were rarefied or extrapolated to equal sample coverage of 80.8% (the lowest sample coverage among samples, Supplementary Figure S2). The standard error and 95% confidence intervals of the Hill numbers were estimated by 1000 bootstrap resampling. The rarefaction, extrapolation, and resampling of Hill numbers follow the description by Chao et al. (2013).

Functional diversity was calculated from species abundance and functional dissimilarities (Gower distances) between species based on nematode functional traits (buccal morphology, tail shape, and life history). We adopt a Hill-number based approach to estimate the effective number of equally abundant and equally distinct functional groups in an assemblage (Chao et al., 2019). In this new approach, Chao et al. (2019) introduce a positive parameter τ to the ordinary Hill numbers to indicate a threshold of functional distinctiveness between any two species (based on Gower distance in functional trait space). Any two species with functional dissimilarity ≥τ are considered functionally distinct. When τ approaches 0, every species is functionally distinct; therefore, the functional diversity reduces to the ordinary Hill numbers. Following the recommendation by Chao et al. (2019), we set the τ = Rao’s quadratic entropy (Q) of the pooled assemblages, or the abundance-weighted mean functional dissimilarity (based on Gower distance) between any two taxa (Botta-Dukát, 2005). Similarly, we computed the functional diversity of order q = 0, 1, 2 to estimate the functional richness, effective number of abundant functional group, and effective number of highly abundant functional group, respectively. Unfortunately, the functional diversities were not standardized by sample coverage. The rarefied and extrapolated functional diversity is still under development and not available at the time of writing this paper (Chao, personal communication).

Trophic diversity (TD) was determined from four trophic groups based on the buccal morphologies (1A, 1B, 2A, 2B) (Wieser, 1953) and calculated as the sum of square of the relative abundance of each trophic group (Heip et al., 1985; Gambi et al., 2003; Gambi and Danovaro, 2016). Maturity index (MI) was computed as the mean individual genus c-p value (from colonizers to persisters) weighted by their relative abundance (Bongers et al., 1991).

We scaled the relative abundance of each nematode species (from 100 randomly selected individuals) by the total nematode abundance in each sample. The species abundance in each sample was square-root transformed, and converted to Bray-Curtis dissimilarity between samples. The same matrix was subjected to the agglomerative hierarchical clustering based on the group average (=UPGMA) method and Non-metric Multi-dimensional Scaling (nMDS). The species which contributed most to the average Bray-Curtis dissimilarity in each cluster were identified by the Similarity Percentage (SIMPER) routine. These characteristic species were projected on to the nMDS plot by the abundance-weighted averages of the ordination scores. Distance-based Redundancy Analysis (dbRDA) was used to select the subset of environmental variables which best explained (with the smallest AIC) the nematode species composition. The selected variables were projected as vectors on to the same nMDS ordination with the length of vectors indicating their correlations with the nMDS ordination scores and direction of vectors indicating the direction of increasing environmental values.

We used multiple-site Sørensen dissimilarity to examine nematode β-diversity among the four sites within the GPSC and adjacent slope transects, as well as between the canyon and slope site pair across four depth strata (i.e., 200–400 m, 400–600 m, 600–800 m, 800–1,100 m). The dissimilarities within-habitat (canyon or slope) or between-habitat (canyon vs. slope) were respectively partitioned to their turnover and nestedness components (Baselga, 2010, 2012). Average dissimilarity, turnover, and nestedness over the two cruises in 2015 were calculated. We used the turnover dissimilarity to examine the contributions of environmental filtering in structuring the nematode species composition within and between the habitats. The nestedness dissimilarity explored the contribution of source-sink dynamics or local extinction.

Nematode community attributes, such as total abundance, number of species, species diversity, functional diversity, trophic diversity, and maturity index of each subsample were averaged per core (i.e., average of three subsamples) before further statistical tests. We used Generalized Least Squares (GLS) modeling to examine the effects of habitat (canyon vs. slope), depth strata, and sampling time (Cruise 1114 vs. 1126) on the mean values of the community attributes mentioned above. Also, we performed a three-way cross permutation analysis of variance (PERMANOVA) to examine the effects of habitat, depth, and sampling time on the average nematode species composition per core. The number of permutation was set to 999. All statistical tests used α-value = 0.05 and pairwise tests used α-value = 0.05/numbers of tests (i.e., Bonferroni correction).

Statistical analyses used software R version 3.6.1 (R Core Team, 2019). Hill numbers were computed by the “iNEXT” package. Functional diversity was computed by R code available at https://github.com/AnneChao/FunD/blob/master/FunD_Rcode.txt. Multivariate analyses used the “vegan” package and β-diversity partition used the “betapart” package. Generalized Least Squares (GLS) modeling used the “nlme” package. All relevant analyses can be reproduced and the data are deposited in an R data package “nema” available at https://github.com/chihlinwei/nema.

Results

A total of 163 nematode species belonging to 32 families and 105 genera were identified from 3,577 individuals. Among them, 53 and 156 species were found in the GPSC and adjacent slope, respectively. In other words, 46 species co-occurred in both habitats and only seven species were exclusively found in the GPSC. Among the nematode species identified in this study, there was a total of 32 different combinations of functional types (excluding the gutless Astomonema sp.) based on their buccal morphologies, tail shapes, and life-history strategies. Only 15 combinations were found in the GPSC, whereas all functional types existed on the slope. The non-selective deposit feeders (1B) with clavate tail and r-selection strategy (c-p value = 2) flourished in the GPSC (GLS, P < 0.01, Supplementary Table S1). In contrast, the different buccal morphologies, tail shapes, and life-history strategies were more evenly represented in the slope assemblages (Figure 2). The relative contributions of epigrowth feeders (2A), conical and elongated/filiform tails, and persisters (c-p value = 4) were higher on the slope (GLS, P < 0.01, Supplementary Table S1). The relative contributions of deposit feeders (1A, 1B) increased, but the omnivores/predators (2B) decreased with depth (GLS, P < 0.05, Supplementary Table S1 and Figure 2A). In the GPSC, the r-selection strategists contributed more toward the head of the canyon, reflected in a decreasing contribution of persisters (GLS, P < 0.01, Supplementary Table S1). Albeit weakly, on the slope, the relative contributions of persisters declined with depth (GLS, P = 0.026, Supplementary Table S1 and Figure 2C).

Figure 2. The relative abundance of (A) buccal morphology, (B) tail shape, and (C) life history strategy of nematodes of the GPSC and adjacent slope. Buccal cavities classified according to Wieser (1953): selective deposit feeders (1A), non-selective deposit feeders (1B), epigrowth feeders (2A), and omnivores/predators (2B). Tail shapes included clavate (cla), conical (co), elongated/filiform (e/f), and short/round (s/r). Life history denoted as c-p values 1–4 from colonizer to persister.

Nematode abundances were not significantly different between the canyon and slope and among different depth strata (Figure 3 and Table 2). The observed number of species (Figure 3B), Shannon diversity, functional diversity (q = 1), trophic diversity and maturity index were significantly depressed in the GPSC compared to the adjacent slope (Figure 4 and Table 2). Similar patterns (i.e., depressed diversity in the GPSC) were also found for the species richness, Simpson diversity, and functional diversity of order q = 0 and 2 (Supplementary Figures S2, S3 and Supplementary Table S2). Individual linear regressions for the canyon and slope samples show that the species diversity (including species richness, Shannon diversity, and Simpson diversity) both increased with depth; however, the increase was more rapid on the slope than in the canyon (Supplementary Figures S4A–C). On the slope, the number of species increased and trophic diversity declined with depth, but no bathymetric pattern was evident in the canyon (Supplementary Figures S4D,E). These subtle discrepancies in bathymetric patterns between the canyon and slope may contribute to the significant interaction between habitat and depth in the GLS model for the number of species, Shannon diversity, trophic diversity and maturity index (Table 2). Similar interactions were also evident for species richness and Simpson diversity, probably for the same reason (Supplementary Table S2). Nevertheless, the effect size of habitat (in the variables mentioned above) is much stronger than the effect of habitat and depth interactions (see GLS coefficients in Table 2 and Supplementary Table S2), which warrants the observation of depressed nematode diversity (all aspects from species, functional to trophic diversity) and maturity in the GPSC. Among the two sampling times, none of the GLS tests suggest significant cruise effect (P > 0.05, Table 2 and Supplementary Table S2); however, the inconclusive time effect might have been confounded by the significant interaction between habitat and cruise for total abundance and Shannon diversity (P < 0.05, Table 2), as well as for the species richness and Simpson diversity (P < 0.05, Supplementary Table S2). Nevertheless, further pairwise GLS tests (with depth, cruise and their interaction as factors) in each of the canyon and slope habitats confirmed no cruise effect (P > 0.05/2, Bonferroni Correction).

Figure 3. Nematode abundance and number of species along the upper GPSC and adjacent slope. The bar chart and error bars show the mean and standard deviation of the nematode (A) abundance and (B) observed number of species per station during the August and November 2015 sampling.

Table 2. Generalized Least Square (GLS) modeling on mean Shannon diversity, functional diversity (q = 1), trophic diversity, and maturity index of nematodes per station.

Figure 4. Nematode species and functional diversity along the upper GPSC and adjacent slope. The bar chart and error bars show the mean and standard deviation of the nematode (A) species diversity as Hill number of order q = 1 or effective numbers of abundant species, (B) functional diversity of order q = 1 or effective numbers of abundant functional group, (C) trophic diversity, and (D) maturity index per station during the August and November 2015 sampling.

Habitats were clearly separated in the nMDS ordination of nematode assemblages without any overlap (Figure 5). This was also reflected by the significant difference found between nematode species composition of canyon and slope habitats (habitat, P = 0.001) in the PERMANOVA (Table 3). In the meantime, no statistical depth effect (P = 0.135) but marginal cruise effect (P = 0.052) and marginal interaction between habitat and depth were detected (P = 0.061, Table 3). Although not statistically significant, the temporal change in nematode composition is visible in the nMDS ordination (Figure 5), especially for the canyon (i.e., better separation between cruises). Based on the similarity percentage (SIMPER) analysis, Metadesmolaimus sp. 1, Daptonema sp. 4, Axonolaimus sp. 1, Sabatieria sp. 2, and Sabatieria sp. 3 have the highest contribution to the average dissimilarity in the GPSC (Figure 5A and Supplementary Table S3). Sabatieria sp. 1, Elzalia gerlachi, Dorylaimopsis variabilis, Cervonema tenuicaudatum, and Sphaerolaimus sp. 1 were the top five most contributing species on the slope (Supplementary Table S3). In terms of species abundance, Daptonema sp. 2, Axonolaimus sp. 1, Metadesmolaimus sp. 1, Daptonema sp. 4, and Sabatieria sp. 2 dominated in the GPSC, whereas Craspodema sp. 1, Elzalia gerlachi, Paramesacanthion tricuspis, Setosabatieria hilarula, and Sabatieria sp. 1 flourished on the slope. For the most contributing species identified by SIMPER (Supplementary Tables S3, S4), the top five species in the GPSC were non-selective deposit feeders (1B) with clavate tail and r-selection strategy, whereas the highly contributed species on the slope including Dorylaimopsis variabilis and Sphaerolaimus sp. 1 were epigrowth feeders (2A) and omnivores/predators (2B), respectively. The same pattern was also found for the most abundant species. The top five most abundant canyon species were non-selective deposit feeders (1B) with clavate tail and r-selection strategy. In contrast, Craspodema sp. 1 and Paramesacanthion tricuspis from the slope habitats belonged to epigrowth feeder (2A) and omnivore/predator (2B), respectively.

Figure 5. Non-metric multi-dimensional scaling (nMDS) based on squared-root transformed species abundance. The rectangular labels in panel (A) indicate the abundance-weighted averages scores for the top five most contributing species in the GPSC and adjacent slope, respectively. The species contributions were decomposed from the average Bray-Curtis dissimilarity within the GPSC and adjacent slope using SIMPER (Supplementary Tables S3, S4). In other words, the label shows where the high abundance of the five most contributing species may occur on the nMDS ordination. Panel (B) shows the correlations between the best-subset of environmental variables and nMDS ordination axes. The environmental subset was selected by distance-based redundancy analysis (dbRDA) with the smallest AIC. The green labels show mean bottom current velocity (Spd) and duration of bottom current velocity exceeding 20 cm/s (Over20), salinity (Salin), temperature (Temp), C/N ratio (CN), total organic carbon (TOC), and light transmission (Trans).

Table 3. Three-way cross permutation analysis of variance (PERMANOVA) on nematode species composition.

Based on distance-based redundancy analysis (dbRDA), the best subset of environmental factors explaining nematode species composition (adjusted R² = 58.8%, p < 0.05) were variables related to:

1. Internal-tide energy, i.e., mean bottom current velocity (Spd), duration of bottom current velocity exceeding 20 cm/s (Over20), and light transmission (Trans),

2. Food supply, i.e., total organic carbon contents (TOC) and carbon to nitrogen ratio (CN),

3. Water masses, i.e., salinity (Salin) and temperature (Temp).

These variables are mapped onto the same nMDS plot using their correlations with the ordination axes (Figure 5B), showing that the canyon assemblages were characterized by increasing Over20 and Spd and the slope assemblages by increasing TOC and Trans.

Multiple-site Sørensen dissimilarities were 78 and 81% along the GPSC and slope transects, in which the spatial turnover component contributed to 88 and 96% and the nestedness component contributed only 12 and 4% of the total β-diversity, respectively (Figure 6A). Between the canyon and slope habitats, the Sørensen dissimilarities were between 72 and 79% across the four depth strata (Figure 6B). The spatial turnover processes still dominated and accounted for 69–75% of the total β-diversity, but the nestedness patterns also played considerable roles and contributed 25–31% of the between-habitat dissimilarity across the four depth strata.

Figure 6. Partition of multiple-site nematode (Sørensen) dissimilarity into spatial turnover and nestedness components. Panel (A) shows the Sørensen dissimilarity and its turnover and nestedness components among the four sites within the GPSC and adjacent slope transects. Panel (B) shows the Sørensen dissimilarity and its turnover and nestedness components between the canyon and slope site pair across different depth strata. The bar chart and error bars show the mean and standard deviation of the multiple-site dissimilarity and its components across different sampling times.

Shannon diversity, functional diversity (q = 1), trophic diversity, and maturity index were negatively correlated with the mean bottom current velocity and positively correlated with total organic carbon (TOC) contents (Figure 7). Environmental gradients mainly drove these significant correlations, which arose from strong bottom currents and low TOC in the GPSC as well as relatively weak currents and high TOC on the adjacent slope (Figure 7).

Figure 7. Correlations of diversity and maturity indices to proxies of internal tide energy and food supplies. The mean effective number of abundant species (A,B), effective number of abundant functional group (C,D), trophic diversity (E,F), and maturity index (G,H) per station are plotted against the mean bottom current velocity and total organic carbon in the surface sediments. Error bar indicates standard deviation.

Discussion

Study Limitations

Despite that 163 nematode species were uncovered from the GPSC and adjacent slope for the first time, the current study is by-no-means a comprehensive inventory of the nematode species diversity in the region. Given the small area of seafloor sampled (i.e., only 283 cm² of surface sediment from 46 subsamples), the gamma diversity estimate of nematodes is no doubt an underestimation. Thus, the interpretation of our results should only focus on comparisons between the canyon and non-canyon communities rather than the absolute alpha, beta, and gamma diversities. Moreover, like many deep-sea benthic studies, we took only one multicore per station and allocated resources to cover a larger survey area. In any deep-sea survey, there is always a trade-off between whether to sample more replicates to give higher statistical power or more stations to cover a larger sampling area. To resolve such dichotomy in the sampling design, Montagna et al. (2017) partitioned the variability of abundance and diversity of macrofauna among sampling locations, multicorer deployments, and core tubes of a multicorer. Gallucci et al. (2009) also partitioned the variability of total and species abundance and assemblage composition of nematodes among core tubes (i.e., from a multicorer) and within a core tube (i.e., subcores). Both studies found that the smaller-scale heterogeneity (i.e., core tubes for macrofauna or subcores for meiofauna) explained more variability in the nested sampling designs. Therefore, given limited resources in deep-sea researches, it may be more cost-effective to allocate efforts to increase the spatial cover rather than repeat deployments at a single location (Montagna et al., 2017). In this study, we have taken a more conservative approach by using the true replication (i.e., averages of subsamples per core) from two short repeated cruises to ensure robust statistical hypothesis testing. Nevertheless, the interpretation of our results still requires caution. While a conservative approach may be effective in reducing Type I error (i.e., false-positive), the risk of Type II error (i.e., false negative) may increase.

Response of Diversity to Disturbance

Submarine canyons represent remarkably heterogeneous habitats on the otherwise seemingly monotonous, mud-covered seafloor of the continental margins. Due to the funneling effect, the seafloor of the canyon thalweg often receives and accumulates a large amount of terrestrial and oceanic organic matter. The highly diverse microhabitats and the ample food supplies in the canyon likely allow species to coexist and thus sustain high biodiversity and ecosystem functioning (McClain and Barry, 2010; Leduc et al., 2014; Fernandez-Arcaya et al., 2017). Therefore, the majority of the studies have found that the diversity and abundance of meiofauna were significantly higher in the canyons than the adjacent slopes (Baguley et al., 2006; Ingels et al., 2009; Román et al., 2016; Gambi et al., 2019). However, contrasting results or lack of difference between the canyon and slope communities have also been found, especially in large-scale comparative studies (Soltwedel et al., 2005; Garcia et al., 2007; Bianchelli et al., 2010; Leduc et al., 2014; Bianchelli and Danovaro, 2019). For example, lower abundance and diversity of meiofauna have been found in the Congo Canyon than on the adjacent slope, presumably, related to the strong bottom currents and unstable seafloor conditions within the canyon (Van Gaever et al., 2009).

The inhabitants of the Gaoping Submarine Canyon (GPSC) are no stranger to extreme environmental conditions including high sediment loads originating from the small mountain river, extreme monsoonal and typhoon rainfalls triggering subsequent turbidity currents, strong internal-tide energy causing swift bottom currents, and high tectonic activity leading to frequent debris flows (Jan et al., 2008; Chiou et al., 2011; Liu et al., 2013, 2016). These physical disturbances within the GPSC substantially suppressed the vulnerable macrofauna taxa (Liao et al., 2017) and depressed the species and functional diversity of nematode assemblages (this study). Nevertheless, the nematode abundance was not affected by the adverse condition in the GPSC, suggesting rapid recovery of the more resilient nematode species. In fact, the dominant species in the GPSC, including members of Daptonema, Sabatieria, and Metadesmolaimus (belonging to Families Comesomatidae and Xylidae), were eurytopic, opportunistic species with short generation time (c-p value = 2). These r-selection strategists may quickly occupy the empty niches left by eliminated predecessors to maintain the population size after disturbance. In comparison, extremely high biomass but low diversity of nematode assemblages was found in the active Kaikoura Canyon off eastern New Zealand, where the Daptonema and Sabatieria pulchra were the dominant taxa due to their ability to adapt to hypoxia and eutrophic conditions (Leduc et al., 2014). This scenario is also common in the benthic communities under turbidite impact or hypoxia after eutrophication (Tsujimoto et al., 2008, 2020; Vanreusel et al., 2010; Zeppilli et al., 2018).

Characteristic Species and Functional Traits

Nematode assemblages of submarine canyons often share a considerable number of genera with adjacent slope habitats, although the canyon communities may comprise more dominant taxa with a lower species evenness (Vanreusel et al., 2010). In this study, we found that the GPSC and the adjacent slope harbored distinct nematode assemblages. On the slope, the species, functional, and trophic diversity and maturity index of nematodes were significantly higher than those in the GPSC. Except for the 46 common species that occurred in both the canyon and slope, seven species were restricted to the GPSC, and 110 species were only found on the slope. The canyon nematode assemblages were dominated by the non-selective deposit feeding (1B) species with r-selection strategy and clavate tail (e.g., Daptonema, Axonolaimus, Metadesmolaimus, and Sabatieria). These same species also contributed highly to the average dissimilarity within the canyon. The Daptonema and Sabatieria species are commonly found in the disturbed, organic-rich, and polluted sediment (Vanreusel, 1990; Schratzberger and Jennings, 2002; Moreno et al., 2008). The Sabatieria species is also known to adapt to anoxic conditions and thus prevalent in reduced sediment (Jensen, 1984; Soetaert and Heip, 1995; Steyaert et al., 2007). The dominance of r-selection strategists and depressed maturity index in the GPSC indicated that not just a few dominant species have low c-p value (=2), but that the canyon assemblage almost entirely consisted of colonizers. The canyon nematodes may benefit from the colonization strategy so they can quickly recover from the extremely unstable conditions associated with the recurrent disturbances by internal tides or periodic disturbances by turbidity currents (Wang et al., 2008; Chiou et al., 2011; Liu et al., 2013; Hsu et al., 2014). Moreover, the clavate-shape tails have been considered as typical morphology of interstitial inhabitants of high-energy, sandy habitats, in which the long tail might be used as a tool to anchor the body (Riemann, 1974; Thistle et al., 1995). The slope assemblages, on the other hand, were dominated (also highly contributed) by species with more diverse feeding guilds including non-selective deposit feeders (1B, e.g., Sabatieria, Setosabatieria, Cervonema and Elzalia), epigrowth feeders (2A, e.g., Dorylaimopsis and Craspodema), and omnivores/predators (2B, e.g., Sphaerolaimus and Paramesacanthion). Some selective deposit-feeding species (1A, e.g., Quadricoma, Pselionema, and Halalaimus) and omnivores/predators (2B, e.g., Abelbolla) with K-selection strategy (c-p value = 4) also had high contribution to the average assemblage dissimilarity on the slope (Supplementary Table S4), indicating that the slope offers a stable habitat where persister species can thrive. In contrast to the canyon assemblages, the slope nematodes had more diverse feeding guilds, tail shapes, and life strategies, which resulted in more functional groups with higher trophic diversity. The increasing relative abundance of nematodes with conical and elongated/filiform tails on the slope also signals the change of sediment types (Riemann, 1974; Thistle et al., 1995). These tail shapes are well-represented in fine sand and muddy sediments with little interstitial space, where nematodes may burrow or push sediment particles to create their own space (Singh and Ingole, 2016; Armenteros et al., 2019).

Environmental Controls on Diversity and Composition

The canyon assemblages differed from those of the adjacent slope in that they experienced stronger internal tide energy, including higher modeled bottom current velocity (i.e., Spd) and longer duration of currents stronger enough to induce sediment erosion (i.e., Over20). The organic matter might be rapidly buried after episodic turbidity currents such as occurred during the 2011 Tohoku Earthquake (Tsujimoto et al., 2020). The strong internal tides may cause intermittent resuspension of sediment and transport of the resuspended material in the bottom nepheloid layer (Liu et al., 2010). The internal tides may also lead to long-lasting sediment winnowing and prevent the organic-rich particles from settling on the seafloor and thus decrease the total organic carbon (TOC) contents in the sediments (Liao et al., 2017). As a result, the slope assemblages were separated from the canyon by the higher TOC contents in the surface sediment and the clearer bottom water (with higher light transmission). In submarine canyons, frequent physical disturbances and abundant organic inputs may favor opportunistic nematode species and lead to high dominance and low evenness assemblages (Garcia et al., 2007; Ingels et al., 2009; Leduc et al., 2014; Gambi and Danovaro, 2016). In the GPSC, unlike most other canyons, the total organic carbon (TOC) was also in short supplies due to the sorting and winnowing of particles by strong bottom currents (Liao et al., 2017). The contrast between the physically dynamic canyon and more tranquil slope habitat not only drove the species composition of nematode assemblages but also caused the species and functional diversity to decrease with the increasing bottom-current velocity and decreasing TOC contents from the slope into the canyon. Within the individual canyon and slope transects, we found that the observed number of species and the species diversity (i.e., Hill numbers of order q = 0–2) both increased with depth. The rates of increase, however, were different between the canyon and slope, suggesting that the underlying drivers may be different. The head of GPSC is known to experience frequent submarine geohazards (Hsu et al., 2008; Su et al., 2012; Gavey et al., 2017). The bottom currents at the canyon head regularly exceed 1 m/s, but then the internal-tide energy and near-bottom currents weaken away from the head region along the thalweg (Wang et al., 2008; Chiou et al., 2011; Liao et al., 2017). Due to the relaxation of physical disturbance, it may be expected that the nematode diversity increased slightly with depth (or away from the canyon head). Other physical disturbances such as bottom trawling could also remove the organic-rich sediment and cause a reduction of meiofauna diversity and ecosystem function (Puig et al., 2012; Pusceddu et al., 2014; Román et al., 2019). In our study region, the bottom trawl fishery typically targets water depths between 100 and 400 m, while the GPSC is avoided to prevent gear damage. Therefore, we cannot rule out the possibility that the upper slope sites (i.e., GS1, GS2, Figure 1 and Table 1) may have experienced some fishing impacts and thus had low nematode diversity. In the meantime, the increasing nematode diversity with depth on the slope could also be explained by a unimodal diversity-depth relationship (Etter and Grassle, 1992; Boucher and Lambshead, 1995), in which the slope sites (278–848 m) happened to be on the ascending limb of the general unimodal diversity gradient.

Mechanisms Shaping the Nematode Assemblages

Along the GPSC and adjacent slope, the spatial turnover (or species replacement) was the main contributor of β-diversity, suggesting that the environmental filtering, not the nestedness, was the primary mechanism driving the bathymetric variation of species composition. Although the depth range of this study was less than 1,000 m, there were large environmental variations with depth, including increasing pressure, decreasing temperature and dissolved oxygen, as well as the increasing total organic carbon (TOC), total nitrogen (TN) and porosity in the sediment (Liao et al., 2017). Within the GPSC, the internal tide energy also weakened considerably with depth, resulting in declining bottom current velocity, duration of sediment erosion, and increasing water clarity or light transmission toward the deeper regions of the canyon (Liao et al., 2017). These environmental variations may replace (or filter) the poorly fitted species with the better-fitted one to cause a gradual, continued species turnover along the length of the continental margin (McClain and Rex, 2015). Moreover, the physical disturbances by bottom currents, turbidity and debris flows may contribute to the heterogeneity of the canyon seafloor (Liu et al., 2016; Chiang et al., 2020) and thus drive the bathymetric variation in species composition (Zeppilli et al., 2016, 2018). In the GPSC, we found that the selective deposit feeding (1A) and K-selection species tend to increase toward the deeper depth strata, providing some indirect evidence of decreasing physical disturbance toward the deeper canyon. Since the environmental conditions in the submarine canyons are obviously more dynamic than the gentle sloping seafloor, many studies found conspicuous bathymetric difference in the canyon nematodes as opposed to the more homogeneous slope assemblages (Ingels et al., 2009; Leduc et al., 2014; Gambi and Danovaro, 2016; Rosli et al., 2016; Román et al., 2019). Nevertheless, we found that the bathymetric β-diversity was comparable between the GPSC and the adjacent slope.

Strikingly, even at a similar depth, the β-diversity between the pairs of GPSC and slope sites were just as high as the bathymetric β-diversity along the canyon or slope transects. For the canyon-slope β-diversity comparison, the spatial turnover still dominated, but the nestedness also played a considerable role (25–31% of the β-diversity) across the four depth strata. For the spatial turnover, it is apparent that the majority of the slope nematode species did not occur in the canyon. The abrupt change in species composition or rapid species replacement is evident in the nMDS ordination, in which the canyon and slope assemblages were clearly separated without overlap. It may be possible that the environmental stress in the GPSC prevented the colonization of K-selection species or larger individuals, which effectively filtered out ∼71% of the slope species. The seven endemic species from the canyon might be well-adapted to the extreme environmental conditions but unable to compete with the abundant persisters on the slope. Possibly, these species also occur on the slope, but in too low numbers to be detected by our sampling. In terms of the 46 species that occurred in both the canyon and slope, these nematodes might be the tolerant species that adapt to the strong physical disturbance and lower TOC in the GPSC and probably immigrated from the adjacent slope by ocean currents (Palmer, 1988). Despite that the nematodes are generally considered as inferior swimmers (Jensen, 1981), they can be suspended by weak currents and move a notable distance (Boeckner et al., 2009). Active colonization by nematode assemblages with different species composition may also influence the recolonized communities and increase the patchiness in a new habitat (Ullberg and Ólafsson, 2003; Fonsêca-Genevois et al., 2006; Lins et al., 2013). This may have given rise to the nestedness pattern in the canyon, where some of the species appear to be subsets of slope communities (e.g., ∼87% of the canyon species also occurred on the slope). It is also possible that some shared species might be morphologically similar (so identified as the same species) from the same genus but displayed distinct ecological functions and adapted to both the canyon and slope conditions (De Meester et al., 2016). The majority of species could be sensitive species from the slope, which may have been subjected to local extinction and disappearance from the canyon or experienced some degree of dispersal limitation (e.g., washed away by strong currents) as they immigrated into the GPSC to contribute to the nestedness patterns between the canyon and slope (Baselga, 2010, 2012).

Contrast Between Meiofauna and Macrofauna

Although all aspects of the nematode diversity (i.e., species, functional, and trophic diversity) measured in this study appeared depressed in the GPSC, the nematode abundance was not affected by the adverse conditions in the GPSC. In contrast, the macrofauna abundance was significantly depressed in the GPSC, especially toward the canyon head, where the internal-tide energy is the strongest (Liao et al., 2017). Within the macrofauna taxa, the peracarid crustaceans (i.e., amphipods, cumaceans, isopods, and tanaids) appeared most sensitive to environmental perturbations and thus declined or disappeared altogether in the GPSC, whereas the tolerant harpacticoids, nematodes, and polychaetes (macrofaunal or meiofaunal size) were common in both the canyon and the slope (Liao et al., 2017). Nematodes are known to be more resilient to physical disturbances such as sediment resuspension and deposition than the larger macrofauna taxa (Lambshead et al., 2001; Schratzberger and Jennings, 2002; Whomersley et al., 2009). After disturbances, nematodes typically show higher turnover and lower mortality rates than macrofauna in the early successional stages (Schratzberger and Jennings, 2002). The meiofaunal abundance may quickly recover within a few months but the diversity tends to remain low and species composition does not recover to its original state for prolonged time (Kitahashi et al., 2014, 2016; Tsujimoto et al., 2020). Considering that the GPSC may be subject to disturbance by internal tide flushing on a daily basis, only the resilient nematode taxa (e.g., Daptonema, Sabatieria, and Axonolaimus) may quickly recover from the disturbance events and flourish in the canyon. The fast growth rates and short generation times (i.e., weeks to months) of these tolerant nematode taxa may help them to recover rapidly from the recurrent disturbance events (Giere, 2009). Nevertheless, some uncertainty remains. During the fall sampling, the nematode abundance appeared lower in the GPSC than on the adjacent slope, but the habitat effect was not detected when both the summer and fall sampling were taken into consideration. The episodic flooding in the Gaoping River and the subsequent turbidity currents in the GPSC were unpredictable but mostly followed the monsoonal rainfalls and typhoons which mostly occurred between May and October (Liu et al., 2016); therefore, the abundance variations between habitats may have been smeared by the potential seasonal effect or unpredictable sediment flows. A long-term, repeated sampling will be needed to resolve the potential effect of seasonality or episodic events on the nematode assemblages.

Conclusion

The Gaoping Submarine Canyon (GPSC), which is subject to frequent turbidity currents and debris flows and strong near-bottom currents driven by internal tides, is an ideal testbed to study the effect of large-scale disturbance on deep-sea benthic communities. In the summer and fall of 2015, we conducted the first detailed deep-sea nematode diversity investigation in the GPSC and adjacent slope. When considering the average effect across both the sampling campaigns, we found no statistical difference in nematode abundance across the depth transects and between the GPSC and adjacent slope. This suggests that the canyon environment had little effect on the nematode population densities. Most strikingly, we found that the nematode species, functional and trophic diversity and maturity were significantly depressed in the GPSC, corresponding to stronger near-bottom currents and lower total organic carbon (TOC) contents in the surface sediments. We also found that the nematode species diversity increased with depth in the upper GPSC and adjacent slope, consistent with the classical paradigm of unimodal diversity-depth relationship in the deep sea. Between the GPSC and adjacent slope, the distinct species composition was mainly contributed by species replacement due to environmental filtering. The canyon assemblages were mainly driven by stronger internal tide energy, and the slope assemblages were associated with higher total organic matter in the surface sediments from which a higher food supply may be inferred. As a result, the canyon nematode assemblages were dominated by non-selective deposit-feeding species with r-selection strategy and clavate tail (e.g., most species of Comesomatidae and Xylidae), whereas on the slope, the longer life-span species with diverse feeding strategies and tail shapes may co-exist. While environmental filtering was the dominant mechanism, the considerable contribution of nestedness to between-habitat nematode dissimilarity suggests that species immigration and local extinction may also play an important role in shaping the canyon assemblages (e.g., opportunistic species from adjacent slope became dominant in the canyon).

In the age of global climate change, the deep seafloor is expected to experience warming, deoxygenation, acidification, and decline in food supplies (Mora et al., 2013; Sweetman et al., 2017; Breitburg et al., 2018). On the continental slope, the changes in climatic drivers may lead to species range shift (Smith et al., 2012), reduction of suitable habitat (Morato et al., 2020), reduction of abundance and biomass (Jones et al., 2014), decline in biodiversity (Sperling et al., 2016; Brito-Morales et al., 2020), and eventually shift of relative contribution of carbon-processing from the large to small benthic size classes (Jones et al., 2014; Sweetman et al., 2017). River-connected canyons such as GPSC may not only be affected by the marine expressions of climate change, but also by terrestrial climate change impacts due to their close association with fluvial processes. For example, the changes in storm and precipitation patterns may lead to changes in the frequency and intensity of flood events (Knutson et al., 2010; Westra et al., 2014; Bacmeister et al., 2018), which may cascade down to trigger turbidity currents and the associated submarine geohazards. Understanding the impact on species and functional diversity from large-scale natural disturbance in submarine canyons like GPSC will enhance our ability to predict the potential impacts of climate change and other anthropogenic effects on deep-sea ecosystems.

Data Availability Statement

The datasets generated for this study can be found in the https://github.com/chihlinwei/nema.

Author Contributions

J-XL and C-LW designed and conceived the study, to which all authors contributed the ideas and discussion. J-XL and C-LW executed the field sampling. J-XL conducted the laboratory work and drafted the manuscript with substantial inputs from C-LW and MY. J-XL, C-LW, and MY contributed to the data analysis. All authors contributed to the article and approved the submitted version.

Funding

This project was part of an integrated project, Fate of Terrestrial/Non-terrestrial Sediments in High Yield Particle-Export River-sea Systems (FATES-HYPERS) and South China Sea Marine Research – Toward the south (SCSMART): Environment Change and Marine Biodiversity, sponsored by the Ministry of Science and Technology (MOST 108-2611-M-002-001 and 108-2119-M-001-019). J-XL was supported by the NTU and MOST Postdoctoral Fellowship (MOST 108-2811- M-002-622).

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.

The reviewer NG declared a past co-authorship with one of the authors C-LW to the handling editor.

Acknowledgments

We would like to thank the Institute of Oceanography (IO), the National Taiwan University (NTU), and the Ministry of Science and Technology for supporting the field work, analysis, and manuscript preparation. We also thank Anne Chao for providing R code to compute functional diversity and for valuable suggestions to improve the manuscript. We also thank Sen Jan and Ming-Da Chiou for providing the internal tide model. We also thank captain, crew members, and technicians of the R/V Ocean Researcher I, as well as the graduate students who participated in OCEAN 7090 Field Work in Marine Biology. The majority of samples and data were collected by students through the field course.

Supplementary Material

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

References