The role of seagrass vegetation and local environmental conditions in shaping benthic bacterial and macroinvertebrate communities in a tropical coastal lagoon

We investigated the influence of seagrass canopies on the benthic biodiversity of bacteria and macroinvertebrates in a Red Sea tropical lagoon. Changes in abundance, number of taxa and assemblage structure were analyzed in response to seagrass densities (low, SLD; high, SHD; seagrasses with algae, SA), and compared with unvegetated sediments. Biological and environmental variables were examined in these four habitats (hereafter called treatments), both in the underlaying sediments and overlaying waters, at three randomly picked locations in March 2017. Differences between treatments were more apparent in the benthic habitat than in the overlaying waters. The presence of vegetation (more than its cover) and changes in sedimentary features (grain size and metals) at local scales influenced the observed biological patterns, particularly for macroinvertebrates. Of note, the highest percentage of exclusive macroinvertebrate taxa (18% of the gamma diversity) was observed in the SHD treatment peaking in the SA for bacteria. Benthic macroinvertebrates and bacteria shared a generally low number of taxa across treatments and locations; approximately, 25% of the gamma diversity was shared among all treatments and locations for macrofauna, dropping to 11% for bacteria. Given the low overlap in the species distribution across the lagoon, sustaining the connectivity among heterogeneous soft sediment habitats appears to be essential for maintaining regional biodiversity. This study addresses a current scientific gap related to the relative contributions of vegetated and unvegetated habitats to biodiversity in tropical regions.


