The green berries are found in the same organic-rich, intertidal pools of Little Sippewissett salt marsh on Cape Cod (MA, United States) where both multicellular magnetotactic bacteria and pink berries have been previously studied (Seitz et al., 1993; Shapiro et al., 2011; Wilbanks et al., 2014). Though less abundant than the pink berries found in these pools (Figure 1A), the green berries form similar irregular ellipsoid aggregates measuring 1–8 mm in diameter, with an average equivalent spherical diameter of 1.7 mm ± 0.1 mm (standard deviation, Figures 1A,D). Green berries were dense and compact aggregates that were typically observed at the sediment-water interface, but were occasionally found to float at the water surface when suspended by bubbles. Microscopic observation of the green berries revealed abundant coccoid unicellular cyanobacteria 5–7 μm in diameter (which we call GB-CYN1, “green berry cyanobacteria 1”), interspersed with pennate diatoms (Figures 1B,C). Filamentous cyanobacteria were observed occasionally, but were rare compared to the unicellular GB-CYN1 morphotype. A clear, extracellular matrix (putatively exopolysaccharides) coated these aggregates of phototrophic cells, and was colonized by a variety of smaller bacteria (Figures 1B,C). GB-CYN1 exhibited absorption maxima at 620, 660, and 680 nm corresponding to the presence of phycocyanin, allophycocyanin and chlorophyll a, respectively.

Sequencing of 18S rRNA genes from the green berries indicated that the eukaryotic community was predominantly made up of two different pennate diatom species related to Navicula cari strain AT-82.04c (96% sequence identity) and Amphora pediculus strain AT-117.11 (95% sequence identity; Supplementary Figure S1). These same diatom species were also the dominant eukaryotic 18S rRNA gene sequences recovered from pink berry aggregates, though diatoms were more abundant in green berries than in pink berries, as observed by microscopy and the relative abundance of 16S rRNA chloroplast sequences (Wilbanks et al., 2014).

Bacterial 16S rRNA gene sequences amplified from the green berries were dominated by sequences related to either diatom chloroplasts (phylum Bacillariophyta, 18/92 clones, representing 4 OTUs at 97% similarity threshold) or Chroococcales unicellular cyanobacteria (18/92 clones, representing 2 OTUs; Figure 2). Unassembled metagenomic sequence reads assigned to rRNA sequences and protein-coding regions support the observed abundance of Chroococcales (accounting for up to 37% relative sequence abundance), but did not recover comparable proportions of diatom chloroplasts (<5% relative sequence abundance; Figure 2).

The overall bacterial community structure of the green berries was significantly different from coexisting pink berry consortia (unweighted UniFrac analysis, Bonferroni corrected p-value < 0.01). Some abundant taxa from the pink berries co-occurred in the green berries as rare OTUs, such as the purple sulfur bacterial species Thiohalocapsa sp. PB-PSB1 (1/92 clones), and a Winogradskyella species (Flavobacteriales, 2/92 clones) (Wilbanks et al., 2014). The persistence of these distinct, co-occurring pink and green berry consortia suggests that the process of macroscopic aggregation enables niche partitioning between oxygenic (green berries) and anoxygenic (pink berries) phototrophs in these marsh pools.

Most of the non-cyanobacterial sequences in the green berry consortia are related to aerobic and facultatively anaerobic marine heterotrophs from the Bacteroidetes, Alphaproteobacteria, and Gammaproteobacteria (Figure 2). Many of these sequences (e.g., taxa from the Rhodobacterales, Kordiimonadales, Sphingomonadales, and Flavobacterales) were most closely related to environmental 16S rRNA sequences associated with aggregates of oxygenic phototrophs. Examples of such habitats included phytodetrital aggregates (marine snow) collected from euphotic and hadal environments (DeLong et al., 1993; Eloe et al., 2011), and epiphytes of marine macroalgae (Burke et al., 2011; Fernandes et al., 2012). The occurrence of related phylotypes in such environments suggests that taxa may be well adapted to an attached lifestyle, degradation of photosynthate, and the fluctuating oxygen conditions in an aggregate environment.

