Sediment Microbiota as a Proxy of Environmental Health: Discovering Inter- and Intrakingdom Dynamics along the Eastern Mediterranean Continental Shelf

Biodiversity is the keystone of ecosystem functioning and environmental health. At that, microbiota comprise the foremost part of the biodiversity and provide the majority of ecosystem services (9). We obtained a high-spatial and -temporal resolution map of the microbiota along the EMS shelf transect. Integrating data across the three main kingdoms was key to resolving the main drivers of the SM assembly (Fig. 2, Fig. S2). Clustering of samples was highly supported in the integrated data set and much less so when considering each kingdom separately. Interestingly, the contribution of Archaea and Bacteria to the clustering of samples was higher than that of Eukaryota. This phenomenon might be explained by the high sensitivity ascribed to Archaea and Bacteria with respect to environmental conditions (30). Another explanation may be the difference in homogeneity of distribution of organisms within the sediment, with patchier distribution of Eukaryotes, particularly meiofauna or larger organisms, compared to prokaryotes (31, 32).

The continental shelf of Israel was characterized by several key properties, including particle size distribution, total organic carbon (TOC), water content and carbonate content (150 km north to south) for seafloor depths between 10 m and 100 m, attributing direct correspondence between seafloor depth and the measured parameters (28). Previous studies of marine and estuarine sedimentary systems show that the properties of sediment are key to determining the composition of microbiota (15, 32–35). Here, seafloor depth, representing these properties, explained 11% of the variance among samples (Fig. 2). Similarly, a survey of costal sediment of the Mediterranean Sea and French Atlantic Ocean estimated that grain size and TOC contribute 10% of the variance in microbiota composition (33). Second in importance was the core-depth variable, which explained 4.2% of microbiota variance among samples (Fig. 2). This significant effect can be attributed to dramatic shifts in the physicochemical conditions within the sediment (36, 37). Distinct physiologies and lifestyles may characterize microbiota even millimeters apart within the sediment core (37, 38). The factor of season was third in relation to its contribution to microbiota variation, amounting to only 2.5%. In sediment systems more exposed to the atmosphere (e.g., rivers, estuaries [39] coastal sands [40, 41], marine water column at shallow seafloor depth [42]) a large portion of the variation in microbiota composition was attributed to season.

Major variation in the composition of all kingdoms was observed (Fig. 3a). Gammaproteobacteria and Desulfobacteria were the dominant bacterial classes in this study and were often found to be dominant in marine sediments (15, 36, 43, 44). One major observed effect on bacteria was the reduction in RA of Gammaproteobacteria from the shallow to deep sites. A compensative increase in RA was attributed to a large number of highly diverse populations, related to over 150 different classes within 61 phyla. In agreement, overall bacterial diversity significantly increased with seafloor depth (Fig. 3b, Table S5).

In Archaea, similarity in the RA of dominant classes was higher between 10 m and 25 m than between any other pair. This similarity may be explained by higher exposure to physical impacts such as waves or storms, and water mixing in shallow depths, or may reflect a dispersal barrier. The four main classes of Archaea are ubiquitous in marine sediment (15, 45). Nitrosphaeria perform ammonia oxidation, even at low concentrations, which characterizes the open ocean (46, 47) and the EMS, in particular (48). The RA of Nitrososphaeria was dramatically higher at the deep site, which is characterized by small grain size (Fig. 3a, Fig. S1). A study of estuarine sediments found negative relatedness between ammonia-oxidizing Archaea abundance, diversity, and grain size (49). The increase in Nitrosphaeria RA was accompanied by a proportional decrease in RA of the classes Thermoplasmata and Nanoarchaeota. The latter group is largely composed of organisms of reduced genome size that, based on current evidence, are considered to be obligate symbionts of other Archaea (45). Recent metagenomic surveys identified several Nanoarchaeota hosts from the class Thermoprotei of Crenarchaeota (48–52). However, Nanoarchaeota were first observed attached to Archaea of the order Thermoplasmatales (Thermoplasmata) in biofilm (53). It remains to be determined whether the Nanoarchaeota found here maintain a symbiotic lifestyle and, if so, what are their hosts.

