Impacts of Zostera eelgrasses on microbial community structure in San Diego coastal waters

Seagrasses provide several critical marine ecosystem services, including provisioning of food and habitat and acting as an indicator for the overall health of the ecosystems in which they are present (Lamb et al., 2017; Hogarth, 2015). Seagrasses also facilitate ecosystem health by trapping fine sediments and bacteria that are introduced by storm runoff or resuspended from the seafloor. By trapping sediments, eelgrasses can have a direct impact on pollutants, aiding in increased water clarity and purity (Gacia and Duarte, 2000; Lamb et al., 2017). Zostera marina and Zostera pacifica are species of seagrass, commonly referred to as eelgrasses, which naturally occur in the National Marine Fisheries Service (NMFS) West Coast Region. Included in this region are California, Washington, and Oregon (NOAA, 2018). Within California, these eelgrasses may be found along the entire coast from Imperial Beach in Southern California to Humboldt County in the North, but they are concentrated within the San Francisco, Humboldt, and San Diego Bays (Merkel & Associates, Inc., 2014).

Previously, the Z. marina microbiome has been studied by Ettinger et al. (2017) in Bodega Bay, California, where researchers examined the microbial communities found in association with different components of the eelgrass bed, including roots, leaves, and sediments. Using 16S rRNA gene sequencing, they were able to identify significant differences between the microbial communities of each sample type, i.e., between samples taken from the roots vs. leaves vs. sediments. Ettinger et al. (2017) also found that microbial communities differed between samples taken from varying distances to the edge of the seagrass patch, but only among the sediment samples. They did not observe this pattern between samples of the leaf and root microbiomes. The researchers examined how microbial communities differed from sediment found within the eelgrass patches compared to sediment from the edge of the eelgrass patches and outside of the patches; they found significant differences in overall microbial composition, with the most abundant orders being Bacteroidales (highest inside the patch), Flavobacteriales (highest on the patch edge), Desulfobacterales (highest inside the patch), Thiotrichales (highest outside of the patch), Clostridiales (highest inside the patch) and Alteromonadales (highest outside of the patch). Although the exact mechanism for this phenomenon is uncertain, they observed a correlation between microbial community structure, sampling location, C:N ratio, and eelgrass density.

The abilities of Z. marina to reduce pathogen load and alter microbial community structure have also been highlighted in previous work. Enterococcus assays showed Z. marina had a significant impact on bacteria in eelgrass meadow ecosystems, reducing viable Enterococcus, an indicator of fecal contamination (Ortega et al., 2009). Lamb et al. (2017) concluded that Z. marina eelgrass meadows found in the intertidal regions of islands in the Spermonde Archipelago, Indonesia, significantly reduced bacterial and pathogenic load in the surrounding waters to the benefit of both humans and marine organisms. The 16S rRNA gene sequencing and Enterococcus assays suggested a 50% reduction in the abundance of animal and human pathogens when eelgrass meadows were present, along with a two-fold reduction in coral disease in adjacent reefs (Lamb et al., 2017).

Although the specific mechanisms by which seagrasses influence microbial community structure are not well constrained, some seagrasses have been found to exhibit antimicrobial properties both in the ocean and in vitro (Alam et al., 1994; Choi et al., 2009; Mayavu et al., 2009; Kannan et al., 2010). For instance, methanol extracts from the eelgrass Enhalus acoroides exhibited antibacterial behavior on Gram-positive bacteria and the model pathogen Pseudomonas aeruginosa (Alam et al., 1994). Along with antibacterial potential, extracts from Z. marina have also demonstrated antioxidant properties (Choi et al., 2009). These experiments suggest that eelgrasses may directly influence microbial community structure in the coastal marine environment, including reducing the pathogen load.

