The effect of carbon subsidies on marine planktonic niche partitioning and recruitment during biofilm assembly

Chuck Pepe-Ranney, Edward Hall



Biofilms are diverse and complex microbial consortia, and, the biofilm lifestyle is the rule rather than the exception for microbes in many environments. Large and small-scale biofilm architectural features play an important role in their ecology and influence their role in biogeochemical cycles (Battin et al., 2007). Fluid mechanics impact biofilm structure and assembly (Hodl et al., 2011; Besemer et al., 2009; Battin et al., 2003), but it is less clear how other abiotic factors such as resource availability affect biofilm assembly. Aquatic biofilms initiate with seed propagules from the planktonic community (Hodl et al., 2011; McDougald et al., 2011). Thus, resource amendments that influence planktonic communities may also influence the recruitment of microbial populations during biofilm community assembly.

In a crude sense, biofilm and planktonic microbial communities divide into two key groups: oxygenic phototrophs including eukaryotes and cyanobacteria (hereafter “photoautotrophs”), and heterotrophic bacteria and archaea. This dichotomy, admittedly an abstraction (e.g. non-phototrophs can also be autotrophs), can be a powerful paradigm for understanding community shifts across ecosystems of varying trophic state (Cotner et al., 2002). Heterotrophs meet some to all of their organic carbon (C) requirements from photoautotroph produced C while simultaneously competing with photoautotrophs for limiting nutrients such as phosphorous (P) (Bratbak et al., 1985). The presence of external C inputs, such as terrigenous C leaching from the watershed (Jansson et al., 2008; Karlsson et al., 2012) or C exudates derived from macrophytes (Stets et al., 2008a; Stets et al., 2008b), can alleviate heterotroph reliance on photoautotroph derived C and shift the heterotroph-photoautotroph relationship from commensal and competitive to strictly competitive . Therefore, increased C supply should increase the resource space available to heterotrophs and increase competition for mineral nutrients decreasing nutrients available for photoautotrophs (assuming that heterotrophs are superior competitors for limiting nutrients as has been observed ). These dynamics should result in the increase in heterotroph biomass relative to the photoautotroph biomass along a gradient of increasing labile C inputs. We refer to this differential allocation of limiting resources among components of the microbial community as niche partitioning.

While these gross level dynamics have been discussed conceptually (Cotner et al., 2002) and to some extent demonstrated empirically (Stets et al., 2008a), the effects of biomass dynamics on photoautotroph and heterotroph membership and structure has not been directly evaluated in plankton or biofilms. In addition, how changes in planktonic communities propagate to biofilms during community assembly is not well understood. We designed this study to test if C subsidies shift the biomass balance between autotrophs and heterotrophs within the biofilm or its seed pool (i.e. the plankton), and, to measure how changes in biomass pool size alter composition of the plankton and biofilm communities. Specifically, we amended marine mesocosms with varying levels of labile C input and evaluated differences in photoautotroph and heterotrophic bacterial biomass in plankton and biofilm samples along the C gradient. In each treatment we characterized plankton and biofilm community composition by PCR amplifying and DNA sequencing 16S rRNA genes and plastid 23S rRNA genes.

Materials and Methods

Experimental Design

Test tube racks were placed in one smaller (185 L, control) and 3 larger (370 L) flow-through mesocosms. All mesocosms were fed directly with marine water from an inflow source in Great Bay, Woods Hole, MA, approximately 200 m from the shore. Each mesocosm had an adjustable flow rate that resulted in a residence time of approximately 12 h. Irregular variation in inflow rate caused the actual mesocosm flow rate to vary around the targeted flow rate throughout the day. However, regular monitoring ensured that the entire volume of each system was flushed approximately two times per day (i.e. maintained a residence time of \(\sim\)12 h). To provide a surface for biofilm formation, we attached coverslips to glass slides using nail polish and then attached each slide to the test tube racks using office-style binder clips. Twice daily 10 mL of 37 mM KPO\(_{4}\) and 1, 5 and 50 mL of 3.7 M glucose were added to each of 3 mesocosms to achieve target C:P resource amendments of 10, 100 and 500. The goal of the resource amendments was to create a gradient of labile carbon. The same amount of P was added to each treated mesocosm to ensure that response to additions of C were not inhibited by extreme P limitation. The control mesocosm did not receive any C or P amendments.

DOC and Chlorophyll Measurements

To assess the efficacy of the C additions, we sampled each mesocosm twice daily during the first week of the experiment to evaluate dissolved organic C (DOC) content. After the initiation of the experiment we collected plankton on filters regularly to evaluate planktonic Chl a and bacterial abundance. Once it was clear that pool size of each community had been altered (day 8) we filtered plankton onto 0.2 \(\mu\)m filters and harvested coverslips to assess bacterial and algal biofilm community composition (16S and plastid 23S rDNA). In addition, all mesocosms were analyzed for community composition a second time (day 17) to assess how community composition of both the plankton and biofilm communities had changed over time. Control samples were only analyzed for community composition on day 17.

