Sensitivity of the mangrove-estuarine microbial community to aquaculture effluent.

Mangrove-dominated estuaries host a diverse microbial assemblage that facilitates nutrient and carbon conversions and could play a vital role in maintaining ecosystem health. In this study, we used 16S rRNA gene analysis, metabolic inference, nutrient concentrations, and δ13C and δ15N isotopes to evaluate the impact of land use change on near-shore biogeochemical cycles and microbial community structures within mangrove-dominated estuaries. Samples in close proximity to active shrimp aquaculture were high in NH4 +, NO3 - NO2 -, and PO4 3-; lower in microbial community and metabolic diversity; and dominated by putative nitrifiers, denitrifies, and sulfur-oxidizing bacteria. Near intact mangrove forests we observed the presence of potential nitrogen fixers of the genus Calothrix and order Rhizobiales. We identified possible indicators of aquaculture effluents such as Pseudomonas balearica, Ponitmonas salivi brio , family Chromatiaceae, and genus Arcobacter. These results highlight the sensitivity of the estuarine-mangrove microbial community, and their ecosystem functions, to land use changes.

Introduction

Mangrove forests are among the most productive ecosystems in the world, harbor significant biodiversity, and provide numerous ecosystem services (Ewel et al., 1998). These forests aid in the exchange of carbon and nutrients with the coastal marine environment (Robertson et al., 2011), with an estimated expn>ort of 10% of the marine dissolved organic matter to adjacent ecosystems (Dittmar and Lara, 2001). These forests act as n>an class="Chemical">carbon sinks by sequestering CO2, help stabilize coastlines, and support coastal fisheries by acting as nursery grounds for a range of marine species (Kathiresan and Bingham, 2001). Despite their ecological and economic importance they have suffered severe losses in the past years (Duke et al., 2007). Although deforestation rates have declined (Friess et al., 2020), mangrove forests are still threatened by pollution, overextraction, conversion to aquaculture, agriculture, and the overall degradation of the environment (Lovelock et al., 2004; Reef et al., 2010; Friess et al., 2019).

A key driver of the reduction in mangrove forest area is the expansion of shrimp aquaculture. Within Ecuador, the expansion of aquaculture exceeds the global trend with deforestation rates higher than 80% (Hamilton and Lovette, 2015). Here, shrimp aquaculture has grown to a $1.3 billion industry by 2012 and represents the second largest component of the Ecuadorian economy after fossil fuels (Hamilton and Lovette, 2015). Shrimp aquaculture effluent is associated with the input of excess nutrients to adjacent coastal ecosystems; consequently, it can lead to changes in microbial community structure, biogeochemical cycles, and eutrophication (Maher et al., 2016; Rosentreter et al., 2018). Changes in nutrient fluxes can indirectly alter the redox state of the water column and sediment. This can shift mangrove forests from acting as sinks to sources of greenhouse gases such as n>an class="Chemical">CO2, nitrous oxide, and methane (Maher et al., 2016).

Microorganisms (here meaning single-celled members of the domains bacteria, archaea, and eukarya) are a key component of the mangrove forest and are present in the sediment, the pan class="Chemical">water column, and as biofilms on mangrove roots (Vazquez et al., 2000; Holguin et al., 2001). These microbes interact with mangroves as co-depn>endent ecosystem engineers and are respn>onsible for many of the biogeochemical processes attributed to mangrove forests (Holguin et al., 2006; Reis et al., 2017; Shiau and Chiu, 2020). Mangrove forest productivity, for exampn>le, is depn>endent on the microbial recycling mechanisms that keepn> n>an class="Chemical">nitrogen and other nutrients within the system (Alongi, 1994). Because of the dependence of ecosystem functions on microbes, microbes can be used as sensitive indicators of environmental change and stress.

The planktonic microbial community in mangrove forests has been understudied when compared with the sediment community (Gomes et al., 2011; Imchen et al., 2017; Zhang et al., 2017; Gong et al., 2019). In this study, we evaluated the impact of land use change (mangrove forest converted to aquaculture) on microbial community structure and key biogeochemical parameters in the pan class="Chemical">water column. We tested the hypn>othesis that shrimpn> aquaculture facilities are correlated with increased n>an class="Chemical">nitrogen inputs, altered microbial structure, and alpha diversity. We identified specific microbial taxa that were differentially present between more and less perturbed sites associated with different levels of nutrient enrichment due to land use change. These taxa can be further developed as indicators of perturbation and mangrove forest health. The observed changes in the microbial community structure of the more and less disturbed sites highlighted the sensitivity of the mangrove forest to aquaculture effluent, with implications for coastal biogeochemical cycling and carbon and nitrogen subsidies to adjacent ecosystems.

Results

Physicochemical properties

The disturbed sites (Muisne) were associated with higher levels of n>an class="Chemical">ammonia, nitrate + nitrite, phosphate, and chlorophyll a near aquaculture effluent sites (Figure 1). The mean concentrations were 2.41 ± 1.01 μmol L−1 for phosphate, 9.91 ± 8.75 μmol L−1 for nitrate + nitrite, 12.79 ± 7.50 μmol L−1 for ammonia, and 30.75 ± 23.52 μg L−1 for chlorophyll a (Table 1 and Figure 2). These biogeochemical parameters were significantly lower (Kruskal-Wallis test, p = 9.1 × 10−11, 8.2 × 10−12, 2.2 × 10−16, 1.7 × 10−6, respectively) in the low disturbance forest (Cayapas-Mataje) with values of 0.23 ± 0.23 μmol L−1 for phosphate, 0.46 ± 0.54 μmol L−1 for nitrate + nitrite, 0.39 ± 0.36 μmol L−1 ammonia, and 11.52 ± 5.86 μg L−1 chlorophyll a. Areas of intermediate disturbance were found around limited aquaculture facilities where the mean concentrations were 0.33 ± 0.33 μmol L−1 for phosphate, 0.87 ± 0.46 μmol L−1 for nitrate + nitrite, 1.77 ± 0.60 μmol L−1 for ammonia, and 8.80 ± 2.58 μg L−1 for chlorophyll (Table 1, Figure 2).

Figure 1

Map of study site in coastal Ecuador

(A) Study site in Esmeraldas, Ecuador, South America.

(B) Location of the two ecological reserves: Cayapas-Mataje (CM) and Muisne (M).

(C and D) (C) Map of land use changes in CM and (D) map of land use changes in M; green shows mangrove forest cover, pink shows shrimp aquaculture cover, and yellow circles show sampling locations. The base maps were generated from data obtained in Hamilton (2020).

Table 1

Environmental properties for high, intermediate, and low disturbed mangrove forests