Metagenomic data indicate that the orders Rhizobiales and Rhodobacterales of the Alphaproteobacteria are abundant in the green berry consortia. While these groups were rarely detected in the PCR-based 16S rRNA survey, we have previously observed this same PCR bias from the 8F primer during studies of the pink berry consortia (Wilbanks et al., 2014). We find the abundance of these clades in the green berries particularly interesting as they include lineages of marine denitrifying bacteria. For example, pelagic Rhizobiales (e.g., Labrenzia and Roseibium species) have been linked to denitrification when found in association with macroscopic Trichodesmium sp. aggregates in oxic waters bordering oxygen minimum zones (Wyman et al., 2013).

The cyanobacterial 16S rRNA gene sequences from the green berries grouped into two closely related OTUs (97% similarity threshold), GB-CYN1a and GB-CYN1b, that can be confidently placed in the order Chroococcales (Figure 3). The GB-CYN1 monophyletic cluster formed a clade basal to the UCYN-A clade (96% sequence identity to CP001842, “Candidatus Atelocyanobacterium thalassa” isolate ALOHA). Using 29 concatenated single-copy phylogenetic marker genes (Wu et al., 2013) assembled from the metagenome (Supplementary Table S1), we reconstructed a phylogenetic tree that placed the GB-CYN1 within a clade including Crocosphaera watsonii and Cyanothece sp. ATCC 51142 as a sister taxa to “Candidatus Atelocyanobacterium thalassa” isolate ALOHA (Supplementary Figure S2).

A phylogenetic tree inferred from nifH gene sequences reveals that the near full-length nifH gene recovered from the GB-CYN1 metagenomic data affiliated with the UCYN-B clade, and was most closely related to Cyanothece sp. 8801/8802 and Crocosphaera watsonii (Figure 4). We conclude that the observed discordance between 16S rRNA, concatenated, and nifH gene phylogenies involving species such as Cyanothece sp. 8801, Gloeothece sp. KO68DGA, and the cyanobacterial endosymbiont of Rhopalodia gibba is most likely due to lateral gene transfer of the nifH gene. Lateral transfer of nifH has been observed in many other species, including mat-forming filamentous cyanobacteria (Cantera et al., 2004; Bolhuis et al., 2010; Parker, 2012; Gaby and Buckley, 2014; Meyer and Huber, 2014).

A full suite of nitrogenase genes were found in the green berry metagenome and were consistently assigned to GB-CYN1, indicating the metabolic potential for nitrogen fixation characteristic of other members of the UCYN A-C clades (Supplementary Table S1). Diazotrophy in the green berries was detected in whole aggregates by acetylene reduction. In two separate experiments (containing 5 berries each), we measured rates of 11 and 20 nanomoles acetylene reduced per hour per milligram of aggregate dry weight (nmol mg^-1 hr^-1). These rates are comparable, though faster than the rates of 3–6 nmol mg^-1 hr^-1 reported in macroscopic aggregates of filamentous cyanobacteria from Bogue Sound, North Carolina (Paerl and Prufert, 1987). Rates measured from actively growing Cyanothece cultures (Reddy et al., 1993) were two orders of magnitude larger (∼1000 nmol mg^-1 hr^-1) than the rates in the green berries.

The total aggregate elemental composition was analyzed using elemental analyzer isotope ratio mass spectrometry (EA-IRMS). The mean observed C:N ratio in the green berries, 7.1 ± 0.6 (n = 6), falls within range of the Redfield molar ratio (Redfield et al., 1963). This observed C:N ratio is higher than that the ratio of 5.4 ± 0.4 (n = 50) observed in similarly large, anoxic, diazotrophic Nodularia spumigena aggregates from the Baltic Sea (Klawonn et al., 2015). Cultures of Crocosphaera watsonii exhibit wide diel fluctuations in C:N ratios (∼5 at dawn to ∼9 at dusk) as a result of temporal partitioning of carbon (day) and nitrogen (night) fixation activities (Dron et al., 2012, 2013). Our samples, collected in the late afternoon on a 14 h light/10 h dark photoperiod, are comparable to reports of C:N = 7 from C. wastsonii at similar late afternoon times in a 16 h light/8 h dark photoperiod (Dron et al., 2013). Future studies investigating the temporal partitioning for such activities in the green berry aggregate would be informative to elucidate the dynamics of carbon and nitrogen flow in the consortia.

