Metagenomic analysis of intertidal hypersaline microbial mats from Elkhorn Slough, California, grown with and without molybdate

Abstract

Cyanobacterial mats are laminated microbial ecosystems which occur in highly diverse environments and which may provide a possible model for early life on Earth. Their ability to produce hydrogen also makes them of interest from a biotechnological and bioenergy perspective. Samples of an intertidal microbial mat from the Elkhorn Slough estuary in Monterey Bay, California, were transplanted to a greenhouse at NASA Ames Research Center to study a 24-h diel cycle, in the presence or absence of molybdate (which inhibits biohydrogen consumption by sulfate reducers). Here, we present metagenomic analyses of four samples that will be used as references for future metatranscriptomic analyses of this diel time series.

Introduction

Microbial mats are amongst the most diverse microbial ecosystems on Earth, inhabiting some of the most inhospitable environments known, including hypersaline, dry, hot, cold, nutrient poor, and high UV environments. Photosynthetic microbial mats found in intertidal environments are stratified microbial communities. Microbial metabolism under anoxic conditions at night results in the generation of significant amounts of H₂ and organic acids. The high microbial diversity of microbial mats makes possible a highly complex series of metabolic interactions between the microbes, the nature and extent of which are currently under investigation. To address this challenge, we are using a combination of metagenomics, metatranscriptomics, metaproteomics, iTags and naturally collected, as well as culture-based simplified microbial mats to study biogeochemical cycling (H₂ production, N₂ fixation, and fermentation) in mats collected from Elkhorn Slough, Monterey Bay, California. We present here the metagenome data, which will be used as a reference for metatranscriptomic analysis in a later paper.

Site information

Cyanobacterial mats are compact, laminated, and highly structured microbial communities (Fig. 1) that comprise great diversity at both the metabolic and phylogenetic level [1] and typically exist in highly saline environments such as lagoons and salterns. These mats notably have a suite of phototrophic organisms and photosynthetic lifestyles, from the dominant cyanobacterial phototroph Coleofasciculus chthonoplastes (basionym Microcoleus chthonoplastes) to purple sulfur and non-sulfur bacteria, and potentially other anoxygenic phototrophs. During the nighttime portion of the diel cycle, phototrophic organisms release fermentation byproducts which in turn help drive a shift from oxic to anoxic metabolism dominated by hydrogen consumption and sulfate reduction by sulfate reducing bacteria such as Desulfobacteriales [2]. Naturally occurring mats have a documented capacity to produce and liberate fermentation by-products (H2 and acetate primarily) [3, 4] and to consume them [5, 6] depending on the point in the diel cycle. Lastly, nitrogen assimilation is dominated by nitrogen fixation in these mats, typically by several members of the phylum Cyanobacteria such as ESFC-1 and Lyngbya sp. and by sulfate reducing bacteria [7,8,9,10,11]. The mats of Elkhorn Slough are situated in an estuary emptying into Monterey Bay, California and are located in a former salt production pond. The MIMS coding is shown in Table 1.

Fig. 1

a. Photograph of location of cores collected in the field from microbial mats at the Moss Landing Wildlife Area in Elkhorn Slough, Moss Landing, California on 07/11/11. Individual samples collected in core tubes were numbered and could be tracked throughout the diel experiment. b. Experimental apparatus used to incubate microbial mats throughout the diel period from 08/11/11 to 09/11/11. Incubation containers containing cores used for control and molybdate treatments are labeled

Table 1 Study information

Microbial mats like the ones at Elkhorn Slough have long been studied as a model for early life and gained prominence with the discovery that hypersaline mats in Guerrero Negro, Baja California, represented one of the most highly species-diverse microbiomes ever studied [1]. Though not as diverse as the Lyngbya mats of the Guerrero Negro system, the Elkhorn Slough mat system captures a similar distribution of organisms observed in laminated seasonal microbial ecosystems [6, 12]. Several areas of microbial mat physiology research are on-going at the Elkhorn Slough site. The site has been used to isolate a novel nitrogen fixer [9] and to show that the majority of fixation is attributable to a Lyngbya sp. [10], and to identify the dominant SRB ( Desulfobacterales ) in the ecosystem [2]. Additionally, the site has been investigated for hydrogen cycling. Burow and colleagues [5], showed that hydrogen flux likely originates from the fermentation of photosynthate. This system has also been subjected to metatranscriptomics and metaproteomics analyses [12, 13].