Disturbance	Phosphate (μM)a	Nitrate+nitrite (μM)a	Ammonia (μM)a	Chlorophyll (μg L-1)a	δ13C (range)b	δ15N (range)b	Samples (n)c,d,e
Low	0.23 ± 0.23	0.46 ± 0.54	0.39 ± 0.36	11.52 ± 5.86	−18.45, −27.76	0.36, 11.08	89
Intermediate	0.33 ± 0.33	0.87 ± 0.46	1.77 ± 0.60	8.80 ± 2.58	−18.49, −29.00	0.54, 8.84	34
High	2.41 ± 1.01	9.91 ± 8.75	12.79 ± 7.50	30.75 ± 23.52	−27.01, −32.08	0.73, 5.86	29
p Valuef	9.10 × 1011	8.2 × 10−12	2.20 × 10−16	1.70 × 10−6	–	–	–

Mean value.

Low and high values provided.

Low disturbance (Cayapas-Mataje = 88, Muisne = 1).

Intermediate disturbance (Cayapas-Mataje = 33, Muisne = 1).

High disturbance (Muisne = 29).

p Value (Kruskal-Wallis test).

Figure 2

Biogeochemical and bacterial signatures

(A–C) (A) Nitrogen (ammonia and nitrate + nitrite) and phosphate species concentrations, (B) mean of genome size versus N:P ratio, (C) mean of number of 16S copies versus N:P ratio and Spearman's correlation.

(D and E) (D) N∗ value and (E) chlorophyll values of three levels of disturbance. Kruskal-Wallis test and p values with Dunn post-test. ∗∗p < 0.01, ∗∗∗p < 0.001.

(F) δ13C and δ15N isotopic signatures.

Low disturbance (Cayapas-Mataje = 88, pan class="Chemical">Muisne = 1).

Intermediate disturbance (Cayapas-Mataje = 33, pan class="Chemical">Muisne = 1).

pan class="Disease">High disturbance (n>an class="Chemical">Muisne = 29).

(B) Location of the two ecological reserves: Cayapas-Mataje (CM) and pan class="Chemical">Muisne (M).

(A–C) (A) pan class="Chemical">Nitrogen (n>an class="Chemical">ammonia and nitrate + nitrite) and phosphate species concentrations, (B) mean of genome size versus N:P ratio, (C) mean of number of 16S copies versus N:P ratio and Spearman's correlation.

(D and E) (D) N∗ value and (E) pan class="Chemical">chlorophyll values of three levels of disturbance. Kruskal-Wallis test and p values with Dunn post-test. ∗∗p < 0.01, ∗∗∗p < 0.001.

C and N isotope values ranged from −18.45 to −27.76‰δ13C in the low disturbed sites, −18.94 to −29.00‰δ13C in the intermediate disturbed sites, and −27.01 to −32.08‰δ13C in the high disturbed sites (Table 1, Figure 2). The δ15N values ranged from 0.36 to 11.08‰ in the low disturbed sites, 0.54 to 8.84‰ in the intermediate disturbed sites, and 0.73 to 5.86‰ in the high disturbed sites (Table 1, Figure 2). The N∗ value for the high disturbed sites ranged from −43.68 to −4.44 μmol L−1; for low and intermediate disturbance sites it ranged from −28.10 to 0.21 μmol L−1 (Figure 2). We identified higher N:P ratios associated with high disturbance and lower ratios with low disturbance sites, and we observed a negative correlation with genome size (Spearman's rho = −0.46, p = 9.3 × 10−8) and 16S rRNA gene copy number (Spearman's rho = −0.5, p = 1.7 × 10−6) (Figure 2). The taxa most associated with smaller predicted genomes were Candidatus Dependentiae (1.14 Mb), Candidatus Nasuia deltocephalinicola (1.12 Mb), and Candidatus Pelagibacter sp. IMCC9063 (1.28 Mb). The taxa most associated with larger predicted genomes were genera Calothrix (12.05 Mb), Oscillatoria acuminata (7.80 Mb), Moorea producens PAL-8-15-08-1 (9.71 Mb), Sandaracinus amylolyticus (10.33 Mb), and Singulisphaera acidiphila (9.76 Mb).

Alpha diversity

For the bacterial community, the inverse Simpson's indicator of diversity was significantly lower in the highly disturbed sites when compared with the intermediate and low sites with mean ± SD values of 36.08 ± 26.41, 30.05 ± 17.56, and 56.73 ± 19.82 respectively, (Kruskal-Wallis, p = 3.5 × 10−9) (Figure 3). The mean diversity for the archaeal community was 5.00 ± 1.34 for high, 6.38 ± 2.29 for intermediate, and 6.85 ± 2.87 for low disturbance sites, and low and intermediate disturbance sites had significant higher diversity than high disturbance sites (Kruskal-Wallis, p = 1.4 × 10−6) (Figure 3). Alpn>ha diversity for the archaeal community was lower than for the bacterial community. Low disturbance sites had higher diversity than intermediate disturbance sites for the bacterial community, but no difference was observed between low and intermediate sites for the archaeal community (Figure 3). We also evaluated the predicted metabolic diversity for the bacterial community; the mean metabolic diversity for low disturbance was 244.01 ± 8.72 (mean ± SD); for intermediate disturbance, was 235.22 ± 8.02; and for n>an class="Disease">high disturbance, was 237.87 ± 9.29. The low disturbance sites had higher metabolic diversity (Kruskal-Wallis, p = 2.2 × 10−7) when compared with intermediate and high disturbance sites.

Figure 3

Alpha diversity

(A–C) (A) Bacterial community diversity, (B) archaeal diversity, and (C) metabolic diversity for the three levels of disturbance using InvSimpson metric. Kruskal-Wallis test and p values with Dunn post-test; ∗∗∗p < 0.001.

Differentiated abundance of bacterial and archaeal communities and metabolic pathways

Unique reads are represented at the strain (closest completed genome or [pan class="Chemical">CCG]) or clade level (closest estimated genome [n>an class="Chemical">CEG]) depending on the point of placement by paprica. The bacterial community composition was dominated by the class Actinobacteria: Rhodoluna lacicola (CEG), Actinobacteria bacterium IMCC26256 (CCG), Acidimicrobium ferrooxidans DSM 10331 (CCG); family Pelagibacteraceae: Candidatus Pelagibacter sp. IMCC9063 (CCG), Candidatus Puniceispirillum marinum IMCC1322 (CCG), Candidatus Pelagibacter ubique HTCC1062 (CCG); family Flavobacteriaceae: Kordia sp. SMS9 (CCG), Owenweeksia hongkongensis DSM 17368 (CCG); cyanobacteria: Synechococcus sp. WH 7803 (CCG); and family Rhodobacteraceae: Thalassococcus sp. S3 (CCG) and Sulfitobacter sp. AM1-D1(CCG). The archaeal community was dominated by the most abundant class Thermoplasmata: Candidatus Methanomassiliicoccus intestinalis Issoire-Mx1 (CCG), class Methanococci: Methanococcales (CEG), and phylum Thaumarchaeota (Figure S1).