Oxygen microsensors were used to characterize the balance between respiration and photosynthesis in the green berry aggregates. Aggregates examined were relatively symmetric ellipsoids of similar size with an equivalent spherical diameter (ESD) of 1.7 ± 0.1 mm (n = 5, average ± standard deviation). Photosynthesis produced supersaturated oxygen concentrations within the aggregates: 380 μM O₂ with illumination at 170 μE m^-2 s^-1 (one lamp) and 520 μM O₂ at 320 μE m^-2 s^-1 (two lamps; Figure 5). During these experiments, bubbles were not observed on aggregate surfaces and the aggregates never floated. Oxygen production rates per aggregate (net photosynthesis) were calculated from these profiles as 13 and 31 nmol O₂ per hour at 170 μE m^-2 s^-1 and 320 μE m^-2 s^-1, respectively (Table 1). We did not rigorously determine the saturating light intensity; however, we observed that illumination with a third lamp failed to stimulate increased oxygen production beyond that with two lamps (320 μE m^-2 s^-1), suggesting that the saturating light intensity lies in the range of 170 – 320 μE m^-2 s^-1 (data not shown).

The green berries were anoxic in darkness with a dark respiration rate of 19 nmol per aggregate per hour (Figure 5B). This rate, equivalent to a volumetric rate of 6.8 μmol cm^-3 hr^-1, is well above the theoretical threshold for respiration rates capable of causing oxygen diffusional limitation from surrounding seawater (5.0 μmol cm^-3 hr^-1 for ESD = 1.7 mm; calculated after Ploug et al., 1997). In light-dark shifts, we observed a rapid response where the aggregate core transitioned from steady state supersaturation to full anoxia in 5–8 min (Figure 5C).

Volumetric gross photosynthetic rates were calculated by two methods: (1) from depth microprofiles via the sum of net photosynthesis and dark respiration and (2) via the light-dark shift technique (Revsbech et al., 1981) performed at a single point in the aggregate center (Table 1). At both light intensities examined, rates calculated via the light-dark shift method were found to be 4.5 μmol cm^-3 hr^-1, lower than those from depth microprofiles. While this difference could arise from biological variability between aggregates, we suspect that the light-dark shift rates measured at the aggregate core were lower than those we might have measured closer to the aggregate surface. Future depth integrated studies of photosynthetic rates will help to clarify this difference and allow better characterization of respiratory activity in the light.

Comparing the green berries’ dark respiration and gross photosynthesis to other photosynthetic mats and aggregates, we find them similar to the high rates measured for large (>1 mm), filamentous aggregates of the heterocystous cyanobacterium, Nodularia spumigena from the Baltic Sea (Ploug et al., 2011). Indeed, our estimates of carbon fixation (Table 1) are close to prediction of 349 ng C per aggregate per hour calculated using Ploug et al.’s (2011) regression of volume to gross photosynthesis from a 2009 Nodularia bloom. The green berry dark respiration rate (19 nmol agg^-1 h^-1) was similar, though slightly higher than that observed for Nodularia aggregates of similar diameter; however, the ratio of dark respiration to gross photosynthesis of 0.38 (for an aggregate volume of 2.8 mm³) was identical to that observed for 7 mm³ aggregates found late in the bloom (Ploug et al., 2011). In large Nodularia aggregates, dark anoxia is associated with active dissimilatory nitrogen cycling, including denitrification, dissimilatory reduction of nitrate to ammonia (DNRA), and significant rates of cryptic nitrification (Klawonn et al., 2015).

