Effects of experimental nutrient enrichment and eutrophication on microbial community structure and function in “marine lakes”

Diverse microbial communities play pivotal roles in the biogeochemical cycling of carbon (C) and nitrogen (N) and are an integral part of all ecosystems. However, human activities have widely modified the global C and N cycles (Rockström et al., 2009), affecting microbially mediated biogeochemical processes in multiple ways (Falkowski et al., 2008; Gruber and Galloway, 2008). Alteration of the N cycle is driven primarily by application of N-rich agricultural fertilizer, as only approximately 50% of applied fertilizer is eventually harvested in crops (Gruber and Galloway, 2008). The remainder is lost to air, soil, and water, where it has multiple unintended consequences (Galloway et al., 2003). In particular, N loading can increase primary production rates in aquatic ecosystems, leading to the elevated rates of organic matter supply and C flux known as eutrophication (Nixon, 1995). Microbes and their diverse enzyme systems subsequently breakdown and recycle organic material. Because even small differences in the balance between primary production and respiration can have substantial implications for the global C cycle (Duarte et al., 2013; Ducklow and Doney, 2013; Williams et al., 2013), determining the degree to which organic carbon is respired back to carbon dioxide is central to our understanding of the C cycle and climate change (Jiao et al., 2024). In aquatic ecosystems, eutrophication can also drive oxygen loss and hypoxia when microbial respiration consumes limited dissolved oxygen (DO) in water out of contact with the atmosphere (Breitburg et al., 2018). Microbial responses are therefore embedded in and integral to multiple aspects of nutrient enrichment and eutrophication—from the response of phytoplankton to nutrient loading in aquatic ecosystems, to the breakdown of organic C in all ecosystems, to the linked consequences for multiple biogeochemical cycles.

Despite their potential importance, our understanding of microbial community and functional responses to increases in N and C availability and flux is limited. Decades of research on nutrient limitation of phytoplankton growth and photosynthetic rates suggest differences across taxa that are still poorly resolved (Moore et al., 2013). Yet the response of microbial heterotrophs that respire organic C is more uncertain: there is still debate about whether photosynthesis and respiration are balanced in the open ocean (Duarte et al., 2013; Ducklow and Doney, 2013; Williams et al., 2013), and evidence that coastal ecosystems have switched from C source to C sink following anthropogenic perturbation (Bauer et al., 2013; Regnier et al., 2013). Notably, N enrichment can increase respiration or bacterial production directly by an order of magnitude (McAndrew et al., 2007)—an overlooked aspect of the N cycle that forms a strong link to the C cycle. In addition, experiments that supply N alongside C frequently elicit a stronger ecological or functional response than addition of C alone, even when labile forms of organic C are added (Carlson et al., 2002; Mills et al., 2008; Martínez-García et al., 2010). This strong effect of combined N and C indicates that microbial heterotrophs are highly responsive to combined nutrient enrichment and eutrophication, and that this combination may result in a community-mediated “priming” effect that elevates overall function (Carlson et al., 2002; Mills et al., 2008; Martínez-García et al., 2010). Yet how different groups of microbes interact, compete for N and C, and drive community responses is not well-known, as few studies have assessed how concurrent changes interact in additive or synergistic ways (Jaiswal et al., 2021). Consistent responses of particular clades to enrichment and eutrophication, regardless of starting conditions, could reflect a keystone role in ocean biogeochemistry (Pedler et al., 2014). Such groups may be disproportionately important if they are dominant under enriched and eutrophied conditions when rates of biogeochemical cycling are rapid.

We examined the strength and consistency of short-term responses in microbial community structure and function through experimental nutrient enrichment (N additions), eutrophication (C additions), and combined enrichment and eutrophication (N+C additions) of distinct starting communities. Increased N frequently leads to eutrophication, and so these are often considered to be synonymous (Nixon, 1995). However, distinguishing the effects of N enrichment versus eutrophication—more accurately defined as an increase in labile organic matter availability (Nixon, 1995)—is important, as N and C may drive distinct responses. Distinguishing these effects is particularly important when considering autotrophic and heterotrophic responses together, as N may drive dual increases in both primary production and respiration at the ecosystem level, while C is more likely to affect respiration alone. Quantifying the effects of combined of N and C is also essential, as responses may be nonlinear (Allison and Martiny, 2008; Wilson et al., 2017). We conducted experimental additions of N, C, and combined N+C in “marine lakes” (bodies of seawater surrounded by land; Figure 1A), for which the functional responses were previously published (Wilson et al., 2017). These ecosystems provide an ideal comparative framework, as Palau’s marine lakes range from well-mixed, oligotrophic, holomictic lakes to eutrophic, productive, and stratified meromictic lakes (Hamner et al., 1982; Hamner and Hamner, 1998; Colin, 2009; Wilson et al., 2017; Wilson et al., 2019). As such, this study system provides levels of biogeochemical variation similar to what one might see from nutrient-rich, productive coastal waters to the open ocean, and an opportunity to test how different starting conditions affect the responses of microbial community structure and function to nutrient enrichment and eutrophication (Table 1).

