The study was conducted in a monospecific meadow of H. stipulacea, located on the western shores of the GoA, Northern Red Sea, Israel (29.497664°N, 34.912737°E). At this site, the shoreline is almost pristine, with very little coastal infrastructure or anthropic pressure (Mejia et al., 2016). This H. stipulacea meadow is considered healthy and harbors a rich and diverse community of corals (Winters et al., 2016). During the study (October 2013), the meadow had an area of about 61,900 m² commencing at 4–5 m and extending to depths greater than 50 m (Sharon et al., 2011; Winters et al., 2016).

Sampling took place in October 2013, with average sea surface water temperature of 24°C and salinity of 40.70 PSU (NMP, 2014). Sampling was conducted using SCUBA along 50 m transect lines placed, parallel to the shore, at 4, 9, 18, and 28 m (see Figure 1). The sampling depths, selected based on recent mapping of the meadow (Winters et al., 2016), cover the different parts of the meadow, as follows: at 4 m, near the coral reef and the shallow edge of the meadow, at 9 m along the flat part of the sea bottom, before the drop to deeper waters, at 18 m, along the actual slope of the meadow, and at 28 m around the area where the slope straightens out again. To characterize the environment, three sediment samples were collected at random sites along the transect line with minicorers (3.0 cm diameter and 50 ml volume) and used for grain size (granulometry) measurements and Total Organic Carbon (%TOC) content analyses. The Photosynthetic Active Radiation (PAR) in the water column was measured at noon, from 0 to 28 m, using the 2π-quantum sensor of a diving-pulse amplitude modulated (PAM) fluorometer (Walz, Germany), with 3–4 light measurements made every meter. The diffuse attenuation coefficient Kd (PAR) was then calculated based on measurements made at 1 and 28 m.

At each depth, samples were taken at three sites, randomly selected, within a 5 m range, at the beginning, middle and end of the transect line. Intact plants (with rhizomes and roots) were collected at each site, separately for morphometry (five plants), biochemical analyses (eight plants), and microbial community analyses (eight plants; Figure 1). For the latter, the eight plants were pooled to obtain two replicates, separating underwater aboveground (leaves) and belowground (rhizomes and roots) compartments. Each replicate of the aboveground compartment was composed of at least 32 leaves, and the belowground of about 16 cm of rhizomes with roots. Four plants per replicate were pooled to obtain a representative measure of bacteriome composition at the ecosystem level, while still permitting assessment of within- and between-depth bacteriome variability through the inspection of two independent replicates per depth below- and aboveground. All samples for morphometry, biochemical and microbial community analyses were then transported to the lab and kept cold (<10°C), in the shade until further processing (within 2 h).

To determine granulometric composition, sediment samples were pretreated to remove organic matter with 20 ml of 15% hydrogen peroxide (H₂O₂) and stirred for 15 min; the solution was rinsed out once bubbling ceased. The samples were then freeze dried for 72 h. Particle size distribution in the 0.02–2000 μm range was measured using a Beckman Coulter LS 13 320 Laser diffraction particle size analyzer. Sediment particles >2000 μm were weighted separately using an analytical scale. Granulometric composition was evaluated by calculating the percent frequency of seven class sizes: <63, 63–125, 125–250, 250–500, 500–1000, 1000–2000, and >2000 μm.

Sediment samples were also used for determining the %TOC in the sediments. This was measured using a Skalar Primacs SLC TOC Analyser (Skalar Analytical BV, Netherlands). For this, dry samples were grounded and split into two fractions. The first was burned at 1050°C for total carbon measurement and the second was titrated with 20% phosphoric acid to detect the inorganic carbon fraction. The subtraction yielded the %TOC.

Leaf morphometrics were measured using ImageJ software (version 1.47; Abramoff et al., 2004). Leaves from each depth, were scanned (Canon Lide 110 scanner) and images were used to measure leaf length (mm), width (mm) and area (mm²); at least 10 leaves from each plant were used for these morphometric measurements. Total phenols were extracted from the leaves (200 mg fresh weight) and rhizomes (200 mg fresh weight) of each plant in duplicate to increase confidence and quantified according to Migliore et al. (2007). Photosynthetic pigments (Chlorophyll a and Chlorophyll b, total Carotenoids) were extracted in duplicate from leaves (250 mg) of each plant, according to Wellburn (1994), modified by Rotini et al. (2013). Both phenols and photosynthetic pigment contents are expressed as mg/g of fresh weight (FW).