Our metagenomic data analysis predicts genes in GB-CYN1 involved in the Calvin cycle, TCA cycle, and photosystems I and II (Supplementary Table S1). Sequence identities for these genes in GB-CYN1 and published Chroococcales genomes were high, ranging from 75 to 99% (average of 85% for assembled contigs, Supplementary Table S1). The observation of photosystem II genes in GB-CYN1 indicates the metabolic potential for oxygenic photosynthesis, unlike UCYN-A which lacks photosystem II, RuBisCo and the TCA cycle. These findings confirm that GB-CYN1 resembles free-living Chroococcales of the clades B and C (Welsh et al., 2008; Bandyopadhyay et al., 2011; Bench et al., 2011, 2013) rather than the metabolically streamlined, endosymbiotic UCYN-A clade (Zehr et al., 2008; Tripp et al., 2010; Thompson et al., 2012; Hagino et al., 2013).

Related Chroococcales species are known to produce copious quantities of exopolysaccharides (EPS) (Pereira et al., 2009; Sohm et al., 2011; Mota et al., 2013), and EPS produced by Crocosphaera watsonii has been linked to the formation of transparent exopolymeric particles in pelagic environments (Passow, 2002; Webb et al., 2009; Sohm et al., 2011). A cassette of genes required for EPS production has been identified by comparative genomic analysis of related Crocosphaera watsonii strains (Webb et al., 2009; Bench et al., 2011, 2013); however, we did not recover homologs of these genes in our metagenomic sequence (either assembled contigs or sequence reads). We conclude that the absence of these sequences is most likely due to our incomplete sampling of the GB-CYN1 genome, though it could also indicate either an unknown pathway for EPS synthesis in GB-CYN1, or an alternate source of the green berry exopolymeric matrix (e.g., diatoms).

Consistent with the dynamic oxygen conditions within the green berries, we found sequences in the metagenomic data that suggest the metabolic potential for aerobic and anaerobic respiration, and anaerobic fermentative pathways assigned to several different phyla of bacteria (data not shown). Abundant fixed nitrogen from diazotrophy and transient anoxia presents an ecological opportunity for dissimilatory nitrogen metabolism, a process that could cause concomitant nitrogen fixation and loss over rapid spatiotemporal scales. Indeed, previous studies of diazotrophic Nodularia spumigena aggregates measured significant rates of both nitrification and denitrification (Klawonn et al., 2015).

We specifically investigated potential marker genes for dissimilatory nitrogen metabolism within the green berries. Metagenomic sequences homologous to the nitrite reductase gene nirK were found in six unassembled reads and on a corresponding 753 bp long contig (Supplementary Table S2). A single read homologous to nitrous oxide reductase nosZ was also found suggesting the presence of a denitrification pathway (Supplementary Table S2). The top database matches to these nirK and nosZ sequences belonged to marine phytoplankton epiphytes from the Alphaproteobacteria (nirK) and Flavobacterales (nosZ). In the case of the nirK metagenomic sequences, the best match was to Roseibium sp. TrichSKD4, an alphaproteobacterial species isolated from a nitrogen-fixing Trichodesmium aggregate in the Atlantic Ocean (Mann and Barbeau, 2014). Homologs to the napA periplasmic nitrate reductase and the nirB nitrite reductase were identified from six and five unassembled reads, respectively, and were most similar to database sequences from marine heterotrophs in the Alphaproteobacteria (Rhodobacterales, Rhizobiales) and Gammaproteobacteria (Oceanospirillales; Supplementary Table S2). We did not detect sequences supporting the presence of nitrifying bacteria or archaea, though given our limited sequence depth, this could be a function of missing data.