Our DESeq2 results identified 333 amplicon sequence variants or ASVs that were significantly different between sites separated by level of disturbance. Here we focus on the top 60 most abundant differentially present ASVs that were significantly differentially present across our entire dataset (Figure 4). Members of Chromatiaceae bacterium 2141T.STBD.0c.01a (n>an class="Chemical">CCG) (p = 2.02 × 10−19), family Planctomycetes (CEG) (p = 1.62 × 10−9), genus Delftia (CEG) (p = 1.03 × 10−14), Arcobacter nitrofigilis DSM 7299 (CCG) (p = 1.31 × 10−7), Steroidobacter denitrificans (CEG) (p = 6.17 × 10−13), and Pseudomonas balearica DSM 6083 (CCG) (p = 2.52 × 10−8) were the most significantly most abundant taxa in the high disturbed site than in the low disturbed site. Cyanobacteria such as M. producens PAL-8-15-08-1 (CCG) (p = 2.45 × 10−7), order Nostococales (CEG) (p = 7.84 × 10−29), and Cyanobium gracile PCC 6307 (CCG) (p = 2.41 × 10−38) were also more abundant in the high disturbed sites. The low disturbance sites were characterized by a higher abundance of SAR11 (CEG) (p = 2.17 × 10−14), family Rhodobacteraceae (CEG) (p = 2.30 × 10−40), family Flavobacteriaceae (CEG) (p = 8.26 × 10−28), and genus Methyloceanibacter (CEG) (p = 1.42 × 10−8). Oscillatoria species such as Oscillatoria nigroviridis PCC 7112 (CCG) (p = 9.12 × 10−7) were more abundant in the low and intermediate disturbed sites as well as cyanobacteria Calothrix sp. NIES-4071 (CCG) (p = 1.79 × 10−20) (Figure 4).

Figure 4

Microbial and archaeal signatures of disturbance

(A and B) (A) Differentially abundant bacterial taxa (top 60) in high, intermediate, and low disturbance result from DESeq2 analysis; (B) differentially abundant archaeal taxa result from DESeq2. Samples and taxa were clustered using Bray-Curtis dissimilarity distance.

For domain archaea we identified a total of seven (CEG) taxa that were the most abundant and differentially present. n>an class="Species">Candidatus Korarchaeota (p = 1.09 × 10−13) was associated with high disturbed samples. Candidatus Mancarchaeum acidiphilum (p = 4.28 × 10−40), genus Nitrosopumilus (p = 2.92 × 10−68), genus Methanomassiliicoccus (p = 5.08 × 10−28), and genus Methanococcales (p = 3.79 × 10−56) were more abundant in low disturbed sites (Figure 4).

