ABSTRACT
Nitrogen flux into the coastal environment via submarine groundwater discharge may be modulated by microbial processes such as denitrification, but the spatial scales at which microbial communities act and vary are not well understood. In this study, we examined the denitrifying community within the beach aquifer at Huntington Beach, California, where high-nitrate groundwater is a persistent feature. Nitrite reductase-encoding gene fragments (nirK and nirS), responsible for the key step in the denitrification pathway, were PCR amplified, cloned, and sequenced from DNAs extracted from aquifer sediments collected along a cross-shore transect, where groundwater ranged in salinity from 8 to 34 practical salinity units and in nitrate concentration from 0.5 to 330 μM. We found taxonomically rich and novel communities, with all nirK clones exhibiting <85% identity and nirS clones exhibiting <92% identity at the amino acid level to those of cultivated denitrifiers and other environmental clones in the database. Unique communities were found at each site, despite being located within 40 m of each other, suggesting that the spatial scale at which denitrifier diversity and community composition vary is small. Statistical analyses of nir sequences using the Monte Carlo-based program ∫-Libshuff confirmed that some populations were indeed distinct, although further sequencing would be required to fully characterize the highly diverse denitrifying communities at this site.
Denitrification, the microbial dissimilatory reduction of nitrate and nitrite to gaseous products (NO, N2O, and N2) under suboxic conditions, has long been considered the main biological loss term for fixed nitrogen from ecosystems. In coastal and estuarine sediments, this process can remove more than half of the inorganic nitrogen inputs from terrestrial systems (27). One source of nitrogen for coastal environments is submarine groundwater discharge from unconfined aquifers, which is increasingly recognized as an important source of nutrients (7, 16). Groundwater below the beach face at Huntington Beach (HB), California, is enriched in nitrate, with reported nitrate concentrations in excess of 400 μM (3). The source of this nitrate is unknown, but potential sources include infiltration of fertilizer-contaminated surface water, microbial nitrification, or leaking subsurface sewage lines. Nitrate concentrations in the groundwater pool decrease 2 orders of magnitude in the cross-shore direction as the nitrate-rich brackish groundwater mixes with saline, nitrate-depleted seawater (Fig. 1). This mixing establishes ideal gradients in both nitrate and salinity over which to examine the geochemical parameters that may be controlling microbial diversity, especially the diversity of denitrifying bacteria.
Sediment sampling stations located along a cross-shore transect at Huntington Beach, California.
Because denitrifying microorganisms span a wide range of taxonomic groups, including over 50 different genera (32), 16S rRNA-based approaches are not generally useful for investigating natural denitrifying communities. Nitrite reduction to nitric oxide is the rate-limiting step in denitrification and is catalyzed by the metalloenzyme nitrite reductase, which occurs in two functionally equivalent but structurally different forms, i.e., the copper-containing NirK and cytochrome cd1 NirS enzymes. In recent years, the functional genes encoding these dissimilatory nitrite reductases, nirK and nirS, have been used extensively and effectively to elucidate denitrifier community structures in marine, estuarine, and groundwater environments (6, 18, 20, 30).
For this study, we examined how the diversity and composition of denitrifying bacteria vary across environmental gradients in both nitrate and salinity within an unconfined beach aquifer. Specifically, one of our primary objectives was to determine whether diversity positively correlated with the concentration of the terminal electron acceptor, nitrate. Because changes in nitrate and salinity over the study site occur over a short distance (∼10 m) relative to those often examined in traditional estuarine systems, we hypothesized that changes in the denitrifying community would exhibit a similarly steep gradient. We compared sequence data from our site to the growing database of environmental nirS and nirK sequences to determine if this unusual coastal environment harbors unique denitrifying genotypes. The Monte Carlo-based ∫-Libshuff program, previously only applied to 16S rRNA gene sequences (26), was used to determine whether patterns observed among functional denitrification gene clones in our samples were indicative of distinct microbial communities.
MATERIALS AND METHODS
Sample collection.Sediment and water samples were collected from five stations along a cross-shore transect on 22 and 23 June 2004 at Huntington State Beach, California (33°38.34′N, 117°58.59′W). Well points were installed at stations C, D, and E, and sediment samples were collected at the interface between the vadose zone and the water table just prior to well point installation (Fig. 1). For the remaining two stations (A and B), well points were not necessary and sediments were taken from just below the beach face as soon as the water table was reached, or in the case of the surf zone sample, surface sediments were collected. Five cubic centimeters of sediment was collected for each sample, using a modified polypropylene syringe. The well was flushed for approximately 2 min, and water samples were collected through the well point using a peristaltic pump and filtered (0.2 μm) into acid-washed containers for nitrate, nitrite, and ammonium analyses. Salinity and temperature were measured using a YSI 85 sonde. All samples were immediately stored on ice and flash frozen in liquid nitrogen within 30 min of collection. Samples were transported on dry ice before being transferred to −80°C for storage.
Nutrient analyses.An Alpkem autoanalyzer was used to measure nitrate, nitrite, and ammonia concentrations in undiluted water samples. We measured [NO3−-N], [NO2−-N], and [NH4+-N] (techniques P/N 000623 and P/N 000156; O/I Analytical). Our detection limit for each was 40 μg/liter.
DNA extraction, PCR amplification, cloning, and sequencing.DNAs were extracted from 0.5 g of sediment for each station, using a FastDNA spin kit for soil (Qbiogene). Purified DNA extracts were used as templates in 50-μl PCR mixtures. For nirS, the nirS1F and nirS6R primers and thermocycling conditions described by Braker et al. (5) were used, with a modified step-down protocol. After an initial denaturing step at 95°C for 5 min, 10 cycles of 95°C for 30 s, 54°C for 30 s, and 72°C for 1 min were carried out, with a 0.5°C step down in annealing temperature each cycle. This was followed by 25 cycles of 95°C for 30 s, 49°C for 30 s, and 72°C for 1 min and a final extension at 72°C for 7 min. The following reaction chemistry was used: 25 μl 2× PCR Premix E (Epicenter), a 0.5 μM concentration of each primer, and 1.25 U Taq polymerase. For nirK, a modified version of the previously described forward primer copper 583F (30) (5′-TCATGGTGCTGCCGCGYGANGG) was used along with the reverse primer (nirK5R or CuNIR4) (5, 9), with the following reaction chemistry: 5 μl 10× PCR buffer, 2 mM MgCl2, 2 μM forward primer, 1 μM reverse primer, a 200 μM concentration of each deoxynucleoside triphosphate, and 1.25 U Taq polymerase. Amplification products were visualized by electrophoresis in 1.5% (wt/vol) ethidium bromide-stained agarose gels. Prior to cloning, triplicate PCR products were pooled (to minimize PCR artifacts), and the products were separated by electrophoresis in 1.5% (wt/vol) ethidium bromide-stained agarose gels. Bands of the predicted product sizes (∼430 and ∼890 bp for nirK and nirS, respectively) were gel purified using a QIAquick gel extraction kit (QIAGEN).
Products were cloned using a TOPO-TA (Invitrogen) cloning kit with the pCR2.1 vector and TOP10 competent cells. Insert-bearing white clones were transferred to 96-well plates containing LB broth (with 50 μg/ml kanamycin) and PCR screened directly for inserts, using T7 and M13R vector primers. Forty-eight clones were randomly selected for sequencing with vector primers, using an ABI 3100 or ABI 3730 capillary sequencer (PE Applied Biosystems). The total number of sequences obtained from each clone library varied due to differences in the length and quality of the sequencing reads.
Sequence and phylogenetic analysis.Nucleotide sequences were aligned and edited using Sequencher v4.2 (Gene Codes) and translated using MacClade v4.06 (21), and the nucleotide alignments were manually corrected based on the resulting translation. Final amino acid alignments of sequences in the present study with database sequences were generated using ClustalX v1.81 (29). Neighbor-joining trees were produced for an ∼140 (nirK)- or ∼280 (nirS)-amino-acid segment by using PAUP 4.0b10 (Sinauer Associates). PAUP was also used to generate a Jukes-Cantor-corrected nucleic acid distance matrix for each site, which was used in subsequent community analyses. Bootstrapping (1,000 replicates) was used to estimate the reproducibility of the trees. The nearest database matches to Huntington Beach clones were determined using BLAST queries (http://www.ncbi.nlm.nih.gov) of the protein database.
Community richness and composition analysis.Operational taxonomic units (OTUs) for the purposes of community analysis were defined by a 5% difference in nucleic acid sequences, as determined using the furthest neighbor algorithm in DOTUR (25). The choice of a 5% cutoff represents what we feel is an appropriate balance between emphasizing small but population-defining genetic differences (highlighted by a small OTU cutoff) and the likelihood of functional differences in enzyme activity (highlighted by large OTU cutoffs) (15, 18). Rarefaction, richness, and diversity statistics were also calculated using DOTUR, including the nonparametric richness estimators ACE and Chao1 and the Shannon diversity index.
Nonparametric pairwise site similarity indices (Chao's estimated abundance-based Sørenson's index, or Labd) (10) were calculated using EstimateS (12). This method has been shown to alleviate some of the negative bias inherent in the classical Sørenson's index and is robust for data sets with unequal sample sizes such as the clone libraries presented in this study.
Statistical analysis of phylogenies.To assess whether observed differences in microbial community composition represented statistically significant community differences, well-aligned subsets of each gene fragment were chosen for analysis, using the ∫-Libshuff (26) program with 10,000 randomizations and a distance interval (D) of 0.01 on PAUP-generated Jukes-Cantor pairwise distance matrices. The software uses Monte Carlo methods to calculate the integral form of the Cramér-von Mise statistic by constructing random communities from the entire data set and comparing the coverage of the random communities to the observed coverage of experimental libraries. Significant P values (P < 0.0026; α = 0.05) were evaluated after correcting for multiple pairwise comparisons based on 20 tests [k(k − 1); k = 5 libraries]. In cases where coverage of the homologous library is significantly different from that of a randomly constructed heterologous library, we can be reasonably certain that the clone libraries represent separate populations (28).
Statistical analyses of environmental variables.Spearman rank, multivariate, and stepwise linear regression analyses of environmental variables (nitrate, ammonia, salinity, and temperature) against the richness estimators Chao1 and ACE were computed using MATLAB (MathWorks, Inc.) on both raw and log-transformed data.
Nucleotide sequence accession numbers.The nirK and nirS sequences reported in this study have been deposited in GenBank under accession numbers DQ159668 to DQ159860 (nirK) and DQ159473 to DQ159667 (nirS).
RESULTS AND DISCUSSION
Environmental gradients across the study site.Nitrate concentrations ranged from 0.5 μM at station A to 331 μM at station D, indicating a pronounced nitrate gradient across the study site (Table 1). Ammonium concentrations ranged from a low of 9.2 μM at station A to a maximum of 160 μM at station C. Although the ammonium concentrations were high compared to those in typical estuarine and marine water columns, they are not atypical of porewaters (8, 11, 17), which have reported concentrations in excess of 1 mM. Salinity increased in the cross-shore direction as nutrient concentrations decreased and ranged from 8 practical salinity units (psu) at station E to 33 psu at station A.
Locations and geochemical properties of microbial community sampling sites along a cross-shore transect at Huntington Beach, Orange County, California
nirK and nirS phylogeny.In order to explore the denitrifier diversity and community composition across environmental gradients within the study site, both nirK and nirS gene fragments were amplified from sediment DNA extracts obtained from all five stations. Clone libraries were generated for each station, and a total of 193 nirK and 195 nirS sequences were obtained (Table 2). Phylogenetic analyses of deduced amino acid sequences for both nirK (Fig. 2) and nirS (Fig. 3) gene fragments revealed distinct communities, with little overlap between stations.
Neighbor-joining phylogenetic tree of nirK gene products (partial, 140 amino acids), with Neisseria gonorrhoeae AniA as the outgroup. Bootstrap values (≥50%) for 1,000 replicates are shown at the branch points. Clones from the present study are shown in color by station (station A, blue; station B, aqua; station C, green; station D, orange; station E, red). Database sequences are shown in black, with GenBank accession numbers in parentheses. Roman numerals refer to clusters discussed in the text.
Neighbor-joining phylogenetic tree of nirS gene products with Paracoccus denitrificans as the outgroup. Bootstrap values (≥50%) for 1,000 replicates are shown at the branch points. Clones from the present study are shown in color by station (station A, blue; station B, aqua; station C, green; station D, orange; station E, red). Database sequences are shown in black, with GenBank accession numbers in parentheses. Roman numerals refer to clusters discussed in the text.
Diversity and predicted richness of nirK and nirS gene fragments from five subsurface samples along a nitrate and salinity gradient at Huntington Beach, California, as estimated by the Shannon diversity index and Chao1 and ACE richness estimators computed using DOTUR
nirK. Four distinct clades were apparent in the NirK phylogeny, with each clade corresponding to a particular salinity/nutrient regimen. Groups I and II, composed of >70% of the sequences from stations A and B, may represent a low-nitrate/high-salinity clade. Similarly, group IV is composed predominantly of clones from stations D and E and thus may represent a high-nitrate/low-salinity group. These groupings suggest that sequence differences in the nirK denitrification gene may reflect physiological adaptations to particular nutrient or salinity conditions. The remaining sequences fall into an intermediate Nitrosomonas-like cluster (group III) containing clones from every station. This cosmopolitan cluster may represent either nirK-containing organisms able to exploit many nutrient/salinity regimens or organisms that have been transported between stations by tidal mixing.
There was little similarity between nirK genotypes discovered in this study, cultured denitrifiers, or environmental database sequences (Fig. 2). The extensive number of unique nirK sequences recovered in this study may be partially attributable to the modified (i.e., more degenerate) forward primer, which was redesigned to amplify a broader range of nirK sequences, as previously published primers may be biased towards well-conserved genes (5). The database sequences with the nearest matches to any of the HB nirK sequences were from engineered systems and contaminated sites, rather than other marine sediment sequences from the Pacific Northwest and Mexico (4, 20). While this may reflect subsurface transport of organisms associated with sewage at our field site, it may also be an artifact of the relatively few marine nirK sequences in the database. The nearest BLAST hits were between sequences from station E and a Tokyo bioreactor clone (BAD17995, 84% amino acid identity) or an activated sludge clone (BAD00126, 85% amino acid identity). The closest match between a sequence from this study and the only other groundwater-derived nirK database sequences (30) was 77% identical at the amino acid level. Groundwater inputs to the surf zone below the beach face at this site have been correlated with periods of increased fecal indicator bacteria (3) and have been the focus of numerous studies investigating potential sewage contamination (2, 19). Therefore, it is not entirely surprising that the microbial community at Huntington Beach shows some similarity to wastewater treatment system communities.
nirS.The station-specific clustering present in the NirK phylogeny was not as apparent in the NirS tree. Two small clusters (Fig. 3, groups I and II), composed primarily (>80%) of clones from stations D and E, may correspond to genotypes specifically associated with (and/or adapted to) high-nitrate/low-salinity conditions. Group V is composed predominantly of clones from stations A and B and is the only other group with sequences specific to a particular environment. The remaining clusters each contain sequences from all stations. The more cosmopolitan nature of nirS sequences compared to nirK sequences suggests that either the nirS gene or, more likely, the organisms that carry the nirS gene do not exhibit as much environmental specialization, taken here to mean a specific association with specific biogeochemical conditions. Alternatively, nirS-containing organisms may rely on other metabolic pathways at this site; therefore, any selection for environment-specific nirS genes is not apparent.
In further contrast to the results for nirK, the nearest database matches to nirS sequences found in this study were primarily from marine and estuarine clone libraries (6, 22). However, the database of nirS sequences from marine and estuarine environments is much larger than that for nirK. Sequences from all clusters were, at most, 91% identical at the amino acid level to environmental clones from Puget Sound sediments (6) and 92% identical to sequences from the coastal Arabian Sea oxygen minimum zone (18). Little similarity was observed between HB groundwater nirS genotypes and the only other previously published groundwater sequences from Tennessee (30) (78% amino acid identity for the closest BLAST match), which is not surprising given the large differences in nitrate concentration and salinity between the two field sites.
nirK and nirS richness.To compare the relative nir-based richness between stations, rarefaction analysis was performed with the sequences from each station, using a 5% cutoff at the DNA level to define an OTU. The richness of both nirK and nirS OTUs was high at all stations, displaying steep rarefaction curves (Fig. 4). However, nirS richness was greater than nirK richness at each station, as illustrated by the rarefaction curves, two different nonparametric estimators of species richness, and the traditional Shannon diversity index (Table 2). According to the Chao1 estimator, richness was lowest at the marine, low-nitrate station A for both nirS and nirK; however, the ACE estimator indicated that richness was lowest at station E for nirS and station A for nirK. These findings are in agreement with those recently reported for nir diversity in a wastewater treatment plant, where diversity was inversely correlated with salinity (31). These results are also in agreement with broader trends in other functional genes involved in nitrification and denitrification, where the highest diversity is often observed at low-salinity, high-nutrient sites (1, 15, 23). Detailed comparisons among stations for these richness estimates must be made with some caution, however, as the estimates have not fully converged at these sample sizes.
Rarefaction curves for nirK (top) and nirS (bottom) gene fragments, as determined by sequence analysis. A 5% difference in nucleic acid sequence alignment was used to determine OTUs.
We were unable to find a clear relationship between measures of richness and any combination of measured environmental parameters by using either multivariate or stepwise linear regression. Previous studies have had similar difficulty relating nir richness and diversity to environmental parameters using both direct and multivariate analyses (30), indicating that our understanding of the factors controlling diversity in these functional genes is still inadequate. It does, however, seem likely that community composition is controlled by environmental attributes that are more difficult to quantify, such as microscale spatial heterogeneity and temporal variability of geochemical parameters. Nonequilibrium conditions, one of the hallmarks of coastal aquifers, have long been postulated to maintain high levels of diversity in ecological systems (13) and likely play an important role in shaping the microbial diversity observed in this study.
Community-based library comparisons.OTU-based similarity indices reinforce the discrete communities evident in the phylogenetic trees (Table 3). For nirK and nirS, there were no OTUs that were represented at every station. Stations nearest each other were the most similar. The two marine sites (stations A and B) shared no OTUs with the brackish groundwater site (station E) for nirK and only a single OTU for nirS. Sørenson's indices (Labd) ranged from 0.00 to 0.92 for nirK and from 0.04 to 0.57 for nirS at the five stations. While the lower Labd values for nirS (relative to nirK) may seem to contradict the results discussed earlier regarding the distinct nutrient/salinity clades observed in the nirK phylogenetic tree, these relatively low Labd values for nirS are heavily influenced by the fact that most OTUs are represented by a single sequence in our clone libraries. A previous study of nitrate- and uranium-contaminated groundwater (30) also found a higher similarity between sites for nirK than for nirS, using restriction fragment length polymorphism data and the traditional Sørenson's index (L = 0.19 to 0.49 for nirK and 0.04 to 0.25 for nirS). Reports of meter-scale similarity indices for nitrous oxide reductase (nosZ) genes, as determined by terminal restriction fragment length polymorphism analysis, were higher, on average, than those observed here (L = 0.67 to 0.76) (24), suggesting that nirS and nirK may exhibit greater degrees of environmental specialization or perhaps are simply more diverse. However, it is difficult to directly compare similarities computed with the classic Sørenson index to the nonparametric Labd values presented here.
Estimated abundance-based Sørenson's indices of similarity (Labd) between stations based on OTUs for nirK (above the diagonal) and nirS (below the diagonal) clone librariesa
To determine if the patterns observed in this study were the result of significant differences in the underlying microbial communities, we compared the coverage of our clone libraries to a hypothetical coverage generated from multiple randomizations of our sequences by calculating population P values (Table 4). If clone libraries from each station are different from each other, as determined by significant P values, then we can be reasonably certain that the patterns displayed by the clone libraries are representative of distinct populations. This method has the advantage of comparing communities on the individual sequence level rather than relying on an arbitrary definition of an OTU or species but has previously only been applied to 16S rRNA gene libraries (14, 26). If libraries X and Y are drawn from separate populations (shown in bold in Table 4), then comparisons of both X versus Y and Y versus X result in significant P values. In cases where only library X versus library Y produces a significant P value, but library Y versus library X does not, library Y is considered a subset of library X (shown in bold italics in Table 4).
Population P values for comparison of nir clone libraries determined using the Cramér-von Mises statistic, as implemented in ∫-LIBSHUFF
For nirK, the results of this analysis suggest that, as seen in the phylogenetic tree, the sequences in most libraries are station specific. Libraries from stations A and B are derived from a single population; libraries C and D are also from a single population which is distinct from that in libraries A and B. The station A library is distinct from those from stations C, D, and E. Though not statistically significant in all cases, small P values for comparisons between stations C, D, and E appear to represent libraries from distinct populations, providing support for the high degree of environmental specificity that appears to be associated with nirK-containing denitrifiers.
As suggested by the nirS phylogenetic tree, ∫-Libshuff analyses indicated that populations at each station exhibit higher degrees of overlap for this gene. The libraries from stations A and B represent the same population and are distinct from those from both stations D and E, while the library at station C is a subset of the libraries from all other stations. However, the latter may be an artifact of the smaller sample size at station C, as the estimated coverage for this clone library was only 0.33 (Table 2).
Conclusion.We demonstrated here that distinct populations can be resolved over small spatial scales with relatively small clone libraries by using iterative statistical methods. The results highlight that geochemistry plays an important role in shaping the bacterial assemblages at a given location; even in connected aquifers, steep gradients in environmental parameters can result in equally steep gradients (i.e., shifts) in community composition. This suggests that diversity in these enzymes results from adaptation to a particular environment rather than just random microheterogeneity alone. The links between the members of the community that are present and those that are active, as well as their relationship to denitrification rates, remain topics for further investigation. Further work at this site will also explore gradients in ammonia-oxidizing communities to determine whether other key players in the coastal nitrogen cycle exhibit similar site-specific patterns.
ACKNOWLEDGMENTS
We thank N. B. Handler for help in the field. P. B. Eckburg and N. J. Nidzieko provided invaluable help with statistical analyses. P. Jewett graciously assisted with nutrient analyses. G. D. O'Mullan, D. P. Keymer, J. M. Beman, and three anonymous reviewers provided helpful comments on the manuscript.
A.E.S. was supported by an NSF graduate research fellowship. A.B.B. was supported by the Clare Boothe Luce Professorship and the Powell Foundation. C.A.F. was supported in part by NSF grant MCB-0433804.
FOOTNOTES
- Received 11 August 2005.
- Accepted 28 December 2005.
- American Society for Microbiology
