Marine microbial communities of the Great Barrier Reef lagoon are influenced by riverine floodwaters and seasonal weather events

Sampling design

Sampling was performed in the Wet Tropics Region of the GBR (Fig. 1), a well-studied coastal area which is regularly exposed to river runoff and flood events (Devlin & Schaffelke, 2009; Schroeder et al., 2012; Turner et al., 2013). The sites surveyed were located on a transect following a gradient of river exposure, from the highly-exposed Tully River mouth (TT1) to the fringing coral reefs of Dunk Island (TT3) that are seasonally reached by flood plumes, and the TT4 off-shore location, rarely exposed to river water. Russel (RI) and Fitzroy islands (FI) were additional reef sites with limited exposure to the waters of the Johnstone and Russel rivers, respectively, and a consistently higher coral health index than Dunk Island (Thompson et al., 2014). All sampling sites were classified based on their proximity to the nearest influencing river mouth: ‘plume’ for <20 km downstream and ‘marine’ if >20km (Table S1). To characterize the influence of riverine microorganisms on marine communities, a freshwater site located 12 km upstream of the river mouth was also selected (TR).

The seven sites were surveyed between January 2011 and October 2013 in the dry season (June), just prior to the wet season (October), and at the end of the wet season (March). At each site and sampling date, a single 2 L seawater sample was taken from a depth of 5 m (just below surface for the TR site, which was a very shallow river bed), passed through a 0.22 µm Sterivex filter, which was stored at −20 °C until further processing. All samples were collected under the auspices of the general permit (G12/35236.1) granted by the Great Barrier Reef Marine Park Authority to the Australian Institute of Marine Science.

Environmental conditions