Polychaeta was the dominant group of eukaryotes in the sediment (Fig. 3c) as in other marine and estuary sediments (54). The RA of Polychaeta was highly variable among sites, ranging from 6% (25 m) and 51% (100 m) with no direct correspondence to seafloor depth. This group plays a major role in the function of benthic communities, responsible for bioturbation of the sediment and burial of organic matter by their burrowing and feeding activity (55, 56). Therefore, their activity radically alters niche properties and, as a result, engineers microbiota assembly and structure (37, 57). Like the Polychaeta, patchy distribution characterized other dominant Eukaryotes. Those include Dinophyceae (Alveolata), Malacostraca, Ostracoda (Arthropoda), and Eurotatoria (Rotifera). These groups include many pelagic species that accumulate in the sediment following their death and sinking (58, 59). Thus, their patchy pattern may be related to abiotic and biotic processes in the water column.

Biological indicators are considered preferable to hydrological or chemical surveys, as biological entities transmit an integrated and cumulative measure of ecosystem variations and/or disturbances, rather than a direct snapshot. Many marine studies have focused on the water column and have demonstrated the ability of environmental genomic surveys to detect environmental shifts and disturbances, even in real-time (60–62). However, biotic and abiotic factors vary temporally (e.g., day versus night) and along the water column (63, 64), and episodic events (e.g., contaminant spill) may dissipate rapidly in the water body, leaving no signature (65). Therefore, frequent water sampling is required to identify the disturbance or standardize sampling to represent specific aspects. In comparison, the sediment is more stable over time, is accumulative in its nature (66), and can register past events, even episodic ones (67). Considering overall effort, sediment monitoring may better serve local -to regional-scale regulatory needs.

Two approaches were examined to demonstrate the applicability of the SM model: (i) the specific markers and (ii) the whole community as a marker.

Many studies attempt to identify specific microbial markers of environmental status based on limited (even a single snapshot) sampling events. However, in order to obtain robust markers, sampling must be repeated to reach stability and assess reproducibility (65, 68). The robustness of our SM model was demonstrated by high reproducibility of specific markers for seafloor depth between the main samples and a set of test samples collected over the study period (Fig. 4).

The marker approach can be further expanded to incorporate not only a set of markers, but also cooccurrence patterns. The cooccurrence network that included seafloor-depth microbial markers demonstrated both intra- and interkingdom connections (Fig. 5). For example, a strong positive correlation between a nematode marker (Euk_41) and a Cytophagales population (Bac_406) was demonstrated (Fig. 5, cluster IV). These two groups were found to be codominant in a previous study of marine costal sediments (13). Cooccurrence of species may reflect shared ecological filtration (i.e., selection) or codependency (e.g., symbiosis). In addition, stochastic processes, such as dispersal limitation, are also manifested by high cooccurrence, as was demonstrated for archaeal communities in coastal sediment (69) and is indicated in the network here. For this reason, examination of cooccurrence networks is appealing, particularly in cross-kingdom data sets. Previous studies have used microbial cooccurrence networks to resolve between different biotic or abiotic conditions (70–73). Hence, analysis of cooccurrence patterns may be useful in detecting changes and perturbations. An additional set of markers was identified by direct analysis of correlation between seafloor depth and RA of different taxa (Fig. 7). This approach may be applicable for gradients in environmental conditions. In this setup, the relevant underlining gradients are sediment grain size and TOC. In order to improve the sensitivity and reproducibility of markers, continued sampling and model refinement are required. Furthermore, specificity of the markers and reliability for biomonitoring may be improved by application of more advanced computational methods, such as machine learning, that were shown to outperform more traditional methodologies, including indicator value and random forest analysis (74).