Microbial sample processing in the laboratory and sequencing data analyses followed the procedures reported in Mejia et al. (2016). Microbial cell pellets were retrieved separately from leaves and rhizomes-roots collected at each depth, followed by the extraction of microbial metagenomic DNA with the Power Soil^® DNA isolation kit (Mo Bio, Carlsbad, CA, USA) according to the manufacturer’s instructions. Pure DNA extracts were sent to Molecular Research LP¹ (MR DNA, Shallowater, TX, USA; Dowd et al., 2008) for 16S rRNA gene PCR amplification and 454-pyrosequencing, using Roche GS FLX Titanium technology. The universal primers Com1 (forward, 5′-CAGCAGCCGCGGTAATAC-3′) and Com2 (reverse, 5′-CCGTCAATTCCTTTGAGTTT-3′) were used for the subsequent PCR amplification of bacterial-specific fragments of 407 bp, encompassing the V4 and V5 regions of the 16S rRNA gene, which are highly variable phylogenetically (Schmalenberger et al., 2000). The raw 454-pyrosequencing sequences were processed with the open source software MOTHUR² (accessed: August 2015), according to the 454 Standard Operating Procedure (SOP) pipeline (Schloss et al., 2009). Within MOTHUR, sequences were depleted of their barcodes, primer sequences and homopolymers, and those <200 bp and with ambiguous base calls were removed. Sequences were then dereplicated, aligned against the Greengenes core-set template alignment, screened and filtered to make them overlap in the same region. Putative chimeras were identified with the “uchime” algorithm and removed. Operational Taxonomic Units (OTUs) were defined by clustering at 3% divergence (97% similarity), using the average neighbor-clustering algorithm. A final matrix (or “complete dataset”) of the OTUs abundance by sample and their taxonomic identification was built in MOTHUR and used in downstream analyses. Before OTU picking, sequences diverging by one or two nucleotides only from a more abundant sequence were considered as if they were representative reads of the most abundant sequence. The most abundant read within each OTU was then chosen and used as representative OTU sequence for taxonomic assignment. OTU representative sequences were classified with naïve Bayesian classifier in MOTHUR using the RDP taxonomy reference database-trainset 10_082014. “Unknown” sequences and those identified as “Mitochondria” were removed, while chloroplast reads were kept to enable further insights into microalgal diversity in the samples. To standardize differences in sampling effort among samples, the complete OTUs dataset was normalized by random subsampling to the common depth of 2,764 sequences (hereafter called “normalized dataset”), which was the lowest number of sequences produced by the samples in the study; this allows for adequate comparisons at community level. The complete set of raw sequences obtained in this study was deposited in GenBank at the Sequence Read Archive (SRA) under the BioProject no. PRJNA317018.

Bacterial community analyses were performed in MOTHUR. Rarefaction curves were first built using the complete dataset to review the overall depth of sequencing of the study. Then each sample library was reduced to 5528 sequences (and a second rarefaction curve was built) to assess the bacterial richness coverage of the microenvironments analyzed. Comparative, alpha and beta diversity analyses were undertaken using the normalized dataset (2,764 sequence reads per sample) as follows: (i) Shannon index (H’) was calculated to estimate bacterial diversity; (ii) Venn diagrams were built to reveal the specific and shared bacterial OTUs by sample type, after merging the OTU profiles of the replicate samples; and iii) stacked bar plots were built to visualize bacterial community composition at high taxonomic ranks. The SIMilarity PERcentage (SIMPER) test was used to identify the contribution of the different taxa to the dissimilarity, using Bray-Curtis distances, in bacterial composition between samples. Analyses used Hellinger-transformed OTU data (square root of OTU relative abundances).

For environmental variables, multivariate non-parametric analysis of similarities (ANOSIM) on Euclidean distances, treating the size classes as variables, was used to test whether significant differences in granulometric composition among depths occurred. Differences in sediment %TOC and water column PAR were evaluated using ANOVA and post hoc Tukey’s pairwise tests.