Green berries were sampled in June and July of 2010 (for all studies except O₂ microsensor work) or July and August 2014 (for oxygen microsensors) from a single intertidal pool formed in the Little Sippewissett Salt Marsh, Falmouth, MA United States (41°34′33.01″N, 70°38′21.24″W). These aggregates, however, have also been observed at other locations throughout Little and Great Sippewissett marshes. Berries were collected from the sediment-water interface by sieving (1 mm mesh size) and were washed three times in 0.2 μm filter sterilized marsh water. Imaging of the berries was conducted using a Zeiss Axio IMAGER MZ epifluorescence microscope equipped with a color camera (AxioCam HRc, Zeiss) and a Zeiss LSM 710 spectral confocal scanning laser microscope (Carl Zeiss, Oberkochen, Germany).

DNA extraction, 16S/18S rRNA PCR, and clone libraries were performed as described previously (Wilbanks et al., 2014). Roche GC Titanium 454 metagenomic sequencing was performed as described previously (Wilbanks et al., 2014), generating 100 megabases (Mbp) of sequence data after quality filtering (275,263 reads of average length = 365 bp). Metagenomic data was assembled using Newbler 2.3 (Margulies et al., 2005), generating a 4.8 Mbp assembly in 6,113 contigs of N50 equal to 941 bp and a maximum contig size of 7,845 bp. Data are publicly available at MG-RAST (reads: mgm4454183.3; assembly: mgm4454167.3). 16S/18S rRNA gene sequence from clone libraries available in NCBI GenBank (MF360994-MF361083; MF372423-MF372506). Metagenome data is associated with NCBI BioProject number PRJNA390846 and raw sequence read data is available in the NCBI SRA (SRR5710041). For the assembled metagenome, a Whole Genome Shotgun project has been deposited at DDBJ/ENA/GenBank under the accession NJIW00000000. The version described in this paper is version NJIW01000000.

The MG-RAST 3.3 pipeline was used to classify unassembled metagenomic sequence reads as ribosomal RNA and protein coding sequences using the M5RNA and M5NR databases to provide a diversity description shown in Figure 2 (Meyer et al., 2008). Metagenomic sequence data (both assembled sequences and unassembled reads) was mined using MG-RAST. Functions of interest were mined using the hierarchical function assignment and were assigned to GB-CYN1 when the best hits (amino acid similarity search) were to sequenced Chroococcales genomes. Marker genes were further investigated via profile HMM (Eddy, 2011) and BLASTX searches (Altschul et al., 1990).

PCR amplified SSU rRNA gene sequences were aligned to the SILVA 115 database (Quast et al., 2013) using SINA (Pruesse et al., 2012) and curated using ARB (Ludwig et al., 2004). 16S rRNA phylogenies were inferred using the GTRGAMMA rate approximation. For functional genes of interest that were found to have frameshift sequencing errors (nifH, nosZ), the sequences were corrected using Framebot (Wang et al., 2013). The metagenomic nifH amino acid sequence was aligned to references sequences from the nifH database available from the Zehr research group (Zehr and Jenkins, 2003) using the ARB software package (Ludwig et al., 2004). ProtTest3 was used to select the fixed WAG model of amino acid evolution with an inverse gamma rate approximation (Darriba et al., 2011) for the nifH phylogeny, though similar topologies were recovered using related models.

Phylogenetic reconstruction for both 16S rRNA and nifH genes were conducted as follows: maximum likelihood phylogeny was constructed with RAxML 7.2.8 (Stamatakis, 2006) with 1000 rapid bootstrap inferences and Bayesian phylogeny with MrBayes 3.3 (Ronquist et al., 2012). For Bayesian tree inferences, MCMC was run with default parameters and convergence was assessed when the reported average standard deviation of split frequencies fell below 0.01.

Hidden Markov model (HMM) profiles of 40 phylogenetic marker genes (phyeco markers) for Bacteria and Archaea (Wu et al., 2013) were used to search the green berry metagenomic peptide sequence database (from six frame translated assembled nucleotide sequence data) using the trusted cutoffs. For each of the 34 markers that have green berry hits, green berry peptide sequences were aligned with all the bacterial and archaeal references sequences by hmmalign in HMMER3 (Eddy, 2011). A maximum likelihood tree was built by Fastree2 (Price et al., 2010) for each alignments, and the alignments and trees were examined. As a result, 33 green berry peptide sequences were selected for further analysis because they branched with Cyanobacteria with good alignments.