While providing a robust mean for classification of samples, the biomarker or specific marker approach is limited in the sense that it absolutely relies on prior knowledge and assumes a steady state. Natural habitats change over time, with or without perturbations, affecting the validity of any marker. Furthermore, specific markers may not be affected by one or another type of perturbation, disabling the detection of an incident or a process. A whole-community approach may therefore have an advantage, as it enables identification of changes over time (30) and relies on measurement of similarities rather than a priori determined criteria. Whole-community examination of SM by NMDS analysis (Fig. 6) demonstrated the applicability of this approach. Samples collected at the test sites were ordinated in a manner expected based on their seafloor depth. Naturally, detection of disturbances is a key goal for monitoring efforts. Samples collected from HM, a highly disturbed site, were markedly divergent from model sites and across the three kingdoms (Fig. 6). Marinas are a common perturbation to the marine environment, typically physically confounded, reducing water circulation, intensely and chronically impacted by diverse anthropogenic activities. Among common pollutants in marinas are organic matter, particularly oil and oil products (75, 76). The intensity and diversity of pollutants was reflected in SM composition (Fig. 8, Table S7). Ciliophora and Dinophyceae, previously associated with organic pollution and oil pollution (77, 78), were higher at the marina than the model site. Similarly, Gammaproteobacteria and Bacteroidetes that increased in RA at the marina were previously identified as dominant bacteria in a chronically oil-contaminated lagoon (79). Anaerolinea (Chloroflexi) was denoted for its contribution to tetrabromobisphenol A dissipation in mangroves (80). Asgardarchaeota, identified as a marker of HM, was previously reported to be associated with petroleum hydrocarbons and polychlorinated biphenyls (81).

Conclusions. This study examined the suitability of marine SM as a bioindicator for environmental health in monitoring programs of the EMS. The main factor determining microbiota composition in EMS costal shelf sediment was the seafloor depth, which in this context, represents basic sediment characteristics: grain size and TOC. This enabled a robust and sensitive model that is applicable for sites outside the monitoring stations and that detects disturbance. Identification of robust markers was feasible. However, further study is required in order to validate them, understand their ecological role, and elucidate their mechanisms that determine niche preference.

Sediment was sampled twice a year (winter and summer) for 3 years (2017 to 2020) along a four-seafloor-depth transect (10, 25, 45, and 100 m) near Sdot-Yam (SY), an undisturbed site (Fig. 1, Table S1). These four sites are termed model sites in this article. To further expand and establish the model, additional sites were selected, Ashdod and Michmoret (70 m and 38 m seafloor depth, respectively), designated the test sites and were characterized by their low anthropogenic impact. Additionally, sediments were taken at Herzliya Marina (HM) and the Orot Rabin power plant in Hadera (HPP) (3 m and 6 m seafloor depth, respectively) since these sites experience high rates of anthropogenic impact and are designated disturbed sites (Fig. 1, Table S1). For model sites, the total organic carbon (percent TOC) and grain size were measured, as previously described (25, 82) (Fig. S1). For further description of the study sites, see File S1.

Two cores (6-cm diameter, 20 cm long) from each seafloor depth were hand-sampled during SCUBA dives from 10- to 45-m depths or sampled using a box corer from 100-m depth. The samples were chilled until transfer to the laboratory. There, each core was split to 1-cm slices, up to 10-cm core-depth, with sterile tools. The slices were kept in separate tubes at −20°C until further processing. Table S2 lists the samples used in this study.

Samples were processed directly following each sampling event. DNA was extracted using the DNeasy PowerSoil kit (Qiagen, Valencia, CA, USA) following the manufacturer’s instructions. Appropriate controls were included in all parts of sample preparations (no input reactions). For Bacteria and Eukaryota, primers were used for PCR amplification of small-subunit (SSU) rRNA gene fragments as described in The Earth Microbiome Project (83, 84). For Archaea, the primer pair used was published in Takahashi et al. (85). Please see File S1 for more details. The raw sequence data are available in the NCBI SRA database under BioProject accession PRJNA910589.

Sequence data were processed using the DADA2 pipeline (86). For each sequencing run, a separate analysis was conducted for quality trimming, error model estimation, sequence error correction, and amplicon sequence variant (ASVs) inference and quantification. The ASVs and count tables of all runs were merged, and suspected chimeras were removed. For each ASV, taxonomy was inferred by alignment to the Silva nonredundant small-subunit rRNA database (v138) using the DADA2 command assignTaxonomy, setting the minimum bootstrap value to 80%. For Eukaryota ASV sequences, taxonomic inference was achieved by the last common ancestor method (LCA) in MEGAN (MEGAN6, community addition) using as input results of a BLASTn search against the NCBI nucleotide database with the 50 best hits (LCA parameters: minimum score, 100, E-value, <10⁻⁹; a coverage, >80%). Table S3a to c present the raw count tables and ASV inferred taxonomy.