For plant descriptors, significant differences in the mean leaf surface area among depths were evaluated using Kruskal–Wallis and Mann–Whitney pairwise tests, since data did not satisfy the homoscedasticity assumptions (Levene’s test). Differences in photosynthetic pigments and total phenols were evaluated using ANOVA and post hoc Tukey’s pairwise tests.

For bacterial community diversity, the non-parametric test Mann-Whitney was used to evaluate pairwise differences between depths. The normalized samples vs. OTUs contingency table and corresponding metadata were used as input for multivariate ordination analysis with the software package Canoco for Windows 4.5 (Microcomputer Power, Ithaca, NY, USA), following the procedures of Costa et al. (2006) and Keller-Costa et al. (2014). To determine the species distribution model (linear or unimodal) that fits best to our data. Detrending Correspondence Analysis (DCA) was employed to assess the extent of variation (“lengths of gradient” as implemented in Canoco) within the OTU dataset. A linear model was considered the best fit, and thus unconstrained ordination via Principal Components Analysis (PCA), was used to explore possible correlations between bacterial community structures, plant ecophysiological descriptors and environmental variables. Analyses used Hellinger-transformed OTU data (square root of OTU relative abundances).

Pearson’s correlation coefficient was used to test for relationships between environmental variables and plant descriptors using log-transformed data, and also with bacterial communities (square root transformed data). Statistical significance was defined at p < 0.05 (Supplementary Table S1).

The granulometric composition of the sediments changed significantly along the depth gradient (ANOSIM, R = 0.53, p < 0.01; see Supplementary Figure S1), from gravel/coral rubble at 4 and 9 m (>40% of sediment grains sizes >2000 μm) to coarse sand at 18 and 28 m (>50% of sediment grains between sizes 250 and 2000 μm). The %TOC showed significant differences among depths (ANOVA, F = 9.155, p < 0.01). Significantly higher values were found at 9 m (0.48%) and 18 m (0.51%), compared to 4 m (0.10%) and 28 m (0.14%) (post hoc Tukey’s test, p < 0.05, see Supplementary Figure S1).

The photosynthetic active radiation (PAR) at noon, decreased significantly with depth (ANOVA, 259 F = 681.5, p < 0.01 and post hoc Tukey tests, p < 0.001). PAR was 182.3 ± 8.3, 258 87.7 ± 3.5, 38.7 ± 3.5, and 16.7 ± 1.2 μmol photon m^-2 s^-1 at 4, 9, 18, and 28 m depth, respectively.

The mean leaf surface area, width, and length (Table 1, Leaf morphometrics) showed significant differences among depths. Mean leaf area increased significantly with depth (Kruskal–Wallis, H = 39.14, p < 0.01) and showed a significant negative correlation with PAR (R = –0.751, p < 0.01).

Photosynthetic pigment content (Figure 2A) increased with depth, as did the ratio of Chl _a+b/Car, while the ratio Chl a/b decreased (Figure 2B). Significant differences along the depth gradient were found in pigments and in the ratios Chl _a+b/Car and Chl a/b (ANOVA, p < 0.001) in particular, between 28 m and all the other depths (post hoc Tukey tests, p < 0.001). Photosynthetic pigment contents and ratios showed significant correlations with PAR, granulometry and plant descriptors (See Supplementary Table S1). The mean total phenol content (Figure 2C) showed significant differences only in leaves at 28 m where the lowest value was found (post hoc Tukey test, p < 0.001), while rhizomes did not show significant differences along the entire depth gradient.

The raw output of the 454-pyrosequencing resulted in 107,782 bacterial 16S rRNA gene sequences, comprising all samples. The processing of sequences with the MOTHUR software resulted in 84,892 sequences, assigned to a total of 6,793 OTUs, including those occurring only once or twice. On average, each individual sample was represented by 5,306 sequences, ranging from 2,764 to 9,001 sequences. The rarefaction curves analysis applied to the complete dataset (see Supplementary Figure S2A) was useful to demonstrate that an asymptote in OTU richness was not fully reached, even for those microenvironments characterized with a high sequencing effort (e.g., aboveground compartments at 4 m, and belowground compartment 28 m). Sectioning each curve at a common sequence depth of 5528 gene reads per sample type (see Supplementary Figure S2B), showed that in both leaves and rhizomes, communities sampled at 9 m were the least rich. Within each depth, belowground communities displayed, overall, higher bacterial richness values than the corresponding aboveground communities.