In San Diego Bay, untreated sewage and storm runoff can influence microbial community structure and introduce pathogenic bacteria. San Diego coastal waters from Imperial Beach to as far north as Coronado are periodically exposed to sewage from the transboundary Tijuana River. During rain events, the Comisión Internacional de Límites y Agua (CILA) in Tijuana shuts down its pump station sewage diverter, resulting in the subsequent discharge of sewage into the Tijuana River estuary (U.S. Department of Justice, 2018). This water contains a variety of sewage and trash, including microbial contaminants, and can result in the closure of public beach access (U.S. Department of Justice, 2018; San Diego County, 2018). San Diego Bay is also periodically exposed to toxic stormwater runoff from Chollas Creek (Schiff et al., 2003; US Department of Justice, 2018). Furthermore, the bay has a history of locally sourced industrial pollution. In the decades preceding 1960, San Diego Bay was exposed to a constant influx of raw sewage and untreated discharge from a variety of industries, including fish canneries, kelp processing facilities, aircraft manufacturing plants, commercial shipyards, and multiple U.S. Naval installations (Fairey et al., 1998). More recently, heavy urban development in San Diego Bay and cross-border sewage outflow have raised concerns about water quality (San Diego County, 2018).

In this study we used 16S rRNA gene sequencing and flow cytometry to test whether local eelgrass meadows influence the structure of the marine microbial community in San Diego Bay, and whether this influence includes a reduced pathogen load. We examined the structure of bacterial and archaeal assemblages from water samples collected from paired sites where eelgrass beds are either present or absent. Sample pairs were collected from two areas, within and just outside of the San Diego Bay.

Sites likely to contain eelgrass were identified from the 2014 San Diego Bay Eelgrass Inventory (Merkel & Associates, Inc., 2014), the most recent survey available. The presence of eelgrass was confirmed visually; however, mapping the extant of eelgrass meadows was beyond the scope of this study. Water samples at eelgrass-present and eelgrass-absent sites were collected from a small boat via a peristaltic pump on 3 July 2017 (Figure 1, Table S1). Water was collected from either 10 cm (estimated) off the seabed, or in the mid-water column as indicated in Table S1. In the latter case the point of sample collection was less than 2 m above the seabed. For DNA samples, approximately 1 L of water was filtered through a sterile 0.2 μm Supor filter (Pall), and the filter was immediately placed on dry ice with longer-term storage at –80°C. Flow cytometry samples were collected into sterile 15-ml falcon tubes and immediately placed on ice. Chlorophyll a concentration, yield, and phytoplankton photosynthetic efficiency (F_v/F₀) were measured with an AquaFlash handheld active fluorometer (Turner Designs) following the manufacturer’s instructions. Temperature, salinity, and turbidity were measured with a YSI ProDss (Xylem).

DNA was extracted with the DNEasy PowerWater DNA extraction kit (Qiagen) within two weeks of sampling. Extracted DNA was quantified using the Qubit HS DNA quantification kit (Invitrogen) and then quality-checked by gel electrophoresis and PCR amplification of the 16S rRNA gene using primers 515F and 806R (Walters et al., 2015). High quality extracted DNA was submitted to the Argonne National Laboratory sequencing center for amplification and library preparation with the same primer set, then paired-end sequenced on the Illumina MiSeq platform. Sequence reads and accompanying metadata were uploaded to NCBI SRA under BioProject PRJNA483963.

Samples for flow cytometry were stored in the dark at 4°C, and fixed to a final concentration of 4% formaldehyde within 12 hours of sampling. Fixed samples were stained with SybrGreen I (Molecular Probes), spiked with a known concentration of 123-count beads (ThermoFisher), and quantified on a Sysmex CyFlow Space flow cytometer. High nucleic acid (HNA), low nucleic acid (LNA), and overall bacterial cell counts were determined from populations identified from the forward scatter (FSC) and 488 nm (FL1) parameters using a self-organizing map (SOM) following Bowman et al. (2017).