These sequences were further filtered to include only those 29 phyeco gene families with a single copy, cyanobacterial green berry hit. Single-copied reference sequences from 126 Cyanobacteria genomes of the 29 markers and the green berry cyanobacterial metagenomic sequences were retrieved from the alignments built in the previous step and were concatenated into a large alignment. A maximum likelihood tree was built using PHYML3.0 with the LG substitution model (Guindon et al., 2010). Tree topology and branch lengths were optimized by the program and SH statistics was used for branch support estimation.

Whole-aggregate in vivo acetylene reduction assays were conducted on two 30 mL serum bottles of five washed green berries each. Acetylene, generated by the hydration of calcium carbide in an evacuated 150 mL serum-vial, was added to 10% of the headspace in a 30 mL serum-vial containing the berries in 5 mL of anoxic 0.2 μm filtered in situ marsh water under an 90% N₂-CO₂ atmosphere (90:10). Experiments were incubated on a 14 h light, 10 h dark cycle with full spectrum illumination at 30°C for 2 days. Acetylene and ethylene were quantified using a Varian 2400 series gas chromatograph (Varian Instrument Group, Walnut Creek, CA, United States) with an H₂ flame ionization detector, as described previously (Lobo and Zinder, 1988). Rates were calculated over the 2 days incubation period with the assumption that nitrogen fixation was restricted to the dark hours of incubation, as previously described for unicellular cyanobacteria (Dron et al., 2013).

Six green berries sampled from the acetylene reduction experiment (control conditions without acetylene added) in late afternoon were dried at 55°C and crimped in aluminum foil. Samples were analyzed using a Europa 20-20 elemental analyzer interfaced with a Europa ANCA-SL mass analyzer (Sercon Limited, Cheshire, United Kingdom) under contract with the Stable Isotope Lab at the Marine Biological Laboratory’s Ecosystem Center (Woods Hole, MA, United States).

Individual aggregates were placed on an agar plate and covered with filter-sterilized seawater collected from the marsh pool. An airstream was directed at the water so that slight ripples were visible on the surface. The water had a salinity of 3.5 psu and the temperature during measurements was either 24 or 27°C, which was taken into account for calibration. We used Clark-type oxygen microsensors (tip diameter 50 μm, 90% response time < 1 s; Unisense A/S, Aarhus, Denmark), and calibrated before and after measurements in air-saturated (by bubbling) and anaerobic (0.1 M acetate, 0.1 M NaOH) 3.5% NaCl. The micro-profiling apparatus and software Sensor Trace Suite was also provided from Unisense. Starting at the water surface, microprofiles through individual aggregates were measured in the light, and light-dark shifts. Time-lapse light-dark-shift recordings included sensor tips placed directly above the aggregate (in the water phase); and inserting the tip into the central core of the aggregate. Recordings of the oxygen signal were taken every second. Light sources were 65 W halogen lamps and experiments were conducted with either 170 μmol photons m^-2 s^-1 (one lamp), or 320 μmol photons m^-2 s^-1 (two lamps). Dark conditions were realized by switching off the lamps, removing them from the table, and carefully placing a carton box over the entire profiling setup to avoid residual light from the room. Background light intensities under the box were < 1 μmol m^-2 s^-1.

Theoretical limits of oxygen and DIC flux and whole aggregate O₂ flux calculations were calculated from depth concentration profiles according to Ploug et al. (1997), with a diffusion coefficient for O₂ in 3.5% saline water of 2.175 × 10^-5 cm² s^-1 at 24°C and 2.3535 × 10^-5 cm² s^-1 at 27°C. Inside the aggregate, the apparent diffusivity of O₂ was assumed to be 0.95 (Ploug and Passow, 2007). Carbon fixation was estimated based on a photosynthetic quotient of 1.2 (Masotti et al., 2007). The diffusion of oxygen in agar (0.2–2%) was not found to be different than in water over a wide range of salinities (Revsbech, 1989).