Water samples collected from 2011 to 2013 were processed according to the long-term GBR Reef Rescue Marine Monitoring Program (Thompson et al., 2013) to assess temperature, salinity, bottom depth and water chemistry: concentrations of suspended solids (SS), particulate organic carbon (POC), particulate phosphorus (PP), particulate nitrogen (PN), dissolved inorganic nitrogen (DIN), silica (Si) and chlorophyll a (CHLA). In addition, diuron concentration was determined by collecting 1 L of water in pre-washed bottles, and storing the water at 4 °C until processing by solid phase extraction liquid chromatography tandem mass spectrometry (SPE-LC-MS/MS) at Queensland Health and Forensic Scientific Services, Coopers Plains, Australia. All these measured environmental parameters were deposited as NCBI BioSamples (accession # PRJNA276058). In addition, meteorological parameters were acquired from public resources: Tully River discharge and water temperature (Department of Natural Resources and Mines, http://watermonitoring.derm.qld.gov.au/host.htm), and solar exposure and rainfall (Bureau of Meteorology, http://www.bom.gov.au/climate/data/index.shtml). An average of these meteorological parameters was calculated for the seven days preceding each sampling date.

A river exposure index was calculated using a hydrodynamic model (http://www.bom.gov.au/environment/activities/coastal-info.shtml), based primarily on the Sparse Hydrodynamic Ocean Code (SHOC) hydrodynamic model (http://www.emg.cmar.csiro.au/www/en/emg/software/EMS/hydrodynamics.html). SHOC is a general purpose model, applicable on spatial scales ranging from estuaries to regional ocean domains (Herzfeld, 2006). We used outputs from the regional application of SHOC to the GBR using a horizontal spatial resolution of about 4 km, with a model grid size of 180 × 600 with 48 vertical layers with 1 m resolution at the surface. In this context, conservative tracers were used in this study to simulate the transport of unique tracers ‘released’ from different rivers. In mathematical terms, a conservative substance is represented without terms of sinks or sources of mass in the transport equation. This means that the change in concentration values happens due to physical processes (advection and diffusion). Furthermore, it was stipulated that the tracers used do not affect the hydrodynamics. This technique enables the identification of marine regions influenced by individual catchments, and provides insight into the mixing and retention of river water along various regions in a given domain (Brinkman et al., 2002; Brinkman et al., 2014; Luick et al., 2007). Model simulations of the 3-dimensional distributions of passive tracers were analyzed to produce weekly estimates of cumulative exposure to tracers above a threshold of 1% of the source concentration. An exposure index was calculated that integrates the tracer concentrations above this threshold, based on a cumulative measurement of the exposure concentration and duration of exposure related to individual river sources, and expressed as Concentration x Days (conc.d). For every location in the model domain, the cumulative exposure index was calculated as: ${\displaystyle \mathrm{Conc.Days}=\sum _{t=0}^T{\mathrm{Conc}}_{\mathit{exceed}}\times t\mathrm{where}{\mathrm{Conc}}_{\mathrm{exceed}}=\left\{\begin{array}{c}\mathrm{Conc}\left(t\right)-{\mathrm{Conc}}_{\mathit{thresh}},\mathrm{if}\mathrm{Conc}\left(t\right)>{\mathrm{Conc}}_{\mathrm{thresh}}\hfill \\ 0,\mathrm{if}\mathrm{Conc}\left(t\right)\le {\mathrm{Conc}}_{\mathrm{thresh}}\hfill \end{array}\right.}$ and Conc_thresh is defined as 1% of the source concentration, Conc(t) represents the time-varying tracer concentration, and t is the time in days from the beginning of the wet season to the end (01 November–31 March). Cumulative exposure was calculated for each grid point in the model domain. Using this representation, the exposure index integrates both concentration above a defined threshold and the duration of exposure. For example, an exposure of 20 days at a concentration of 1% above the threshold would produce an index value of 0.2, which is equivalent to 10 days exposure at 2% above the concentration threshold. This index provides a consistent approach to assess relative differences in exposure of inshore GBR waters to inputs from various rivers. For each of the wet seasons simulated by the model, spatial maps of river exposure indices were calculated for the target rivers: Herbert, Tully, Murray, Johnstone, Mulgrave and Russel rivers (Wet Tropics catchment), Burdekin and Haughton rivers (Burdekin catchment, affecting the south of the Wet Tropics catchment).

16S rRNA gene amplicon sequencing

DNA was extracted from each Sterivex filter using a modified method from Suzuki et al. (2004). In brief, the filters were thawed on ice with Invitrogen’s P1 buffer with lysozyme at a final concentration of 2 mg/mL, and incubated for 30 min at 37 °C, while rotating at 10 rpm. Proteinase K (0.75 mg/mL final concentration) and 10% sodium dodecyl sulfate (1% final concentration) were added and the sample was incubated, with rotation, at 55 °C for 2 h. DNA was extracted using phenol:choloroform:isoamyl alcohol (25:24:1; pH 8.0) followed by an overnight ethanol precipitation and purified using a MO BIO PowerClean DNA Clean-Up kit (Carlsbad, CA, USA).

Amplicons were generated by PCR-amplifying the V6–V8 variable regions of the 16S rRNA gene using the pyroLSSU926F and pyroLSSU926F universal primers as described in Dove et al. (2013). The resulting DNA amplicons were sequenced on a Roche-454 GS-FLX instrument at the Australian Centre for Ecogenomics and deposited in the NCBI Short Read Archive (accession # PRJNA276058).

Bioinformatic processing

Amplicon reads were processed using Hitman (https://github.com/fangly/hitman), a bioinformatic workflow based around the UPARSE methodology (Edgar, 2013). In brief, Hitman: (1) joins read pairs with PEAR (Zhang et al., 2014), but keeps the forward read when pairs cannot be joined; (2) truncates the 3′ end of sequences at the first residue below a threshold quality value (Q) using TRIMMOMATIC (Bolger, Lohse & Usadel, 2014); (3) trims the 3′ end of all sequences to a target length (L) using TRIMMOMATIC, discarding all smaller sequences, (4) removes sequences exceeding the maximum number of expected errors (E) using USEARCH’s fastq˙filter (Edgar & Flyvbjerg, 2015); (5) uses USEARCH’s cluster˙otus to form operational taxonomic units (OTUs) from high-fidelity sequences (stringent quality processing in steps 2 and 4) that are sorted by decreasing abundance, occur at least twice in the dataset and meet a minimum percentage of similarity (O); (6) discards chimeric OTUs using USEARCH’s cluster˙otus in a reference-independent, and using UCHIME (Edgar et al., 2011) based on a reference database (C); (7) assigns regular-fidelity sequences (less stringent quality processing in steps 2 and (4) to each OTU using USEARCH’s usearch˙global (Edgar, 2010); (8) formats the results in BIOM format using Bio-Community’s bc˙convert˙files (Angly, Fields & Tyson, 2014); (9) gives a taxonomic assignment to each OTU by globally aligning their representative sequences against a database (T) of reference sequences trimmed to the target region (keeping only the best-matching alignment with a minimum required identity percentage (I) using USEARCH’s usearch˙ global; (10) removes OTUs belonging to specific taxa (W) using Bio-Community’s bc˙manage˙ samples; (11) rarefies the microbial profiles at the given depth (D) with Bio-Community’s bc˙accumulate assuming an infinite number of bootstrap replicates; and (12) corrects gene-copy number bias using CopyRighter (Angly et al., 2014).

In this study, Hitman was run using the following parameters: L = 250 bp, Q = 7 (16 for HiFi sequences), E = 3.0 expected errors (0.5 for HiFi sequences), O = 97% identity (species-level), C = GOLD database (Bernal, Ear & Kyrpides, 2001), T = merged Silva (Quast et al., 2012) and Greengenes (McDonald et al., 2012) databases (https://github.com/fangly/merge˙gg˙silva), I = 95% identity (genus-level), W = “Eukaryota* *Chloroplast*” and D =279 for Bacteria & Archaea (100 for Eukaryotes). In addition, rarefaction curves were generated using Bio-Community (Angly, Fields & Tyson, 2014).

Statistical analysis

All statistical analyses were carried out using the R language (R Core Development Team, 2015)). Comparisons of diversity between groups of samples were carried out using the non-parametric, unilateral Mann–Whitney U test (wilcox.test() function). Principal coordinates analysis (PCoA) and PERMANOVA were performed using the capscale() and adonis() functions of the vegan package (Dixon, 2003). The indicspecies package (Cáceres & Legendre, 2009) was used to determine indicator species with the multipatt() function. Redundancy analysis (RDA) model selection was based on the AIC (Akaike information criterion) and calculated by ordistep() in vegan. Pearson correlations between environment variables were computed using rcorr() from the Hmisc package. The functions fa.parallel(), fa(), target.rot() and fa.diagram() from the psych package were used to conduct exploratory factor analysis (EFA), i.e., to identify groups of co-varying variables. EFA was performed on several subsets of the data including different environmental parameters and the results were summarized.

Sampling and environmental context

Seven inshore GBR sites exposed to different levels of river runoff from the Wet Tropics catchments were surveyed over three years for water chemistry assessments and determination of microbial community structure. The classification of these sites as plume or marine sites was based on their distance to the nearest influencing river mouth (Table S1) and matched well with their river exposure index as calculated by oceanographic modeling (Fig. 1); sites <20 km from a river mouth were more highly exposed to riverine water (>20 conc.d) than sites >20 km away (<15 conc.d).

The weather in the Tully catchment from 2010 to 2013 followed the expected seasonal dynamics, dry and cool conditions between the months of May and October, hot and humid with most of the annual rainfall from November to April (Fig. 2). However, the 2010–11 wet season was marked by extreme weather (Table S2) and the landfall of category 5 tropical cyclone Yasi (3 February 2011), that significantly affected coral reefs (Perry et al., 2014). Annual river discharges reached a record high (Schaffelke et al., 2011), causing elevated exposure to nutrients, PSII herbicides and sediments across most inshore GBR regions (Devlin et al., 2012b; Kennedy et al., 2012; Perry et al., 2014). The following wet season (2011–12) represented a return to typical weather conditions, with river discharge close to the long-term median (Schaffelke et al., 2012a; Wallace et al., 2014). However, coral reef recovery was delayed until 2014, when the coral health index reached pre-2011 levels again (Thompson et al., 2014).

Microbial diversity

The archaeal and bacterial microbial profiles obtained by Roche-454 sequencing of 16S rRNA gene amplicons were rarefied (279 counts per sample) to allow comparison of microbial profiles (Fig. S1). A ten times higher sequencing effort would have been needed to sample nearly all the OTUs present in these aquatic samples (richness; Fig. S1A). Nevertheless, the chosen rarefaction depth recovered the vast majority of archaeal and bacterial diversity (Shannon–Wiener index; Fig. S1B). The microbial diversity of the rarefied profiles was calculated, with a median richness of ∼65 OTUs in the river, and ∼90–100 at the marine and plume sites (Fig. 3A). At plume sites, richness (Fig. 3A) and evenness (Fig. 3B) were higher in the dry season than in the wet season, which corresponded to a larger overall diversity (Shannon–Wiener index) in the dry season (Fig. 3C) (Mann–Whitney U test; p < 0.05).

Establishing a microbial monitoring baseline for inshore GBR

Taxonomic assignments and gene copy number correction were conducted to produce accurate estimates of microbial relative abundance (Angly et al., 2014) (Fig. S2) that yield an understanding of the prevalence of microbial taxa in the GBR lagoon waters. At a coarse taxonomic level, the bacterial orders Sphingobacteriales, Burkholderiales, and Xanthomonadales dominated the TR river site (Fig. 4), while the archaeal order E2 and bacterial orders Rickettsiales and Synechococcales were prevalent at the plume and marine sites (Fig. 4). On average, Rickettsiales and E2 had higher relative abundance at the marine sites (29.6 and 17.7% respectively) compared to the plume sites (20.4 and 10.4% respectively). This distribution was not constant over time, and notably, the plume sites were characterized by a large fraction of Burkholderiales in the January–March 2011 period, i.e. the extreme wet season 2010–11. Flavobacteriales were found in all sites, riverine, plume and marine. Burkholderiales are Betaproteobacteria commonly found in rivers (Cottrell et al., 2005; Ghai et al., 2011; Liu et al., 2012). Sphingobacteriales have been reported at riverine locations affected by waste water treatment effluent (Drury, Rosi-Marshall & Kelly, 2013), where they may degrade complex compounds such as herbicides and antibiotics (Kämpfer, 2011). Synechococcales such as those found in the plume and marine sites include the Synechococcus and Prochlorococcus genera, which represent the main photosynthetic bacteria in oceanic waters (Partensky, Blanchot & Vaulot, 1999; Partensky, Hess & Vaulot, 1999), and the small heterotrophic Rickettsiales are also commonly reported in the ocean (Morris et al., 2002; Carlson et al., 2008). More specifically, both orders were previously identified in coral reef waters (Kelly et al., 2012; Lu et al., 2015). However, Thermoplasmata E2 Archaea have only occasionally been reported in temperate surface seawater (Massana et al., 1997; Pernthaler et al., 2002). These Archaea are motile photo-heterotrophic cells focused on the degradation of protein and lipids (Iverson et al., 2012) and are found in coral mucus, suggesting a potential role in maintaining coral health (Kellogg, 2004).

Genomic methods such as 16S rRNA gene amplicon sequencing are effective for marine monitoring (Bourlat et al., 2013). The present taxonomic characterization of the microbial communities in inshore GBR waters represents a baseline against which future microbiological studies can be compared. This baseline may prove valuable for assessing future change in this reef ecosystem, be it further degradation or recovery.

Geographical and temporal distribution of OTUs

Non-constrained ordination (PCoA) was applied to get a precise account of the dynamics and distribution of specific microbial taxa in inshore GBR. At the finer operational taxonomic unit (OTU) level (Fig. 5A), significant spatial (distance to river mouth and site type; PERMANOVA, p < 0.05) and temporal effects (wet or dry; PERMANOVA, p < 0.05) were confirmed. Site type (riverine, plume or marine) had the strongest effect (27.8% explained variance along PC1, Fig. 5A), followed by seasons (wet or dry) (10.0% explained variance along PC2, Fig. 5A).

Marine and plume samples were characterized by indicator OTUs, i.e. OTUs characteristic of the marine and plume areas (indicator species analysis; p < 0.05). These indicator OTUs belonged to the Rickettsiales (e.g., Pelagibacteraceae), Synechococcales (e.g., Prochlorococcus), Acidimicrobiales (e.g., OCS155), and the Archaeal order E2 (e.g., Marine Group II) (Fig. 5A), which is largely consistent with the trends observed at the order level (Fig. 4). The main difference between the OTU- and order-level analyses (Figs. 5 and 4 respectively) was the detection of OCS155, an OTU indicator of marine and plume sites (Fig. 5A). OCS155 belongs to the Acidimicrobiales, a group of likely planktonic, free-living microorganisms with a photo-heterotrophic lifestyle (Mizuno, Rodriguez-Valera & Ghai, 2015) found in all tropical and temperate photic areas (Ghai et al., 2013).

Indicator OTUs for the river samples were also identified (indicator species analysis; p < 0.05), including members of the orders Burkholderiales (family Comamonadaceae, e.g. Limnohabitans, Acidovorax and Rhodoferax; family Oxalobacteraceae), and Sphingobacteriales (family Flexibacteraceae, e.g., Arcicella; family Chitinophagaceae) (Fig. 5A). These results are congruent with those at the order level (Fig. 4). The presence of riverine OTUs belonging to the order Burkholderiales and Sphingobacteriales at some plume sites surveyed in the wet months (Figs. 4 and 5A), when river discharge is high, is a sign that these sites are affected by riverine water effluent.

When repeating the PCoA with the marine samples only, i.e. restricting the analysis to samples unaffected by river output, a clear partition between dry and wet seasons (Fig. 5B) was seen (PERMANOVA, p < 0.05). All samples clustered in one of two season-specific groups arranged along PC1, except for one sample collected at TT4 in December 2011. Several indicator OTUs were identified (indicator species analysis; p < 0.05) (Table S3), e.g., Pelagibacteraceae (3 OTUs), Pirellulaceae (2 OTUs) and Marine group II (2 OTUs), which were specific of the dry season. Generally though, Marine group II OTUs seemed to be more abundant in the wet season (Figs. 5A and 5B). Previous studies have shown seasonality in near-shore microbial communities (Treusch et al., 2009; Gilbert et al., 2009; Gilbert et al., 2012). In particular, Pelagibacter is known to exhibit seasonality (Alonso-Sáez et al., 2007; Carlson et al., 2008; Eiler et al., 2009; Fortunato et al., 2013). Overall, the microbial communities of the GBR lagoon seem to respond to seasonal influence, although it is not as pronounced as the influence of geographical location, an observation that was also made in a previous investigation of river to ocean gradient (Fortunato et al., 2012).

In addition to a seasonal effects, GBR lagoon samples were susceptible to short term effects of potentially high magnitude. Rapid community changes were evident when looking at the Tully transect samples (TT1–4) collected in the wake of Cyclone Yasi (03 Feb 2011) (Fig. S3). The influence of Tully river input was marked, with elevated levels of Sphingobacteriales and Burkholderiales recorded between 13 and 17 February. While these changes were especially pronounced at the Tully River mouth, even the offshore site TT4 experienced analogous changes in this period. Especially large changes in community structure were observed between consecutive days, on February 12–13 (Fig. 2). These findings add to previous investigations of soil and gut microbiota, which have shown that community composition can change within a few days (Schmidt et al., 2007; Michelland et al., 2011; David et al., 2014). While baselines may be established using e.g., monthly sampling, it is clear that future research should adopt an intensive daily sampling regimen to better characterize such transient changes in estuarine and marine microbial communities.

Environmental drivers of coastal microbial dynamics

A range of measurements describing the environmental context of the microbial samples were collected in this study, including temperature, salinity, bottom depth and water chemistry variables (concentrations of suspended solids, silica, chlorophyll a, the herbicide diuron, dissolved and particulate form of carbon, nitrogen and phosphorus), and complemented with weather condition data (Tully River discharge, local rainfall and solar exposure). Many of these environmental parameters were highly correlated (Pearson test; 48.6% of the pairs with r > 0.3), requiring the need for EFA to extract groups of uncorrelated parameters. Following this procedure (Fig. S4), four independent factors were identified. Factor MR1 included local rainfall and Tully River discharge. Factor MR2 included chlorophyll a, suspended solids and other particulates (POC, PN, PP). Factor MR3 included DIN and salinity and can be interpreted as mixing with freshwater. Factor MR4 combined water temperatures and solar exposure. Although larger discharge (factor MR1) can lead to increased suspended solids (factor MR2) and salinity decrease (factor MR3), factor MR1, MR2 and MR3 were not directly correlated because salinity and the amount of suspended solids depend not only on site location (proximity to the river mouth) but also on the action of wind and waves, which homogenize the water column and resuspend sediments (Larcombe et al., 1995; Orpin & Ridd, 2012; Fabricius et al., 2014). A water quality index based on comparison to water quality guidelines (GBRMPA, 2010) and ranging from −1 for poor quality to +1 for very good quality was previously introduced (Thompson et al., 2013). This index aggregates scores given to four indicator concentrations (suspended solids, chlorophyll a, particulate nitrogen, particulate phosphorus), which are all part of the suspended material factor (factor MR2) identified by EFA in the present study. Diuron has a strong association with sediments (Stork, Bennett & Bell, 2008; Balakrishnan, Takeda & Sakugawa, 2012; Xu et al., 2013) and is also included in this factor. Overall, this suggests that factor MR2 can be interpreted as the quality of the water.

Constrained ordinations (RDA) were carried out to study the relative importance of the four independent factors (MR1–4) on microbial community structure. A single environmental parameter was chosen to represent each factor prior to conducting RDA with model selection: rainfall for factor MR1, water quality index for MR2, salinity for MR3, and water temperature for MR4 (Fig. 6). These four environmental parameters were significantly associated with microbial community composition (PERMANOVA, p < 0.05). Microbial community structure was affected by temperature (Fig. 6A), consistent with previous findings in the Western English Channel (Gilbert et al., 2009), which could explain the seasonality identified in marine sites. Further, seasonality can be attributed to the higher rainfall typical of the wet season (Fig. 6A). However, large changes could be attributed to decreases in water quality index and salinity, especially for plume sites (Fig. 6A), with negative consequences for corals (Fabricius et al., 2014). These changes can be interpreted as the effect of riverine floodwater, which was especially pronounced for sites located near the Tully river mouth, and helps explain the geographical pattern seen in microbial communities. The effects of salinity and water quality index were mostly collinear in Fig. 6A, despite the fact that they represent two separate factors (MR2 and MR3). Repeating the RDA analysis with the influence of salinity removed (conditioning term, Fig. 6B), water quality index was still identified as a significant driver of microbial communities (PERMANOVA, p < 0.05). This suggests that, as a whole, nutrient, organic compounds and herbicides such as diuron brought by riverine water into the ocean have an effect on microbial communities.

RDA analysis was also conducted to identify the factors driving the large, transient changes detected by microbial profiling at the plume sites on 12–13 February 2011 (Fig. S5). Being only nine days after the landfall of tropical cyclone Yasi, this period was characterized by extreme weather, e.g., high river discharge (12 Feb: 37.1 ML/d; 13 Feb: 36.2 ML/d; 7–13 Feb mean: 27.5 ML/d). While river discharge was similar on the two days, the February 13 microbial communities were characterized by elevated levels of suspended solids and lower salinity (PERMANOVA, p < 0.05), suggesting that other environmental factors may contribute to microbial community structure. We speculate that larger wind or waves were present on February 13, which would have decreased the vertical stratification of the water column, resulting in lower salinity and higher concentration of suspended solids. This example illustrates that salinity and suspended solids at inshore locations may be a function of both river flow and wind (Schaffelke et al., 2012b; Fabricius et al., 2014). The 12 February was also marked by an increase in the eukaryotic to prokaryotic reads ratio (EPR) (Fig. S6), with a maximum EPR of 1.45 reached on February 13 at TT1. The most abundant eukaryal taxon in this sample was assigned to the hydrozoan genus Merona, that includes very small organisms that can spawn eggs (Schuchert, 2004). These data could be interpreted as the potential spawning of hydrozoan or the displacement from their usual habitat (e.g., the benthic zone) and they illustrate that extreme weather may dramatically change environmental conditions, thereby affecting all microbial kingdoms (Archaea, Bacteria and Eukaryota).