Material and methods
Study site. The present study was conducted at the Al Qadimah lagoon (22° 22′ 39.3″ N 39°, 07′ 47.2″ E; approximate area of 13.9 km 2 ), located in the central part of the western coast of the Saudi Arabian Red Sea (18°-25°N), where the highest seagrass diversity is observed 55 . The lagoon is dominated by shallow subtidal sands (average depth 2.19 m) with a maximum depth of approximately 10 m in the channel that connects the lagoon to the open sea. This lagoon does not experience direct anthropogenic disturbances typical of other coastal lagoons (e.g. sewage discharges, fisheries practices, habitat destruction for coastal development). Riverine inputs are also non-existent. However, it is located between two urbanized areas developed over the last two decades (KAUST, King Abdullah University of Science and Technology, 7,000 inhabitants; KAEC, King Abdullah Economic City, < 5,000 inhabitants but it is expected to reach 50,000 over the next years). Despite the lack of direct anthropogenic influences, recent studies 56 quantified moderate levels of metals in the seagrass meadows.
The bottom of the lagoon, particularly in the inner sheltered areas, is covered by seagrass meadows often consisting of at least two species. The areas are mostly dominated by Cymodocea rotundata, Cymodocea serrulata and Enhalus acoroides, sometimes with algae interspersed. In general, multispecies meadows are dominant. Therefore, the selection of seagrass areas was not limited by the species composition. Previous studies in tropical areas did not find an effect of the vegetation type on densities, species richness or even composition of epifaunal and infaunal communities 19 .
Sampling design and strategy. The  www.nature.com/scientificreports/ microbial and macrobenthic benthic communities to habitat fragmentation. All sampling sites were surveyed within a depth range of 1.6-2.8 m.
In each location, four different treatments were established. These are representative of what is commonly observed in the lagoon but also of different alternative states associated with seagrass fragmentation. The treatments were established based on the seagrass density: unvegetated (UnV) with estimated percentage cover of 0-5%; seagrass low density (SLD) with estimated percentage cover of 10-25%; high density (SHD) with estimated percentage cover of 50-75%; and seagrass and algae (SA) mixed meadows. In the field, triplicate quadrats (50 × 50 cm) were placed at random within each treatment in order to obtain the percentage of seagrass and algae cover (%) and to measure shoot density (m 2 ).
A combination of sediment and water samples was collected at each location and treatment, following a block sampling design with two factors: Location (fixed, three levels: A, B and C) and Treatment (fixed, four levels: UnV, SA, SLD, SHD).
Water samples. At each sampling site, CTD casts were carried out for conductivity, temperature and depth (CTD) with a multiparameter probe (OCEAN SEVEN 316 Plus and 305 Plus CTD). Using Niskin bottles (5 L), water samples were collected at each location. Divers, using 1 L bottles, also collected samples a few centimeters above the bottom while avoiding sediment disturbance (overlaying water samples). For nutrients, triplicate 50 mL of water (both in the water column and overlaying water) were filtered on the boat (Isopore membrane filters, 0.2 μm GTTP) and then stored in the dark and cool until they were frozen at − 50 °C in the lab. For chlorophyll a (Chl a), water samples were later filtered in the lab through 47 mm GF/F filters, wrapped in aluminum foil and preserved at − 80 °C until extraction of the pigments. For dissolved organic carbon (DOC) and total dissolved nitrogen (TDN) concentrations, water was filtered through pre-combusted GF/F filters, acidified with phosphoric acid 87% down to pH 1-2 and stored at 4 °C until analysis. For fluorescent dissolved organic matter (FDOM) analysis, water was filtered through pre-combusted GF/F filters, kept at 4 °C and analyzed within 6 h of collection.
Sediment samples. In each location and treatment, three replicate core samples were taken (12 cm internal diameter, 20 cm depth) for the study of macroinvertebrates. Each sample was composed of three cores (total area of 0.03 m 2 ); 36 replicates in total were collected. In order to avoid edge effects typical of seagrass meadows 57 , samples were whenever possible taken from the innermost area of each treatment. In the field, sediment was washed over a 1 mm sieve, and samples were preserved in ethanol (96%). For benthic bacterial analyses, 5 g of www.nature.com/scientificreports/ sediment was taken in duplicates. These 5 g samples were placed in RNA later, and frozen until further laboratory sampling. Additional sediment cores were collected to quantify associated environmental variables, specifically, total organic carbon (TOC), grain-particle size, and nutrients in pore-water ( ) were quantified using a Continuous Flow Analyzer (SEAL AutoAnalyser 3 with XY2/3 Sampler). To quantify nutrients in pore water samples, sediment was centrifuged for 30 min at 3,000 rpm to separate the pore water from the solid fraction.
Chl a quantification was undertaken using the Acetone (10 mL at 90%) extraction method. Briefly, samples were left in acetone for 24 h in low light conditions to minimize degradation. A Turner Trilogy fluorometer (Turner Designs) was used to quantify the Chl a content using an acidic module. The degradation of the Chl a to phaeophytin was accomplished by acidifying the sample with 60 µL of 0.1 N HCl. DOC and TDN measurements from the water column were performed by high temperature catalytic oxidation (HTCO) using a Shimadzu TOC-L. To monitor the accuracy of the estimated carbon, we used reference material of deep-sea carbon (42-45 μmol C L −1 and 31-33 μmol N L −1 ) and low carbon water (1-2 μmol C L −1 ). Dissolved organic nitrogen (DON) concentrations were calculated after subtracting the dissolved inorganic nitrogen (DIN) to the TDN (DON = TDN − DIN), where DIN (μmol N L −1 ) = [NO 3 − ] + [NO 2 − ], and DOM C:N ratios were calculated by dividing DOC by DON. FDOM analysis were performed using a HORIBA Jobin Yvon AquaLog spectrofluorometer with a 1 cm path length quartz cuvette. The three-dimensional fluorescence excitation emission matrices (EEMs) were recorded following the same steps as in Calleja et al. 58 . The EEMs were modeled together with other FDOM samples (> 400) collected at the same lagoon within a time series study conducted between 2015 and 2017 (unpublished data). A four-component model was validated using split-half validation and random initialization 59 and the following peaks were identified: peak C1 at Ex/Em 252/ 397 nm (humic-like peak A 60 ), peak C2 at Ex/Em 252/467 nm (a combination of humic-like peaks M and C 60 peak C3 at Ex/Em 303/337 (protein-like peak T 60 and attributed to Tryptophane), and peak C4 at Ex/Em 270/312 nm (protein-like peak B 60 and attributed to Tyrosine). The maximum fluorescence (Fmax) of each component is reported in Raman units (RU). Autotrophic and heterotrophic picoplankton abundance in the water column samples was measured by flow cytometry using standard protocols described in detail in Gasol and Morán 61 and Silva et al. 62 . Briefly, 1.8 mL samples were fixed with 1% paraformaldehyde + 0.05% glutaraldehyde, deep frozen in liquid nitrogen and stored at − 80 °C until analysis. After samples were thawed, 0.4-0.6 mL aliquots were run through a FACSCanto II flow cytometer (BD Sciences) to determine the abundance of cyanobacteria and heterotrophic prokaryotes. Only Synechococcus were detected among the latter since Prochlorococcus avoid shallow, enclosed environments such as the lagoon 62 . The widespread groups of low and high nucleic acid content of heterotrophic prokaryotes were distinguished based on the relative green fluorescence after staining the samples with SybrGreen (Molecular Probes). Flow rates were determined empirically to obtain cell abundances.
Sediment samples. These samples were dried at 60 °C until constant weight prior to their homogenization in an automatic mill grinder (model PM200 from Retsch) where the samples were milled for 4 min at 420 rpm. Afterward, the fine powder was collected and prepared for analyses. Total carbon and nitrogen were determined for the powder samples. Inorganic carbon and nitrogen fractions were also calculated after loss on ignition (450 °C for 3 h). Analysis was performed in a CHNS-O "Thermo Scientific" organic elemental analyzer (model Flash 2000), with soil standard 2G NCS as reference material. Particulate and total sediment organic carbon (POC and TOC, respectively) and nitrogen (PON and TON, respectively) were calculated by the difference between the total and inorganic fractions. For trace element analysis (Al, Ba, Cd, Co, Cr, Cu, Fe, Pb, Mg, Mn, Ni, Se, V, and Zn), samples were dried at 40 °C and 0.2 g was analyzed by flame atomic absorption spectrophotometer AAS (Perkin Elmer, Model 2380A Spectrophotometer) following the EPA 2007 methods, as detailed in 56 . Particle size analysis was conducted by wet separation of the silt and clay fractions from the sandy fractions. Replicates of sediment samples were sieved through a 63 μm mesh (silt and clay fraction), and the retained sediment fraction was dried at 80 °C for 24-48 h. The dried sample material was sieved by using a column of sieves. Fractions were classified as: gravel (> 2 mm), very coarse sand (> 1 mm and < 2 mm), coarse sand (> 500 μm and < 1 mm), medium sand (> 250 μm and < 500 μm), fine sand (> 125 μm and < 250 μm), very fine sand (> 63 μm and < 125 μm), and siltclay (fines, < 63 μm). Seagrass material (roots and leaves) were separated and biomass measured as dried weight.
Benthic macrofauna. Samples were sorted and organisms identified to the highest taxonomic separation level (50% of OTUs at the species level) and counted. The lack of taxonomic knowledge for most of the soft-sediment invertebrates in the region, together with the preservation of samples in ethanol precluded the identification of organisms to the same taxonomic level.
Benthic bacteria. Samples were processed using the PowerMax Soil DNA isolation kit (Qiagen) as per the manufacturer's instructions with the exception of step 3 where the bead beating was replaced by the addition of 0.4 mg mL −1 of Proteinase K and an overnight incubation in a shaking incubator at 56°C 63 . The bacterial component was amplified with PCR by targeting the v3 and v4 regions of the 16S rRNA gene using the primers (341F: 5′-CCT ACG GGNGGC WGC AG-3′ and 805R: 5′-GAC TAC HVGGG TAT CTA ATC C-3′ 64

Data visualization and statistical analyses. Environmental data. Principal Component Analysis
(PCA) was used to create an ordination plot of the samples based on the environmental variables that were measured. Prior to the analysis, the correlation between variables was assessed (Pearson R). Highly correlated variables (> 0.7) were excluded, which primarily related to the concentration of metals and grain size. When correlation between two variables was higher than 0.7 only one of the variables was chosen based on the following criteria: (1) the one that was less correlated to other variables; (2) the explanatory importance/ potential of the variable for the specific case; (3) the extension of the gradient (i.e., variables with a larger gradient (range) were preferred to small gradients). Correlations between 0.5 and 0.7 were also subjected to the same procedure to further reduce the number of correlated variables in the analysis. Finally, 10 variables were used in the PCA analysis (fines, TOC, DOC, Chl a, TDN, NO 2 − pw , PO 4 3− pw , NO 2 − and NO 3 − ). Information regarding the multiple correlations between variables is provided in Supplementary Material (S1). Variables that did not present a normal distribution were transformed (log e or Box-Cox transformation). Environmental data was normalized by scaling each variable to zero mean and unit variance before the analyses using function "decostand" from R package "vegan" 65 .
Macrofauna. Several univariate metrics were calculated including density (ind. m −2 ), and the average number of taxa (i.e., number of taxa × 0.03 m −2 ). Two-way ANOVAS were conducted for the factors Treatment and Location after checking for normality and homogeneity of the data. Normality was checked using Shapiro-Wilk normality test and by using quantile-quantile plots. Homogeneity of variance was tested using Levene's test. Whenever normality or homogeneity of variance was not verified the data was transformed (log e ) to meet ANOVA assumptions. Post-hoc tests (Tukey Honest significance test) for ANOVA significant factors were used to ascertain for differences between groups of samples. Community structure was investigated by ordination techniques (Principal Coordinate Analysis), based on the Bray-Curtis dissimilarity coefficient applied to untransformed abundance data as the dataset is characterized by low abundance values and there is no need to reduce the weight of the high abundance taxa. Permutational multivariate analysis of variance (PERMANOVA), for the two-factors in analysis (i.e. Treatment and Location) was applied to assess for significant differences in the associated fauna. When significant differences were detected, pair-wise tests were applied. Distance-based redundancy analysis (dbRDA) was also applied in order to understand the relationship between each environmental variable and benthic community structure (given by the direction and length of vectors for each variable). A set of 16 variables was selected from the original 37 to reduce collinearity problems using the criteria used for the PCA but only removing the variables that were highly correlated in each domain. The set comprised grain size (fines and coarse sand), dissolved inorganic nutrients ( , trace elements (Cu, Pb, As), productivity proxies (Chl a) and other nutrients (TOC, DOC, TDN) and seagrass density (density of seagrass shoots). Due to differences in the scales and units of environmental variables, environmental data was normalized by scaling each variable to zero mean and unit variance before the analyses. The model used for the dbRDA was manually built to maximize the explanatory power of the variables, starting from a model with no explanatory variables and adding or removing variables according to their capability to further explain the variability of the community data (i.e., evaluating all single terms that can be added or dropped from the constrained ordination model using "add1.cca" and "drop1.cca" functions; R package "vegan" 65 ). The selection of the variables was based on statistics (F and AIC for the single terms) but also considering the ecological processes (several processes instead of a set of variables related to one process). The significance of each variable retained in the model was assessed using marginal and sequential tests.
Benthic bacteria. The "dada2" package 71 was used to process the raw sequencing reads within R 66 . The raw reads were truncated based on quality and then filtered with a maxEE of 2 for forward reads and 6 for reverse. Error rates were calculated for the filtered reads and subsequently paired ends were merged with a maximum mismatch of 1 and a minimum overlap of 10. The function removeBimeraDenovo was used to remove chimeras and taxonomy was assigned to each representative sequence against the Silva v128 database 72 within the "dada2" package resulting in a table of operational taxonomic units (OTUs) against sample, using the "phyloseq" package 73 . Any sequences assigned as eukaryotes and Archaea, as well as chloroplasts and mitochondria were subsequently removed bioinformatically. This removed 6.2% of OTUs but only 0.48% of reads. Samples were rarefied to an even depth (10,000 reads) for comparison. Samples not meeting this criterion were removed from analysis. After the creation of the OTU table, the same statistical analytical routines were undertaken as described for the macrofauna. Multivariate analysis was based on the log-transformed data to down weight the weighting of taxa with a high number of sequences. A two-way ANOVA was conducted to test for differences regarding Treatment and Location (see previous section "Macrofauna"). Due to the unbalanced design a Type III ANOVA was used on transformed values (log e ). www.nature.com/scientificreports/

Responses of taxa to habitat drivers.
A forward step-wise regression model using the Bayesian Information Criterion (BIC) exit criterion was used to determine the environmental parameters that explained benthic taxa differences between treatments. The taxa identified as important contributors for the biological patterns based on Similarity Percentage (SIMPER) analysis (Primer v7) were analyzed using step-wise regression. Regression models that account for a large number of potentially relevant variables may lead to over-fitting and poor prediction performance. Hence, environmental variables excluding those that were highly correlated were carried forward in the regression. Linear regression models were then applied to identify the magnitude and direction effect of the important environmental variables identified from the stepwise regression that significantly contributed to benthic abundance data.  www.nature.com/scientificreports/ vegetated areas presented sediment with higher content of phosphate. Chlorophyll a and TOC (and variables correlated to these) separate vegetated samples from Location C from those locations A and B. Within locations A and B, nutrients in the porewater (silicates, phosphates and nitrates), TOC and Chl a presented higher values in vegetated compared with the unvegetated treatments.

DOC, FDOM and autotrophic and heterotrophic bacterioplankton in the water column.
DOC concentration and C:N ratios averaged (± SD) 108 ± 4 µmol C L −1 and 13.7 ± 1.4, respectively, and did not show significant differences between locations or treatments. Although changes were small, C:N ratios significantly increased with NO 3− concentrations (p < 0.05), indicating that NO 3− depleted samples (due to less production and/or higher consumption) presented DOM richer in nitrogen and vice-versa, i.e., samples with higher NO 3− in the water presented DOM richer in carbon. The fluorescent characteristics of DOM identified two humic-like components (C1 and C2), and two protein-like components (C3 and C4). C1 and C2 were always more abundant (range 0.032-0.045 R.U.) than C3 and C4 (range 0.014-0.021 R.U.). The humic fraction ranged between 66 and 69% and did not show significant differences between locations or treatments. However, when correlating the FDOM data to environmental parameters we observed that all components except C4 significantly increased with salinity, despite the limited variability (41.7-42.2). The correlation was particularly strong for the humic components (C1 p < 0.001, r = 0.89; C2: p < 0.001, r = 0.85) and weaker but still significant for the Tryptophanelike peak (C3: p < 0.05 r = 0.70). Humic components also showed a positive and significant relationship with Chl a values (C1: p < 0.001, r = 0.77, C2: p < 0.001, r = 0.80) although no relationship was found with the protein-like components. Autotrophic bacterioplankton was represented by Synechococcus, which ranged from 0.55 to 9.33 × 10 4 cells mL −1 . Heterotrophic planktonic bacteria were one order of magnitude more abundant (0.58 to 3.22 × 10 5 cells mL −1 ) than autotrophs in the pelagic environment. A trend of highest abundance of bacterioplankton in the SLD and lowest in the SHD, with similar intermediate values in unvegetated and mixed meadows was found. Indeed, the abundance of autotrophic and heterotrophic bacteria were strongly correlated (r = 0.93, p < 0.001, n = 10). While differences between the four treatments were not significant, abundances of both autotrophic and heterotrophic planktonic bacteria at location B were significantly lower (ANOVA, p < 0.05) than at the other two locations.
Overall, benthic macrofaunal density (log e transformed data) and number of taxa (log e transformed data) did not vary significantly among treatments (ANOVA; Density: F = 1.558, p = 0.226; Number of OTUs: F = 0.937, p = 0.440) but there was a significant effect of Location (ANOVA; Density: F = 7.201, p = 0.004; Number of OTUs: F = 14.491, p < 0.001). Both metrics were significantly higher in location C, compared to locations A and B (Fig. 3A, B). Regarding the number of macrofauna taxa (log e transformed data), significant differences were also observed between locations A and B (Tukey post-hoc test, p = 0.044). For bacteria, the number of OTUs was generally higher in locations A and B than in location C (Fig. 3C). The observed number of taxa showed, however, a significant interaction between Location and Treatment (F = 2.546, p = 0.043). Post-hoc tests showed significant differences between locations C and both A and B in all treatments, except SHD (no significant differences). In locations A and B no significant differences were detected between treatments whereas in location C, values in UnV and SHD were significantly different (Tukey post-hoc test, p = 0.012) (results in Supplementary Material, Table S1).
A total of 18 macrofaunal taxa representing 23% of gamma diversity (i.e. all taxa identified) were shared among treatments across the three locations. Among the shared taxa, the most abundant was the sipunculid Phascolion strombus strombus (20.2% of total abundance). Other main contributors to total abundance were the opportunistic polychaetes Capitellidae und. (7.8%), the bivalves Barbatia foliata (7.8%) and Cardiolucina semperiana (6.1%), as well as polychaetes from the family Eunicidae (6%). SHD showed the highest percentage of exclusive taxa (18% of gamma-diversity), twice as much as the unvegetated treatment and 3.5 times more than the mixed meadow (Fig. 4). Regarding the factor Location, 20 taxa were present at the three locations, with location C showing the highest percentage of exclusive taxa (i.e., 24; 31% of gamma diversity) (Fig. 4a, b).
A total of 1,454 bacteria OTUs were shared amongst all treatments, representing 10.5% of gamma diversity (yet accounting for 78.6% of the relative abundance, Fig. 4c, d). The majority of the abundance in these shared OTUs was classified as belonging to the families Desulfobacteraceae and Syntrophobacteraceae. The mixed meadow (SA) treatment showed the highest number of exclusive taxa (19.3%). Location C showed the lowest percentage of exclusive taxa (8%) with the highest percentage being observed in location B (35%). A total of 62.1% of gamma diversity was observed in location B. This was higher than location A (53.7% of gamma diversity) and location C (23.8% of gamma diversity).
The principal component analysis (PCoA) showed no clear pattern in the structure of macrofauna assemblages among the treatments with vegetation (i.e., SA, SLD, and SHD) (Fig. 5). However, samples from the unvegetated treatment tended to group separately from other treatments. PERMANOVA analysis showed significant and independent differences in the structure of macrofauna both for Treatment (Pseudo-F = 3.1521, p = 0.001) and Location (Pseudo-F = 2.5591, p = 0.001) factors ( www.nature.com/scientificreports/ that assemblages inhabiting unvegetated areas were significantly different from those associated with vegetation (UnV ≠ SA = SLD = SHD). Analysis also showed that location C was significantly different from A and B (p < 0.01). www.nature.com/scientificreports/ Multivariate patterns of benthic bacterial communities also showed that communities associated with location C were separated along the first axis from those sampled in locations A and B. Samples from the unvegetated habitats were also mainly separated from the remaining treatments on the lower half of the plot (Fig. 5). PERMANOVA analysis identified a significant interaction between Treatment and Location (Pseudo-F = 1.291, p < 0.05). Pair-wise comparisons showed significant differences between locations B and C for the SHD treatment. Also, unvegetated assemblages in location C were significantly different from those associated with the high density (SHD) ( Table 1). environmental drivers of benthic communities. For macrofauna, the dbRDA (Fig. 6a) indicated that the percentage of fine particles in the sediment, seagrass shoot density, total organic carbon and the concentration of Cu (as proxy for most of other trace elements) were the environmental variables significantly driving community composition. The first two axes explained approximately 18% of the total variation. For bacteria (Fig. 6b), along with those variables significant for the macrofauna, nitrate in the pore water (NO 3 − pw ), coarse sand and phosphate (PO 4 3− ) were also significant. However, a lower percentage of variation was explained (15.5%) compared to macrofauna. The general pattern observed in the PCoA analysis for both communities is also observed here, i.e. the separation of the unvegetated samples from the remaining sites. In both plots, the discrimination between these two communities was mainly due to fine particles in the sediment, concentration of Cu, TOC and shoot density. Nutrients, namely nitrate in the pore water and phosphate seem to be responsible for some separation of the samples from location C regarding the bacterial community composition.
Responses of taxa to habitat drivers. Taxa selected for regression modeling showed species-specific responses to a combination of seagrass properties, temperature, depth, salinity and grain-size. Taxa that were associated with seagrass habitats, and more specifically that were positively associated with higher seagrass biomass, included the polychaete families Eunicidae, Terebellidae, Capitellidae and Syllidae ( Table 2). The brittle star, Amphiura sp., was also negatively associated with seagrass beds (notably seagrass shoot density). Conversely, the bivalve Cardiolucina semperiana was negatively associated with higher seagrass biomass ( Table 2). www.nature.com/scientificreports/ Eunicidae and Cardiolucina semperiana were also negatively associated with higher TOC which is not surprising given TOC was slightly elevated in the seagrass beds which often trap and retain organic matter. Although not very variable, DOC (103-114 µmol C L −1 ) in the overlaying waters followed an opposite pattern, with higher concentrations in the unvegetated samples. A number of macrofaunal taxa were found to be sensitive to increasing metal concentrations in the sediment, with slightly higher concentrations of metals recorded at the unvegetated sites which were also associated with a higher percentage of fines (< 63 µm). For example, the polychaetes Capitellidae, Terebellidae, Eunicidae, and Syllidae were all negatively associated with increasing Ba, Cd and Pb levels ( Table 2). For bacterial communities, genera associated with higher seagrass density included Pelobacter and Tistlia, whereas Maricurvus, the environmental clade genus of the family Spirochaetaceae M2PT2-76, and Candidatus Omnitrophus showed the opposite response (see Table 2). Yet, M2PT2-76 was associated with  www.nature.com/scientificreports/ higher seagrass biomass. This clade showed significant responses with most of the environmental parameters tested. In particular, it was negatively associated with Ba, Cd, and Pb (Table 2).

Discussion
This study provides evidence of the role of seagrass vegetation and local scale environmental processes in influencing benthic community patterns. Here, we have considered both the communities of macroinvertebrates (based on morphological characters) and bacteria (based on 16S rDNA amplicon sequencing approaches) present in the sediment of a Red Sea coastal lagoon. The inclusion of the two endpoints in the size spectrum of benthic heterotrophs (i.e. macroinvertebrates and bacteria) provides a sound basis for the evaluation of the role of seagrass cover as a proxy for habitat fragmentation in an era of escalating degradation of these sensitive habitats 20 . Recent studies have highlighted the utility of benthic bacteria for the assessment of environmental changes using molecular technologies 74 . Assessing responses of both macrofaunal and the microbiome associated with seagrasses may therefore increase our knowledge of the responses of these communities to habitat change.
Vegetation and environmental conditions drive benthic patterns. We observed that the presence/ absence of vegetation (more than its cover) along with changes in local scale sedimentary features (notably grain size, organics and metals) did not significantly influence the amount and characteristics of dissolved organic material in the water column nor did these conditions affect the abundance of free-living planktonic bacteria, including autotrophic and heterotrophic species. We believe that this could be explained by the constant exchange of water between the different sampled locations within the lagoon where the study was conducted. However, the presence/absence of vegetation and different sedimentary features did influence the observed patterns in benthic communities, particularly for macroinvertebrates. Some of the biological parameters analyzed were positively influenced by seagrass density. The high density treatment supported the highest number of exclusive macroinvertebrate taxa (18% of the gamma diversity), twice as much as the exclusive taxa observed in unvegetated areas and more than three times the levels recorded in mixed meadows. These results highlight the role of vegetation, most likely through its secondary effects on local hydrodynamics, food availability and quality, as well as protection from predators, in determining benthic patterns. However, this role was not always consistent across the biological variables analyzed. For example, we did not find unequivocal evidence of a positive correlation between seagrass canopy and the diversity, density and community structure of either bacterial or macrobenthic communities, as responses were site-and taxa-specific. Also, we did not find significant differences between seagrass and mixed meadows containing seagrass and algae. There is a considerable amount of literature assessing the effects of seagrass canopy on associated assemblages, particularly in temperate regions [75][76][77][78] . Less research has investigated differences in benthic assemblages www.nature.com/scientificreports/ between vegetated and unvegetated substrates within coastal lagoons 79 , and particularly in sub-tropical and tropical regions. Seagrass meadows have traditionally been considered as biodiversity hotspots 80 . However, researchers worldwide have reported a broad range of ecological responses of benthic fauna to seagrass canopy with a number of studies that do not support the paradigm of increased diversity and abundance in seagrass beds compared to unvegetated sediments [81][82][83][84][85] . In the present study, depending on the variable analyzed we found evidence of the positive role of vegetation on benthic communities but responses were inconsistent denoting the influence of local scale processes reflective of differences in the environmental characteristics. For example, the combination of characteristics such as the nature of the meadow (i.e., multispecific versus monospecific; species composition 28 , shoot density and biomass 86    www.nature.com/scientificreports/ of the meadow 87 , location of sample collection (i.e., "edge-effect", the increase in diversity or abundance in the borders of the meadow), environmental variables within the meadow 88,89 , seasonality 90 and light regulation 91 , create specific environmental and biological conditions at the local scale that might change the direction and magnitude of the observed responses. The variability in the patterns across the three sites does not allow us to disentangle the relative contribution of each factor in explaining the benthic patterns. For example, in location C, higher similarities were noted among different treatments (e.g., unvegetated versus vegetated), than between the same treatments in different locations (e.g., unvegetated location A versus unvegetated location C). This was particularly evident for benthic bacterial communities and has been previously reported for invertebrates 92 .

Phaeodactylibacter
Spatial patterns of diversity. Our results provide mixed evidence of the positive effect of seagrass cover on macrobenthic faunal density, as a primary influence of local environmental conditions over cover was apparent. Nevertheless, our ordination analysis showed that both macrofauna and bacteria communities in unvegetated sediments tended to cluster differently from those with the presence of seagrasses. Indeed, there was some evidence to support the theory of higher diversity in seagrass beds, with a higher number of exclusive macroinvertebrate taxa being recorded in the highest density seagrass treatment (but not benthic bacteria that peaked at mixed meadows). There are several plausible reasons that can contribute to the high variability in biological patterns. One of them is the temporal dynamics of the seagrass canopy. As demonstrated for temperate regions, significant positive effects of vegetation on diversity and biomass of macrobenthic communities may only be detected when seagrass cover reaches its maximum 90 . Therefore, a sound understanding of the seasonal dynamics of seagrass meadows is critical for a thorough assessment of the effect of seagrass cover in shaping benthic communities. Even at low latitude regions, seasonal differences in seagrass growth are noticeable 93 . Regarding benthic bacteria, biodiversity assessed by means of the number of exclusive taxa does not seem to be promoted by the presence of seagrasses as both unvegetated and mixed meadows harbored higher number of species compared to seagrass treatments. A second potential explanation is the focus of the present study on sediment samples and the associated communities rather than those more dependent on above sediment niches such as leaves. The method used for collecting the samples (hand corer) is more appropriate for infaunal organisms. Macroinvertebrates living in sediments have been suggested to be less affected by fragmentation than large motile fauna, which usually display epibenthic behavior and are more associated with the aboveground biomass. Frost et al. 27 suggested that infaunal species are less likely to show species-specific responses resulting in small shifts in community composition rather than clear changes in density as found for epifauna. Even though differences between unvegetated and vegetated treatments were less evident for bacteria, in one of the locations, significant changes were still observed between the two most contrasting treatments (unvegetated and SHD). One of the reasons may be related to the responses of bacteria to different below ground niches associated with the seagrasses 94 . The strategy used in the present study to establish the different treatments was based on the above ground density but a posteriori we realized that below ground biomass of the root system was not always consistent with the pre-established treatment. This might have resulted in high variability in the communities and the lack of significant changes related to treatment. However, Holmer et al. 95 showed that the stable isotope signals in bacteria inhabiting unvegetated patches was similar to those of adjacent seagrass meadows, indicating that seagrass detritus and exudates may still contribute substantially to available organic matter even in distant, unvegetated patches.
The present results suggest that DOM supplied by seagrass meadows to the water column may fuel the pelagic microbial communities of unvegetated areas nearby (DOM values were consistently above 100 µmol C L −1 and higher than in shallow nearby waters influenced by urban settlement 62 ), ultimately attenuating the difference in abundance of free-living planktonic bacteria between treatments. Organic matter detritus from the vegetated meadows may also become a source of carbon for sediment bacteria inhabiting unvegetated patches nearby, as shown in Holmer et al. 95 . Yet, this was not quantified in the present study. Humic DOM components also showed a positive and significant relationship with Chl a values but not with the protein-like components, suggesting that the humic material from freshly produced DOM accumulated while the protein material was being rapidly consumed by heterotrophic prokaryotes. Heterotrophic bacterioplankton was, indeed, one order of magnitude more abundant than their autotrophic counterparts. Although mean values only reached 3 × 10 5 cells mL −1 in the SLD treatment (minimum values), these low surface abundances of heterotrophic bacteria seem the norm in pelagic environments from the central Red Sea, whether shallow 62 or open waters (unpublished data). This emphasizes the need to complement the structural aspects of the community reported here with functional analysis in the future.
Biological responses to environmental drivers. The overall changes in community composition between treatments of varying cover was attributed to the individual habitat preferences of selected species based on regression analysis. Overall, seagrass biomass had a positive effect on the abundance of dominant benthic taxa while fine particles (i.e., > 63 μm fraction) and metals had a primarily negative effect. Previously, researchers have found a positive effect of fragmentation on the density and species richness when compared to homogenous bottoms 96 . Dense seagrass meadows do not necessarily support high benthic diversity and abundance. The broad range of ecological traits observed in macroinvertebrates may contribute to the inconsistency in patterns that have been reported in the literature regarding the relationship between canopy and associated communities. For example, large tube-dwelling and burrowing species can be affected by high shoot density (and root-rhizome mat) that will limit their capacity to penetrate into the sediment 97 . However, complex and well-developed rhizome structures may not affect small polychaetes, also tube-dwellers, which do not require big areas in the sediment to build their tubes 19 . Seagrass bed cover can also limit the movement of mobile organisms 82 and ultimately affect their ability to escape from predators. In the present study, shoot biomass was most often associated with Scientific RepoRtS | (2020) 10:13550 | https://doi.org/10.1038/s41598-020-70318-1 www.nature.com/scientificreports/ increases in the abundance of some macrobenthic taxa, whose traits can be favored by the presence of the seagrasses. Based on the regression analysis, we observed that certain taxa were preferentially associated with seagrass habitats. Taxa that were positively associated with higher seagrass biomass included polychaetes from the families Eunicidae and Syllidae (mostly mobile and carnivores) as well as Terebellidae and Capitellidae (sessile or of limited mobility, mostly suspension feeders/detritivores). The increase in the trophic trait suspension feeders/ detritivores in areas of higher seagrass was not surprising as flow rates within seagrass meadows have been demonstrated to be reduced resulting in higher rates of deposition of organic matter providing food source for these organisms [98][99][100] . The analysis of the exclusive species associated with each treatment also showed a prevalence of gammarid amphipods associated with high density of seagrasses. These are preferential prey items for fish 101,102 and higher abundance of these crustaceans has previously been reported for seagrass beds when compared to unvegetated sediments possibly due to the protection that seagrass beds provide 83,103 . Interestingly, the bacterial genera Tistlia and Pelobacter, both known for nitrogen fixation [104][105][106] , were also associated with higher seagrass cover. Both genera may form a symbiotic relationship with the seagrass rhizome by providing a supply of fixed nitrogen to the growing seagrass and this relationship has previously been shown for Pelobacter 105 . In contrast, the bivalve Cardiolucina semperiana was negatively associated with both seagrass biomass and TOC possibly due to seagrass roots limiting its population density.

conclusion
The present findings emphasize the spatially dynamic nature and complex interplay of ecological and environmental processes driving benthic communities in a tropical coastal lagoon with heterogeneous bottoms. It is known that habitat loss or fragmentation can affect the resilience and persistence of vegetated habitats 107 and in turn influence faunal connectivity among habitats, species diversity and ecological interactions 108,109 . Variables found to be important in predicting the response of benthic communities were related to processes working across multiple spatial scales. Conservation strategies should prioritize connectivity between habitats to maintain natural migration of species 110 and to preserve species that are restricted to mature habitats 111 and references therein . Conservation measures should thus consider both vegetated and unvegetated habitats, which will also allow for the evaluation of potential changes associated with anthropogenic and natural impacts given the dynamic nature of the ecological linkages across the seascape 112,113 . Considering that only approximately 25% of invertebrate taxa were shared among treatments or locations, with 11% shared OTUs for benthic bacteria sampled, and that high density treatments supported the highest number of exclusive taxa for macroinvertebrates but not bacteria (highest at unvegetated and mixed meadows), sustaining ecological connectivity across the seascape is therefore critical in maintaining the biodiversity at regional scales.