Metagenomic analysis of intertidal hypersaline microbial mats from Elkhorn Slough, California, grown with and without molybdate
Standards in Genomic Sciences volume 12, Article number: 67 (2017)
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.
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 H2 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 (H2 production, N2 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.
Cyanobacterial mats are compact, laminated, and highly structured microbial communities (Fig. 1) that comprise great diversity at both the metabolic and phylogenetic level  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 . 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.
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 . 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  and to show that the majority of fixation is attributable to a Lyngbya sp. , and to identify the dominant SRB ( Desulfobacterales ) in the ecosystem . Additionally, the site has been investigated for hydrogen cycling. Burow and colleagues , 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 , 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.
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.
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  (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 .
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.
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.
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
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.
Trimmed, screened, paired-end Illumina reads were assembled using SOAPdenovo v1.05  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  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) . Read depths were estimated based on read mapping with BWA . These sequences are currently available to the public at IMG/M and the JGI genome portals. Metagenome statistics are summarized in Table 5.
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 . 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 . 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 . Hits with high similarity were kept; all other parameters were set to default values. Ribosomal RNA genes were predicted using hmmsearch  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  and PILERCR . 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) ,Metagene (v. 1.0) , Prodigal (V2.50: November, 2010)  and FragGenescan (v. 1.16)  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 , the COG database using the NCBI developed PSSMs , the Pfam database  using hmmsearch. Assignment to KEGG Ortholog protein families was performed using the algorithm described in . Annotation parameters are summarized in Table 6.
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.
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 , which uses the phylogeny described as part of the genomic encyclopedia of Bacteria and Archaea (GEBA) project . 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).
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.
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.
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.
Basic local alignment search tool
Clusters of orthologous groups
Integrated Microbial Genomes
Sulfate reducing bacteria
Ley RE, Harris JK, Wilcox J, Spear JR, Miller SR, Bebout BM, et al. Unexpected diversity and complexity of the Guerrero negro hypersaline microbial mat. Appl Environ Microbiol. 2006;72:3685–95.
Burow LC, Woebken D, Marshall IP, Singer SW, Pett-Ridge J, Prufert-Bebout L, Spormann AM, Bebout BM, Weber PK, Hoehler TM. Identification of Desulfobacterales as primary hydrogenotrophs in a complex microbial mat community. Geobiology. 2014;12:221–30.
Hoehler TM, Bebout BM, Des Marais DJ. The role of microbial mats in the production of reduced gases on the early earth. Nature. 2001;412:324–7.
Skyring GW, Lynch RM, Smith GD. Quantitative relationships between carbon, hydrogen, and sulfur metabolism in cyanobacterial mats. In: Cohen Y, Rosenberg E, editors. Microbial Mats: physiological ecology of benthic microbial communities. Washington, DC: American society for Microbiology; 1989. p. 170–9.
Burow LC, Woebken D, Bebout BM, McMurdie PJ, Singer SW, Pett-Ridge J, et al. Hydrogen production in photosynthetic microbial mats in the Elkhorn slough estuary, Monterey Bay. ISME J. 2012;6:863–74.
Lee JZ, Burow LC, Woebken D, Everroad RC, Kubo MD, Spormann AM, et al. Fermentation couples Chloroflexi and sulfate reducing bacteria to Cyanobacteria in hypersaline microbial mats. Front Microbiol. 2014;5:61.
Bebout BM, Paerl HW, Bauer JE, Canfield DE, Des Marais DJ. Nitrogen cycling in microbial mat communities: the quantitative importance of N-fixation and other sources of N for primary productivity. In: Stal LJ, Caumette P, editors. Microbial Mats NATO ASI series. Berlin; Heidelberg: Springer; 1994. p. 265–71.
Omoregie EO, Crumbliss LL, Bebout BM, Zehr JP. Determination of nitrogen-fixing phylotypes in lyngbya sp. and microcoleus chthonoplastes cyanobacterial mats from Guerrero negro, Baja California, Mexico. Appl Environ Microbiol. 2004;70:2119–28.
Woebken D, Burow LC, Prufert-Bebout L, Bebout BM, Hoehler TM, Pett-Ridge J, et al. Identification of a novel cyanobacterial group as active diazotrophs in a coastal microbial mat using NanoSIMS analysis. ISME J. 2012;6:1427–39.
Woebken D, Burow LC, Behnam F, Mayali X, Schintlmeister A, Fleming ED, et al. Revisiting N2 fixation in Guerrero negro intertidal microbial mats with a functional single-cell approach. ISME J. 2015;9:485–96.
Steppe TF, Paerl HW. Potential N2 fixation by sulfate-reducing bacteria in a marine intertidal microbial mat. Aquat Microb Ecol. 2002;28:1–12.
Burow LC, Woebken D, Marshall IP, Lindquist EA, Bebout BM, Prufert-Bebout L, et al. Anoxic carbon flux in photosynthetic microbial mats as revealed by metatranscriptomics. ISME J. 2013;7:817–29.
Stuart RK, Mayali X, Lee JZ, Everroad RC, Hwang M, Bebout BM, Weber PK, Pett-Ridge J, Thelen MP. Cyanobacterial reuse of extracellular organic carbon in microbial mats. SME J. 2016;10:1240–51.
Bebout BM, Carpenter SP, Des Marais DJ, Discipulo M, Embaye T, Garcia-Pichel F, et al. Long-term manipulations of intact microbial mat communities in a greenhouse collaboratory: simulating earth's present and past field environments. Astrobiology. 2002;2:383–402.
Field D, Garrity G, Gray T, Morrison N, Selengut J, Sterk P, et al. The minimum information about a genome sequence (MIGS) specification. Nat Biotechnol. 2008;26:541–7.
SOAPdenovo v1.05. http://soap.genomics.org.cn/soapdenovo.html
Chaisson M, Pevzner P. Short read fragment assembly of bacterial genomes. Genome Res. 2007;18:324–30.
Aligner BW. (BWA). http://bio-bwa.sourceforge.net/
Edgar RC. Search and clustering orders of magnitude faster than BLAST. Bioinformatics. 2010;26:2460–1.
Lowe TM, Eddy SR. tRNAscan-SE: a program for improved detection of transfer RNA genes in genomic sequence. Nucleic Acids Res. 1997;25:955–64.
Altschul SF, Madden TL, Schäffer AA, Zhang J, Zhang Z, Miller W, Lipman DJ, Gapped BLAST. PSI-BLAST: a new generation of protein database search programs. Nucleic Acids Res. 1997;25:3389–402.
Eddy SR, Accelerated Profile HMM. Searches. PLoS Comput Biol. 2011;7:e1002195.
Bland C, Ramsey TL, Sabree F, Lowe M, Brown K, Kyrpides NC, Hugenholtz PCRISPR. Recognition tool (CRT): a tool for automatic detection of clustered regularly interspaced palindromic repeats. BMC Bioinformatics. 2007;8:209.
Edgar RC. PILER-CR: fast and accurate identification of CRISPR repeats. BMC Bioinformatics. 2007;8:18.
Besemer J, Borodovsky M. GeneMark: web software for gene finding in prokaryotes, eukaryotes and viruses. Nucleic Acids Res. 2005;33:W451–4.
Noguchi H, Park J, Takagi T. Meta gene: prokaryotic gene finding from environmental genome shotgun sequences. Nucleic Acids Res. 2006;34:5623–30.
Hyatt D, Chen GL, Locascio PF, Land ML, Larimer FW, Hauser LJ. Prodigal: prokaryotic gene recognition and translation initiation site identification. BMC Bioinformatics. 2010;11:119.
Rho M, Tang H, Ye Y. FragGeneScan: predicting genes in short and error-prone reads. Nucleic Acids Res. 2010;38:e191.
Tatusov RL, Fedorova ND, Jackson JD, Jacobs AR, Kiryutin B, Koonin EV, et al. The COG database: an updated version includes eukaryotes. BMC Bioinformatics. 2003;4:41.
Finn RD, Mistry J, Tate J, Coggill P, Heger A, Pollington JE, et al. The Pfam protein families database. Nucleic Acids Res. 2010;38:D211–22.
Mao X, Cai T, Olyarchuk JG, Wei L. Automated genome annotation and pathway identification using the KEGG Orthology (KO) as a controlled vocabulary. Bioinformatics. 2005;21:3787–93.
Markowitz VM, Chen IM, Chu K, Szeto E, Palaniappan K, Grechkin Y, et al. IMG/M: the integrated metagenome data management and comparative analysis system. Nucleic Acids Res. 2012;40:D123–9.
Wu D, Hugenholtz P, Mavromatis K, Pukall R, Dalin E, Ivanova NN, et al. A phylogeny-driven genomic encyclopaedia of bacteria and Archaea. Nature. 2009;462:1056–60.
We thank Jeff Cann, Associate Wildlife Biologist, Central Region, California Department of Fish and Wildlife, for coordinating our access to the Moss Landing Wildlife Area.
This research was supported by the U.S. Department of Energy Office of Science, Office of Biological and Environmental Research Genomic Science program under the LLNL Biofuels SFA, FWP SCW1039, and by JGI Community Sequencing Program award #701. Work at LLNL was performed under the auspices of the U.S. Department of Energy under Contract DE-AC52-07NA27344. Work at LBNL, the National Energy Research Scientific Computing Center (NERSC), and the DOE Joint Genome Institute (JGI) was performed under the auspices of the U.S. Department of Energy Office of Science under Contract No. DE-AC02-05CH11231.
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
About this article
Cite this article
D’haeseleer, P., Lee, J.Z., Prufert-Bebout, L. et al. Metagenomic analysis of intertidal hypersaline microbial mats from Elkhorn Slough, California, grown with and without molybdate. Stand in Genomic Sci 12, 67 (2017). https://doi.org/10.1186/s40793-017-0279-6