Samples for dissolved organic C (DOC) analysis were collected in acid washed 50 mL falcon tubes after filtration through a 0.2 \(\mu\)m polycarbonate membrane filter (Millipore GTTP GTTP02500, Sigma Aldrich P9199) attached to a 60 mL syringe. Syringes and filters were first flushed multiple times with the control sample to prevent leaching of C from the syringe or the filter into the sample. Samples were then frozen and analyzed for organic C content with a Shimadzu 500 TOC analyzer (Wetzel et al., 2000). Biomass of all biofilm samples were measured by difference in pre-(without biofilm) and post-(with biofilm) weighed GF/F filters after oven drying overnight at 60 \(\,^{\circ}\)C.

For Chl a analysis, we collected plankton on GF/F filters (Whatman, Sigma Aldrich Cat. # Z242489) by filtering between 500 mL and 1L from the water column of each mesocosm for each treatment. For biofilm samples, all biofilm was gently removed from the complete area of each coverslip (3 coverslips for each treatment per sampling event) before being placed in a test tube for extraction with 90-95% acetone for \(\sim\)32 hours at -20 \(\,^{\circ}\)C and analyzed immediately using a Turner 10-AU fluorometer (Wetzel et al., 2000).

Bacterial abundance of the planktonic samples was analyzed using DAPI staining and direct visualization on a Zeis Axio epifluorescence microscope after the methods of Porter et al. (1980). Briefly, 1-3 mL of water was filtered from three separate water column samples through a 0.2 \(\mu\)m black polycarbonate membrane filter and post stained with a combination of DAPI and Citifluor mountant media (Ted Pella Redding, Ca) to a final concentration of 1 \(\mu\)L mL\(^{-1}\).

DNA extraction

For plankton, cells were collected by filtering between 20 30 mL of water onto a 0.2 \(\mu\)m pore-size polycarbonate filter (Whatman Nucleopore 28417598, Sigma-Aldrich cat# WHA110656). For biofilm communities, biomass from the entire coverslip area of three separate slides was collected and combined in an Eppendorf tube by gently scraping the slip surface with an ethanol rinsed and flamed razor blade. DNA from both the filter and the biofilm was extracted using a Mobio Power Soil DNA isolation kit (MoBio Cat. # 12888).


Samples were amplified for 454 sequencing using a forward and reverse fusion primer. The forward primer was constructed with (5-3) the Roche A linker, an 8-10 bp barcode, and the forward gene specific primer sequence. The reverse fusion primer was constructed with (5-3) a biotin molecule, the Roche B linker and the reverse gene specific primer sequence. The gene specific primer pair for bacterial SSU rRNA genes was 27F/519R (Lane, 1991). The primer pair p23SrV_f1/p23SrV_r1 was used to target 23S rRNA genes on plastid genomes (Sherwood et al., 2007). Amplifications were performed in 25 \(\mu\)L reactions with Qiagen HotStar Taq master mix (Qiagen Inc, Valencia, California), 1 \(\mu\)L of each 5 uM primer, and 1 \(\mu\)L of template. Reactions were performed on ABI Veriti thermocyclers (Applied Biosytems, Carlsbad, California) under the following thermal profile: 95 \(^{\circ}\)C for 5 min, then 35 cycles of 94\(^{\circ}\)C for 30 sec, 54\(^{\circ}\)C for 40 sec, 72\(^{\circ}\)C for 1 min, followed by one cycle of 72\(^{\circ}\)C for 10 min and 4 \(^{\circ}\)C hold. Amplification products were visualized with eGels (Life Technologies, Grand Island, New York). Products were then pooled equimolar and each pool was cleaned with Diffinity RapidTip (Diffinity Genomics, West Henrietta, New York), and size selected using Agencourt AMPure XP (BeckmanCoulter, Indianapolis, Indiana) following Roche 454 protocols (454 Life Sciences, Branford, Connecticut). Size selected pools were then quantified and 150 ng of DNA was hybridized to Dynabeads M-270 (Life Technologies) to create single stranded DNA following Roche 454 protocols (454 Life Sciences). Single stranded DNA was diluted and used in emPCR reactions, which were performed and subsequently enriched. Sequencing followed established manufacture protocols (454 Life Sciences).

Sequence quality control

The 16S rRNA gene and plastid 23S rRNA gene sequence collections were demultiplexed and sequences with sample barcodes not matching expected barcodes were discarded. We used the maximum expected error metric (Edgar, 2013) calculated from sequence quality scores to cull poor quality sequences from the dataset. Specifically, we discarded any sequence with a maximum expected error count greater than 1 after truncating the sequence to 175 nt. The forward primer and barcode was trimmed from the remaining reads. We checked that all primer trimmed, error screened, and truncated sequences were derived from the targeted region of the LSU or SSU rRNA gene (23S and 16S sequences, respectively) by aligning the reads to Silva LSU or SSU rRNA gene alignment (“Ref” collection, release 115) with the Mothur (Schloss et al., 2009) NAST-algorithm (DeSantis et al., 2006) aligner and inspecting the alignment coordinates. Reads falling outside the expected alignment coordinates were culled from the dataset. Remaining reads were trimmed to consistent alignment coordinates such that all reads began and ended at the same position in the SSU rRNA gene and screened for chimeras with UChime in “denovo” mode (Edgar et al., 2011) via the Mothur UChime wrapper. 19,978 of 56,322 16S rRNA gene sequencing reads and 44,719 or 78,695 plastid 23S rRNA gene sequencing reads passed quality control.

Taxonomic annotations

Sequences were taxonomically classified using the UClust (Edgar, 2010) based classifier in the QIIME package (Caporaso et al., 2010) with the Greengenes database and taxonomic nomenclature (version ”gg_13_5” provided by QIIME developers, 97% OTU representative sequences and corresponding taxonomic annotations, (McDonald et al., 2012)) for 16S reads or the Silva LSU database (“Ref” set, version 115, EMBL taxonomic annotations, (Quast et al., 2013)) for the 23S reads as reference. We used the default parameters for the algorithm (i.e. minimum consensus of 51% at any rank, minimum sequence identity for hits at 90% and the maximum accepted hits value was set to 3).

Sequence clustering

Reads were clustered into OTUs following the UParse pipeline. Specifically USearch (version 7.0.1001) was used to establish cluster centroids at a 97% sequence identity level from the quality controlled data and to map quality controlled reads to the centroids. The initial centroid establishment algorithm incorporates a quality control step wherein potentially chimeric reads are not allowed to become cluster seeds. Additionally, we discarded singleton reads. Eighty-eight and 98% of quality controlled reads could be mapped back to our cluster seeds at a 97% identity cutoff for the 16S and 23S sequences, respectively.

Alpha and Beta diversity analyses

Alpha diversity calculations were made using PyCogent Python bioinformatics modules (Knight et al., 2007). Rarefaction curves show average OTU counts from 25 re-samplings at intervals of 10 sequences for each sample. Beta diversity analyses were made using Phyloseq (McMurdie et al., 2014) and its dependencies (Oksanen et al., 2013). A sparsity threshold of 25% was used for ordination of both plastid 23S and bacterial 16S libraries. Additionally, we discarded any OTUs from the plastid 23S rRNA gene data that could not be annotated as belonging in the Eukaryota or cyanobacteria for differential abundance, ordination and Adonis analyses. Cyanobacterial DNA sequences were removed from 16S rRNA gene sequence collections for ordination, Adonis and differential abundance analyses. We operated under the assumption that non-cyanobacterial bacteria are predominantly heterotrophs in our mesocosm setup and refer to non-cyanobacterial bacteria as “heterotrophs” in the manuscript (this abstraction is useful however there are likely exceptions – i.e. autotrophs among the non-cyanobacterial bacteria). All DNA sequence based results were visualized using GGPlot2 (Wickham, 2009). Adonis tests and principal coordinate ordinations were performed using the Bray-Curtis similarity measure for pairwise library comparisons. Adonis tests employed the default value for number of permutations (999) (“adonis” function in Vegan R package, Oksanen et al. (2013)). Principal coordinates of OTUs were found by averaging site principal coordinate values for each OTU with OTU relative abundance values (within sites) as weights. The principal coordinate OTU weighted averages were then expanded to match the site-wise variances (Oksanen et al., 2013).

Identifying enriched OTUs

We used an RNA-Seq differential expression statistical framework to find OTUs enriched in the given sample classes (R package DESeq2 developed by Love et al. (2014)) (for review of RNA-Seq differential expression statistics applied to microbiome OTU count data see McMurdie et al. (2014)). We use the term differential abundance coined by McMurdie et al. (2014) to denote OTUs that have different relative abundance across sample classes. We were particularly interested in two sample classes: 1) lifestyle (biofilm or planktonic) and, 2) high C (C:P = 500) versus not high C (C:P = 10, C:P = 100 and C:P = control). A differentially abundant OTU is enriched on one side of a sample class (e.g. enriched in the biofilm versus the plankton). Differential abundance could mark an enrichment of the OTU towards either side of the sample class and the direction of the enrichment is apparent in the sign (positive or negative) of the enrichment. Differential abundance was moderated (see Love et al. (2014)) such that the fold change OTU of an OTU across two categories of a sample class can be used to rank the enrichment of OTUs that span a wide range of base abundance. The DESeq2 RNA-Seq statistical framework has been shown to improve power and specificity when identifying differentially abun