Previous studies in other systems have used gradients in space and time to relate natural microbial community variations to changes in community respiration (CR), net community production (NCP), or C export (Guidi et al., 2016; Wang et al., 2018); others have relied on experimental manipulations to link changes in function to community dynamics (McCarren et al., 2010; Robidart et al., 2019). However, combining experimental manipulations with natural variation to link community structure and function is rare, so we conducted factorial manipulations of N, C, and N+C in multiple lakes. The selected marine lakes contain well-known marine bacterial clades (e.g., Synechococcus, SAR11, SAR86, SAR324; Giovannoni and Vergin, 2012) but vary systematically in terms of their microbial ecology (Meyerhof et al., 2016; Rapacciuolo et al., 2019; Figure 1B). We sequenced and analyzed 16S rRNA genes to examine changes in microbial community structure, and measured rates of CR, NCP, and gross primary production (GPP; by difference) using in situ incubations (Wilson et al., 2017) to examine functional responses. We expected to find specific subsets of the microbial community correlated with functional responses to N and/or C additions, with the initial community dictating the potential taxonomic and functional change. We found tightly coupled changes in microbial community structure and function in combined N+C additions that exceeded the individual effects of N and C additions, that were consistent despite initial community differences, and that were driven by multiple networks of marine microbes.

Palau’s marine lakes formed 6,000–12,000 years ago (Rapacciuolo et al., 2019), when rising sea levels at the end of the last glaciation flooded inland basins with seawater, and marine organisms colonized the lakes (Hamner et al., 1982). Variations across marine lakes in shape, size, depth, and especially connectivity to the surrounding sea create a spectrum of starting conditions and communities (Rapacciuolo et al., 2019), ranging from well-mixed, low-productivity, holomictic lakes inhabited by many marine fauna to stratified, productive, meromictic lakes with subsurface anoxic water, such that their vertical structure resembles hypoxic zones or oxygen minimum zones (Colin, 2009; Meyerhof et al., 2016; Rapacciuolo et al., 2019). In this study, C and N additions were conducted across a gradient of seven of these lakes (Figure 1A), including: two large holomictic lakes, Heliofungia (HLO) and Mekeald (MLN), that resemble the surrounding ocean; two smaller holomictic lakes, Ngeruktabel (NLN) and Uet era Ngchas (ULN), which are more productive than lakes MLN and HLO; two permanently stratified meromictic lakes, Ngermeuangel (NLK) and Ongeim’l Tketau (OTM, aka Jellyfish Lake), that are anoxic at depth; and a transitional lake, T Lake (TLN), that stratifies and mixes every few years, but was stratified at the time of sampling. Stratification in the meromictic lakes and transitional lake is driven primarily by differences in salinity—and therefore density—with depth (Hamner et al., 1982; Hamner and Hamner, 1998). Surface water salinities are typically 25 ppt in OTM and NLK and increase to 32–33 ppt at depth (Hamner et al., 1982; Hamner and Hamner, 1998; Meyerhof et al., 2016), while salinities in all holomictic lakes are 32–35 ppt throughout the water column (Colin, 2009).

N and C additions were conducted for six lakes, while only N additions were conducted in lake HLO (for specific dates see Wilson et al., 2017). Water samples were collected at 1 m depth, and all but HLO experiments included unamended controls, N additions, C additions, and combined N+C additions—all with triplicate light and triplicate dark bottles for each control and treatment (Figure 1A). Incubations took place in 300 mL glass biological oxygen demand (BOD) bottles. Control bottles consisted of water collected in bottles and incubated in situ without any amendment. N additions were in the form of $NH4+$ (5 μM NH₄Cl), C additions were in the form of sucrose (100 μM C₁₂H₂₂O₁₁), and N+C additions included both at the above concentrations. Although these concentrations are high for many open ocean systems, they are more typical in coastal systems, including marine lakes (Hamner et al., 1982; Meyerhof et al., 2016; Wilson et al., 2017). For example, combined dissolved inorganic nitrogen (DIN) concentrations were all in the micromolar range at the onset of our experiments (Table 1), and previous dissolved organic carbon measurements in lake OTM were in the 100 μM range (Orem et al., 1991). $NH4+$ is also a more widely usable form of DIN, while sucrose is a labile C substrate that has been used in previous experiments in the open ocean (Kuparinen et al., 2011; Babbin et al., 2014), particularly in joint N and C manipulation experiments conducted by Babbin et al. (2014) in the eastern tropical North Pacific Ocean. Sucrose is central to terrestrial plant metabolism, and marine cyanobacteria and algae also produce sucrose (Salerno and Curatti, 2003; Abushattal et al., 2020). Moreover, recent analytical advances demonstrate that sucrose is particularly prevalent in tropical marine ecosystems (Sogin et al., 2019). For all lakes, bottles were incubated at 1 m depth on a floating array attached to a surface buoy (tied to the sides of the lake), such that all incubations took place under in situ conditions.