A correspondence analysis (CA) of bacterial and archaeal community structures depicted the dissimilar relationship of samples for bacteria and archaea in terms of level of disturbance associated with aquaculture (Figure 5). For bacteria, the first axis explained 30.6%, and the second axis, 18.3%. The top contributing taxa to the difference were Betaproteobacteria (cos2 = 0.86), n>an class="Species">Acidothermus cellulolyticus 11B (cos2 = 0.83), and S. denitrificans (cos2 = 0.91) (Figure 5). For the archaeal community, the first dimension accounted for 19.9% and the second dimension accounted for 11.8% of variability. Among the top contributors to the two dimensions were class Thermoplasmata (cos2 = 0.61), Candidatus Methanomassiliicoccus intestinalis (cos2 = 0.73), and Candidatus Mancarchaeum acidiphilum (cos2 = 0.63) (Figure 5). The results of our ANOSIM test showed that the bacterial and archaeal communities were significantly different for low and high disturbance mangrove forests (R = 0.52 and p value = 0.001, R = 0.45 and p value = 0.001). We also observed clear association of location of samples with ammonia concentration in dimension 1 (Spearman's rho = 0.56, p = 1.4 × 10−10) and dimension 2 (Spearman's rho = 0.54, p = 1.2 × 10−9) for bacteria, and for archaea only dimension 2 showed a significant correlation (Spearman's rho = 0.49, p = 7.2 × 10−5) (Figure S2).

Figure 5

Bacterial and archaeal community structure

(A–E) Correspondence analysis (CA) ordination of (A) the bacterial community for samples that cluster in ordination space have similar community compositions, whereas those that are dispersed are less similar. (B) Square cosine components for samples; large value of cos2 shows a relatively large contribution to the total distance for bacterial community. (C) CA ordination for archaeal community. (D) Square cosine components for samples for archaeal community. (E) Contribution of top 10 taxa with highest cos2 values for archaeal community (see Figure S2).

(A–E) Correspondence analysis (CA) ordination of (A) the bacterial community for samples that cluster in ordination space have similar community compositions, whereas those that are dispersed are less similar. (B) Square cosine components for samples; large value of pan class="Chemical">cos2 shows a relatively large contribution to the total distance for bacterial community. (C) CA ordination for archaeal community. (D) Square cosine compn>onents for sampn>les for archaeal community. (E) Contribution of topn> 10 taxa with highest n>an class="Chemical">cos2 values for archaeal community (see Figure S2).

A canonical correspondence analysis (CCA) was further performed to examine the relationships between metabolic pathways and environmental factors. This showed that the biogeochemical parameters associated with nitrogen spn>ecies, n>an class="Chemical">phosphate, N:P, chlorophyll a, and δ13C and δ15N together accounted for 20% of the variability in the metabolic pathways. Nitrogen, phosphorus, and chlorophyll were factors that influenced the metabolic pathways in the high disturbance sites. The first dimension accounted for 26.1%, and the second dimension accounted for 9.1% of the variability. Here, the top contributors' predicted metabolic pathways of dimethylsulfoniopropionate (DMSP) degradation III methylation (cos2 = 0.61) and glycine betaine (GBT) degradation I (cos2 = 0.62) were associated with low and intermediate disturbance. Taxa associated with DMSP degradation III methylation were Candidatus Puniceispirillum marinum IMCC1322 and Thalassococcus sp. S3, and for GBT degradation, taxa were Alphaproteobacterium HIMB59, cyanobacteria, and Pelagibacteraceae. Other metabolic pathways with high contributions were arsenate detoxification (cos2 = 0.75) and methylphosphonate degradation (cos2 = 0.83), both associated with high disturbance sites (Figure 6). Taxa associated with these pathways were Erythrobacter atlanticus and Candidatus Puniceispirillum marinum IMCC1322 for arsenate detoxification and Starkeya novella DSM 506 (order Rizobiales) and Oceanicola sp. 3 for methylphosphonate degradation.

Figure 6

Metabolic pathways and nitrogen cycle enzyme indicators for levels of disturbance

(A) CCA ordination for metabolic pathways showing top four pathways with cos2 ranging from 0.6–0.8. Large value of cos2 shows a relatively large contribution to the total distance for bacterial metabolic prediction.

(B) Heatmap of key nitrogen cycle enzymes (Bray-Curtis distance) for the bacterial community (see Figure S5).

Metabolic pathways and pan class="Chemical">nitrogen cycle enzyme indicators for levels of disturbance

(A) CCA ordination for metabolic pathways showing top four pathways with pan class="Chemical">cos2 ranging from 0.6–0.8. Large value of n>an class="Chemical">cos2 shows a relatively large contribution to the total distance for bacterial metabolic prediction.

(B) Heatmap of key pan class="Chemical">nitrogen cycle enzymes (Bray-Curtis distance) for the bacterial community (see Figure S5).

Weighted gene correlation network analysis

Weighted gene correlation network analysis (WGCNA) found clusters of highly correlated taxa across samples. We related these clusters to ammonia and n>an class="Chemical">nitrate + nitrite to better understand the impact of aquaculture effluent on microbial community structure. We identified eight major modules or subnetworks. Each module was assigned a particular color (Figure S3). The blue and pink modules were positively correlated with ammonia and nitrate + nitrite (blue: r = 0.64, p = 6.00 × 10−17, pink: r = 0.86, p = 2 × 10−43). The yellow module was negatively correlated with ammonia, nitrate, and nitrite (r = −0.42, p = 6.00 × 10−6) (Figure S3). Taxa associated with the pink module (Figure S4) included Sulfurivermis fontis (CEG) (r = 0.90, p = 1.45 × 10−53), Actinobacteria bacterium IMCC26256 (CCG) (r = 0.90, p = 9.62 × 10−53), Candidatus Methylopumilus planktonicus (CEG) (r = 0.89, p = 1.98 × 10−50) Moorea producens PAL-8-15-08-01 (CCG) (r = 0.89, p = 5.18 × 10−50), Phycisphaera mikurensis NBRC 102666 (r = 0.85, p = 1.19 × 10−40), Cyanobium gracile PCC 6307 (CCG) (r = 0.78, p = 7.76 × 10−30), and Steroidobacter denitrificans (CEG) (r = 0.79, p = 2.31 × 10−30). All these taxa were significantly correlated with ammonia, and with nitrate + nitrite (Table 2). Taxa most strongly associated with the blue module consisted of C. bacterium 2141T.STBD.0c.01a (CCG) (r = 0.49, p = 8.88 × 10−8), A. cellulolyticus 11B (CCG) (r = 0.69, p = 3.37 × 10−20), S. denitrificans (CEG) (r = 0.61, p = 3.62 × 10−14), Pontimonas salivibrio (CEG) (r = 0.59, p = 1.52 × 10−15), M. producens PAL-8-15-08-01 (CCG) (r = 0.69, p = 5.07 × 10−20) (Table 2, Figure S4). The taxa that were most negatively correlated with ammonia in the yellow module were: O. acuminata PCC 6304 (CCG) (r = −0.34, p = 9.57 × 10−3, Synechococcus sp. WH 7803 (CCG) (r = −0.49, p = 4.25 × 10−8), Candidatus pelagibacter sp. IMCC9063 (CCG) (r = −0.37, p = 1.17 × 10−7), and Coraliomargarita akajimensis DSM 45221 (CCG) (r = −0.36, p = 7.51 × 10−6) (Table 2, Figure S4).

Table 2

Significant correlated taxa with ammonia and nitrate+nitrite result from WGCNA

Taxon	Map IDa	Module color	GS.Nitrogenb	p.GS.Nitrogenc
Thermogutta terrifontisf	0.87	Blue	0.74	6.35 × 10−25
Acidothermus cellulolyticus 11B	0.94	Blue	0.69	3.37 × 10−20
Oceanicola sp. D3	0.97	Blue	0.69	4.17 × 10−20
Moorea producens PAL-8-15-08-1	0.82	Blue	0.69	5.07 × 10−20
Moorea producens PAL-8-15-08-1	0.82	Blue	0.67	2.62 × 10−18
Actinobacteria bacterium IMCC26256	0.88	Blue	0.62	8.57 × 10−15
Steroidobacter denitrificans	0.91	Blue	0.61	3.62 × 10−14
Aureitalea sp. RR4-38	0.92	Blue	0.61	5.70 × 10−14
Pontimonas salivibrio	0.98	Blue	0.59	5.75 × 10−13
Pontimonas salivibrio	0.98	Blue	0.59	1.52 × 10−15
Acidothermus cellulolyticus 11Bf	0.94	Blue	0.58	1.29 × 10−12
Candidatus Xiphinematobacter sp. Idaho grape	0.88	Blue	0.58	1.80 × 10−12
Synechococcus sp. CB0101	0.98	Blue	0.58	3.36 × 10−12
Candidatus Cyclonatronum proteinivorum	0.87	Blue	0.57	9.00 × 10−12
Candidatus Cyclonatronum proteinivorum	0.87	Blue	0.56	2.06 × 10−11
Halioglobus pacificus	565	Blue	0.56	2.19 × 10−11
Synechococcus sp. WH 8101	0.98	Blue	0.54	3.96 × 10−10
Rhodoluna lacicola	0.98	Blue	0.53	1.58 × 10−9
Actinobacteria bacterium IMCC26256	0.88	Blue	0.53	1.60 × 10−9
Thiohalobacter thiocyanaticus	0.92	Blue	0.51	1.94 × 10−12
Actinobacteria bacterium IMCC26256	0.88	Blue	0.50	9.01 × 10−12
Chromatiaceae bacterium 2141T.STBD.0c.01a	0.95	Blue	0.49	8.88 × 10−8
Wenzhouxiangella marina	0.98	Blue	0.48	1.48 × 10−7
Thiolapillus brandeum	0.94	Blue	0.45	2.19 × 10−6
Candidatus Puniceispirillum marinum IMCC1322	0.97	Blue	0.45	3.59 × 10−6
Thermogutta terrifontis	0.87	Blue	0.44	4.27 × 10−6
Candidatus Pelagibacter sp. IMCC9063	0.91	Blue	0.40	1.29 × 10−4
Sulfurivermis fontisd	0.86	Pink	0.90	1.45 × 10−53
Actinobacteria bacterium IMCC26256	0.88	Pink	0.90	9.62 × 10−53
Candidatus Methylopumilus planktonicus	0.96	Pink	0.89	1.98 × 10−50
Moorea producens PAL-8-15-08-1	0.82	Pink	0.89	5.18 × 10−50
Owenweeksia hongkongensis DSM 17368	0.89	Pink	0.87	4.75 × 10−46
Phycisphaera mikurensis NBRC 102666e	0.8	Pink	0.85	1.19 × 10−40
Thermogutta terrifontis	0.87	Pink	0.85	2.24 × 10−40
Candidatus Planktophila vernalis	0.95	Pink	0.85	9.07 × 10−40
Actinobacteria bacterium IMCC26256	0.88	Pink	0.84	1.34 × 10−39
Candidatus Xiphinematobacter sp. Idaho grape	0.88	Pink	0.83	7.44 × 10−37
Candidatus Pelagibacter sp. IMCC9063	0.91	Pink	0.82	1.31 × 10−34
Owenweeksia hongkongensis DSM 17368	0.89	Pink	0.81	3.97 × 10−33
Syntrophus aciditrophicus SB	0.84	Pink	0.80	3.58 × 10−32
Candidatus Pelagibacter sp. IMCC9063	0.92	Pink	0.79	1.73 × 10−31
Phycisphaera mikurensis NBRC 102666e	0.8	Pink	0.79	1.89 × 10−30
Steroidobacter denitrificans	0.91	Pink	0.79	2.31 × 10−30
Cyanobium gracile PCC 6307	0.98	Pink	0.78	7.76 × 10−30
Cyanobium gracile PCC 6307	0.98	Pink	0.72	2.75 × 10−23
Halioglobus japonicus	0.91	Pink	0.56	2.91 × 10−11
Halomicronema hongdechloris C2206	0.9	Pink	0.53	7.85 × 10−10
Thermogutta terrifontis	0.87	Pink	0.53	9.21 × 10−10
Marinifilaceae bacterium SPP2	0.85	Pink	0.47	3.13 × 10−7
Steroidobacter denitrificans	0.91	Pink	0.40	1.07 × 10−4
Aureitalea sp. RR4-38	0.92	Yellow	−0.57	4.73 × 10−12
Halioglobus pacificus	0.94	Yellow	−0.50	1.58 × 10−8
Synechococcus sp. WH 7803	0.99	Yellow	−0.49	4.25 × 10−8
Candidatus Methylopumilus planktonicus	0.96	Yellow	−0.47	4.55 × 10−7
Flavobacteriaceae bacterium	0.91	Yellow	−0.45	2.36 × 10−6
Thiolapillus brandeum	0.94	Yellow	−0.42	3.22 × 10−5
Owenweeksia hongkongensis DSM 17368	0.89	Yellow	−0.39	2.45 × 10−4
Thermogutta terrifontis	0.87	Yellow	−0.38	7.07 × 10−4
Candidatus Pelagibacter sp. IMCC9063	0.92	Yellow	−0.37	1.17 × 10−3
Coraliomargarita akajimensis DSM 45221	0.89	Yellow	−0.36	2.84 × 10−3
Oscillatoria acuminata PCC 6304	0.81	Yellow	−0.34	6.15 × 10−3
Acidothermus cellulolyticus 11Bf	0.94	Yellow	−0.34	9.57 × 10−3

Map ID phylogenetic classification. Value = 1 represents a perfect placement on the tree.

GS = Pearson correlation to ammonia and nitrate + nitrite.

p.GS = p-adjusted value (Bonferroni correction) for correlation to ammonia and nitrate+nitrite.

Represents presence of nitrogenase enzyme EC.1.18.61.

Represents presence of nitrate reductase enzyme EC.1.7.99.4.

Represents presence of nitrate reductase enzyme EC.1.7.1.4.

Significant correlated taxa with ammonia and n>an class="Chemical">nitrate+nitrite result from WGCNA

GS = Pearson correlation to n>an class="Chemical">ammonia and nitrate + nitrite.

p.GS = p-adjusted value (Bonferroni correction) for correlation to n>an class="Chemical">ammonia and nitrate+nitrite.

Represents presence of pan class="Chemical">nitrogenase enzyme EC.1.18.61.

Represents presence of pan class="Chemical">nitrate reductase enzyme n>an class="CellLine">EC.1.7.99.4.

Represents presence of pan class="Chemical">nitrate reductase enzyme n>an class="CellLine">EC.1.7.1.4.

We further explored taxa that correlated with salinity to better understand the impact of tide on the microbial community structure. The red module was positively correlated with salinity (r = 0.47, p = 7.00 × 10−8) (Figure S3). Taxa associated with the red module included Acidimicrobium ferrooxidans DSM 10331 (n>an class="Chemical">CCG) (r = 0.53, p = 8.19 × 10−10), Haliglobus japonicus (CEG) (r = 0.52, p = 2.66 × 10−9), Synechococcus sp. CC9605 (CCG) (r = 0.50, p = 2.28 × 10−8), Candidatus Puniceispirillum marinum IMCC1322 (CCG)) (r = 0.47, p = 4.56 × 10−7), and Prochlorococcus marinus str. MIT 9301 (r = 0.40, p = 1.98 × 10−4) (Table 3).

Table 3

Significant correlated taxa with salinity result from WGCNA

Taxon	Map IDa	Module color	GS.Salinityb	p.GS.Salinityc
Acidimicrobium ferrooxidans DSM 10331	0.94	Red	0.53	8.19 × 10−10
Kordia sp. SMS9	0.91	Red	0.53	1.49 × 10−9
Halioglobus japonicus	0.93	Red	0.52	2.66 × 10−9
Synechococcus sp. CC9605	1.00	Red	0.50	2.28 × 10−8
Synechococcus sp. RCC307	1.00	Red	0.47	4.36 × 10−7
Candidatus Puniceispirillum marinum IMCC1322	0.98	Red	0.47	4.56 × 10−7
Salipiger profundus	0.82	Red	0.46	8.37 × 10−7
Halioglobus pacificus	0.95	Red	0.46	1.02 × 10−6
Roseovarius mucosus	0.96	Red	0.44	7.04 × 10−6
Acidimicrobium ferrooxidans DSM 10331	0.82	Red	0.41	9.08 × 10−5
Candidatus Pelagibacter ubique HTCC1062	1.00	Red	0.40	1.10 × 10−4
Owenweeksia hongkongensis DSM 17368	0.89	Red	0.40	1.64 × 10−4
Prochlorococcus marinus str. MIT 9301	1.00	Red	0.40	1.98 × 10−4
Synechococcus sp. KORDI-100	1.00	Red	0.39	2.07 × 10−4
Sulfurivermis fontis	0.87	Red	0.39	3.93 × 10−4

Map ID phylogenetic classification. Value = 1 represents a perfect placement on the tree.

GS = Pearson correlation to salinity.

p.GS = p-adjusted value (Bonferroni correction) for correlation to salinity.

pan class="Disease">GS = Pearson correlation to salinity.

p.pan class="Disease">GS = p-adjusted value (Bonferroni correction) for correlation to salinity.

We analyzed the Hellinger-transformed enzyme level output from paprica to better understand the enzymatic potential of those CEG and n>an class="Chemical">CCG that were correlated with ammonia. We found 35 enzymes associated with the nitrogen cycle (Figure 6). The enzyme nitrogenase EC 1.18.6.1 had a mean value of 0.23 ± 0.05 for the low disturbed sites, significantly higher than that in the intermediate (0.17 ± 0.06) and high (0.17 ± 0.05) disturbance sites (p = 1.2 × 10−10) (Figure S5). Nitrate reductase EC 1.7.99.4 had a mean value of 0.23 ± 0.05 for low disturbance site, 0.45 ± 0.14 for intermediate disturbance site, and 0.44 ± 0.13 for high disturbance site, and the nitrate reductase value was significantly higher in the high disturbance sites (p = 8.7 × 10−14). The same was observed with nitrate reductase NADH EC 1.7.1.4 with a mean of 0.13 ± 0.05 for low disturbed sites, 0.23 ± 0.14 for intermediate disturbance, and 0.23 ± 0.13 for highly disturbed sites (p = 2 × 10−15) (Figures 6 andS5). The taxa that were associated with nitrogenase were Methylocella silvestris BL2, genus Calothrix, and Synechococcus sp. CC9605. For nitrate reductase members of the Betaproteobacteria, Desulfococcus oleovorans Hxd3, and P. mikurensis NBRC 102666 were found to contribute to enzyme abundance. The taxa that were associated with nitrate reductase NADH were A. cellulolyticus 11B and members of the Rhodobacteraceae (Table 4).

Table 4

Number of enzymes copies for nitrogenase and nitrate reductase enzymes and top 10 associated taxa

Taxon	Nitrogenase EC.1.18.6.1	Nitrate reductase EC.1.7.99.4	Nitrate reductase NADH EC.1.7.1.4
Methylocella silvestris BL2	16,984	4,246	0
Genus Calothrix	15,186	3,796	0
Synechococcus sp. CC9605	38,937	0	0
Family Rhodobacteraceae	0	0	1,273
Oscillatoria nigroviridis PCC 7112	0	5,418	5,418
Acidothermus cellulolyticus 11B	0	0	23,205
Phycisphaera mikurensis NBRC 102666	0	8,288	0
Rhodopirellula baltica SH 1	0	10,956	21,912
Desulfococcus oleovorans Hxd3	0	4,773	0
Class Betaproteobacteria	0	17,102	0

Number of enzymes copies for pan class="Chemical">nitrogenase and n>an class="Chemical">nitrate reductase enzymes and top 10 associated taxa

Discussion

Mangrove forests are experiencing a high degree of perturbation through nutrient enrichment, pollution, and deforestation. Shrimp aquaculture effluent in particular is associated with the input of excess nutrients to mangrove forests. In this study we found that shrimp aquaculture effluent is associated with changes in microbial community structure with likely consequences for biogeochemical cycles and mangrove forest health. Previous work suggests that for intensive shrimp farming, 2.22 km2 of mangrove forest is required to remove effluent from one pond of 0.01 km2, whereas 0.20 km2 is required for less-intensive farming from one pond of 0.01 km2 (Robertson and Phillips, 1995). As of 2014 in the Muisne region there were 20.47 km2 of shrimpn> farms and 12.06 km2 of mangrove forests, indicative of an intensive farming system. Cayapn>as-Mataje had 11.04 km2 of shrimpn> aquaculture farms and 302.05 km2 of mangrove forest, suggesting less intensive farming (Figure 1) (Hamilton, 2020). As the areal extent of shrimpn> aquaculture increases so does the volume of the effluent, elevating the flux of n>an class="Chemical">ammonia and nitrate to the surrounding ecosystem. Based on our observations we found that microbial communities in mangrove forests are significantly altered by this perturbation.

The bacterial communities in our mangrove systems were characterized by members of the Pelagibacteraceae, Flavobacteriaceae, Rhodobacteraceae, n>an class="Species">Actinobacteria, and cyanobacteria (Figure 4, Figure S1). The archaeal community was dominated by members of the Thermoplasmata, Thaumarchaeota, and Methanococcales (Figure 4, Figure S1). This was in accordance with other studies that have identified Rhodobacteraceae, SAR86 clade, Actinobacteria, and Flavobacteriaceae, and Thaumarchaeota as the most abundant taxonomic groups (Dhal et al., 2020). Rhodobacteraceae has been found to be dominant in mangrove-dominated estuaries, and members of this family are associated with marine phytoplankton blooms where they play a role in transformations of derived phytoplankton organic matter (Ghosh et al., 2010; Simon et al., 2017). The presence of Actinobacteria has been documented previously in mangrove ecosystems (Azman et al., 2015; Gong et al., 2019), and it has been suggested that they could play a role in carbon cycling by decomposing the plant biomass including refractory lignins (Scott et al., 2010). Thaumarchaeota are the most abundant archaea in the surface ocean (Santoro et al., 2015), and Thermoplasmata have been found in mangrove ecosystems (Zhang et al., 2019). Both these groups play an important role in the nitrogen cycle by carrying out the oxidation of ammonia in nitrification (Santoro et al., 2015; Zhang et al., 2019).

Both bacterial and archaeal communities were less diverse at our more disturbed sites. This pattern extended to predicted metabolic diversity (Figure 3). We hypothesize that this reduction in diversity could cause reductions in ecosystem functions. This has been observed in previous mangrove forest studies, for example, where lower microbial diversity was associated with a reduction in microbial productivity in sites with high levels of deforestation, sewage, and fishing activities (Carugati et al., 2018).

Our results showed differences in biogeochemical parameters between sites at varying levels of disturbance (Figure 2). In particular, nitrogen was a driver of the microbial community structure leading to segregation into three clusters of disturbance in the CA analysis based on our ANOSIM test and significantly correlated with n>an class="Chemical">ammonia concentrations (Figure 5, Figure S2). We note that the variance explained by the first and second dimensions in our ordination analyses is relatively low (30.6% and 19.9% for the bacterial and archaeal communities, respectively). We attribute this to the complexity associated with the mangrove ecosystem and the large number of physical, chemical, and biological factors that could impact changes in the microbial community.

We found a strong connection between N:P ratio and genome size among planktonic bacteria across study sites (Figure 2). Generally, smaller predicted genomes and lower 16S rRNA gene copy number was associated with higher N:P ratios, whereas larger predicted genomes and higher 16S rRNA gene copy number was associated with lower N:P ratios. The differences in genome sizes between communities associated with different levels of disturbance suggest differing ecological strategies. Studies suggest that generalists possess larger genomes in contrast to the smaller genomes in more specialized microbes (Sriswasdi et al., 2017; Willis and Woodhouse, 2020). This falls from the generalist requirement for a larger gene repertoire to boost activity in multiple environmental conditions and to cope with different stressors associated with a broad physicochemical niche (such as low levels of nitrogen and tidal fluctuations in mangrove-dominated estuaries). The low disturbance sites showed a higher metabolic diversity and larger genomes, which we interpn>ret as a more generalist microbial community. Taxa with larger genomes included Planctomycetes such as n>an class="Species">Singulisphaera acidiphila. This taxon has been found in other wetland ecosystems (Kulichevskaya et al., 2008; Dedysh and Ivanova, 2019), and it has been shown to play an important role in degradation of plant-derived polymers such as pectin and xylan (Dedysh and Ivanova, 2019). The S. acidiphila genome encodes several dozen proteins that do not belong to any of the currently carbohydrate-active enzymes, but the enzymes display a distant relationship to glycosyltransferases and carbohydrate esterases, suggesting that this taxon has a diverse glycolytic and carbohydrate metabolic potential (Dedysh and Ivanova, 2019). Other taxa included Sandaracinus amylolyticus. This taxon has been found in association with plant residues (Mohr et al., 2012), in coral ecosystems (Rubio-Portillo et al., 2016), and it is known to survive in poor nutrient conditions by developing desiccation-resistant spores (Mohr et al., 2012).

We also observed larger genomes in cyanobacteria including members of the genus pan class="Species">Calothrix, n>an class="Disease">genus Oscillatoria, and M. producens PAL-8-15-08-1. Cyanobacteria are known to have large genomes with low coding density and a high level of gene duplication; it has been proposed that the large non-protein-coding sequences contribute to the genome expansion and metabolic flexibility observed in diazotrophs (nitrogen fixers) that are associated with nitrogen-limited environments (Sargent et al., 2016). The high diversity of cyanobacteria observed in mangrove ecosystems suggests that they play a key role in the ecosystem. Relevant functions associated with cyanobacteria include nitrogen and carbon fixation and the production of herbivory-defense molecules and plant growth-promoting substances (Alvarenga et al., 2015).

In the disturbed sites the parasite C. pan class="Species">Dependentiae accounted for much of the decrease in genome size. Studies have found that C. n>an class="Species">Dependentiae infects a wide range of protists, including heterotrophs and phytoplankton (Deeg et al., 2019). Other studies have shown that C. Dependentiae is associated with free-living ameba, suggesting that it could be an endosymbiont (Delafont et al., 2015). C. Dependentiae has very limited metabolic capability, lacks complete biosynthetic pathways for various essential cellular building blocks, and has protein motifs to facilitate eukaryotic host interactions (Yeoh et al., 2016; Deeg et al., 2019). C. Nasuia deltocephalinicola was also identified as having a small genome. C. Nasuia deltocephalinicola is an obligate symbiont of plant phloem-feeding pest insects, and its main role is to provide essential amino acids that the host can neither synthesize nor obtain in sufficient quantities from a plant diet (Bennett and Moran, 2013).

The increase in the concentration of nitrogen spn>ecies associated with aquaculture effluent could further select for spn>ecialist microbes with reduced metabolic potential and lower diversity in n>an class="Chemical">nitrogen-processing enzymes. Our nitrogen isotope values are consistent with this, showing reduced variability for the highly disturbed sites (Table 1, Figure 2). Previous work in mangrove systems (Bernardino et al., 2018) associated reduced isotopic variability with a loss of trophic diversity. Higher variability in stable carbon isotopes has also been observed in salt marshes due to contribution dominated by allochthonous material derived from the phytoplankton community (Boschker et al., 1999). The larger variation in the isotopic signal observed in the low and intermediate disturbance sites suggests that these pristine systems contain a more diverse trophic food web as result of a wide range of metabolic and fixation pathways, and environmental conditions in the mangrove-estuarine ecosystem (Boschker and Middelburg, 2002).

The low N:P ratios we observed in the low disturbance sites suggest that the system is N limited. Pristine mangrove forests tend to be N limited, although nutrients are not uniformly distributed within the mangrove ecosystem and they can switch from N to P limitation. It has been shown that mangrove trees within fringe and tidally exposed zones tend to be N limited (Feller et al., 2003). One way mangrove trees cope with N limitation is through associations with diazotrophs that play a crucial role in N cycling within the mangrove forest (Holguin et al., 1992). Here we showed that the biological pan class="Chemical">nitrogen fixation signal, confirmed by the N∗ value (the linear combination of n>an class="Chemical">nitrate and phosphate that eliminates the effect of nitrification; thus, the remaining variability can be explained by nitrogen fixation and denitrification) (Gruber and Sarmiento, 1997) and nitrogenase EC.1.18.6.1 abundance, were higher at low disturbance sites in contrast to high disturbance sites (Figures 2, 6, and S5). The microbial denitrification signal was further confirmed by negative N∗ values in the highest disturbance sites (Figure 2) (Gruber and Sarmiento, 1997). Because excess nitrate is being introduced into the system via aquaculture effluent, we expect denitrification rates to be high. Conversely, the lowest disturbance sites have a slight positive N∗ consistent with our identification of putative diazotrophs such as genus Calothrix, genus Oscillatoria, and taxa of the order Rhizobiales such as M. silvestris (Essien et al., 2008; Liu et al., 2019).

GBT degradation I was one of the major pathways contributing to the differences observed between low and high disturbed sites (Figures 6 and S3). n>an class="Chemical">GBT is an important source of nitrogen in oligotrophic systems, acts as an organic osmolyte, and plays an important role in phytoplankton-bacteria interactions (Becker et al., 2019; Jones et al., 2019; Zecher et al., 2020). The intertidal coastal mangrove ecosystem experiences daily fluctuations in a range of environmental conditions, including water levels and salinity. Organisms living in this dynamic environment cope with changing environmental conditions by synthesizing a range of organic and inorganic osmolytes including GBT. The results from WGCNA showed that Pelagibacteraceae taxa correlated with salinity (Table 3) and primary contributors of the GBT degradation I pathway. This suggests that osmolyte production is an important adaptation to salinity intrusions from oceanic waters into the mangrove environment, and GBT could be an additional pool of organic N within this system.

Shrimp aquaculture impacts the water quality in adjacent mangrove forests by creating eutrophic conditions that can lead to anoxia. Eutrophic conditions were evident through high levels of nutrients and chlorophyll a (Figure 2, Table 1). Although we did not measure oxygen concentrations, we observed taxa indicative of hypoxic or anoxic conditions. These included purple sulfur bacteria (PSB), such as family Chromaticeae, and sulfur-oxidizing bacteria (SOB), such as genus Sulfurivermis (Figure 4, Table 2). PSB use sulfide, elemental sulfur, and thiosulfate as electron donors in anoxygenic photosynthesis and have been shown to play an important role in regime shifts from oxygenated to anoxic conditions (Diao et al., 2018). PSB flourish in micro-aerobic conditions oxidizing sulfide into sulfate (Diao et al., 2018). As the oxygen influx is reduced below a critical threshold, sulfate-reducing bacteria (SRB) and PSB can take over and outcompete the SOB. This suggests a more anoxic regime in the high disturbance site, allowing for PSB groups and SRB to become more abundant.

Based on our WGCNA analysis we also found nitrate-reducing bacteria (NRB)—indicative of reduced n>an class="Chemical">oxygen availability—that strongly correlated with the level of ammonia, nitrate, and nitrite. Putative NRB taxa included P. mikurensis and S. denitrificans (Table 2). In addition, we also saw a microbial signature associated with dissimilatory nitrate reduction to ammonium (DNRA) with the presence of genus Acidothermus, and anaerobic ammonium oxidation (annamox) with the presence of Planctomycetes (Thermogutta terrifontis) (Table 2); the presence of the genes involved in these pathways (denitrification, annamox, DNRA) were inferred by paprica, although further work is needed to confirm the presence and activity of these enzymes. Overall, as nitrate and ammonia inputs increased with aquaculture effluent the relative abundance of NRB increased.

We identified specific microbes that can be used as sensitive indicators of aquaculture impacts. These included P. balearica (Figure 4), which has been associated with other contaminated wetland systems, suggesting that this taxon could be a potential bio-indicator of a disturbed mangrove ecosystem (Salvà-Serra et al., 2017). Similar studies have also identified aquaculture effluent as a source of pathogens to the coastal ecosystem (Garren et al., 2009). In the disturbed site we saw the presence of members of the genus Arcobacter (Figure 4). These bacteria have been identified in coral systems expn>osed to aquaculture effluents, and have been associated with feces (n>an class="Species">human, porcine, and bovine) and with sewage-contaminated waters (Garren et al., 2009). PSB taxa such as family Chromaticeae have also been shown to be potential bio-indicators for anthropogenic contamination associated with other agriculture effluent systems (Mohd-Nor et al., 2018). P. salivibrio in the order Micrococcales was elevated at the disturbed sites. This taxon has been isolated from high-salinity systems and aquaculture farms (Jang et al., 2013); high salinity levels have been associated with shrimp aquaculture effluent due to high evaporation in the ponds (Barraza-Guardado et al., 2013). Previous studies have shown that taxa in the Micrococcales order are part of the core microbiome signal in shrimp ponds (Chen et al., 2017). Thus, P. salivibrio is a sensible indicator of shrimp aquaculture effluent. Further work is needed to establish robust spatiotemporal baselines of microbial indicators of aquaculture to effectively monitor biogeochemical changes and health of the mangrove forests.

Aquaculture could impact the health of mangrove ecosystems involving the direct loss of mangrove forests, effluent associated with high levels of nutrients, and the development of anoxic and pan class="Chemical">sulfuric n>an class="Chemical">water conditions (Robertson and Phillips, 1995). Aquaculture effluent released into mangrove forests may be sequestered and processed by bacteria. However, processing efficiency could change with increasing input. High organic loadings, for example, may shift the balance from aerobic to anaerobic systems (Lønborg et al., 2020). Anaerobic systems are less efficient in nutrient cycling. The signals of SOB, SRB, denitrifiers, and potential pathogenic taxa associated with the perturbed site suggest that aquaculture effluent is playing a role in shifting the microbial community to a more pathogenic and less nutrient efficient community that could impact the health of the mangrove forest.

Conclusion

In this study, we showed the impacts of aquaculture effluent on the microbial community structure in mangrove forests and identified microbial signals associated with NRB, pan class="Species">PSB, and n>an class="Chemical">SRB taxa that could have impacts in nutrient cycling. The high level of nutrients in the perturbed sites were associated with changes in microbial community structure that could impact ecosystem functions. In the low disturbance sites, we saw that the presence of Calothrix species and nitrogen fixers could be important in increasing nitrogen inventories via nitrogen fixation. Denitrification reduces excess inorganic nitrogen concentration, and in the highly disturbed sites we saw the presence of NRB-associated microbes. Nutrient cycling in mangrove habitats is a balance between nutrient inputs, availability, and internal cycling, and the changes in microbial community structure we see in disturbed sites could be indicators of biogeochemical changes. The results of the study highlight the sensitivity of the mangrove-estuarine microbial community to aquaculture effluent, and the impacts of land use changes could be amplified by climate change such as changing precipitation patterns, heat, and rising sea level with severe consequences for the ecosystem.

Limitations of the study

Our analysis was based on comparison between sites of low, intermediate, and pan class="Disease">high disturbance in two mangrove systems in coastal Ecuador. Although n>an class="Chemical">ammonia concentration is a good proxy for disturbance from shrimp aquaculture effluent, quantification of land use changes, and the hydrological connections between aquaculture facilities and our sampling sites was beyond the scope of the current work. We considered salinity, macronutrient concentrations, and isotopes in our analysis, but anticipate that other variables not considered here are contributing to differences in microbial community structure. These include physical processes such as tides and hydrology. The complexity of these environments is evident in our CCA and CA analyses, which capture a relatively small amount of variability in the first two dimensions (see discussion section). Other limitations of note are typical of microbial community structure analyses. These include primer bias and a dependence on relative rather than absolute abundance.

Resource availability

Lead contact

Further information and requests for resources should be directed to and will be fulfilled by the lead contact, pan class="Chemical">Natalia Erazo (n>an class="Chemical">nerazo@ucsd.edu).

Materials availability

This study did not generate new unique reagents.

Data and code availability

The data that support the findinpan class="Disease">gs of this study and sequences were submitted to the NCBI sequence read archive (SRA) under BioProject ID: PRJNA633714. Code for analysis is available on github repn>ository: httpn>s://github.com/galud27.

Methods

All methods can be found in the accompanying Transparent methods supplemental file.