Out of a total of 218 samples taken from the Sdot-Yam site, we obtained full high-quality results, including Archaea, Bacteria, and Eukaryota representation for 156 samples (Table S2). This set of samples, here termed main samples, served for the analysis of microbiota composition and structure. Samples from the Sdot-Yam site for which only partial data were obtained (i.e., not all three kingdoms) are here termed test samples. In order to take advantage of information resulting from the three data sets per sample (i.e., Archaea, Bacteria, and Eukaryota), we chose the similarity network fusion (SNF) approach as implemented in the SNF tool (87) For SNF analysis, each of the data sets was filtered to retain only ASVs detected in at least 10 of the 156 samples. For each data set, counts were log-ratio-transformed, and a distance matrix was calculated. The R package SNFtools (v2.3.1) was used for the calculation of affinity matrices from the distance matrices (k = 15, sigma = 0.5) and the fusion network. The number of clusters supported by the fused network was calculated by the “rotation cost best” method, and a subsequent concordance matrix was calculated. Silhouette analysis was also performed to support the chosen number of clusters (R package cluster [v2.1.3]). For comparison, networks were calculated from each of the data sets independently, using the same methodology. Additionally, using the fused network (Fig. 2), the ASVs from each kingdom were ranked by their normalized mutual information (NMI). The top 10 ASVs for each kingdom at seafloor depth were recorded as markers of the network clusters.

To examine the contribution of the environmental factors (seafloor depth, core depth, and season) to variation in microbiota composition in the main samples, permutational analysis of variance (PERMANOVA) was performed (command adonis in the R package vegan [v2.5.7]) based on Euclidean distances on log-ratio-transformed counts. Additionally, for model and auxiliary samples in each data set, nonmetric multidimensional scaling (NMDS) analysis was performed. NMDS was based on Bray-Curtis dissimilarities, calculated after cumulative sum square normalization of raw count data.

For each of the three data sets, raw counts were subsampled to the minimum library size, and then a Shannon index was calculated (diversity function in the R package vegan). In order to test differences in Shannon index as a function of the environmental variables, aligned rank-transformed (ART) analysis of variance (ANOVA) was performed (88). As some of the season or the core-depth levels were missing in the model main sample data set, two models were tested: (i) Shannon ~ seafloor-depth × core-depth and (ii) Shannon ~ seafloor-depth × season. Differences were considered significant for a P value of <0.05.

Linear discriminant analysis effect size (LEfSe) analysis was chosen to calculate differential abundance and identify putative markers. This method is effective in determining which features, in this case ASVs, are most likely to explain observed differences among factor levels (89). LEfSe was performed using the online Galaxy module (http://huttenhower.sph.harvard.edu/galaxy) for three sets of data: (i) main samples, (ii) test samples (in total, 62 samples; 28 Archaea, 39 Bacteria, and 47 Eukaryota), and (iii) disturbance set-HM versus main samples from 10-m seafloor depth. In all cases, Cumulative Sum Scaling (CSS) normalized count data were used as input. LEfSe was applied to each of the three sets separately. The factors seafloor depth (as main class) and core depth (as subclass) were used for detection of putative markers for seafloor depth. Test parameters were set to default, but the LDA effect size threshold was set to 1.25 for the main and test sets of samples and to 2.7 for the disturbance set of samples. For LEfSe results of the disturbance set, cladograms linking ASV markers based on taxonomic relatedness were drawn using the R package igraph (v1.2.0).

For the set of identified markers for seafloor depth (above), cooccurrence with additional ASVs from all three kingdoms was examined. For this purpose, Spearman correlations were calculated between this set of ASVs and all other ASVs. Significant (false-discovery rate [FDR]-adjusted P < 0.05) and strong positive correlations (Spearman rho ≥ 0.75) were selected and combined to provide a list of associated ASVs. Then, Spearman correlations among all ASVs in the combined list were calculated. Only significant correlations (P < 0.05, Spearman rho ≥ 0.75) were maintained. An affinity matrix was then calculated based on this data set and used to create an association network using the R package igraph.

Spearman correlations were calculated using count data from all samples excluding the disturbance set, summed to genera, and log-transformed. Spearman correlations were calculated and tested using the command rcorr in the R package Hmisc (v4.6.0). Correlations were considered significant with the Benjamini-Hochberg FDR adjusted P value was <0.05.

The data sets generated and analyzed during the current study are available in the NCBI SRA repository, under accession number PRJNA910589.