Triplicate light and dark bottle incubations were used previously to measure NCP and CR (Wilson et al., 2017). Briefly, beginning and end-point DO measurements were made in 300 mL glass BOD bottles filled to at least two times overflowing via slow laminar flow in order to exclude bubbles. All incubations were conducted for 24 hours so that changes in CR with time of day would be integrated and eliminated (Williams, 2000). Temperature and DO measurements were collected in the field using a YSI ProODO (YSI Incorporated, Yellow Springs, OH, USA). NCP and CR were calculated from changes in DO in light (NCP) and dark (CR) bottles over 24 hours.

For DNA extraction and microbial community analysis, 250 mL of water from each lake and experimental treatment were filtered through 0.22 µm filters (Millipore, Darmstadt, Germany) using a peristaltic pump. Samples were frozen in Sucrose-Tris-EDTA (STE) buffer until DNA was extracted following Beman et al. (2012). DNA concentrations were measured using the Quant-iT PicoGreen dsDNA Assay Kit and the manufacturer’s protocol (Invitrogen Corporation, Carlsbad, CA, USA) on a Stratagene MX 3005P (Agilent Technologies, Inc., Santa Clara, CA, USA). Samples were diluted to a common concentration and sent to Argonne National Laboratory (Lemont, IL, USA) for 16S rRNA sequencing using Illumina MiSeq with the universal primers 515F-Y (5′-GTGYCAGCMGCCGCGGTAA) and 926R (5′-CCGYCAATTYMTTTRAGTTT) according to Earth Microbiome protocols; these primers have been shown to be effective for analysis of marine microbial communities (Parada et al., 2016). Although considerable microbial functional diversity can be missed by 16S rRNA gene sequencing, we used this approach as a way to examine overall microbial community responses to experimental treatments.