Metagenome sequencing information

Metagenome project history

Building on previous work examining gene expression patterns associated with fermentation pathways in microbial mat systems [12], a 24-h study of Elkhorn Slough, CA microbial mats was conducted in 2011. Briefly, field-collected mats were incubated at NASA Ames in seawater media and repeatedly sampled over one diel cycle. In addition, to understand gene expression across the diel cycle, DNA and RNA were extracted from molybdate and control samples for metagenome and metatranscriptome sequencing. Study information is summarized in Table 1.

Sample information

To understand the variation in gene expression associated with the daytime oxygenic phototrophic and nighttime fermentation regimes in hypersaline microbial mats, a contiguous mat piece was sampled at regular intervals over a 24-h diel period. Additionally, to understand the impact of sulfate reduction on biohydrogen consumption and impacts on community-structure, molybdate was added as an inhibitor to a parallel experiment. Contiguous mat samples were incubated and sampled at regular intervals throughout a 24-h period (8 time points). Four metagenome samples (two time points 12 h apart, from mats with and without molybdate added to the overlying water) and 13 metratranscriptomes (including nine time points for the control time series, four for the molybdate time series, and duplicates for most time points) were sequenced using Illumina technology.

Sample preparation

Microbial mats used in the experiment were collected using 3 in. acrylic core tubes on the morning of 07/11/11 and transported to Ames Research Center (about one hour by car). The mats were collected from a single contiguous section of mat (Fig. 1a) and were not covered with water at the time of collection (low tide). The microbial mats were immediately transferred to temperature controlled water baths on a rooftop facility [14] (Fig. 1b) containing either seawater or seawater amended with 30 mM (final concentration) sodium molybdate to inhibit the activities of sulfate reducing bacteria. The seawater used was obtained from the boat launch in the Moss Landing harbor at the time of collection of the mats. Two replicate containers each were used for mat incubations: 1) seawater alone and 2) seawater with molybdate water baths.

Mat samples for metagenomic analysis were subsampled from the acrylic core tubes using smaller metal coring tubes (having an area of 1.15 cm2, and a depth of 0.5 cm) on 09/11/11 at 01:30 h and 13:30 h (PST), corresponding to the 2nd and 6th time point in the larger diel time series (one control and one molybdate sample at each time point). Samples were placed in liquid nitrogen immediately after collection and, after flash freezing, were stored in a − 80 °C freezer for later extraction.

The four samples, and resulting metagenomes presented here will be referred to by a 4-character code: CD2A (Control, DNA, time point 2, replicate A), CD6A (Control, DNA, time point 6, replicate A), MD2A (Molybdate, DNA, time point 2, replicate A), MD6A (Molybdate, DNA, time point 6, replicate A). Sample information is provided in Table 2 as per minimal information standards [15].

Table 2 Sample information

DNA extraction

