http://orcid.org/0000-0003-2488-8609
Figures
Abstract
The presence of large densities of the piston-pumping polychaete Lanice conchilega can have important consequences for the functioning of marine sediments. It is considered both an allogenic and an autogenic ecosystem engineer, affecting spatial and temporal biogeochemical gradients (oxygen concentrations, oxygen penetration depth and nutrient concentrations) and physical properties (grain size) of marine sediments, which could affect functional properties of sediment-inhabiting microbial communities. Here we investigated whether density-dependent effects of L. conchilega affected horizontal (m-scale) and vertical (cm-scale) patterns in the distribution, diversity and composition of the typical nosZ gene in the active denitrifying organisms. This gene plays a major role in N2O reduction in coastal ecosystems as the last step completing the denitrification pathway. We showed that both vertical and horizontal composition and richness of nosZ gene were indeed significantly affected when large densities of the bio-irrigator were present. This could be directly related to allogenic ecosystem engineering effects on the environment, reflected in increased oxygen penetration depth and oxygen concentrations in the upper cm of the sediment in high densities of L. conchilega. A higher diversity (Shannon diversity and inverse Simpson) of nosZ observed in patches with high L. conchilega densities (3,185–3,440 ind. m-2) at deeper sediment layers could suggest a downward transport of NO3− to deeper layers resulting from bio-irrigation as well. Hence, our results show the effect of L. conchilega bio-irrigation activity on denitrifying organisms in L. conchilega reefs.
Citation: Yazdani Foshtomi M, Leliaert F, Derycke S, Willems A, Vincx M, Vanaverbeke J (2018) The effect of bio-irrigation by the polychaete Lanice conchilega on active denitrifiers: Distribution, diversity and composition of nosZ gene. PLoS ONE 13(2): e0192391. https://doi.org/10.1371/journal.pone.0192391
Editor: Kay C. Vopel, Auckland University of Technology, NEW ZEALAND
Received: June 16, 2017; Accepted: January 23, 2018; Published: February 6, 2018
Copyright: © 2018 Yazdani Foshtomi et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: The nosZ data set (dereplicated, chimera-corrected, and Framebot-filtered nucleotide sequences) is available in the MG-RAST database (https://metagenomics.anl.gov) under accession number “mgm4762338.3”. The nosZ nucleotide sequences (502 sequences), corresponding to unique amino acid transcripts (nosZ-UATs), are available in EMBL’s European Nucleotide Archive (ENA, http://www.ebi.ac.uk/ena) (GenBank/DDBJ) under accession numbers “LT908519-LT908520, LT908570-LT908579, LT908638-LT908727, LT908734-LT908833, LT908874-LT909023, and LT909024-LT909173” under the project number “PRJEB22661”. All other relevant data are within the paper and its Supporting Information files.
Funding: Funded by Fonds voor Wetenschappelijk Onderzoek Vlaanderen. http://www.fwo.be/. Grant number G0033.11 to MV, JV, AW. Special Research Fund of Ghent University. https://www.ugent.be/nl/onderzoek/financiering/bof. Grant number 01G01911 to MV and AW. Ministry of Science, Research and Technology, Iran. https://www.msrt.ir/en. Grant number 24-88200022 to MYF. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: The authors have declared that no competing interests exist.
Introduction
Denitrification is a key process in the biogeochemical cycling of nitrogen. It is a primary loss mechanism for nitrogen in the nitrogen budget of coastal ecosystems [1]. It is a four-step respiratory process, in which nitrate (NO3−) is reduced sequentially to nitrite (NO2−), nitric oxide (NO), nitrous oxide (N2O), and nitrogen gas (N2), and is mediated by a taxonomically diverse group of microorganisms, mainly bacteria [2]. Denitrification counteracts eutrophication by removing N2, but it can contribute to global warming and ozone depletion due to the release of NO and N2O [3–5].
The distribution of denitrifying organisms and denitrification rates in marine ecosystems is affected by different environmental factors including organic matter content, concentration gradients of dissolved inorganic nitrogen (DIN: NO3−, NO2− and NH4+) and the availability of oxygen [6,7]. Stratification and vertical distribution of denitrifying organisms also results from the vertical distribution of these environmental factors in the sediment [7,8].
Denitrifiers are facultative anaerobes and denitrification is often coupled with nitrification across oxic/anoxic interfaces in the sediment [9]. With increasing sediment depth, oxygen concentrations decrease and denitrifying communities shift from facultative anaerobes to strict anaerobes [8]. However, maximum denitrification rates do not necessarily occur at the lowest O2 concentration in deep layers [8,10, 11] due to low/absent NO3− and NO2− concentrations. The NOx− concentration may therefore act together with O2 to control denitrification rates and the diversity and structure of the denitrifying community [8,10,12]. In addition, the effects of species-specific adaptation and variations in substrate threshold among denitrifying organisms need to be considered as well. For example, while oxygen has an inhibitory effect on the denitrification enzymes of Agrobacterium tumefaciens [13], Pseudomonas stutzeri is able to denitrify under aerobic conditions [14].
In general, microbial community patterns in marine ecosystems can vary from large (km or m) [15, 16, 17] to small scales (cm or less) [12,17,18]. In large scales the different patterns are due to the wide range of environmental conditions including changes in median grain size and organic matter concentration [15,16,17]. At small scales vertically in the sediment, microbial variations can be related to the different factors such as sediment types [12,19], the steepness of redox gradients [17], substrate availability [20] as well as activity and sediment mixing effects of larger fauna [17,21]. The boundary between oxic and anoxic conditions, for example, along the vertical gradient within the sediment migrates on spatial scales of micrometers to centimeters and on time scales of seconds to hours [22]. Denitrifiers have the ability to adapt to oscillations of oxygen concentrations during the rapid transition from oxic to anoxic conditions [23,24]. Such oxygen oscillations and different oxygen regimes have effects on the structure of the denitrifying community [25].
Functional diversity and abundance of larger organisms affect the biogeochemical cycles and variations in composition and diversity of microbial communities [16], creating small-scale heterogeneity in the sediment as microniches [20,26]. Building and irrigating of burrows by bioturbating organisms like worms, for example, can provide a unique environment for micro-organisms inhabiting the burrow wall [27,28] or enhance total microbial metabolism in adjacent surrounding sediments [29]. However, the impact of bioturbation depends on individual functional effects of burrowing organisms [20,28,30].
Lanice conchilega is a tube-building worm (Polychaeta) that glues sediment grains together to form tubes reaching 10–30 cm vertically into the sediment [31,32] and extending 2–3 cm into the water [33]. It is present throughout the North Sea and can be found in dense aggregations [34] (over 3000 ind. m-2 [35]), referred to as biogenic reefs [36]. This polychaete is widely distributed in soft bottom environments and acts as a bio-engineer both in intertidal and subtidal areas [37,38]. L. conchilega manifests both autogenic and allogenic ecosystem engineering properties via, respectively, its own physical structures and by transforming living/non-living materials from one physical state to another [39]. The autogenic engineering effect provides new habitat for associated species by increasing bed stability [36] and trapping organic matter from the water column [40]. It substantially affects the structure and abundance of associated communities and food-web properties [40–42]. The allogenic engineering capacity of L. conchilega is reflected in its piston-pumping activity [31] causing water and solute exchanges between tubes and the overlying water (bio-irrigation). This stimulates nutrient fluxes, mineralisation and denitrification processes in coastal sediments (as deduced from lab experiments) and seems to be density-dependent [43]. However, the effect of L. conchilega on the microbial communities mediating these processes remains unexplored.
The capacity for denitrification can be deduced from functional genes involved in denitrification pathways, which reflect the distribution and function of denitrifying organisms in the environment [1]. The detection of a gene in the environment (DNA-based methods), however, does not imply that the corresponding activity is present [44,45]. Instead, the expression of genes (RNA-based methods) and the detection of enzymes (protein-based methods) allow investigation of the direct link between the composition and density of the transcribed denitrification genes and rate of denitrification [44].
In this study, we investigated the expression of the nosZ gene encoding the enzyme which catalyses the final step of denitrification (conversion of N2O to N2 gas; [2]). Microbial taxa possess divergent nosZ clusters (typical and atypical nosZ) with genes that are related yet evolutionarily distinct from each other [4]. Typical nosZ genes, which occur more commonly in bacteria with a complete denitrification pathway (co-occurrence of denitrifying genes; [4,46]), play the major role in N2O reduction in coastal ecosystems [25]. Atypical nosZ, in contrast, seems to occur frequently in non-denitrifying organisms that do not possess any other denitrification genes and is potentially important in terrestrial environments [4,46].
Here, we investigated how distribution, diversity and composition of the typical nosZ gene, as well as the biogeochemical environment (oscillation in oxygen concentration, sediment grain size and chlorophyll a) are affected by the bio-irrigation activity of L. conchilega vertically at small (cm) scales along the sediment depth profiles and also at larger horizontal scales (m) at different natural L. conchilega densities within the reef. The observed diversity patterns in nosZ were then statistically linked to environmental variables to objectively assess the link between faunal-mediated heterogeneity and typical nosZ gene diversity.
To our knowledge, the current study is novel as most of the previous studies were limited to the bioturbation and bioirrigation impact of larger organisms on the rate of denitrification [47–51] and only a few studies have investigated this effect on the denitrifying functional genes [52–54] and still the active community investigation was neglected. As such, the current study can improve our understanding of the effects of bio-irrigation on diversity of nosZ gene in the active denitrifying community in marine sediments at small spatial scales.
Materials and methods
Study site and sampling
Sampling was done in October 2014 at the intertidal zone of the sandy seashore of Boulogne-sur-Mer, along the northern part of the English Channel (50° 44.01’ N, 01° 35.15’ E; Northern France; Fig 1). Sediments of the L. conchilega reef were collected by core (Plexiglas, 78.5 cm2 surface area; height: 25 cm) from a limited area of the western reef zone (about 4000 m2) located higher on the beach and exposed at every low tide [38]. No specific permits for sampling and ethics requirements were needed since our research was approved by the FWO-research proposal (G.0033.11).
Sediment samples were collected from L. conchilega patches in an area of almost 100 m2 within the reef. Tidal height and sediment composition were similar for all samples. Three replicate cores were taken from three different areas of the patches: (i) an area with an average density of L. conchilega individuals (based on counts of fringed tubes, [55]) of 25–27 tubes per core surface (3185–3440 ind. m-2) (hereafter referred to “high L. conchilega treatment” or “H”), (ii) an area with lower L. conchilega densities (5 tubes per core surface; 637 ind. m-2) located on the edge of the patches (“low L. conchilega treatment” or “L”), and (iii) sediments without L. conchilega between the reef patches (“control treatment” or “C”).
We used intact sediment cores from the reef. However, L. conchilega aggregations have been shown to support a high diversity and abundance of macrofauna in comparison with the adjacent sediments without L. conchilega [41]. Therefore, differences observed in the nosZ community might also result from activities of the associated community which is different from the control situation. In our sampling site, the amphipods from the genus Urothoe and the polychaetes Eumida sanguinea, Pygospio elegans, Heteromastus filiformis and Capitella sp. were the dominant macrofaunal species after L. conchilega [41]. Some of these species, i.e. the surficial modifier Urothoe and the fixed tube builder Pygospio elegans, are small (max 1.5 cm; [32,56,57]) and cannot affect the ecosystem as strongly as the bio-irrigator L. conchilega, making up 20–25% of total macrofaunal density in the sampling month [41]. Visual inspection did not show abundance of other large organisms in our sampling cores while Heteromastus filiformis and Capitella sp. show limited movements in the sediment [57]. Eumida sanguinea is a biodiffuser [57] but there is no record of the functional importance of this species for N-cycle processes. In addition, considering the high ability of L. conchilega to build reef structures [36,37], the effects on the denitrifying community in our study are therefore considered to result mainly from the presence and activity of L. conchilega.
The intact sediment cores (3 treatments × 3 replicates) were transported to the lab and submerged, uncapped without creating suspension of the sediment in tanks containing continuously aerated seawater at an in situ temperature-controlled room (17 °C in October). Every core was aerated separately by a stream of fine bubbles in the overlying water. In order to create water circulation inside the cores, Teflon-coated magnets were inserted at appropriate distance from the sediment surface and rotated by a central magnet in the tanks at a speed below the resuspension limit.
Within two days after sampling, vertical profiles of sediment oxygen concentration were measured (three profiles per core) in the ambient sediment at 1 to 2 cm distance from tubes, using Unisense oxygen microsensors (type ox100) in vertical increments of 250 μm. Bio-irrigation activities of L. conchilega individuals were assessed as well, by logging changes of sedimentary oxygen concentrations every second for 30–35 min at 1.5 mm sediment depth, at 1 cm distance from the tube, using the same Unisense microsensors.
At the end (three days after sampling), the upper sediment layers were sliced in 0.5 cm intervals (0–0.5 [d1], 0.5–1 [d2], 1–1.5 cm [d3]) to investigate microbial communities in the oxic and anoxic zones based on the measured oxygen profiles. This is the smallest possible scale without disturbing the sediment in the adjacent layers. A fourth layer (2.5–3 cm [d4]) was collected as well because L. conchilega was reported to ventilate its tube up to 2.4 cm tube length [31]. All sediment slices were homogenized before collecting subsamples for further analyses.
The sediment was subsampled by taking 3 to 5 mL using a cut-off syringe for the analysis of labile organic matter (chl-a), % clay and silt content and median grain size (MGS). Around 4 g of sediment (wet weight) was subsampled using sterilised tools for microbial analyses and stored in sterile 15-mL falcon tubes. Microbial and chl-a samples were immediately frozen at -80 °C, whereas the samples for grain size were dried at 60 °C before analysis.
Chl-a was determined by HPLC (Gilson, Middleton, Wisconsin, USA) analysis according to [58] and sediment granulometry by laser diffraction (Malvern Instruments, Malvern, UK).
RNA extraction and nosZ sequencing
RNA was extracted from 4 g homogenised sediment (wet weight) from each depth layer (36 samples = 3 treatments × 3 replicates × 4 layers) using the RNA Power Soil Total RNA isolation kit (MO BIO Laboratories). Integrity and purity of RNA extractions were checked following [16]. The DNA-free RNA samples were reverse transcribed into cDNA using Omniscript Reverse Transcriptase Kit (Qiagen) according to the instructions of the manufacturer and using 10 μM of Random Hexamer primers (Applied Biosystems) and 5 μL RNA template per total volume of a reaction (20 μL).
RT-PCR amplification of nosZ gene fragments (267 bp) was performed in three technical replicates for each sample using specific primer sets (nosZ-2F: 5’-CGCRACGGCAASAAGGTSMSSGT-3’ / nosZ-2R: 5’-CAKRTGCAKSGCRTGGCA GAA-3’) [3] for the typical nosZ gene cluster [4,46]. The forward primer contained the 5’ Illumina adaptor, forward primer pad and linker and the nosZ-2F primer. The reverse primer consisted of the reverse complement of the 3’ Illumina adapter, 12-bp multiplex identifiers (MIDs), the reverse primer pad and linker and the nosZ-2R primer. MIDs in the reverse primer were used to identify the different samples (see S1 Table for adaptor, MID, pad, linker and primer sequences for the forward and reverse data; [59]). A control (no template) was also included to ensure that no contamination occurred in the laboratory.
Amplifications were performed in volumes of 25 μL containing 5 μL of target cDNA, 5 μL of 5xKAPA HiFi Buffer (containing 2 mM MgCl2 at 1x; Kapa Biosystems, Boston, Massachusetts, United States), 0.4 μg/μL BSA, 0.2 mM dNTPs, 0.4 μM of each primer, and 1.5 U KAPA HiFi HotStart DNA Polymerase (1 U/μL) (Kapa Biosystems). The ‘touchdown’ PCR conditions were as follows: after an initial denaturation step at 95 °C for 5 min, 20 cycles were performed consisting of three steps: denaturation (98 °C, 20 sec), annealing (70 °C, 1 min, decreasing 0.5 °C cycle-1 to 60 °C) and extension (72 °C, 1 min) followed by 15 additional cycles in which the annealing temperature was 60 °C [60]. The final elongation step was performed at 72 °C for 10 min.
PCR amplicons of the three technical replicates of each sample were combined, purified using E-Gel (Invitrogen, Life technologies) and measured with a Qubit fluorometer (Life Technologies). Samples were then pooled in equimolar amounts and loaded on a Bioanalyzer 2100 (Agilent Technologies) to check the presence of the single peak. Purified pooled libraries were submitted to Genomics Core (Center for Human Genetics UZ—K.U. Leuven) for 150×2-cycle paired-end sequencing on an Illumina Miseq platform.
Sequence analyses
Demultiplexing of the data was carried out by Genomics Core using the Illumina standard procedure by applying the bcl2fastq tool of Illumina which allows for 1 mismatch between the barcodes. Raw demultiplexed Illumina fastq files were then quality-filtered and merged using Pear v0.9.5 (Paired-End read mergeR; [61]). Reads shorter than 100 bp or longer than 400 bp and low-quality reads (scores <25) were removed from the output. Any reads containing uncalled bases and singletons were also discarded. The final file was checked for quality control with Fastqc v0.11.3 (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) [62].
We eliminated potential chimeras using the USEARCH package [63,64]. Non-target reads (non-nosZ reads) were filtered out using FrameBot included in the Fungene pipeline (http://fungene.cme.msu.edu/). We retrieved a reference set of nosZ protein sequences of 163 different species with high scores (>919) and 98% coverage with a hidden Markov model (HMM) from Fungene. FrameBot compared each member of these reference sequences to the query nucleotide sequence in both forward and reverse directions. In addition, FrameBot corrected insertion and deletion errors, and translated DNA sequences to frameshift-corrected protein sequences [64,65].
nosZ diversity was determined at the amino acid level. Therefore, a complete identity threshold (100%) of amino acid sequences was applied in USEARCH. These units (unique AA nosZ sequences) were considered to be operational taxonomic units (OTUs).
Prior to further analysis, OTUs with <0.005% relative abundance in all samples, which together accounted for less than 2% of the total reads and mostly represented singletons and doubletons, were discarded [66,67].
Maximum likelihood phylogeny of unique nosZ AA sequences was inferred using RAxML v8.2.6 on the CIPRES Science Gateway [68], under a JTT+I+G model as determined using the lowest AIC criterion in ProtTest version 2.4 [69]. Node confidence was determined using 200 bootstrap replicates. In order to construct a phylogeny with reference sequences, we performed protein BLAST searches against the NCBI non-redundant protein database (http://blast.ncbi.nlm.nih.gov/Blast.cgi) of 25 OTU representatives selected randomly across the phylogenetic tree (S2 Fig). For each blastp search, the best 100 blast hits were retained. To reduce the size of the tree, identical sequences were removed as well as sequences differing in 1 or 2 AA positions. A final alignment of 208 nosZ sequences from the database was used to construct a ML tree using RAxML as described above.
Data analyses
Environmental data analyses.
We used permutational multivariate ANOVA (PERMANOVA; Primer v6.1.10., Primer-E Ltd., Plymouth, United Kingdom with the PERMANOVA + add-on package; [70]) to perform a two-way fixed factor model design with “treatment” and “depth” as factors. “Replicate” was nested as a random factor in “treatment” [71]. We investigated whether significant differences in environmental factors (MGS, % silt+clay [mud] content and chl-a concentration) among treatments and at different depths were observed. To test for the differences in maximum oxygen penetration depth (max OPD) and oxygen profiles (both depth and time profiles) between levels of L. conchilega density, PERMANOVA with a one-way fixed factor model design (“Treatment” as a factor) was performed. Prior to the PERMANOVA analyses for oxygen profiles, a multivariate data matrix was constructed in which each depth horizon (for depth profiles) and each second (for time profiles) was considered a “variable” [72]. Pairwise tests were performed for significant terms using Monte Carlo p-values when the number of possible permutations was restricted [70].
Euclidean distance was used to calculate the resemblance matrices in PERMANOVA for the univariate environmental data [70] and also for the multivariate data matrix corresponding to oxygen depth and time profiles [72]. Homogeneity of multivariate dispersion (‘variance’) was tested with PERMDISP for any of the significant terms in PERMANOVA to verify any dispersion effects on the significant results from the PERMANOVA tests. For oxygen depth profiles, SIMPER analysis was used to determine which depth was responsible for the identified differences.
Microbial community analyses and diversity estimates.
To correct for technical bias related to read number variations among the samples, the data was normalised [73]. All samples were normalised by taking the proportions of each OTU, multiplying it with the minimum sample size (1022 reads) and rounding to the nearest integer [74,75].
To calculate alpha diversity and to assess how the sampled community of nosZ genes reflects a “true” diversity per sample [76], we made random sub-samples of the non-normalised data matrix to the minimum number of reads (1022) in R version 3.3.1. This was repeated 1000 times and the average diversity (richness, Shannon-Wiener [log e] and inverse Simpson) of each sample was calculated from the estimated diversity of each trial.
Rarefaction curves were also calculated in R for the number of observed OTUs vs. the total number of reads to investigate the variation of OTU numbers in our samples. Generalised UniFrac distances (α = 0.5) [77] were calculated on the composition of nosZ from normalised data with the GUniFrac package in R [78]. PERMANOVA was conducted on these UniFrac distances using the Adonis package in R. PERMDISP and pairwise difference tests were also performed in R.
Euclidean distance was used to calculate the resemblance matrices in PERMANOVA (Primer v6.1.10) for the diversity indices (estimated as average values).
PERMANOVA was performed as a two-way fixed factor model design (“treatment” and “depth” as factors while “replicate” was nested in “treatment”) to test significant differences among factors in composition and diversity indices of nosZ.
Principal coordinates analysis (PCoA) plots were generated using the Ade4 package in R to visualize variation in composition of nosZ among treatments and depths.
In addition, we investigated whether differences among factors were caused by differences in non-abundant OTUs (<1% relative abundance), by constructing a dataset with abundant OTUs (defined as OTUs with >1% relative abundance in at least one treatment-depth combination). This resulted in a second dataset with 21 abundant OTUs out of the total 502 OTUs and with 439624 reads. Statistical analyses on the dataset with the abundant OTUs were performed as described above.
Multivariate data corresponding to nosZ sequences were reduced to binary data (OTUs were counted as present or absent) in Excel [79]. Jaccard similarity index and Venn diagrams were calculated and drawn based on the number of shared and unique OTUs among treatments per depth layer [79,80].
Multiple linear regressions
Multiple linear regressions were performed to identify the variables contributing significantly to the observed variations in the diversity indices of nosZ (richness, Shannon diversity and Inverse Simpson as response variables) in our study. Multivariate data related to oxygen profiles, both depth and time profiles, were also included in the regression analyses but as two univariate variables: mean values of oxygen concentration per depth layer (“OX”) and coefficient of oxygen variations over time (“CV”), respectively. To avoid multi-collinearity among independent variables, the regression analyses were performed after removing highly correlated variables based on (i) the correlation coefficient (lrl ≥ 0.85) among numerical independent variables, and (ii) graphical exploratory techniques (boxplots) displaying the distribution of each numerical variable at different depths and treatments (categorical independent variables). No collinearity was observed among numerical variables while “depth” and “treatment” were collinear with chl-a and max OPD, respectively. Six explanatory variables were ultimately entered into the models: chl-a concentration, MGS, % silt+clay [mud] content, CV, max OPD and OX. All assumptions [81] were checked after model selection. An initial linear regression analysis showed violation of homogeneity of residuals in all models. We therefore applied linear regression with generalised least-squares (GLS) using restricted maximum-likelihood (REML) estimation. This incorporates variance–covariates to model the variance structure. To find the optimal random structure, the full linear regression model was compared to the GLS models of specific variance structures using lowest Akaike information criteria (AIC). This procedure resulted in the use of a variance structure that allowed for different variances per stratum for “max OPD” (varIdent function, R package nlme). The optimal fixed components in the final model were obtained by applying a backward selection using the likelihood ratio test obtained by ML estimation. A graphical model validation was applied to check for homogeneity and normality in the final models [82].
Results
Sediment environmental factors
The granulometric variables, including MGS and mud content (%), were not significantly influenced by interactive and single effects of treatments (high, low and no (= control) L. conchilega) and depths (PERMANOVA, p > 0.05). chl-a concentration was significantly affected by depth (PERMANOVA, pseudo-F = 12.455, p = 0.000). Highest chl-a concentrations were observed at the top layer (0–0.5 cm) (pairwise test, all 0.001 < p < 0.01) followed by the second layer which was significantly higher than the deepest layer (2.5–3 cm) (pairwise test, p = 0.014; Table 1).
Oxygen depth profiles (Fig 2) were significantly affected by the treatment effect (PERMANOVA, pseudo-F = 5.33, p = 0.006) (S2 Table). Pair-wise tests (S2 Table) revealed significant differences in vertical oxygen distribution between control and both L. conchilega treatments (p = 0.002 and p = 0.039 for high and low L. conchilega density, respectively). There were no significant differences between both L. conchilega treatments. SIMPER analysis showed that the difference in oxygen profiles between the control treatment and high L. conchilega treatment was attributed to 1–4 mm sediment depth whereas this difference was between 1–3.25 mm for the control and low L. conchilega treatments.
Control (A), low (B) and high L. conchilega (C) treatments. Nine profiles were measured per treatment (shown with different symbols and colours per core replicate). The zero depth in the profiles represents the sediment–water interface.
Oxygen concentrations over time (over a period of 30–35 min) measured at 1 cm distance from the tube and at 1.5 mm sediment depth remained high in the high L. conchilega (93–235 μmol L-1) and in the low L. conchilega treatment (60–203 μmol L-1). The lowest values were observed in the control treatment (80–157 μmol L-1) (Fig 3). However, oxygen time profiles were not significantly affected by treatment (p > 0.05; S2 Table).
Measured in 30–35 min at 1.5 mm depth in control (A), low (B) and high L. conchilega (C) treatments. Three profiles were measured per treatment. “CV” indicates coefficient of oxygen variations over time.
Maximum OPD was significantly affected by treatment (pseudo-F = 12.19; p = 0.000) and pair-wise tests showed a significantly higher max OPD in the high L. conchilega density treatment (5.17 ± 0.44 mm; mean ± se) compared to the low (3.69 ± 0.17 mm; p(MC) = 0.008) and control (3.30 ± 0.12 mm; p(MC) = 0.001) treatments (S2 Table; Fig 2).
PERMDISP analyses showed some heterogeneity of multivariate dispersion (F = 3.16, p = 0.028) in the analysis of maximum OPD but this was not the case for the analysis of the entire oxygen depth profiles (p > 0.05).
Composition of nosZ and diversity indices
In total, 545,742 nosZ reads were assigned to 502 OTUs. The current sampling effort was sufficient in most samples to reflect the variation of the observed OTUs (Fig A in S1 Fig). Only samples of the Hd4 treatment in all three replicates (H1d4, H2d4 and H3d4) did not reach an asymptotic rarefaction curve. Instead, read numbers of the high L. conchilega treatment decreased with increasing depth (lowest read numbers in Hd4; see S3 Table). This pattern was consistent across replicate samples, suggesting that the low read numbers are not a technical artifact but probably a reflection of the lower read numbers of OTU in this treatment.
Only 21 of 502 OTUs were found to be abundant (>1% relative abundance in at least one treatment-depth combination). Of these abundant OTUs, the two dominant OTUs 1 and 2 combined constituted 50–77% of the community where OTU 1 alone accounted for 44–71% of the relative abundance in different samples (Fig 4).
Abundant OTUs were defined as OTUs with >1% relative abundance in at least one treatment-depth combination (values are averages of replicates) and non-abundant OTUs defined as OTUs with <1%. “H” indicates high L. conchilega treatment. “L”: low L. conchilega treatment, “C”: control treatment, “d1”: 0–0.5 cm depth, “d2”: 0.5–1 cm, “d3”: 1–1.5 cm, “d4”: 2.5–3 cm.
Numbers of reads and total numbers of OTUs per sample as well as the numbers of abundant OTUs (>1% relative abundance) in treatment-depth combinations are provided in S3 Table.
The results of the protein blast indicate that among the abundant OTUs, OTUs 1, 2, 4 and 10 were most closely (>97% AA identity) related to nosZ sequences from Pseudomonas stutzeri constituting 76% of the community. The other abundant OTUs were less than 95% similar to protein sequences available in sequence databases. In addition, a blast search was performed of all 502 nosZ OTUs against NCBI’s non-redundant protein database (nr) with an e-value threshold of 1E-20 and taking into account only the first blast hit for each OTU. These results showed that 380 OTUs (which corresponds to 93.3% of all reads) were most closely related to Pseudomonas stutzeri.
The maximum likelihood tree of nosZ sequences retrieved from our samples showed considerable phylogenetic diversity (S2 Fig). No clear phylogenetic signal could be observed with respect to treatment, i.e. samples from the three treatments (control, low and high L. conchilega) contained sequences from across the phylogenetic tree (S2 Fig). Phylogenetic analysis of representative OTUs together with related sequences derived from blastp searches showed that there were three major clusters belonging to the Alpha, Beta and Gamma classes of the Proteobacteria and abundant as well as non-abundant OTUs are dispersed among these clusters (S3 Fig).
PERMANOVA based on Generalized Unifrac distance showed that the composition of nosZ gene with and without including non-abundant OTUs was significantly affected by the single effect of treatment and depth, regardless of the dataset used (Table 2).
Pairwise tests (S4 Table) indicated that the composition of nosZ in the treatment with high density of L. conchilega was significantly different from the control treatment (p = 0.036). However, no difference in nosZ gene composition was observed between the low and high L. conchilega treatments, and also between low L. conchilega and control treatments (p > 0.05).
In the vertical depth distribution, pairwise tests revealed that the composition of nosZ at the top layer (0–0.5 cm) was significantly different from the other layers (all pairwise test p < 0.05; S4 Table).
The highly significant PERMDISP results in both analyses with all OTUs (factor “treatment”: F = 10.07, p = 0.000; factor “depth”: F = 8.95, p = 0.000) and only abundant OTUs (factor “treatment”: F = 10.93, p = 0.000; factor “depth”: F = 7.56, p = 0.000) indicated a dispersion effect as well. Visualization by PCoA plot (Fig 5 and S4 Fig) showed that dispersion increased with increasing L. conchilega densities and with increasing depth in the sediment.
A: per treatment; B: per depth. Each point represents a sample. Three treatments (High: high L. conchilega treatment, Low: low L. conchilega treatment, Control) and four depth layers (0–0.5, 0.5–1, 1–1.5 and 2.5–3 cm) are shown.
The Venn diagram and Jaccard similarity index (Fig 6) indicated the highest similarity in OTUs among treatments at the top layer (0–0.5 cm) with more than 91% similarity. With increasing depth, the number of shared OTUs was generally reduced between the high L. conchilega density and the two other treatments, reaching the lowest similarity (<50%) at the deepest sediment layer (2.5–3 cm). However, 231 out of 247 (93.5%) and 234 out of 247 (94.7%) of OTUs in Hd4 were shared with Cd4 and Ld4, respectively (Fig 6). Among different depths per treatment, the lowest value (47%) of the Jaccard similarity index was observed when comparing the top layer (0–0.5 cm) and the 2.5–3 cm layer of the high L. conchilega treatment (Table 3).
Treatments: high L. conchilega treatment (H), low L. conchilega treatment (L), Control (C); depth: 0–0.5 (d1), 0.5–1 (d2), 1–1.5 (d3) and 2.5–3 cm (d4).
PERMANOVA did not show any significant effect on Shannon diversity and inverse Simpson for the factors “depth”, “treatment” or the interaction (all p > 0.05) although higher values were observed in the high L. conchilega treatment at highest depth (S5 Table; Fig 7). On the other hand, OTU richness was affected by the interaction effect “treatment × depth” (PERMANOVA, pseudo-F = 5.45, p = 0.043; S5 Table). Pairwise tests (S6 Table; Fig 7) did not reveal significant differences in OTU richness among treatments at the top layer (0–0.5 cm) (p(MC) > 0.05). In deeper sediment layers (2.5–3 cm), lower values were detected in the high L. conchilega treatment (pairwise tests, p(MC) = 0.015 and p(MC) = 0.016 compared with the control and low L. conchilega treatments). In addition, richness was significantly lower in the high L. conchilega compared to the control at the depth 0.5–1 cm (pairwise test, p(MC) = 0.026).
Obtained from 1000 sub-samples of the data matrix to the minimum number of reads at three treatments (high L. conchilega, low L. conchilega and control treatments) and four depth layers (0–0.5 [d1], 0.5–1 [d2], 1–1.5 [d3] and 2.5–3 [d4] cm). The solid horizontal line shows the median. Box is drawn between the quartiles (the middle 50% of the data). Dotted lines extend to the minimum and maximum values.
We generally found a decreasing trend in OTU richness in the high L. conchilega treatment with increasing depth (Fig 7). These changes were significant between the top (0–0.5 cm) and deepest (2.5–3 cm) layer (pairwise test, p(MC) = 0.022; S6 Table).
Statistical modeling to link abiotic sediment characteristics with the diversity indices of nosZ gene in L. conchilega reefs
Maximum oxygen penetration depth (max OPD; collinear with “treatment”) was an important variable incorporated in all diversity models. However, it was not significant in the model for the inverse Simpson index (p > 0.05). We found significant negative relation (β = -8.10x10-3, p < 0.01) between max OPD and OTU richness; and a positive relation between this variable and Shannon diversity (β = 1.96 × 10−4, p < 0.01) (Table 4). The average oxygen concentrations per depth layer (OX) and % mud content were further important explanatory variables in the model of Shannon diversity. Shannon diversity was negatively related with both OX (β = -3.19 × 10−3) and % mud content (β = -0.06). Apart from max OPD, chl-a concentration contributed significantly (β = 7.24, p < 0.001) to the model of richness.
Discussion
In the current study, we investigated the denitrifying community in L. conchilega aggregations at the intertidal zone by focusing on expression of the gene encoding nitrous-oxide reductase (nosZ), the enzyme catalysing the final step of denitrification. We investigated typical nosZ genes, which play the major role in N2O reduction in coastal ecosystems [25]. The denitrifiers with a complete pathway were identified as the dominant community [25].
In comparison with previous studies on functional genes for denitrification (nirS, nirK and nosZ) in marine sediments [12,18,83,84], our results show a diverse active nosZ community (502 OTUs) in the L. conchilega reef (area of L. conchilega aggregations plus zones without L. conchilega as a control). This is what can be expected in intertidal sediments where regular oxygen oscillations support a high genetic potential for denitrification [25].
The highest phylogenetic diversity was found in the Alphaproteobacteria and a high percentage of the sequences were closely related to Pseudomonas stutzeri (class Gammaproteobacteria). This is in agreement with previous studies in shelf sediments [25,85] and shallow coastal sediments [12,86,87]. In addition, the high relative abundance of Gammaproteobacteria-related phylotypes in our RNA-derived libraries suggests high levels of metabolic activity within this group [86].
Horizontal and vertical patterns of the active nosZ community in L. conchilega aggregations
Our study revealed differences in the composition of nosZ at both horizontal (m) and vertical (cm) scales between high L. conchilega and control treatments, and between top oxic (0–0.5 cm) and deeper layers. These differences in composition were mainly due to changes in the relative contributions of abundant OTUs to the community. The same pattern was also observed in the nosZ community structure in the high L. conchilega treatment at a depth of 2.5–3 cm (Fig 6). Statistical modelling of diversity measures suggested that biogeochemical heterogeneity (oxygen concentrations and oxygen penetration depth) created by the activity of high L. conchilega aggregations is the driving factor behind the observed differences, next to chl-a and mud content (%) (Table 4). This suggests a structuring effect of L. conchilega on the expression of nosZ genes in our study site. However, MGS and % mud were not affected by treatment nor depth, and chl-a was only affected by depth in the sediment. This is in agreement with the previous study in the same sampling location in Boulogne-sur-Mer and at the same sampling month in October [41]. However, temporal fluctuations in L. conchilega population density can significantly affect environmental properties [88]. Therefore, autogenic effects of L. conchilega on the denitrifying communities cannot be completely excluded. In our case study, however, there was a mainly allogenic effect of L. conchilega on the environment.
Allogenic effect of L. conchilega
Lanice conchilega bio-irrigates its tube through piston-pumping activity for 1.5 min generally every 4 min, thereby transporting about 3 mmol O2 m-2 d-1 into the sediment [31]. Lanice conchilega introduces oxygen-rich water in layers where oxygen is absent [31]. This intermittent ventilation results in temporal changes in oxygen concentrations in the upper layers of the sediment [31] and a deeper OPD especially close to the L. conchilega tubes [31,43]. Our results reveal a density-dependent effect of L. conchilega bio-irrigation on the OPD (Fig 2) as well as a significantly higher oxygen content at the top sediment layer (mainly 1–4 mm) in high and low L. conchilega treatments compared to the control sediment (Fig 2). In addition, the periodic piston-pumping activity resulted in oxygen variations over time. The oscillations were generally higher in both L. conchilega treatments compared to the control treatment (Fig 3), although there is no statistically significant difference. This could be caused by the fact that for some measurements, variations in oxygen concentration (the values of CV) were comparable with those in the control treatment. This can be explained by merging pulses from piston pumping activity of adjacent animals or by merging individual pulses in one animal resulting from a high frequency of occurrence of this pumping activity [31].
Generally low richness and high diversity (Shannon and inverse Simpson) of nosZ was observed in the high L. conchilega treatment compared with the control treatment at depth with most prominent differences at the deepest sediment layer (2.5–3 cm). At this depth, a significant difference between the high and low L. conchilega treatments was also observed. In a vertical sediment depth profile, nitrate coexists with oxygen in the top layer. In the deeper layers where oxygen is low or absent, the nitrate concentration is an important factor determining similarity or difference in anaerobic denitrifying communities [89,90]. As anoxic conditions continue at depth, nitrate concentrations decrease [8] while denitrification progresses. This trend coincides with the decrease in diversity of denitrifying organisms and the rate of their activity with increasing sediment depth [8,91,92]. The differences observed in diversity indices of nosZ between treatments at depth may thus be explained by the different availability of NO3− between treatments in the low oxygen or anaerobic deeper layers. The higher availability of NO3− in the high L. conchilega treatment can be by directly pumping NO3− [47] and/or by bringing O2 to deeper layers [31] increasing the surface for coupled nitrification–denitrification in the sediment along the tube [29] and also by more intensive irrigation activity in this treatment, compared to the low L. conchilega and control treatments. The most prominent differences between the high L. conchilega and the other two treatments were at the deepest sediment layer (2.5–3 cm). These differences are not due to adaptation of certain OTUs specific to the high L. conchilega treatment: 93.5% and 94.7% of OTUs in Hd4 were shared with Cd4 and Ld4, respectively (Fig 6). Instead, observed differences between treatments are due to reduction in the number of non-abundant and rare OTUs in Hd4 compared with Cd4 and Ld4.
Lanice conchilega ventilates its tube in each pumping activity by a volume of water equivalent to 2.4 cm tube lengths [31]. Therefore, the most suitable condition for denitrification (concomitant high availability of NO3− and little or no oxygen in the surrounding sediment) was met at 2.5–3 cm sediment. In the upper layers, high richness values resulted from the presence of a larger number of non-abundant and rare OTUs mainly belonging to the Alphaproteobacteria. This highly diverse pool may provide a background reservoir for suboptimal denitrification conditions increasing niche creation [90,93,94]. In addition, considering denitrification is an anaerobic process, the presence of active nosZ community in the top oxic layer in all treatments (Fig 7) may point to co-consuming oxygen and NO3− as electron acceptors (aerobic denitrification [10]) by denitrifying organisms in the sediment of our study area. This process has also been observed before in marine sediments [10].
Finally, the higher variations in composition of nosZ in the L. conchilega treatments compared with the control (Fig 5A) are in accordance with a previous study indicating higher variations in denitrification rates in bioturbated sediment [51]. In bioturbated sediments, nitrogen cycle processes are most likely driven by the formation of microniches surrounding the burrow systems [51]. In our study, this pattern can also be visualised in “depth” with increasing depth (Fig 5B).
Conclusions
This study improves our understanding of the effects of bio-irrigation on the nosZ gene. Our results reveal that large densities of L. conchilega (>3000 ind. m-2) affect composition, structure and diversity indices of the active nosZ community both vertically and horizontally in the sediment. These patterns can be linked with heterogeneity in geochemical environments in the sediment that are affected by allogenic ecosystem engineering effects of L. conchilega. The piston-pumping activity of the polychaete transports oxygen rich waters to the deeper layers of the sediment. This could also provide higher availability of NO3− at depth by increasing the surface for coupled nitrification-denitrification in the sediment along the tube.
Since L. conchilega is present throughout the North Sea, often in dense aggregations, our results provide evidence for the implementation of this polychaete in marine Ecosystem-Based Management. This study illustrates that a future loss of macrobenthic diversity, and especially ecosystem engineers, will have important repercussions for benthic ecosystem functioning.
Supporting information
S1 Fig. Rarefaction curves of the nosZ gene sequeces.
Plotting the number of observed OTUs (unique AA nosZ sequences) as function of the number of sequences screened in samples with >4000 reads (A) and ≤4000 reads (B). “H” indicates high L. conchilega treatment. “L”: low L. conchilega treatment, “C”: control treatment, “d1”: 0–0.5 cm depth, “d2”: 0.5–1 cm, “d3”: 1–1.5 cm, “d4”: 2.5–3 cm.
https://doi.org/10.1371/journal.pone.0192391.s001
(TIF)
S2 Fig. Maximum likelihood phylogeny of nosZ amino acid sequences retrieved from our samples.
10% of OTUs (50 OTUs out of 502) making up around 84% of the total reads are shown. Abundant OTUs (>1% relative abundance in at least one treatment–depth combination) are indicated in red. The heat map (on the right) illustrates the average relative abundance of each OTU per sample. The dominant OTU 1 is shown in red. Numbers at nodes are bootstrap values (values <60% not shown). The scale bar represents 10% sequence divergence (10 mutations per 100 sequence positions).
https://doi.org/10.1371/journal.pone.0192391.s002
(TIF)
S3 Fig. Maximum likelihood phylogeny of nosZ amino acid sequences retrieved from our samples and reference sequences.
To construct a phylogeny with reference sequences, protein BLAST searches of 25 OTU representatives selected randomly across the phylogenetic tree in S2 Fig were performed against the NCBI non-redundant protein database. Abundant OTUs (>1% relative abundance) are shown in red. Numbers at nodes are bootstrap values (values <50% not shown). The scale bar represents 30% sequence divergence (30 mutations per 100 sequence positions).
https://doi.org/10.1371/journal.pone.0192391.s003
(TIF)
S4 Fig. PCoA based on generalized UniFrac distances on the normalised data including only abundant OTUs.
A: per treatment; B: per depth. Each point represents a sample. Three treatments (High: high L. conchilega treatment, Low: low L. conchilega treatment, Control) and four depth layers (0–0.5, 0.5–1, 1–1.5 and 2.5–3 cm).
https://doi.org/10.1371/journal.pone.0192391.s004
(TIF)
S1 Table. Primers for paired-end nosZ sequencing on the Illumina Miseq platform.
Adaptor, Multiple Identifier (MID), pad, linker and primer sequences for the forward and reverse data sets are given.
https://doi.org/10.1371/journal.pone.0192391.s005
(DOCX)
S2 Table. Results from PERMANOVA analysis for differences in oxygen concentrations (depth and time profiles and max OPD) among treatments.
Treatments: high and low L. conchilega densities and control. Analyses were based on Euclidean distance similarity matrix. Where the treatment effect was significant, pairwise tests were performed within treatments. P-values were obtained by permutation and drawn from Monte-Carlo samplings (p(MC)) if the number of unique permutations was less than 100 [70]. Max OPD indicates maximum oxygen penetration depth. Bold values indicate significant differences at p < 0.05.
https://doi.org/10.1371/journal.pone.0192391.s006
(DOCX)
S3 Table. Numbers of reads and total numbers of OTUs per sample as well as the numbers of abundant OTUs (>1% relative abundance) in treatment-depth combinations are provided.
https://doi.org/10.1371/journal.pone.0192391.s007
(DOCX)
S4 Table. Results from PERMANOVA analysis pairwise tests for differences in composition of nosZ influenced by the single effects of “treatment” and “depth”.
Analyses were carried out based on Generalized UniFrac distances (α = 0.5) on the data sets using all OTUs or only abundant OTUs with relative abundance >1%. Treatments: high and low L. conchilega densities and control; Depths: 0–0.5, 0.5–1, 1–1.5 and 2.5–3 cm.
https://doi.org/10.1371/journal.pone.0192391.s008
(DOCX)
S5 Table. Results from PERMANOVA analysis main tests for differences in diversity indices of nosZ among treatments and depths.
Treatments: high and low L. conchilega densities and control; depths: 0–0.5, 0.5–1, 1–1.5 and 2.5–3 cm. Diversity indices (richness, Shannon-Wiener [log e] and inverse Simpson) were calculated from the average values obtained from 1000 sub-samples of the data matrix to the minimum number of reads (1022). Univariate analyses of diversity indices were based on Euclidean distance similarity matrix. P-values obtained by permutation.
https://doi.org/10.1371/journal.pone.0192391.s009
(DOCX)
S6 Table. Results from PERMANOVA analysis pairwise tests for differences in OTU richness influenced by the interaction effect “treatment × depth”.
Three treatments (High: high L. conchilega treatment, Low: low L. conchilega treatment, Control) and four depths (0–0.5, 0.5–1, 1–1.5 and 2.5–3 cm). Richness was calculated from the average values obtained from 1000 sub-samples of the data matrix to the minimum number of reads (1022). Analyses were based on Euclidean distance similarity matrix. P-values obtained by permutation.
https://doi.org/10.1371/journal.pone.0192391.s010
(DOCX)
Acknowledgments
We are very grateful to Bart Beuselinck for the analysis of granulometric variables, to Dirk Van Gansbeke for analysis chlorophyll a. We are very grateful to Mehrshad Taheri and Renata Alves Mamede da Silva for sampling assistance. We also thank Nele De Meester from Marine Biology Department, Ghent University for her assistance and great thanks to Frederiek-Maarten Kerckhof and Ruben Props from Department of Biochemical and Microbial Technology, Ghent University for their useful suggestions in some statistical analyses. We thank the anonymous reviewers for their constructive comments that improved the quality of the manuscript.
References
- 1.
Jenkins BD, Zehr JP. Molecular approaches to the nitrogen cycle. In: Capone DG, Bronk DA, Mulholland MR, Carpenter EJ, editors. Nitrogen in the marine environment. 2nd ed. Amsterdam: Elsevier; 2008. pp 1303–1344.
- 2.
Canfield DE, Thamdrup B, Kristensen E. Aquatic geomicrobiology: Advances in marine biology. California, USA: Elsevier Academic Press; 2005. Vol 48.
- 3. Knowles R. Denitrification. Microbiol Rev. 1982; 46: 43–70. pmid:7045624
- 4. Sanford RA, Wagner DD, Wu Q, Chee-Sanfordd JC, Thomasc SH, Cruz-Garcíac C, et al. Unexpected Nondenitrifier nitrous oxide reductase gene diversity and abundance in soils. PNAS. 2012; 109: 19709–19714. pmid:23150571
- 5.
Kohn RA. Nitrogen emissions from animal agricultural systems and strategies to protect the environment. In: Malik PK, Bhatta R, Takahashi J, Kohn RA, Prasad CS, editors. Livestock production and climate change. CABI Climate Change Series 6. Wallingford, UK: CABI Nosworthy Way; 2015. pp 61–72.
- 6. Deutsch B, Forster S, Wilhelm M, Dippner JW, Voss M. Denitrification in sediments as a major nitrogen sink in the Baltic Sea: an extrapolation using sediment characteristics. Biogeosciences. 2010; 7: 3259–3271.
- 7. Huang S, Chen C, Yang X, Wu Q, Zhang R. Distribution of typical denitrifying functional genes and diversity of the nirS-encoding bacterial community related to environmental characteristics of river sediments. Biogeosciences. 2011; 8: 5251–5280.
- 8. Tiquia SM, Masson SA, Devol A. Vertical distribution of nitrite reductase genes (nirS) in continental margin sediments of the Gulf of Mexico. FEMS Microbiol Ecol. 2006; 58: 464–475. pmid:17117989
- 9. Zehr J, Ward BB. Nitrogen cycling in the ocean. New perspective on processes and paradigms. Appl Environ Microbiol. 2002; 68: 1015–1024. pmid:11872445
- 10. Gao H, Schreiber F, Collins G, Jensen MM, Kostka JE, Lavik G, et al. Aerobic denitrification in permeable Wadden Sea sediments. ISME J. 2010; 4: 417–426. pmid:20010631
- 11. Li F, Yang R, Ti C, Lang M, Kimura SD, Yan X. Denitrification characteristics of pond sediments in a Chinese agricultural watershed. J Soil Sci Plant Nutr. 2010; 56: 66–71.
- 12. Zheng G, Wang X, Hannides AK, Sansone FJ, Wang GY. Impact of redox-stratification on the diversity and distribution of bacterial communities in sandy reef sediments in a microcosm. Chin J Oceanol Limn. 2011; 29: 1209–1223.
- 13. Kampschreur MJ, Kleerebezem R, Picioreanu C, Bakken L, Bergaust L, de Vries S, et al. Metabolic modeling of denitrification in Agrobacterium tumefaciens: a tool to study inhibiting and activating compounds for the denitrification pathway. Front Microbiol. 2012; 3: 370. pmid:23087683
- 14. Su JJ, Liu BY, Liu CY. Comparison of aerobic denitrification under high oxygen atmosphere by Thiosphaera pantotropha ATCC 35512 and Pseudomonas stutzeri SU2 newly isolated from the activated sludge of a piggery wastewater treatment system. J Appl Microbiol. 2001; 90 (3): 457–462. pmid:11298242
- 15. Scala DJ, Kerkhof LJ. Horizontal heterogeneity of denitrifying bacterial communities in marine sediments by terminal restriction fragment length polymorphism analysis. Appl Environ Microbiol. 2000; 66: 1980–1986. pmid:10788370
- 16. Yazdani Foshtomi M, Braeckman U, Derycke S, Sapp M, Van Gansbeke D, Sabbe K, et al. The link between microbial diversity and nitrogen cycling in marine sediments is modulated by macrofaunal bioturbation. PLoS ONE. 2015; 10(6): e0130116. pmid:26102286
- 17. Braker G, Ayala-del-Rio HL, Devol AH, Fesefeldt A, Tiedje JM. Community structure of denitrifiers, Bacteria, and Archaea along redox gradients in pacific northwest marine sediments by terminal restriction fragment length polymorphism analysis of amplified nitrite reductase (nirS) and 16S rRNA genes. Appl Environ Microbiol. 2001; 67: 1893–1901. pmid:11282647
- 18. Falk S, Hannig M, Gliesche C, Wardenga R, Köster M, Jürgens K, et al. nirS-containing denitrifier communities in the water column and sediment of the Baltic Sea. Biogeosciences. 2007; 4: 255–268.
- 19. Franco MA, De Mesel I, Demba Diallo M, Van der Gucht K, Van Gansbeke D, van Rijswijk P, et al. Effect of phytoplankton bloom deposition on benthic bacterial communities in two contrasting sediments in the southern North Sea. Aquat Microb Ecol. 2007; 48: 241–254.
- 20. Bertics VJ, Ziebis W. Biodiversity of benthic microbial communities in bioturbated coastal sediments is controlled by geochemical microniches. ISME J. 2009; 3(11): 1269–1285. pmid:19458658
- 21. Birchenough SNR, Parker RE, McManus E, Barry J. Combining bioturbation and redox metrics: Potential tools for assessing seabed function. Ecol Indic. 2012; 12: 8–16.
- 22. Joye SB, Hollibaugh JT. Influence of sulfide inhibition of nitrification on nitrogen regeneration in sediments. Science. 1995; 270: 623–625.
- 23. Bergaust L, Shapleigh J, Frostegård Å, Bakken L. Transcription and activities of NOx reductases in Agrobacterium tumefaciens: the influence of nitrate, nitrite and oxygen availability. Environ Microbiol. 2008; 10(11): 3070–3081. pmid:18312398
- 24. Liu B, Mao Y, Bergaust L, Bakken LR, Frostegård Å. Strains in the genus Thauera exhibit remarkably different denitrification regulatory phenotypes. Environ Microbiol. 2013; 15: 2816–2828. pmid:23663391
- 25. Wittorf L, Bonilla-Rosso G, Jones CM, Backman O, Hulth S, Hallin S. Habitat Partitioning of Marine Benthic Denitrifier Communities in Response to Oxygen Availability. Environ Microbiol. 2016; 8(4): 486–492.
- 26. Kristensen E. Organic matter diagenesis at the oxic/anoxic interface in coastal marine sediments, with emphasis on the role of burrowing animals. Hydrobiologia. 2000; 426: 1–24.
- 27. Laverock B, Gilbert JA, Tait K, Osborn AM, Widdicombe S. Bioturbation: impact on the marine nitrogen cycle. Biochem Soc Trans. 2011; 39: 315–320. pmid:21265795
- 28. Laverock B, Tait K, Gilbert JA, Osborn AM, Widdicombe S. Impacts of bioturbation on temporal variation in bacterial and archaeal nitrogen-cycling gene abundance in coastal sediments. Environ Microbiol Rep. 2014; 6(1): 113–121. pmid:24596269
- 29.
Kristensen E, Kostka JE. Macrofaunal burrows and irrigation in marine sediment: microbiological and biogeochemical interactions. In: Kristensen E, Haese RR, Kostka JE, editors. Interactions between macro- and microorganisms in marine sediments. Washington: American Geophysical Union; 2005. pp. 125–158.
- 30. Papaspyrou S, Gregersen T, Kristensen E, Christensen B, Cox RP. Microbial reaction rates and bacterial communities in sediment surrounding burrows of two nereidid polychaetes (Nereis diversicolor and N. virens). Mar Biol. 2006; 148: 541–550.
- 31. Forster S, Graf G. Impact of irrigation on oxygen flux into the sediment: intermittent pumping by Callianassa subterranea and “piston-pumping” by Lanice conchilega. Mar Biol. 1995; 123 (2): 335–346.
- 32.
Degraer S, Wittoeck J, Appeltans W, Cooreman K, Deprez T, Hillewaert H, et al. The macrobenthos atlas of the Belgian Part of the North Sea. Brussels, Belgium: Belgian Science Policy; 2006.
- 33. Heuers J. A model on the distribution and abundance of the tube building polychaete Lanice conchilega (Pallas, 1766) in the intertidal of the Wadden Sea. Verh Ges Ökol. 1998; 28: 207–215.
- 34. Van Hoey G, Guilini K, Rabaut M, Vincx M, Degraer S. Ecological implications of the presence of the tube-building polychaete Lanice conchilega on soft-bottom benthic ecosystems. Mar Biol. 2008; 154: 1009–1019.
- 35.
Read G, Bellan G. Lanice conchilega. 2012. In: Read G, Fauchald K, editors. 2017. World Polychaeta database. World Register of Marine Species.http://www.marinespecies.org/aphia.php?p=taxdetails&id=131495
- 36. Rabaut M, Vincx M, Degraer S. Do Lanice conchilega (sandmason) aggregations classify as reefs? Quantifying habitat modifying effects. Helgol Mar Res. 2009; 63(1): 37–46.
- 37. Rabaut M, Guilini K, Van Hoey G, Vincx M, Degraer S. A bio-engineered soft-bottom environment: the impact of Lanice conchilega on the benthic species-specific densities and community structure. Estuar Coast Shelf Sci. 2007; 75(4): 525–536.
- 38. Rabaut M, Braeckman U, Hendrickx F, Vincx M, Degraer S. Experimental beam-trawling in Lanice conchilega reefs: impact on the associated fauna. Fish Res. 2008; 90(1–3): 209–216.
- 39. Godet L, Toupoint N, Olivier F, Fournier J, Retiere C. Considering the functional value of common marine species as a conservation stake: the case of sandmason worm Lanice conchilega (Pallas 1766) (Annelida, Polychaeta) beds. Ambio. 2008; 37: 347–355. pmid:18828280
- 40. De Smet B, Braeckman U, Soetaert K, Vincx M, Vanaverbeke J. Predator effects on the feeding and bioirrigation activity of ecosystem-engineered Lanice conchilega reefs. J Exp Mar Bio Ecol. 2016; 475: 31–37.
- 41. De Smet B, D’Hondt A-S, Verhelst P, Fournier J, Godet L, Desroy N, et al. Biogenic reefs affect multiple components of intertidal soft-bottom benthic assemblages: the Lanice conchilega case study. Estuar Coast Shelf Sci. 2015; 152: 44–55.
- 42. Callaway R. Tube worms promote community change. Mar Ecol Prog Ser. 2006; 308: 49–60.
- 43. Braeckman U, Provoost P, Gribsholt B, Van Gansbeke D, Middelburg JJ, Soetaert K. Role of macrofauna functional traits and density in biogeochemical fluxes and bioturbation. Mar Ecol Prog Ser. 2010; 399: 173–186.
- 44. Philippot L, Hallin S. Finding the missing link between diversity and activity using denitrifying bacteria as a model functional community. Curr Opin Microbiol. 2005; 8: 234–239. pmid:15939345
- 45. Wang Z, Zhang X-X, Lu X, Liu B, Li Y, Long C, Li A. Abundance and diversity of bacterial nitrifiers and denitrifiers and their functional genes in tannery wastewater treatment plants revealed by high-throughput sequencing. PLoS ONE. 2014; 9(11): e113603. pmid:25420093
- 46. Graf DRH, Jones CM, Hallin S. Intergenomic comparisons highlight modularity of the denitrification pathway and underpin the importance of community structure for N2O emissions. PLoS ONE. 2014; 9(12): e114118. pmid:25436772
- 47. Gilbert F, Stora G, Bonin P, Le Dréau Y, Mille G, Bertrand J-C. Hydrocarbon influence on denitrification in bioturbated Mediterranean coastal sediments. Hydrobiologia. 1997; 345: 67–77.
- 48. Howe RL, Rees AP, Widdicombe S. The impact of two species of bioturbating shrimp (Callianassa subterranea and Upogebia deltaura) on sediment denitrification. J Mar Biol Assoc UK. 2004; 84: 629–632.
- 49. Webb AP, Eyre BD. Effect of natural populations of burrowing thalassinidean shrimp on sediment irrigation, benthic metabolism, nutrient fluxes and denitrification. Mar Ecol Prog Ser. 2004; 268: 205–220.
- 50. Nizzoli D, Welsh DT, Fano EA, Viaroli P. Impact of clam and mussel farming on benthic metabolism and nitrogen cycling, with emphasis on nitrate reduction pathways. Mar Ecol Prog Ser. 2006; 315: 151–165.
- 51. Bertics VJ, Sohm JA, Magnabosco C, Ziebis W. Denitrification and nitrogen fixation dynamics in the area surrounding an individual ghost shrimp (Neotrypaea californiensis) burrow system. Appl Environ Microbiol. 2012; 78(11): 3864–3872. pmid:22447588
- 52. Fernandes SO, Michotey VD, Guasco S, Bonin PC, Loka Bharathi PA. Denitrification prevails over anammox in tropical mangrove sediments (Goa, India). Mar Environ Res. 2012; 74: 9–19. pmid:22197479
- 53. Laverock B, Tait K, Gilbert JA, Osborn AM, Widdicombe S. Impacts of bioturbation on temporal variation in bacterial and archaeal nitrogen-cycling gene abundance in coastal sediments. Environ Microbiol Rep. 2014; 6(1): 113–121. pmid:24596269
- 54. Poulsen M, Kofoed MV, Larsen LH, Schramm A, Stief P. Chironomus plumosus larvae increase fluxes of denitrification productsand diversity of nitrate-reducing bacteria in freshwater sediment. Syst Appl Microbiol. 2014; 37: 51–59. pmid:24054696
- 55. Van Hoey G, Vincx M, Degraer S. Some recommendations for an accurate estimation of Lanice conchilega density based on tube counts. Helgol Mar Res. 2006; 60: 317–321.
- 56.
Marine Species Identification Portal, Developed by ETI BioInformatics. http://species-identification.org/
- 57. Queirós AM, Birchenough SNR, Bremner J, Godbold JA, Parker RE, Romero-Ramirez A, et al. A bioturbation classification of European marine infaunal invertebrates. Ecol Evol. 2013; 3: 3958–3985. pmid:24198953
- 58.
Wright SW, Jeffrey SW. High-resolution HPLC system for chlorophylls and carotenoids of marine phytoplankton. In: Jeffrey SW, Mantoura RFC, Wright SW, editors. Phytoplankton pigments in oceanography: guidelines to modern methods. Paris: UNESCO; 1997. pp 327–341.
- 59. Caporaso JG, Lauber CL, Walters WA, Berg-Lyons D, Huntley J, Fierer N, et al. Ultra-high-throughput microbial community analysis on the Illumina HiSeq and MiSeq platforms. ISME J. 2012; 6: 1621–1624. pmid:22402401
- 60. Bai Y, Shi Q, Wen D, Li Z, Jefferson WA, Feng C, et al. Bacterial communities in the sediments of Dianchi Lake, a partitioned eutrophic waterbody in China. PLoS ONE. 2012; 7(5): e37796. pmid:22666393
- 61. Zhang Y, Xie X, Jiao N, Hsiao SS-Y, Kao S-J. Diversity and distribution of amoA-type nitrifying and nirS-type denitrifying microbial communities in the Yangtze River Estuary. Biogeosciences. 2014; 11: 2131–2145.
- 62.
Marconi M, Rodriguez-Romero J, Sesma A, Wilkinson MD. Bioinformatics tools for next-generation RNA sequencing analysis. In: Sesma A, von der Haar T, editors. Fungal RNA biology. Switzerland: Springer; 2014. pp 371–391.
- 63. Edgar RC, Haas BJ, Clemente JC, Quince C, Knight R. UCHIME improves sensitivity and speed of chimera detection. Bioinformatics. 2011; 27(16): 2194–2200. pmid:21700674
- 64. Fish JA, Chai B, Wang Q, Sun Y, Brown CT, Tiedje JM, et al. FunGene: the functional gene pipeline and repository. Front Microbiol. 2013; 4: 291. pmid:24101916
- 65. Wang Y-F, Feng Y-Y, Ma X, Gu J-D. Seasonal dynamics of ammonia/ammonium-oxidizing prokaryotes in oxic and anoxic wetland sediments of subtropical coastal mangrove. Appl Microbiol Biotechnol. 2013; 97(17): 7919–7934. pmid:23099914
- 66. Bokulich NA, Joseph CML, Allen G, Benson AK, Mills DA. Next-Generation Sequencing Reveals Significant Bacterial Diversity of Botrytized Wine. PLoS ONE. 2012; 7(5): e36357. pmid:22563494
- 67. Kartzinel TR, Chen PA, Coverdale TC, Erickson DL, Kress WJ, Kuzmina ML, et al. DNA metabarcoding illuminates dietary niche partitioning by African large herbivores. PNAS. 2015; 112: 8019–8024. pmid:26034267
- 68.
Miller MA, Pfeiffer W, Schwartz T. Creating the CIPRES science gateway for inference of large phylogenetic trees. In: Proceedings of the Gateway Computing Environments Workshop (GCE). New Orleans, LA, USA: IEEE; 2010.
- 69. Abascal F, Zardoya R, Posad D. ProtTest: selection of best-fit models of protein evolution. Bioinformatics. 2005; 21 (9): 2104–2105. pmid:15647292
- 70.
Anderson MJ, Gorley RN, Clarke KR. PERMANOVA + for PRIMER: Guide to software and statistical methods. Plymouth, United Kingdom: PRIMER-E; 2008.
- 71. Braeckman U, Provoost P, Moens T, Soetaert K, Middelburg JJ, Vincx M, et al. Biological vs. physical mixing effects on benthic food web dynamics. PloS ONE. 2011; 6(3): e18078. pmid:21455308
- 72. Widdicombe S, Beesley A, Berge JA, Dashfield SL, McNeill CL, Needham HR, Øxnevad S. Impact of elevated levels of CO2 on animal mediated ecosystem function: The modification of sediment nutrient fluxes by burrowing urchins. Mar Pollut Bull. 2013; 73: 416–427. pmid:23218873
- 73. Dillies MA, Rau A, Aubert J, Hennequet-Antier C, Jeanmougin M, Servant N, et al. A comprehensive evaluation of normalization methods for Illumina high-throughput RNA sequencing data analysis. Brief Bioinform. 2012; 14(6): 671–683. pmid:22988256
- 74. McMurdie PJ, Holmes S. Waste not, Want not: Why rarefying microbiome data is inadmissible. PLoS Comput Biol. 2014; 10(4): e1003531. pmid:24699258
- 75. Props R, Kerckhof FM, Rubbens P, De Vrieze J, Hernandez Sanabria E, Waegeman W, et al. Absolute quantification of microbial taxon abundances. ISME J. 2016; 1–4.
- 76. Hughes JB, Hellmann JJ, Ricketts TH, Bohannan BJM. Minireview: counting the uncountable: statistical approaches to estimating microbial diversity. Appl Environ Microbiol. 2001; 67: 4399–4406.
- 77. Chen J, Bittinger K, Charlson ES, Hoffmann C, Lewis J, Wu GD, et al. Associating microbiome composition with environmental covariates using generalized UniFrac distances. Bioinformatics. 2012; 28: 2106–2113. pmid:22711789
- 78.
Team RDC. R: A language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing; 2008.
- 79. Fulthorpe RR, Roesch LFW, Riva A, Triplett EW. Distantly sampled soils carry few species in common. ISME J. 2008; 2: 901–910. pmid:18528413
- 80. Schloss PD, Handelsman J. Introducing SONS, a tool for operational taxonomic unit-based comparisons of microbial community memberships and structures. Appl Environ Microbiol. 2006; 72: 6773–6779. pmid:17021230
- 81. Zuur AF, Ieno EN, Elphick CS. A protocol for data exploration to avoid common statistical problems. Methods Ecol Evol. 2010; 1: 3–14.
- 82.
Zuur AF, Ieno EN, Walker NJ, Saveliev AA, Smith GM. Mixed effects models and extensions in ecology with R. New York: Springer; 2009.
- 83. Scala DJ, Kerkhof LJ. Diversity of nitrous oxide reductase (nosZ) genes in continental shelf sediments. Appl Environ Microbiol. 1999; 65: 1681–1687. pmid:10103268
- 84. Decleyre H, Heylen K, Sabbe K, Tytgat B, Deforce D, Van Nieuwerburgh F, et al. A doubling of Microphytobenthos biomass coincides with a tenfold increase in denitrifier and total bacterial abundances in intertidal sediments of a temperate estuary. PLoS ONE. 2015; 10(5): e0126583. pmid:25961719
- 85. Hunter EM, Mills HJ, Kostka JE. Microbial community diversity associated with carbon and nitrogen cycling in permeable shelf sediments. Appl Environ Microbiol. 2006; 72: 5689–5701. pmid:16957183
- 86. Mills HJ, Hunter E, Humphrys M, Kerkhof L, McGuinness L, Huettel M, Kostka JE. Characterization of nitrifying, denitrifying, and overall bacterial communities in permeable marine sediments of the northeastern Gulf of Mexico. Appl Environ Microbiol. 2008; 74(14): 4440–4453. pmid:18487394
- 87. Gihring TM, Humphrys M, Mills HJ, Huettel M, Kostka JE. Identification of phytodetritus-degrading microbial communities in sublittoral Gulf of Mexico sands. Limnol Oceanogr. 2009; 54: 1073–1083.
- 88. Alves RMS, Vanaverbeke J, Bouma TJ, Guarini J-M, Vincx M, Van Colen C. Effects of temporal fluctuation in population processes of intertidal Lanice conchilega (Pallas, 1766) aggregations on its ecosystem engineering. Estuar Coast Shelf Sci. 2017; 188: 88–98.
- 89. Liu X, Tiquia SM, Holguin G, Wu L, Nold SC, Devol AH, et al. Molecular diversity of denitrifying genes in continental margin sediments within the oxygen-deficient zone off the Pacific Coast of Mexico. Appl Environ Microbiol. 2003; 69(6): 3549–3560. pmid:12788762
- 90. Jayakumar A, O’Mullan GD, Naqvi SWA, Ward BB. Denitrifying bacterial community composition changes associated with stages of denitrification in oxygen minimum zones. Microb Ecol. 2009; 58: 350–362. pmid:19238477
- 91. Middelburg JJ, Levin LA. Coastal hypoxia and sediment biogeochemistry. Biogeosciences. 2009; 6: 1273–1293.
- 92. Jäntti H, Hietanen S. The effects of hypoxia on sediment nitrogen cycling in the Baltic Sea. Ambio. 2012; 41: 161–169. pmid:22246635
- 93. Tilman D, Lehman CL, Thomson KT. Plant diversity and ecosystem productivity: theoretical considerations. PNAS. 1997; 94: 1857–1861. pmid:11038606
- 94. Bent SJ, Forney LJ. The tragedy of the uncommon: understanding limitations in the analysis of microbial diversity. ISME J. 2008; 2(7): 689–695. pmid:18463690