Sequence processing and community analysis were conducted in mothur (http://www.mothur.org/) following the approach of Wilson et al. (2018) modified from the mothur MiSeq SOP (Kozich et al., 2013; 2019): as contigs were formed from forward and reverse reads, a delta 1 was allowed between quality scores of mismatched bases; during the screening process, a minimum length of 280 was required (75% of average at the recommendation of the Earth Microbiome Project), a maximum length of 385 was required, and up to 3 ambiguous bases were allowed. After quality screening, 3,648,696 16S rRNA sequences remained, with individual sample library sizes ranging from 4,750 to 149,919 sequences (mean and standard deviation of 61,842 ± 33,413). Sequences were aligned to the SILVA SEED database (version 123) and chimeras were removed using UCHIME (Edgar et al., 2011) in mothur. All samples were subsampled to 18,980 reads, which eliminated two anomalously low-read samples from ULN; however, these two samples were drawn from the N treatment where we saw no functional response to N. Sequences were clustered into operational taxonomic units (OTUs) based on a 97% identity threshold using the furthest-neighbor algorithm in mothur. Taxonomy was determined by classifying the 16S rRNA sequences based on the SILVA (version 123) database in mothur.

Analysis of Variance (ANOVA), Tukey honestly significant difference (HSD) tests, regression, self-organizing maps (SOM), weighted gene correlation network analysis (WGCNA), and permutational analysis of variance (PERMANOVA) were carried out in the R statistical environment (RStudio Version 1.0.136). Significant differences in NCP and CR were determined using ANOVA and Tukey HSD post-hoc tests; data were first tested for normality using the Shapiro-Wilk test and were log-transformed when necessary to achieve normality prior to ANOVA and Tukey HSD post-hoc tests. Regression was performed using the lm function in R.

SOM is a neural network algorithm previously applied to analysis of complex microbial communities (Wehrens and Buydens, 2007; Bowman et al., 2017); it was used to reduce the dimensionality of the community using the kohonen R package (Wehrens and Buydens, 2007). A Hellinger-transformed OTU abundance matrix of the 1,500 most abundant OTUs was used (after 1,500 the frequency of OTUs in any sample was likely to be zero), and samples were assigned to a series of 3 × 2 circular nodes in a hexagonal, non-toroidal grid. A total of 100 iterations were run with an alpha from 0.05 to 0.01. Clusters of nodes were identified using k-means clusters, with k = 4 selected via evaluation of a scree plot of within-clusters sum of squares against different k values and finding the inflection point.

Weighted gene correlation network analysis (WGCNA) was used to relate “subnetworks” of microbial taxa to experimental treatments and functional responses to treatments. WGCNA is a network analysis tool originally created to find and summarize clusters of highly correlated genes, and then relate clusters to one another and external traits (Langfelder and Horvath, 2007; 2008), but it has also been applied to microbial community data (Duran-Pinedo et al., 2011; Aylward et al., 2015; Guidi et al., 2016; Wilson et al., 2018; Wilson et al., 2021). We used WGCNA to assess co-occurrence of OTUs via an adjacency function that (i) factors in the degree of shared neighbors between two OTUs and (ii) is magnified by a power-law function, so that the topology of the graph becomes scale-free. OTUs are sorted into subnetworks (or “modules”), and the co-occurrence profiles for each subnetwork are represented by the first principal component of that subnetwork (aka eigengene) which can be compared to variables of interest. We used a Hellinger-transformed OTU abundance matrix of the 1,500 most abundant OTUs; a signed adjacency measure was calculated for each pair of OTUs by raising a function involving their Pearson correlation coefficient to a soft thresholding value of 16. This thresholding value optimized a scale-free topological overlap measure that considered a weighted pairwise correlation of each pair of OTUs, as well as weighted correlations with other OTUs in the subnetwork. We used a minimum module size of 20 OTUs and merged modules that were more than 70% similar (assessed via clustering of the eigengene of each module), resulting in 16 total subnetworks. Eigengenes for each subnetwork were compared with rate data using regression in R.

We used permutational analysis of variance (PERMANOVA; in R), as well as homogeneity of molecular variance (HOMOVA), Libshuff, and UniFrac (all as implemented in mothur) to test for significant differences in microbial communities. Using the adonis function in the R package vegan (Oksanen et al., 2013), PERMANOVA was conducted on θ values (Yue and Clayton, 2005), with different lakes, experimental treatments, and their interaction as potential explanatory variables. HOMOVA tests whether the overall diversity of two communities is statistically different (Schloss, 2008); Libshuff randomizes sequences and tests for significant differences between observed patterns and null models (Schloss et al., 2004); and UniFrac examines differences in phylogenetic structure (Lozupone and Knight, 2005). Distance matrices were calculated for each experiment in mothur; then we used HOMOVA, Libshuff (form = default), and weighted UniFrac to test for significant differences (1,000 iterations for each test) between treatments and the initial community, as well as communities in corresponding control bottles. Non-metric multidimensional scaling (NMDS) was also conducted in mothur using multiple beta-diversity indices; the θ index (Yue and Clayton, 2005) produced the lowest stress value (0.12) and highest r² value (0.82) in the fewest dimensions (four).

Minimum entropy decomposition (MED; Eren et al., 2015) was performed on all sequences in the subsampled dataset (after removing sequence gaps and padding shorter reads). Maximum nucleotide variation in each node was set to 7, and all other default parameters were used. The network of MED results (Eren et al., 2015) was visualized in Gephi using a Force Atlas projection. Cohen’s d effect sizes were calculated for MED nodes in each of the experimental treatment types compared with controls. Cohen’s d demonstrates the strength of an experimental treatment by dividing the mean difference between treatments by the pooled standard deviation of the two means, and therefore shows the relative size and direction of change for each treatment and rate (Cohen, 2013). MED nodes were also compared with rate data using regression in R.

Finally, once subnetworks that correlated with variables of interest were identified, we searched taxa with module membership >0.70 in each subnetwork (indicating that they are central nodes in the subnetwork) for potential homologs within the TARA Ocean database (http://tara-oceans.mio.osupytheas.fr/ocean-gene-atlas/; Villar et al., 2018).

Sequence data are available in the Sequence Read Archive under accession number PRJNA555354: https://www.ncbi.nlm.nih.gov/sra/PRJNA555354.

Microbial community structure varied across marine lakes and was significantly altered by experimental enrichment with N, eutrophication with C addition, and especially the combined addition of N and C. We first used SOM to identify broad community types across lakes and experimental treatments, and found that marine lake communities fell into 4 SOM taxonomic groupings (between sums of squares/total sums of squares = 0.794). For initial communities, the holomictic lakes (HLO, MLN, NLN, and ULN) grouped together into a single collective community type, whereas each of the three meromictic lakes (NLK, OTM, and TLN) exhibited a distinct community type (Table S1). We observed similar patterns with NMDS: meromictic lakes separated from holomictic lakes along axis 1, and meromictic lakes were distinguished from each other along axis 2 (Figure 2).

Holomictic lake communities were dominated by Synechococcus, SAR11, and marine Flavobacteria (Figure 1B), all of which are typically found throughout wide areas of the ocean (Meyerhof et al., 2016; Rapacciuolo et al., 2019). The meromictic lakes (NLK, OTM, and TLN) differed from holomictic lakes and each other (Figure 2). While many marine microbial OTUs observed in the holomictic lakes were also present in large proportions in meromictic lakes (Figure 1B), other groups were observed only in meromictic lakes, leading to distinct initial communities across lakes (PERMANOVA: lake r² = 0.47, P < 0.001).

We subsequently perturbed these communities through N enrichment and C eutrophication. N, C, and especially N+C frequently induced significant shifts in microbial communities and C cycling rates (Wilson et al., 2017; Figure 2, Table 2). In spite of differences in initial communities, experimental treatments had a significant effect on community composition across lakes (PERMANOVA: treatment r² = 0.18, P < 0.001). This effect was strongest in N+C treatments, with four lakes—the holomictic lakes NLN and ULN, the transitional lake TLN, and the meromictic lake OTM—converging toward similar composition (Figure 2). SOM similarly showed that N+C additions to lakes NLN and ULN produced a new community type, while lake OTM also changed to this community type following C additions to dark bottles and N+C additions (Table S1). These results indicate that N and C addition can drive initially dissimilar communities toward similar composition.

We used three additional and complementary statistical approaches to test whether experimentally driven shifts in microbial communities were significant: (1) UniFrac to test for differences in phylogenetic structure; (2) HOMOVA to test for differences in overall diversity; and (3) Libshuff to test for differences in observed OTU patterns versus null models (see Materials and methods). For all of these tests, we conducted two sets of statistical comparisons (Text S1): (i) tests for significant changes in comparison with initial starting communities (Table S2) and (ii) tests for significant differences between control incubations and experimental additions of N, C, and N+C (Table 2). Although N and C drove changes on their own (Text S1), we found that combined N and C additions consistently altered microbial communities based on all tests (Table 2). N+C additions significantly shifted communities in all cases based on UniFrac, and in all but light bottles in lakes MLN and NLK based on both HOMOVA and Libshuff. Overall, N, C, and N+C additions significantly altered microbial communities in 67% (Libshuff) to 100% (UniFrac) of the experiments. Experimental effects were more sporadic for N or C additions on their own, but were highly consistent in N+C additions.

To examine these patterns in greater detail, we used minimum entropy decomposition (MED) to partition our dataset into “MED nodes” based on Shannon entropy (Figure 3A); MED resolves fine scale differences in microbial communities and is therefore ideal for examining shifts across lakes and between treatments (Eren et al., 2015). Of 179 total MED nodes, 142 showed at least a weak response (Cohen’s d > 0.2) to one or more experimental addition, while 65 MED nodes showed a moderate response (d > 0.5), and 21 showed a strong response (d > 0.8; visualized in Figure 3B). Many nodes responded to more than one treatment group; however, most of the strong responses were to N+C additions (14 in light bottles, 16 in dark bottles). In contrast, there were several moderate responses to N (7 light, 9 dark), and more moderate responses to C (22 light, 13 dark) than to N+C (15 light, 8 dark). We classified consensus sequences for each MED node and found that the strongest responses were typically observed in multiple groups of gamma- and alpha-proteobacteria (Figure 3B; see below). Overall these data demonstrate that N or C addition alone elicited more moderate responses, whereas the combination of N+C strongly affected multiple MED nodes.

Our measures of microbial function are consistent with the finding that N+C additions elicited the strongest microbial community responses, as we previously found that CR and NCP showed sharp changes in response to combined N enrichment and C eutrophication (Table 2; Wilson et al., 2017). N+C co-additions increased CR in all lakes except MLN (and was not added to HLO), leading to sharp changes in NCP, and a switch from net autotrophy to heterotrophy (Wilson et al., 2017). In contrast, individual N and C additions affected rates less frequently (Table 2). Notably, the effects of combined N+C substantially exceeded those of N or C alone, indicating that interactions among heterotrophic bacteria result in consumption of dissolved organic matter in greater amounts than through the summed individual effects of N and C.

With the aim of understanding links between microbial community structure and function, we used weighted gene correlation network analysis (WGCNA; Langfelder and Horvath, 2007; 2008) to identify significant community and rate responses to N+C additions. WGCNA has proven useful for partitioning microbial communities into different subnetworks within the broader community, and then identifying which subnetworks are most strongly related to environmental variations or variations in function (Duran-Pinedo et al., 2011; Aylward et al., 2015; Guidi et al., 2016; Wilson et al., 2018; Wilson et al., 2021). Regardless of whether we examined the entire dataset, or only N+C treatments compared with controls, and regardless of whether we examined absolute or relative changes in rates and subnetwork abundance, three subnetworks correlated significantly and consistently with CR, NCP, or both (Figure 4A and B). Vibrio OTUs comprised the majority of one of these subnetworks (66 of 108 OTUs in subnetwork 13; Figure 4C), and this subnetwork correlated positively with CR (r² = 0.594, P < 0.001; Figure 4B) while inversely with NCP (r² = 0.599, P < 0.001; Figure 4A). Subnetwork 13 responded to individual N and C additions in all lakes, and these community structure responses were strongest where rate responses were strongest. Another subnetwork (subnetwork 11) correlated inversely with NCP (r² = 0.309, P < 0.001; Figure 4A) and contained a mix of gammaproteobacteria (particularly Alteromonas), several groups of alphaproteobacteria (SAR11 and Rhodobacteria), and Flavobacteria (Figure 4C). A third subnetwork (subnetwork 15), composed of Marinomonas gammaproteobacteria, Flavobacteria, and other groups, was related inversely to NCP across the whole dataset (r² = 0.599, P < 0.005; Figure 4).

These data are also consistent with results from MED analyses. In addition to responses to the experimental manipulations identified above, we found that 11 MED nodes correlated positively with CR rates (r² = 0.351–0.564, all P < 0.005). These nodes all correlated inversely with NCP (r² = 0.314–0.719, all P < 0.005), as did two others (r² = 0.251 and 0.3721, both P < 0.05). For both WGCNA and MED analyses, many subnetworks and nodes correlated both positively with CR and negatively with NCP rates, due to the strong effects of the biogeochemical additions. However, cases where subnetworks and nodes correlated with only one rate—rather than both NCP and CR—likely reflect the fact that decreases in NCP in amended light bottles were typically smaller than increases in CR in dark bottles for the same treatment groups.

In addition to the microbial network responses identified above, two subnetworks in the meromictic lake NLK showed responses that were overlaid on the patterns above. Subnetwork 9 was comprised almost entirely of the broadly distributed marine cyanobacteria Synechococcus, while subnetwork 1 was comprised of chloroplast sequences and multiple marine microbial groups (e.g., SAR11, SAR116; Table S3). These subnetworks were present in the other lakes, but were more abundant in lake NLK. In contrast with the other lakes, lake NLK showed a significant increase in NCP and GPP in N additions to light bottles (Table 2; Wilson et al., 2017), indicating that photosynthetic rates increased due to added N. WGCNA results are further consistent with a phytoplankton response to N. This response—partitioned into multiple subnetworks within lake NLK—led to a different response in overall community structure compared with other lakes (Figures 2 and 3).

Finally, previous work using WGCNA to identify relationships between the global ocean microbiome and C export similarly found that many of the taxa identified here were present in a subnetwork that correlated with C export across the ocean, including members of the genera Vibrio and Pseudoalteromonas (Guidi et al., 2016). Patterns in the North Atlantic based on sequencing and oxygen measurements also showed relationships between NCP and numerous gammaproteobacterial taxa, including Vibrio (Wang et al., 2018). Despite differences in approach, scale, and system, our results are notably consistent with these earlier findings, suggesting that several of the specific heterotrophic bacterial groups identified in our experiments may be broadly relevant throughout the surface ocean. We directly assessed this suggestion by comparing 16S rRNA sequences of OTUs in the significant subnetworks identified above to the frequency and distribution of 16S rRNA sequences in the TARA Oceans database (Villar et al., 2018). We found that the majority of OTUs in the significant subnetworks were closely related, and in some cases identical, to homologs found throughout the global oceans (Text S2, Table S4), including Alteromonas macleodii homolog OM-RGC.v1.001583373 (95.7%–99.5% identical), Vibrio campbellii homolog OM-RGC.v1.001516210 (97.6%–100% identical), and Marinomonas sp. MWYL1 homolog OM-RGC.v1.001593297 (95.2%–96.2% identical), which indicates that the bacterial groups responding to our experimental additions are present across a wide range of marine ecosystems (Figure 5).

We examined the effects of N enrichment and C eutrophication on microbial community structure and function through experimental additions of N alone, C alone, and combined N+C to a spectrum of starting communities. We found the strongest responses to combined N+C treatments, although in some experiments significant changes in community structure were also observed for additions of N or C alone. However, individual N or C additions were not typically accompanied by significant changes in NCP and CR rates (Table 2). These results demonstrate that when N or C are added alone, such biogeochemical rates can remain statistically similar in spite of significant changes in microbial community structure. In contrast, combined N+C perturbations frequently drove strong simultaneous responses in communities and rates (Figure 2, Table 2). (The only exception to this pattern was lake MLN, which showed a significant shift in community structure, but rate changes that did not reach statistical significance.) These findings more broadly indicate that microbially mediated biogeochemical processes may be resilient to a single perturbation (i.e., the addition of either a potentially limiting nutrient or labile carbon), but susceptible to multiple, coincident perturbations, and respond in a nonlinear way. While performing multi-factorial experiments is challenging, our results underline the fact that performing complex experiments is essential in order to fully understand and forecast responses in microbial community structure and function to coincident changes (Boyd and Hutchins, 2012; Hutchins and Fu, 2017).

The combination of N+C selected for particular nodes and subnetworks that correlated strongly with NCP and CR rates (Figure 4). Several of the bacterial groups responding to enrichment and eutrophication were closely related to those identified as biogeochemically important for C-cycling in other marine ecosystems—such as the Alteromonadales, Flavobacteriales, and Rhodobacterales (Pedler et al., 2014; Guidi et al., 2016; Wang et al., 2018). Members of the orders Alteromonadales and Rhodobacterales, as well as the genera Vibrio and Marinomonas, typically increased in relative abundance from <1% to >10% of the community, particularly following N+C additions. We further show that the bacteria responding to enrichment and eutrophication in our experiments closely match homologs that are widely distributed throughout the oceans (Figure 5, Table S4). Combined together, our results and earlier work indicate that these and other microbial groups contribute to rapid biogeochemical cycling under N- and C-rich conditions (Figures 3 and 4), such as those found during our experiments, during phytoplankton blooms, and on sinking particles in the ocean (Pedler et al., 2014; Guidi et al., 2016).

However, our systems differ from open ocean locations in that marine lakes receive significant organic matter inputs from land, in addition to in situ production of organic C by phytoplankton and other marine primary producers (Hamner et al., 1982; Hamner and Hamner, 1998; Colin, 2009; Wilson et al., 2017; Wilson et al., 2019). Terrestrial plants produce a wide range of C compounds, including sucrose (Salerno and Curatti, 2003), which we used in C and N+C additions as a labile C substrate employed in previous C/N manipulation experiments in the ocean (Babbin et al., 2014). Yet in spite of the addition of labile C, not all marine lake microbial communities responded—again highlighting the additional importance of N (Table 2). Our data also underscore the well-known concept of dual macronutrient limitation of biological activity by N and phosphorus (P; i.e., the Redfield ratio), as the two lakes that were least responsive to experimental N enrichment and C eutrophication—MLN and NLK—were deficient in P (Table 1). The lack of a strong response in these lakes may be indicative of P limitation of photosynthesis and/or community respiration. Notably, N and C were added at concentrations intending to mimic enrichment and eutrophication (Hamner et al., 1982; Colin, 2009; Meyerhof et al., 2016) and in a ratio (C:N of 20:1) expected to produce a growth and CR rate response (del Giorgio and Cole, 1998; Taylor and Townsend, 2010). Applying this approach to different lakes with different starting conditions and communities allowed us to build a more comprehensive view than studies that (i) have less initial community variation, (ii) rely on observation alone without experimentation, (iii) manipulate only a single variable, and/or (iv) measure only community properties or (meta)genomic profiles without true biogeochemical rate information.

Using this overall approach, we found that N and C additions produced a mosaic of responses in multiple diverse microbial networks. For example, across all lakes, all three subnetworks that correlated with the measured rates were present in similar proportions initially. However, subnetwork 13 responded to N+C additions in all lakes, whereas the other two significant subnetworks varied in their responses between lakes. Subnetwork 11 was particularly important in lake TLN, while subnetwork 15 responded strongly in lakes NLN, ULN, and OTM. MED results show a similar pattern, with some nodes being more influential in different lakes (Figure 3). Community responses were therefore multi-dimensional, with a universal primary response in subnetwork 13 that was accompanied by responses in additional subnetworks. Slight differences among these networks, and the organisms within them, may reflect underlying genomic differences that also affect their distributions throughout the ocean (Figure 5). The fact that subnetworks 11 and 15 are more diverse in composition (Figure 4C) also fits with the idea that although a few bacterial groups can dominate C flow, biogeochemical cycling rates are further enhanced by a more diverse community (Bell et al., 2005; Pedler et al., 2014). This relationship between diversity and rates in our experiments aligns with biodiversity-ecosystem functioning concepts developed for different communities and processes in other ecosystems (Naeem and Li, 1997; Cardinale, 2011; Naeem et al., 2012).

In addition, the meromictic lake NLK showed a significant increase in NCP and GPP in N additions to light bottles (Table 2; Wilson et al., 2017), with a parallel positive response in two subnetworks comprised of phytoplankton (Synechococcus and chloroplast sequences) and additional bacteria. The response of these subnetworks to N resulted in an overall shift in community structure that was distinct from other lakes. This finding indicates that competition for N among different microbial groups—in this case, heterotrophs and phytoplankton—can lead to different outcomes. Another key factor may be the form of N, as different microbial groups have varying abilities to use different forms of inorganic and/or organic N. Finally, decreases in NCP driven by N+C additions were always lower than corresponding increases in CR (Table 2). This quantitative difference indicates that a positive effect on GPP occurs simultaneously with a stronger effect on CR. This response was particularly strong in lake NLK, resulting in further additional layers of subnetworks shifting in response to additions.

N enrichment and C eutrophication are poorly understood and difficult to disentangle. Our results indicate that their collective effects on microbial community structure can exceed their independent effects, resulting in net heterotrophy driven by globally distributed marine bacteria. Although our results may be specific to the types and amounts of N and C substrates that we used in our experimental manipulations, the underlying potential for enrichment and eutrophication to have strong effects on microbial function and biogeochemical cycling appears to be widespread in the ocean. However, our work shows in several ways that the interaction between N, C, and microbial community structure and function is multi-layered. On its own, labile C rarely had strong effects; only combined addition with N resulted in a tipping point for both community structure and function. Future work could determine where this tipping point lies for different concentrations and forms of N and C, and for different microbial communities. In particular, our experimental manipulation of distinct starting communities indicates that multiple subnetworks respond simultaneously to N and C. These and additional subnetworks likely respond to different forms, ratios, and concentrations of N and C in complex ways, regulating the strength and severity of the overall biogeochemical response to enrichment and eutrophication. Given the implications for the oceanic C cycle (Duarte et al., 2013; Ducklow and Doney, 2013; Williams et al., 2013), quantifying their relative activity and balance under a range of conditions—including variations in N and C—is critical, as their collective interactions control the biogeochemical fate of N and C in the changing ocean.

Sequence data are available in the Sequence Read Archive under accession number PRJNA555354: https://www.ncbi.nlm.nih.gov/sra/PRJNA555354.

The supplemental files for this article can be found as follows:

Text S1.pdf; Table S1.Docx; Table S2.Docx; Table S3.Docx; Table S4.Docx.

We thank the Coral Reef Research Foundation, Dr Michael Dawson, and NECO Marine for supporting fieldwork in Palau. Permits were issued by the Bureau of Marine Resources (RE-13-11, RE1409, RE-16-08, RE-15-08, RE-17-15) and Koror State (13-233, 14-05, 15-10, 16-15, 17-07).

This work was supported by the United States National Science Foundation Dimensions of Biodiversity Program (OCE-1241255 to JMB and Michael Dawson) and by the UC Merced Graduate Division (to JMW).

JMB holds equity in Lillianah Technologies, with prior approval of the University of California, Merced. The other authors declare no competing interests, financial or otherwise.

JMW and JMB designed the research; JMW and SSA performed the research; JMW and JMB contributed new analytical tools; JMW and JMB analyzed the data; JMB and JMW wrote the paper; JMW, JMB, and SA edited the paper and approved the submission.

How to cite this article: Wilson, JM, Abboud, SS, Beman, JM. 2024. Effects of experimental nutrient enrichment and eutrophication on microbial community structure and function in “marine lakes.” Elementa: Science of the Anthropocene 12(1). DOI: https://doi.org/10.1525/elementa.2024.00007

Domain Editor-in-Chief: Jody W. Deming, University of Washington, Seattle, WA, USA

Knowledge Domain: Ocean Science