Nucleic acids were extracted from the samples between 2/2/2012 and 24/3/12. For each time point and treatment, the top 2–2.5 mm (photosynthetic layer) of 4 1-cm diameter cores were extracted by initially placing each core in 2 ml tubes containing a mixture of 0.5 ml of RLT buffer (RNeasy Mini Elute Cleanup Kit #74204; Qiagen, Valencia, CA, USA) and 5 μl of 2-mercaptoethanol (cat. # 0482–100) (Amresco, Solon, OH, USA). Samples were homogenized using a rotor-stator homogenizer (Omni International, Kennesaw, GA, USA), followed by the addition of 0.5 mm zirconium beads (OPS Diagnostics, Lebanon, NJ, USA) and then bead-beaten for 40 s using a FastPrep FP120 Cell Disrupter (Qbiogene, Inc., Carlsbad, CA, USA). Samples were spun down and the supernatant for each tube was transferred into a new tube containing an equal volume of phenol:chloroform:isoamyl alcohol (25:24:1) (cat. # 0883–400) (Amresco, Solon, OH, USA). Samples were vortexed, incubated for 5 min at room temperature, and spun down. The supernatant from each tube was transferred to a new tube containing an equal volume of 100% ethanol (Fisher #BP2818, Waltham, MA, USA) and was vortexed. Replicates of supernatant and ethanol mix for each time point and treatment were pooled, run through a QIAmp spin column (QIAmp DNA mini kit #51304, Qiagen, Valencia, CA, USA), and further purified according to the QIAmp DNA mini kit protocol. DNA quality and concentration were measured using a QUBIT fluorometer model Q32857 (Invitrogen, Carlsbad, CA, USA). Samples were submitted to JGI for sequencing.

Library generation

500 ng of genomic DNA (2 μg for sample MD2A) was sheared using the Covaris E210 (Covaris) and size selected using Agencourt Ampure Beads (Beckman Coulter). The DNA fragments were treated with end repair, A-tailing, and adapter ligation using the TruSeq DNA Sample Prep Kit (Illumina) and purified using Agencourt Ampure Beads (Beckman Coulter). The prepared libraries were quantified using KAPA Biosystem’s next-generation sequencing library qPCR kit and run on a Roche LightCycler 480 real-time PCR instrument. The quantified libraries were then prepared for sequencing on the Illumina HiSeq sequencing platform utilizing a TruSeq paired-end cluster kit, v3, and Illumina’s cBot instrument to generate a clustered flowcell for sequencing. The library information is summarized in Table 3.

Table 3 Library information

Sequencing technology

Sequencing of the flowcell was performed on the Illumina HiSeq2000 sequencer using a TruSeq SBS sequencing kit 200 cycles, v3, following a 2 × 150 indexed run recipe. All sequencing was performed by the Joint Genome Institute in Walnut Creek, CA, USA.

Sequence processing, annotation, and data analysis

Sequence processing

Raw Illumina metagenomic reads were screened against Illumina artifacts with a sliding window with a kmer size of 28, step size of 1. Screened reads were trimmed from both ends using a minimum quality cutoff of 3, reads with 3 or more N’s or with average quality score of less than Q20 were removed. In addition, reads with a minimum sequence length of <50 bps were removed. The sequence processing is summarized in Table 4.

Table 4 Sequence processing

Metagenome processing

Trimmed, screened, paired-end Illumina reads were assembled using SOAPdenovo v1.05 [16] at a range of Kmers (85, 89, 93, 97, 101, 105). Default settings for all SOAPdenovo assemblies were used (options "-K 81 -p 32 -R -d 1"). Contigs generated by each assembly (6 total contig sets), were de-replicated using in-house Perl scripts. Contigs were then sorted into two pools based on length. Contigs smaller than 1800 bp were assembled using Newbler [17] in attempt to generate larger contigs (flags: -tr, −rip, −mi 98, −ml 80). All assembled contigs larger than 1800 bp, as well as, the contigs generated from the final Newbler run were combined using minimus 2 (flags: -D MINID = 98 -D OVERLAP = 80) [18]. Read depths were estimated based on read mapping with BWA [19]. These sequences are currently available to the public at IMG/M and the JGI genome portals. Metagenome statistics are summarized in Table 5.

Table 5 Metagenome statistics

Metagenome annotation

Prior to annotation, all sequences were trimmed to remove low quality regions falling below a minimum quality of Q13, and stretches of undetermined sequences at the ends of contigs were removed. Low complexity regions were masked using the dust algorithm from the NCBI toolkit and very similar sequences (similarity >95%) with identical 5′ pentanucleotides were replaced by one representative, typically the longest, using uclust [20]. The gene prediction pipeline included the detection of non-coding RNA genes (tRNA and rRNA) and CRISPRs, followed by prediction of protein coding genes.

Identification of tRNAs was performed using tRNAScan-SE-1.23 [21]. In case of conflicting predictions, the best scoring predictions were selected. Since the program cannot detect fragmented tRNAs at the end of the sequences, we also checked the last 150 nt of the sequences by comparing these to a database of nt sequences of tRNAs identified in the isolate genomes using blastn [22]. Hits with high similarity were kept; all other parameters were set to default values. Ribosomal RNA genes were predicted using hmmsearch [23] with internally developed models for the three types of RNAs for the domains of life. Identification of CRISPR elements was performed using the programs CRT [24] and PILERCR [25]. The predictions from both programs were concatenated and, in case of overlapping predictions, the shorter prediction was removed.

Identification of protein-coding genes was performed using four different gene calling tools, GeneMark (v. 2.8) [26],Metagene (v. 1.0) [27], Prodigal (V2.50: November, 2010) [28] and FragGenescan (v. 1.16) [29] all of which are ab initio gene prediction programs. We typically followed a majority rule based decision scheme to select the gene calls. When there was a tie, we selected genes based on an order of gene callers determined by runs on simulated metagenomic datasets (Genemark > Prodigal > Metagene > FragGeneScan). At the last step, CDS and other feature predictions were consolidated. The regions identified previously as RNA genes and CRISPRs were preferred over protein-coding genes. Functional prediction followed and involved comparison of predicted protein sequences to the public IMG database using the usearch algorithm [20], the COG database using the NCBI developed PSSMs [30], the Pfam database [31] using hmmsearch. Assignment to KEGG Ortholog protein families was performed using the algorithm described in [32]. Annotation parameters are summarized in Table 6.

Table 6 Annotation parameters

Metagenome properties

Metagenomes were sequenced and assembled into 141,229 (CD6A) to 292,231 (MD2A) contigs, covering 90.6 to 173.6Mbp. GC content of the metagenomes ranged from 46% to 52%. These metagenomes include between 206,164 and 399,161 genes each. More than 99% of these are protein coding, and around 40% have some level of function annotation. Metagenome properties are summarized in Table 7.

Table 7 Metagenome properties

Taxonomic diversity

The taxonomic diversity and phylogenetic structure of the metagenomes was determined based on the best BLASTp hits of assembled protein-coding genes with 60% or more identity to protein in the listed phyla, as calculated by the Phylogenetic Distribution of Genes feature in IMG/M. The phylogeny reported is the one in use in IMG/M [33], which uses the phylogeny described as part of the genomic encyclopedia of Bacteria and Archaea (GEBA) project [34]. Taxonomic composition is summarized in Table 8. Gene copies are estimated based on the number of genes in the assembled metagenome, multiplied by the average read depth of each gene. This provides a better estimate for the total number of reads coming from each taxon, which is proportional to the abundance of those taxa in the microbial mats. Across the assembled metagenomes, the fraction of annotated genes (not accounting for gene copies) that are unassigned at the 60% sequence identity level ranges between 64% and 67%, with 7–13% mapping to phylum Bacteroidetes , 8–13% phylum Cyanobacteria , and 9–16% phylum Proteobacteria . However the estimated gene copies show that these samples are in fact dominated by Cyanobacteria sequences (27–49% of estimated gene copies), with smaller contributions from Proteobacteria , Bacteroidetes , and a variety of other bacterial phyla, and only 34–44% unassigned. The majority of cyanobacterial sequences map to Coleofasciculus chthonoplastes (19–39% of the total estimated gene copies) and Lyngbya sp. PCC 8106 (3.5–5.5% of estimated gene copies). Other individual bacterial species that capture a large fraction of estimated gene copies at 60% identity include Erythrobacter sp. NAP1 ( Alphaproteobacteria ; up to 3.6% in MD6A), Allochromatium vinosum ( Gammaproteobacteria ; up to 3.3% in CD6A), and Marivirga tractuosa ( Cytophagia ; up to 2% in MD6A).

Table 8 Taxonomic composition

There are noticeable differences in taxonomic composition among the four metagenomes. For example, the molybdate treated samples MD2A and MD6A contain fewer sequences from phylum Cyanobacteria and more from phylum Bacteroidetes than the control samples. Some of these differences may be due to spatial heterogeneity in the mat from which the samples were collected.

Functional diversity

The distribution of COG functional categories is very similar between the four genomes (Table 9), with Pearson correlation of the log of the number of genes assigned to each category ranging from 0.986 (CD2A vs. CD6A) to 0.999 (CD2A vs. MD6A), suggesting a broad functional similarity between the samples, despite differences in species composition.

Table 9 Functional diversity

Conclusions

We sequenced and assembled metagenomes for four samples of microbial mat from the Elkhorn Slough estuary in Monterey Bay, California, to be used as reference data for a diel metatranscriptomic study in the presence or absence of molybdate. All four metagenomes were dominated by cyanobacterial sequences, primarily Coleofasciculus chthonoplastes . Despite some differences in community composition between the four metagenomes (which may be partly due to spatial heterogeneity in the mat), their functional composition in terms of COG functional categories is quite similar.