Illumina MiSeq reads were demultiplexed using the ‘iu-demultiplex’ command in IlluminaUtils (Eren et al., 2013). Demultiplexed reads were quality-controlled and denoised using the ‘FilterandTrim’ and ‘dada’ commands within the R package dada2 (Callahan et al., 2016), and assembled with the ‘mergePairs’ command. The final, merged reads all had mean quality scores ≥ 30. The R package dada2 is designed to denoise 16S rRNA amplicon data such that the final “unique” reads reflect true biological diversity. The non-redundant fasta files of unique reads produced by dada2 were made redundant using a custom Python script (https://github.com/bowmanlab/seq_data_scripts) and used as input for the paprica pipeline for microbial community structure and metabolic inference (Bowman and Ducklow, 2015). The paprica method for determining microbial community structure differs from most methods in that it relies on the placement of reads on a phylogenetic tree created from the 16S rRNA gene reads from all completed bacterial and archaeal genomes in Genbank. Because the metabolic potential of each phylogenetic edge on the reference tree is known, a reasonable estimate of gene content can be made for the organism of origin for each read. Edge abundance data were normalized to predict 16S rRNA gene copy number prior to statistical analysis.

Phylogenetic placement of some phylotypes of interest yielded relatively low posterior probabilities. The identity of these phylotypes was further investigated with the Ribosomal Database Project Bayesian classifier (Wang et al., 2007) and blastn megablast (Altschul et al., 1990) against the NCBI 16S rRNA gene database.

The count matrix generated from the 16S rRNA gene read data was used in combination with phyloxml-format phylogenetic trees generated by Guppy (Matsen et al., 2011), as implemented by paprica, along with the Archaeopteryx visualization software tool (Han and Zmasek, 2009), to identify the phylogeny of observed microbial taxa. All statistical analysis was carried out in R (R Core Team, 2012) and RStudio (RStudio Team, 2016). The R package vegan (Oksanen et al., 2013) was used to examine microbial community ecology by generating a dendrogram from the count matrix via dissimilarity analysis with hierarchical clustering based on Bray-Curtis dissimilarity index. The R package clustsig and the similarity profile analysis (SIMPROF) tool were used to determine the number of significant sample clusters based on gene content using the ‘hclust’ function, assuming no predetermined clustering (Yoshioka, 2008; Clarke et al., 2008). These results were utilized in conjunction with the gplots package and its ‘heatmap.2’ function to generate all abundance heatmaps (Warnes et al., 2016).

Principal coordinate analysis (PCoA) was performed using the ‘pcoa’ function from the R package ape v5.2 (Paradis and Schliep, 2018), and taxon- (as phylogenetic edges or clades) weighted means for each principal coordinate were calculated with the ‘wascores’ function in package vegan (Oksanen et al., 2018). We used the R package DESeq2 (Love et al., 2014) to identify phylogenetic edges or clades that were differentially present between eelgrass-present and eelgrass-absent sites, and between sites inside and outside of San Diego Bay. DESeq2 performs differential abundance analysis based on the negative binomial/Gamma-Poisson distribution. We applied the default analysis, which uses estimation of size factors with the “median ratio method” described in Anders and Huber (2010), followed by estimation of dispersion. Next, the Wald test for generalized linear model coefficients was used to test for significance of coefficients, taking into account the size factors and dispersion estimates that were previously calculated. The DESeq2 ‘results’ (‘res’) function, with a p-value set to 0.05, allowed for testing of the null hypothesis: sampling site has no effect on which microbial taxa are present based on sampling location or eelgrass presence. This process was repeated for the domain Archaea. After the most abundant, differentially present bacterial taxa were identified, we determined which of these taxa had potential pathogens in their lower-order clades. To determine if specific clades contained bacterial strains of pathogenic bacteria, we utilized the National Institutes of Health (NIH)- National Center for Biotechnology Information (NCBI) Pathogen Detection Isolates Browser, which includes bacterial pathogen genomic sequences sourced from food, hospital patients, and the environment, as well as a literature search. Clades were labeled as potentially pathogenic if the strain itself or members of the taxonomic clade have been previously identified as being pathogenic to animals or humans.

The HNA content, LNA content, and total bacterial abundance as determined by flow cytometry were used to test for a correlation between eelgrass presence/absence and marine bacterial abundance in our samples. A mixed-effects modeling approach was applied to evaluate differences in HNA, LNA, and total bacterial abundance between sites. The lme4 package (Bates et al., 2012) was used to construct linear mixed effects models to determine and understand the interaction between eelgrass presence, bacterial abundance, and HNA and LNA content in the collected samples. For the mixed effects models, the fixed effect was “eelgrass” (present vs. absent). The random effect was an intercept for sample site location, in conjunction with by-location random slopes for the effect of eelgrass presence. Random slope models were generated to observe how the effect of eelgrass presence differed between locations (i.e., inside bay or outside bay). Significance was determined by p-values generated with likelihood ratio tests (ANOVA) against the full mixed effects models.

We determined significant differences in the microbial community structure across sample location and seagrass presence or absence through use of 16s rRNA sequencing, flow cytometry, and statistical analysis. SIMPROF clustering identified nine significant sample clusters from our collection of samples (n = 45), each cluster having similar overall bacterial community composition among its members (Figures 2 and 3).

The taxa identified using the ‘results’ function in DESeq2 as having a p-value < 0.05 were determined to be differentially present between the inner and outer bay samples and between eelgrass-present and absent sites. From our list of 416 identified bacterial clades, 127 of these clades were identified as differentially present between the inner and outer bay sites (Figure 2). From these 127 clades, we focused on the top 30 taxa ranked by abundance (Table 1), which enabled differences in presence and abundance of specific bacterial clades between sampling conditions to be more apparent visually (Figure 3). In outer bay samples, bacterial taxa found in significantly different abundance from inner bay samples include Fluviicola taffensis, Pseudohongiella spirulinae KCTC 32221, Sulfitobacter, Marinobacter, Pelagibacter ubique HTCC 1062, family Fimbriimonadaceae, Formosa sp. Hel3 A148, Cyanobacteria, Thioglobus singularis PS1, and Planktomarina temperata RCA23 were (Figure 2, Table 1). In comparison, bacterial taxa found in significantly different abundance in samples from the inner bay include Martelella endophytica YC6887, Polaribacter vadi, genus Oxynema, Pelagibacter sp. IMCC9063, and order Rickettsiales (Figure 2 and Table 1).

For both the inner and outer bay, there were minor but significant differences in taxa identified between paired sites. Using DESeq2, we identified 13 differentially present bacterial clades between the eelgrass-present and eelgrass-absent sites (p < 0.05) (Table 2). Taxa found in greater abundance in samples collected at eelgrass-present sites include family Rhodobacteraceae, Candidatus Puniceispirillum marinum IMCC1322, Halioglobus pacificus RR3-57, Teredinibacter turnerae T7901, genus Colwellia, Thioglobus singularis PS1, genus Tenacibaculum, Fluviicola taffensis, Synechococcus sp. CC9311, and Phylum Cyanobacteria (Table 2). In comparison, bacterial taxa found in greater abundance at eelgrass-absent sites include Halobacteriovorax marinus SJ, Oscillatoria acuminata PCC 6304, and Candidatus Promineofilum breve Cfx-K (Table 2).

To further examine and explain variation in the microbial community structure of each sample (n = 45) due to the influence of eelgrass presence and sampling location, we conducted a principal coordinate analysis (PCoA). The first principal coordinate (PC1) accounted for 80.7% of the variance between samples, and the second principal coordinate (PC2) accounted for 10.1% of variance (Figure 4). The distribution of samples in PC1 is strongly dependent on sample location (inner vs. outer bay) (Student’s t-test, p approaches 0, n = 45). We therefore interpret PC1 as a good indicator of sample location. Sample distribution along PC1 is largely the result of select bacterial taxa, including Vibrio gazogenes ATCC 43942, Sulfuricurvum kujiense DSM 16994, and genus Bartonella, being elevated in the inner bay samples with positive PC1 values. Sample distribution along PC2 is largely the result of select bacterial taxa, including Altererythrobacter ishigakiensis NBRC 107699, Kangiella sediminilitoris KCTC 23892, and Archangium gephyra DSM 2261, being associated with samples with more positive PC2 values, most noticeably outlier Sample 16.

Of the 63 archaeal clades identified across the 45 samples, we identified five as differentially present between inner and outer bay sampling sites (Table 3). None of the differentially present archaeal clades have known pathogenic members (Table 3).

We then identified five archaea as differentially abundant between eelgrass-present vs. absent sites. As for the archaea that were differentially abundant between inner and outer bay sites, none of these are known pathogens (Table 4).

We then utilized flow cytometry data and determined whether there was a significant (p < 0.05) correlation between eelgrass presence and marine bacterial abundance in our samples. Specifically, we conducted linear regression models with random intercepts and random slopes to examine whether there was a significant correlation between eelgrass presence and bacterial abundance, with random intercepts for the effect of sample site location in conjunction with by-location random slopes for the effect of eelgrass presence. When comparing the random intercept-random slope model to the null model for HNA in a likelihood ratio test (ANOVA), we did not find evidence that eelgrass presence affected HNA abundance (χ² (1) = 2.42, p = 0.1198, n = 45) between seagrass-present and -absent sites. Similarly, when comparing the random intercept-random slope model to the null model for LNA with an ANOVA, we found no significant effect of eelgrass presence on LNA abundance (χ² (1) = 2.81, p = 0.09393, n = 45) between seagrass-present and -absent sites. In terms of overall bacterial abundance, when comparing the random intercept-random slope abundance model to the null model with an ANOVA, we could not conclude that eelgrass presence significantly affected overall bacterial abundance between eelgrass-present vs. -absent sites (χ² (1) = 2.91, p = 0.08799, n = 45).

In this study we utilized 16S rRNA gene sequencing to examine the microbial community structure of sampling areas located within and just outside of the San Diego Bay, at sites where eelgrass beds were either present or absent. For these samples we determined that while eelgrass presence does have some influence on the microbial community structure of surrounding waters, site location (inside vs. outside of San Diego Bay) exerts a stronger influence on overall microbial community structure. These findings offer an examination of how Zostera eelgrass beds impact their ecosystem within and just outside of San Diego Bay. Regardless of eelgrass presence or absence, however, orders that may include pathogenic bacteria were relatively more abundant in samples taken from inside than from outside of the bay. This trend may result from the more quiescent conditions inside the bay, or from a greater concentration of organic carbon or nutrients, for example, from runoff or sewage input to the bay (not addressed in this study).

Among the 30 most abundant and differentially present bacterial clades between the inner and outer Bay sites, order Rickettsiales (Paddock et al., 2003; Nicholson et al., 2010; Kang et al., 2014), family Flavobacteriaceae (Touchon et al., 2011; McBride, 2014; Småge et al., 2016), and order Pseudomonadales were identified as lineages containing possible pathogens (Figure 3 and Table 1). Family Flavobacteriaceae (map id = 0.87, post-prob = 0.72) was found somewhat evenly in samples collected either within or outside of the bay as well as in eelgrass-present or -absent sites (Figure 3). Pathogens of this taxonomic family include fish pathogens Flavobacterium branchiophilum and Flavobacterium columnare and human pathogens such as Elizabethkingia meningoseptica (McBride, 2014). A phylotype classified as order Pseudomonadales (map id = 0.86, post-prob = 1.00) was evenly distributed across all sites, suggesting that eelgrass does not influence the presence of this phylotype at local scales. Some pathogenic phylotypes of this order, as determined with the NCBI Pathogen Detection Isolates Browser, include 14 strains of Acinetobacter baumannii: PHEA-2, BJAB0868, AB31, 6411, B8300, XH856, AP_882, IEC338SC, YMC2010/8/T346, CA16, SSA3, SSA6, USA15, and HUMV-6483.

Order Rickettsiales was found predominantly in samples collected from within San Diego Bay, without regard to eelgrass presence or absence. While phylogenetic placement could not identify a specific strain, this order is of interest as it includes both human and animal pathogens such as Ehrlichia canis, Ehrlichia chaffeensis, Ehrlichia ewingii, Rickettsia typhi, and Rickettsia hellongjiangensis (Paddock et al., 2003; Nicholson et al., 2010; Kang et al., 2014). Although it is unlikely that the recovered 16S rRNA gene reads reflect the presence of these terrestrial, tick-based pathogens, these reads could indicate the presence of alternate parasites or symbionts with a similar phylogenetic signal. This interpretation is supported by the low posterior probabilities obtained for these reads during phylogenetic placement. Members of the ubiquitous SAR11 clade of marine bacteria are often commonly misidentified as Rickettsiales in 16S rRNA gene studies due to similar GC bias (Luo, 2014). While our analysis did associate reads with the sequenced genomes of the SAR11 genus Pelagibacter, it is possible that these Rickettsiales represent a coastal clade with weak 16S rRNA gene similarity to the sequenced Pelagibacter genomes. Further analysis of Rickettsiales reads by blastn megablast suggested close sequence homology with members of the genus Ehrlichia (E-value of 2 × 10^–28, 76% identity across the alignment region); however, this analysis does not preclude misclassification.

In contrast to Lamb et al. (2017), we did not find conclusive evidence that eelgrass beds reduced overall microbial contamination in San Diego coastal waters. This discrepancy could reflect geographical differences; the Indonesian locations where Lamb et al. (2017) conducted their research are exposed to significant levels of microbial contaminants and pollution due to a lack of sanitation systems and porous surface soils that are unable to retain wastewater. The point-source nature of sewage inflow on the Indonesian islands could also account for differences, with a “ring” of eelgrass beds around the islands reducing the transport of pathogens outside the ring. Sewage contamination of the San Diego coastal environment is associated with significant rain events, with the last heavy rainfall prior to sampling coming nearly 6 months earlier on 9 January 2017 (https://www.ncdc.noaa.gov/stormevents). The difference in pathogen load between eelgrass-present and eelgrass-absent sites may be much more pronounced immediately following rain events.

Regardless of any overall reduced pathogen load, we did identify 13 bacterial taxa that were differentially present between eelgrass-present and eelgrass-absent sites. While none of these taxa are known to include human pathogens, some are of interest for other reasons. Teredinibacter turnerae is known to exhibit antibiotic activity against both Gram-positive and Gram-negative bacteria including Bacillus cereus and Staphylococcus sciuri (Trindade-Silva et al., 2009). Halobacteriovorax marinus SJ and other members of genus Halobacteriovorax are predatory bacteria that attack Gram-negative bacteria including those of the Vibrio, E. coli, and Pseudomonas variety (Enos et al., 2018). Genus Tenacibaculum contains a pathogenic phylotype, Tenacibaculum maritimum, which is a known fish pathogen that infects fish kidneys and causes symptoms including skin lesions and fin rot (Avendaño-Herrera et al., 2006). Teredinibacter turnerae and genus Tenacibaculum were significantly more abundant in samples collected from eelgrass-present sites, while Halobacteriovorax marinus SJ was significantly more abundant at eelgrass-absent sites. Each of these phylogenetic placements had relatively high posterior probabilities (0.96, 0.60, and 0.92 respectively), indicating high confidence in the phylogenetic placement. We note that bacteria that are potentially beneficial to human health were present at both eelgrass-present and -absent sites; i.e., the presence of Teredinibacter turnerae at sites with eelgrass and Halobacteriovorax marinus SJ at sites without. Due to their antibiotic and predatory nature against human pathogens, further research into the presence of these bacteria would allow an examination of how they may be altering the local marine microbiome, whether in conjunction with Zostera eelgrass beds or not, and how their characteristics might be leveraged for biotechnological applications.

Eelgrass beds are critical components of the coastal marine ecosystem, with demonstrated influence on microbial community structure. In some environments, however, their influence on microbial community structure may be less significant than other factors, e.g., site location, physical processes, or other external factors not addressed in this experiment. Although we are cautious about extrapolating our findings to other coastal environments, this experiment presents useful information regarding specific microbial community structure of San Diego Bay and its immediately adjacent waters. Follow-on work to identify differences in microbial community structure when pathogen load is elevated, such as after heavy rainfall, and across an expanded selection of sites, will further clarify the influence on eelgrass beds on microbial community structure.

DNA sequences: NCBI SRA: PRJNA483963.

We would like to thank Brett Pickering for assistance with sampling and Jane Teranes for facilitating the SURF REU program.

TR was supported by the Scripps Institution of Oceanography SURF REU program (NSF-OCE 1659793). JSB was supported by the Simons Foundation Early Career Investigator in Marine Microbial Ecology and Evolution fellowship. Additional support came from the Doherty Career Development Fund through Scripps Institution of Oceanography.

The authors have no competing interests to declare.

• Contributed to conception and design: <TR, JSB>

• Contributed to acquisition of data: <TR, JSB, NE>

• Contributed to analysis and interpretation of data: <SJW, JSB>

• Drafted and/or revised the article: <SJW, JSB>

• Approved the submitted version for publication: <SJW, TR, NE, JSB>