At most depths, the Shannon diversity index (H′) values, using the normalized dataset, were slightly lower aboveground than belowground (Table 2). The highest diversity values in both above- and belowground plant parts were observed at 4 m. Diversity values aboveground fluctuated among depths but belowground diversity decreased with depth. Diversity indices between 4 and 28 m were significantly different (Mann–Whitney, p < 0.05).

Halophila stipulacea-associated bacterial communities (below and aboveground) were dominated by the phylum Proteobacteria (Figure 3), which accounted for 72% of the analyzed sequence reads across all samples. Within this phylum, the class Alphaproteobacteria was the most dominant, representing 52% of the whole community composition in the study, being more abundant aboveground (60%) than belowground (45%). Within this class, the order Rhodobacterales (71%) and family Rhodobacteraceae (71%) were the most abundant. The second and third most abundant classes were Gammaproteobacteria (10%) and Deltaproteobacteria (8%), and both were more abundant belowground than aboveground. Surprisingly, only 12% of the Gammaproteobacteria reads could be assigned to up to seven formally recognized orders, with Alteromonadales (8%) being the most abundant, while the remainder of the reads were considered unclassifiable or possessing uncertain placement at the order level. The Deltaproteobacteria community was mainly composed of bacteria belonging to the orders Myxoccocales (33%) and Desulfobacterales (50%), out of seven orders observed. The latter comprised the family Desulfobulbaceae (47%), which was the most abundant among the 12 families observed, as well as the genus Desulfopila (27%), the most abundant of the 32 genera observed. Overall, the 16 phyla detected comprised about 29 classes, 40 orders, 82 families, and over 100 genera. Bacteria not classifiable at the phylum level and chloroplast-derived sequences accounted for 5.2 and 6%, respectively, of the whole community composition.

The SIMPER test (see Supplementary Table S2) showed dissimilarities of 20.91% in the bacterial community structure between plant parts. Community dissimilarities across depths were slightly lower aboveground (range of 15–19%) than belowground (range of 19–24%). The top five taxa contributing to dissimilarities were Rhodobacteraceae (734 OTUs, 16722 seqs), Gammaproteobacteria (711 OTUs, 4240 seqs), Desulfobulbaceae (145 OTUs, 1579 seqs), Roseibium (33 OTUs, 2792 seqs), and Desulfopila (14 OTUs, 896 seqs). The combined sequences of the groups mentioned above accounted for about 60% of the normalized dataset.

The distribution of specific and shared OTUs showed differences by plant part and depth (Figure 4). In general, the total number of OTUs observed aboveground was lower than belowground (2,598 vs. 3,166 OTUs, respectively; using the normalized dataset). Bacterial communities showed a higher number of specific than shared OTUs. The shared OTUs revealed the ‘core bacteriome’ of H. stipulacea by plant part along the gradient (see Supplementary Table S3). These OTUs usually displayed high abundance values, within the range of 100–1000 sequences per OTU and represented 63 and 52% of the total number of sequences observed above and belowground, respectively. The most abundant phylotype found in the core bacteriome, belonged to the phylum Proteobacteria, class Alphaproteobacteria.

Ordination diagrams obtained by Principal Components Analysis (PCA, Figure 5) showed the different forces shaping the bacterial communities along the gradient. Overall, aboveground and belowground samples displayed a clearer and stronger separation by plant part than depth (Figure 5A). The communities at the edges of the meadows, 4 and 28 m, depicted best the simultaneous influence of parts and depths (Figure 5B), with the latter factor exerting greater influence on the structure of belowground communities than on aboveground communities. Both aboveground (Figure 5C) and belowground (Figure 5D) communities showed higher within-depth variability at 9 and 18 m than at 4 and 28 m, where replicate samples (within each depth) were clearly more similar to each other. The PCA diagrams illustrated the grouping of aboveground communities from 28 m which had the largest leaves, highest carotenoids and chlorophyll contents, and lowest phenols (Figure 5C); and of belowground communities with lower TOC values and smaller grain sizes.