Deciphering the community structure and the functional potential of a hypersaline marsh microbial mat community

Abstract

Microbial mats are vertically stratified communities of microorganisms characterized by pronounced physiochemical gradients allowing for high species diversity and a wide range of metabolic capabilities. High Throughput Sequencing has the potential to reveal the biodiversity and function of such ecosystems in the cycling of elements. The present study combines 16S rRNA amplicon sequencing and shotgun metagenomics on a hypersaline marsh in Tristomo bay (Karpathos, Greece). Samples were collected in July 2018 and November 2019 from microbial mats, deeper sediment, aggregates observed in the water overlying the sediment, as well as sediment samples with no apparent layering. Metagenomic samples’ coassembly and binning revealed 250 bacterial and 39 archaeal metagenome-assembled genomes, with completeness estimates higher than 70% and contamination less than 5%. All MAGs had KEGG Orthology terms related to osmoadaptation, with the ‘salt in’ strategy ones being prominent. Halobacteria and Bacteroidetes were the most abundant taxa in the mats. Photosynthesis was most likely performed by purple sulphur and nonsulphur bacteria. All samples had the capacity for sulphate reduction, dissimilatory arsenic reduction, and conversion of pyruvate to oxaloacetate. Overall, both sequencing methodologies resulted in similar taxonomic compositions and revealed that the formation of the microbial mat in this marsh exhibits seasonal variation.

16S rRNA, extreme environments, metagenome-assembled genomes (MAGs), salt crust, sediment, shotgun metagenomics

Introduction

Microbial mats are vertically stratified communities of functional groups of microorganisms embedded in an organic matrix, which may also contain minerals such as silicates and carbonates (Stal 2012, Bolhuis et al. 2014, Prieto-Barajas et al. 2018a). They grow on a solid substrate (e.g. sand) and the vast majority of microbial mats utilize inorganic carbon as carbon source, hence they are autotrophic (Bolhuis et al. 2014). Microbial mats are highly productive ecosystems (Villagrasa et al. 2019) characterized by pronounced physiochemical gradients, which allow for the presence of high species diversity, encompassing a wide range of metabolic capabilities; thus, mats are ideal models to study a whole ecosystem (Al-Thani et al. 2014) and are considered as natural laboratories (Villanueva et al. 2007). These physicochemical gradients provide microenvironments for various microbial functional groups, which exhibit a certain physiology with which they fulfil a specific function (van Gemerden 1993).

Microbial mats comprise an extensive diversity of microorganisms belonging to different species, which are embedded in a matrix of extracellular polymers (EPS) and exchange signals and nutrients, thus enabling a flow of resources and energy for the survival of the overall community (Ruvindy et al. 2016, Prieto-Barajas et al. 2018b). The role of microbial mats has been vital throughout Earth’s history since they produced and released reduced gases, e.g. O₂, H₂, and CH₄, in the early Earth’s atmosphere (Hoehler et al. 2001) and they are regarded as modern analogs of early ecosystems (Krumbein et al. 2003). In addition, they are considered to constitute the first ecosystems, along with stromatolites (Prieto-Barajaset al. 2018b), and probably are the oldest (van Gemerden 1993) and simplest of the self-organized structures that may have first appeared on earth (Guerrero et al. 2002).

Regardless of the vertical structure, marine microbial mats are comprised of four main functional groups: (i) oxygenic phototrophs (CYN; primarily Cyanobacteria), (ii) aerobic heterotrophic bacteria (HET), (iii) sulphate-reducing bacteria (SRB), and (iv) sulphide-oxidizing bacteria (SOB; Visscher and Stolz 2005). It is suggested that archaeal abundance and diversity can be limited in microbial mats (Hugenholtz et al. 1998, Sievert et al. 2000, Robertson et al. 2009, López-López et al. 2013), with halophilic and methanogenic taxa to be the most predominant among them (Bolhuis et al. 2014). However, archaeal taxa such as methanogens, DPANN, and Asgard archaea might play a key role in some cases (Potter et al. 2009, Kelley et al. 2012, Wong et al. 2020) while novel archaeal taxa are still discovered in such environments (Kozubal et al. 2013, Wong et al. 2020). Although microbial mats can have this very pronounced layered redox stratification during the day, they can become anoxic during nighttime (Des Marais 1995). Microbial mats function as a consortium where coupling of biogeochemical cycles and processes occurs (Paerl et al. 2000), allowing the products of the metabolism of one group to be available and used by another (Prieto-Barajas et al. 2018b). In addition, the metabolic rates of mat microorganisms are so high that the community production per unit mass competes with that of rainforests (Jørgensen 2001, Krumbein et al. 2003).

Microbial mats can be distinguished in six categories (Bolhuis et al. 2014, Prieto-Barajas et al. 2018b): (i) intertidal or coastal, (ii) hypersaline, (iii) hot spring, (iv) mats in oligotrophic environments, (v) psychrophile, and (vi) acid microbial mats. Intertidal mats are formed on beaches with low slopes and fine sandy sediments (Stal 2012) and they experience strong salinity fluctuations, large temperature changes (Bolhuis et al. 2014) and irregular floods (Prieto-Barajas et al. 2018b). On the other hand, hypersaline microbial mats are found in natural occurring salt lakes and man-made salterns (Bolhuis et al. 2014) and are exposed to salinities up to the crystallization point of halite (Jørgensen 1994), high temperatures and high solar radiation (Bolhuis et al. 2014).

High Throughput Sequencing (HTS) technologies and methods have been widely used to study real-world microbial communities. They have enabled the study of ecosystems with no prior knowledge of the resident species, uncovering unknown and uncultivated strains (Hedlund et al. 2014). Metabarcoding studies are common, well-established, and less computationally demanding than shotgun metagenomics (Bell et al. 2021). However, taxonomic biases may arise from differential efficiency of PCR primer pairing in different species (van der Loos and Nijland 2021) while the short barcoding sequences may limit the resolution. On the other hand, shotgun metagenomic sequencing enables profiling up to the level of strains by obtaining information from random sampling of virtually all genomic regions (Clooney et al. 2016, Segata 2018, Dávila-Ramos et al. 2019). Therefore, microbiome metabolic functions and entire biochemical pathways that occur in a sample can be explored after processing the metagenomic information (Sharpton 2014). Over the recent years, HTS approaches have been used to study the taxonomic and the functional profiles of the microbial communities present in microbial mats (Chen et al. 2020, Wong et al. 2020, Kindler et al. 2022). Several novel high-level taxa have been discovered, e.g. Zixibacterial order GN15 (Wong et al. 2020), and a better understanding on both their adaptive responses in such environments has been established. On top of that, further insight on the mechanisms governing such assemblages has been gained, e.g. the role of photoheterotrophy (Kindler et al. 2022).

The present study was conducted in the Tristomo marsh in the island of Karpathos (Aegean Sea, Greece; Fig. 1A). The study area is included in the Natura 2000 network (site GR4210003) and also in the Greek catalogue of small island wetlands (Y421KAR001; Fig. 1B). It is a seasonal brackish water marsh formed at the edge of a small plain where a seasonal stream ends, characterized as an intertidal marsh (type H) according to the Ramsar convention. Freshwater enters the marsh from the precipitation and drainage basin, while the wetland interacts mainly with the sea through the waves but also underground (WWF Greece 2022). On the coastal front, the cobbled beach is full of litter accumulated by the waves. Due to the close proximity of the marsh with the sea, it occasionally receives saline water, therefore, could be characterized as intertidal; however, evaporation creates hypersaline conditions and crystallized salt forms an upper layer above the actual microbial mat, something, i.e. observed in hypersaline mats (Fig. 1C and D).

Figure 1.

(A) Location of Karpathos island (red pin) in the south east of the Aegean Sea, (B) satellite image of Tristomo bay and the Tristomo marsh (red pin), (C) overview of the marsh in July 2018, (D) overview of the marsh in July 2018, where orange aggregates floating in the water are shown, and (E) a sediment core from 2018 including salt crust, microbial mat, sediment, and an orange aggregate. (Map data: ©2022 Google Earth).

The aim of the present study was to identify the microbial communities in samples from the hypersaline Tristomo marsh, as well as their functional and metabolic capabilities. Our main goals were to test whether the formation of the microbial mat is characterized by seasonality, and which are the mechanisms of osmoadaptation of the microbial communities in this hypersaline environment, as a response to the increased levels of salinity.

Methods

Sample collection

Samples were collected in July 2018 and November 2019 from the Tristomo marsh (Fig. 1). Details on the sample collection are given in Table 1. Sediment samples were collected using cylindrical sampling corers (internal sampling surface 15.90 cm²; Fig. 1E). In the cases where microbial mat layers were clearly observed (July 2018), the top layers were collected separately from the bottom layer. In addition, microbial aggregates observed floating in the marsh were also collected. In the cases where microbial mat layers were not clearly formed (November 2019), there was no slicing during sample collection. In November, samples were collected from three different locations in the marsh, distinguished by the colour of the sediment’s upper layer (black, purple, and orange).

Table 1.

Details on the sample collection.

Sample	Type	Date	16S rRNA accession number	16S rRNA analysis accession number	Shotgun accession number (lane 1)	Shotgun accession number (lane 2)	Analysis accession number (lane 1)	Analysis accession number (lane 2)	Latitude	Longitude
Elos01	Top sediment layer	8/7/2018	ERR9657902	MGYA00607840	ERR6290772	ERR6290778	MGYA00607850	MGYA00607845	35.81936	27.20984
Elos02	Bottom sediment layer	8/7/2018	ERR9657903	MGYA00607829	ERR6290773	ERR6290779	MGYA00607847	MGYA00607843	35.81936	27.20984
Elos03	Orange aggregate	8/7/2018	ERR9657904	MGYA00607839	ERR6290774	ERR6290780	MGYA00607849	MGYA00607852	35.81936	27.20984
Elos04	Top sediment layer	13/7/2018	ERR9657905	MGYA00607828					35.81936	27.20984
Elos05	Bottom sediment layer	13/7/2018	ERR9657906	MGYA00607838					35.81936	27.20984
Elos06	Orange aggregate	13/7/2018	ERR9657907	MGYA00607837					35.81936	27.20984
Elos07	Pink aggregate	13/7/2018	ERR9657908	MGYA00607836	ERR6290775	ERR6290781	MGYA00607848	MGYA00607851	35.81936	27.20984
Elos08	Top sediment layer	9/7/2018	ERR9657909	MGYA00607835					35.81936	27.20984
Elos09	Bottom sediment layer	9/7/2018	ERR9657910	MGYA00607834					35.81936	27.20984
Elos10	Sediment (combined layers, black upper layer)	12/11/2019	ERR9657898	MGYA00607832	ERR6290776	ERR6290782	MGYA00607846	MGYA00607842	35.81962	27.2102
Elos11	Sediment (combined layers, black upper layer)	12/11/2019	ERR9657899	MGYA00607841					35.81962	27.2102
Elos12	Sediment (combined layers, orange upper layer)	12/11/2019	ERR9657900	MGYA00607831	ERR6290777	ERR6290783	MGYA00607939	MGYA00607844	35.81942	27.21007
Elos13	Sediment (combined layers, purple upper layer)	12/11/2019	ERR9657901	MGYA00607830					35.81943	27.20998

Samples were placed in 50 ml falcon tubes (Sarstedt, Nümbrecht, Germany) and were stored at −20°C, until further processing in the laboratory. Upon return to the laboratory, they were used for molecular analysis, i.e. DNA extractions, as well as for the measurement of the Particulate Organic Carbon (POC) and chloroplast pigments concentration [chlorophyll-a, phaeopigments, and chloroplastic pigment equivalents (CPE)]. For the latter, the samples were processed at the Environmental Chemistry Lab of the IMBBC (HCMR), based on standard techniques (Yentsch and Menzel 1963, Hedges and Stern 1984). Water temperature and dissolved oxygen concentration were measured in the water overlaying the sediments by means of a portable multiparameter (WTW Multi 3420 SET G). Salinity was also measured with the portable multiparameter but after dilution of samples with dH₂O since the initial measurement was out of limits (TetraCon® 925 sensor range: 0–70). Sampling was conducted under authorization from the relevant licensing authority (Directorate General for the Protection and Development of Forests and the Rural Environment, Directorate of Forest Management) of the Ministry of Environment and Energy. Additional authorization was also provided from the Management Agency of Dodecanese Protected Areas.

DNA extraction, PCR amplification, and 16S rRNA sequencing

DNA was extracted as in (Henckel et al. 1999) and (Lueders et al. 2004). Approximately, 0.7 g of wet sediment were added to a 2-ml screw-cap vial, prefilled with ∼0.7 g of 0.1 mm (diameter) zirconia/silica beads (11079101z, BioSpec, USA). The vials were filled with 750 μl of 120 mM NaPO4 buffer (pH 8) and 250 μl TNS solution [500 mM Tris-HCl pH 8, 100 mM NaCl, 10% SDS (w/v)] and placed horizontally in a vortex for 10 min at maximum speed. Immediately after that the vials were centrifuged for 10 min at 20 800 rcf and 4°C and the supernatants were transferred to new 2-ml vials. One volume of phenol/chloroform/isoamylalcohol (P/C/I; 25:24:1; pH 8; Carl-Roth, Karlsruhe, Germany) was added to the aqueous supernatant. Vials were vigorously shaken for 20 s and centrifuged for 5 min at 20 800 rcf and 4°C. Supernatants were transferred to new 2-ml vials, and one volume of chloroform/isoamylalcohol (C/I; 24:1; Carl-Roth) was added. Vials were again vigorously shaken for 20 s and then centrifuged for 5 min at 20 800 rcf and 4°C. Supernatants were transferred to new 2-ml vials and C/I extraction was repeated to successfully remove all phenol remnants. Supernatants were transferred to new 2-ml vials and 1.5 ml of polyethylene glycol [30% (w/v) polyethylene glycol 6000 in 1.6 M NaCl] was added to precipitate nucleic acids and the vials were centrifuged for 90 min at 20 800 rcf and 4°C. Supernatants were discarded and the pellets were washed with 1 ml 70% ethanol (4°C) and centrifuged for 30 min. Supernatants were again discarded, pellets were left for air drying (∼5 min) to remove leftover ethanol and resuspended with 50 μl 10 mM Tris.

PCR amplification, library preparation, and MiSeq sequencing was performed as in (Pavloudi et al. 2017). Briefly, PCR amplification was performed following the two-Step PCR approach and targeting the V3–V4 region of the 16S rRNA gene with the primers 341F (5′-CCTACGGGNGGCWGCAG-3′) (Herlemann et al. 2011, Klindworth et al. 2013) and 805RB (5′-GACTACNVGGGTATCTAATCC-3′) revised for detection of SAR11 bacterioplankton (Apprill et al. 2015, Pavloudi et al. 2017).

The PCR negative control sample (blank) was also sequenced, so that possible contamination during the library preparation could be assessed. The raw sequence reads were processed with PEMA (version 2.1.4; Zafeiropoulos et al. 2020) using VSEARCH for the creation of OTUs. Taxonomic assignment was performed with the SILVA database (version 132; Quast et al. 2013). The detailed parameters of the PEMA processing are given in Supplementary File 1. The phyloseq (version 1.36; McMurdie and Holmes 2013), vegan (version 2.5.7; Oksanen et al. 2020), and ggplot2 (version 3.3.5; Wickham et al. 2016) packages were used in R (version 4.1.1; R Core Team 2021) for the creation of barcharts, for the nMDS and PERMANOVA, variation partitioning analysis, db-RDA and mantel test. The scripts of Steinberger (2020) were used for the simper and the Kruskal–Wallis tests.

Shotgun metagenomics sequencing

A total of six samples were selected for shotgun sequencing (Elos01, Elos02, Elos03, Elos07, Elos10, and Elos12), based on the results of the 16S rRNA amplicon sequencing, their potential for the identification of novel unknown lineages and the representation of the majority of collected sample types. Sample preparation was performed using the Nextera^TM DNA Flex Tagmentation and sequencing was done at two lanes of a HiSeq 4000 (2 × 150 bp) at the Norwegian Sequencing Centre (NSC). All the raw sequence files of this study (both 16S rRNA and shotgun metagenomes) were submitted to the European Nucleotide Archive (ENA; Cummins et al. 2022) with the study accession number PRJEB46254 (available at http://www.ebi.ac.uk/ena/data/view/PRJEB46254). In addition, the analysis of both the amplicon as well as of the shotgun metagenomic data was performed by MGnify’s pipeline version 5.0 and is available at https://www.ebi.ac.uk/metagenomics/studies/MGYS00006059 and https://www.ebi.ac.uk/metagenomics/studies/MGYS00006060, respectively (Mitchell et al. 2020).

Assembly and binning

Since the samples were sequenced in two lanes, the fastq files of each sample were concatenated before proceeding with the analyses. Metagenome raw reads were processed with the MetaWRAP workflow (version 1.3.2; Uritskiy et al. 2018). Reads were trimmed and qualified using Trim Galore (version 0.5.0; Krueger 2022), which is a wrapper around Cutadapt (version 1.18; Martin 2011) and FastQC. The clean reads were concatenated, and their coassembly was implemented through the corresponding metaWRAP module, using MEGAHIT v1.1.3; by the term ‘coassembly’, we denote the analysis of the community across the samples, meaning the reads from all metagenomic samples were handled as a single entity. The quality of the coassembly was evaluated using QUAST (Gurevich et al. 2013). Binning was then performed using the clean reads and the coassembly. The metaWRAP module for binning was performed using MetaBAT 2 (version 2.12.1; Kang et al. 2019: 2) and MaxBin 2 (version 2.2.6; Wu et al. 2016). CheckM (version 1.0.12; Parks et al. 2015) was used by the metaWRAP module to assess the quality of the bins produced by MetaBAT 2 and MaxBin 2. Bins were then consolidated and refined using Binning_refiner (Song and Thomas 2017) as wrapped in the Bin_refinement module of metaWRAP. The Bin_refinement module was invoked with the default values for minimum completion (70%) and maximum contamination (5%). The consolidated bins set was further improved using the reassemble_bins module of metaWRAP. To this end, bwa (version 0.7.17-r1188; Li and Durbin 2009), spades (version v3.13.0; Nurk et al. 2017), and CheckM were used. To estimate bins’ abundances in each sample (in genome copies per million reads), the corresponding metaWRAP module was performed invoking Salmon (version 0.13.1; Patro et al. 2017). The refined bin-set was also used for the blobology module of metaWRAP; taxonomic annotation of the coassembled contigs was performed using megaBLAST and the nt database of NCBI.

The coassembled contigs and the refined bins set were then used as input to Anvi’o (version 7.1; Eren et al. 2015). Bowtie 2 (version 2.3.5; Langmead and Salzberg 2012) was used to build BAM files and mapping and Prodigal (version 2.6.3; Hyatt et al. 2010) for gene prediction. BAM files were also made from the clean reads of each sample. A contigs database was built (using the anvi-gen-contigs-database program) after converting the contigs name as Anvi’o suggests (see contigs-per-bin.sh script) and it was decorated with hits from HMM models (anvi-run-hmms). An anvi profile was then built for each of the samples’ bam file (anvi-profile) and they were merged (anvi-merge) into a single profile. The refined bins along with their corresponding renamed contigs were imported as a collection in the merged profile database (anvi-import-collection). At this point, a first Anvi’o summary was recovered (anvi-summarize). Bins with a redundancy >10% were manually refined and a second summary of the bins set was created; ‘redundancy’ here is referring to the anvi’o measure of ‘how many copies of each single-copy core gene is found within a genome’.

Afterwards, a simper analysis was performed to identify MAGs that could significantly differentiate the sample categories (i.e. mat, aggregates, and sediment).

Taxonomic composition

Based on the returned coassembly from METAWRAP and the clean reads, communities’ taxonomic composition was assessed using Kraken2 (Wood et al. 2019) and the standard Kraken 2 database (NCBI: January 2022). GTDB-Tk (version 1.7.0; Chaumeil et al. 2020) was used to classify genomes with the Genome Taxonomy Database (GTDB, version r202; Parks et al. 2022). GTDB-Tk made use of pplacer (version 1.1.alpha19-0-g807f6f3; Matsen et al. 2010), and FastANI (version 1.32; Jain et al. 2018).

Functional annotation

Functions were predicted at two levels: both at the MAG level, as well as at the sample level. For the functional annotation at the MAG level, using the anvio’ contigs database and the anvi-run-kegg-kofams program, the anvio’ contigs database was annotated with HMM hits from KOfam, a database of KEGG Orthologs (KOs). Likewise, using the anvi-run-ncbi-cogs, NCBI’s Clusters of Orthologous Groups (COGs) based annotations were added. The MAGs that correspond to the refined bins as they were retrieved after the metaWRAP and the anvio refinement steps, were annotated with KEGG modules; manually defined functional units of gene and reaction sets (Kanehisa et al. 2012). MAGs were ‘translated’ to an anvio collection (i.e. a virtual construct storing bins of items in an anvi’o profile database) and this collection was used along with the anvi-estimate-metabolism program to determine which enzymes are present in each MAG and compute the completeness of each metabolic module (scripts can be found under the anvio folder). An nMDS was constructed based on the presence/absence of modules in the MAGs using the jaccard similarity index.

For the functional annotation at the sample level, the clean reads as they were returned by the corresponding metaWRAP module and the DiTing tool (Xue et al. 2021) were used to estimate the contribution of each sample to the biogeochemical cycles incorporated in DiTing. DiTing used MEGAHIT (Li et al. 2015) to build the assembly of each sample separately (so, the coassembly described in the ‘Assembly and binning’ section was not used for this step) and Prodigal (Hyatt et al. 2010) to retrieve the Open Reading Frames (ORFs). KofamScan (Aramaki et al. 2020) was used for the annotation of the ORFs using KEGG ORTHOLOGY terms. The relative abundances of metabolic and biogeochemical functional pathways in each sample were then determined by DiTing.

Extracting molecular functions indicating osmoadaptation

Based on a review on strategies of adaptation of microorganisms to high salt concentrations (Gunde-Cimerman et al. 2018), a list of compatible solutes in halotolerant/halophilic microorganisms, alkali metal–cation transporters and putative osmosensors was retrieved. Subsequently, a list of molecular function (KO) terms related to the processes responsible for all the above was created, available on Github, and their presence in the MAGs was investigated.

In addition, based on a review on photosynthetic light harvesting (Croce and van Amerongen 2014), a list of light harvesting systems was created, and their related KO terms were retrieved. Afterwards, the samples where those systems were present were identified.

MAGs reference phylogenies

Bacteria_76 and Archaea_71 are sets of single copy genes (SCG) that are shared among all bacteria and archaea (Eren et al. 2015). Their intersection leads to a set of 25 SCG; i.e. the maximum number of SCG that 2 taxa potentially share. This set was used to build the phylogenetic tree of the reconstructed MAGs. The anvi-get-sequences-for-hmm-hits program of Anvi’o was used to extract and align the amino acid sequences of each of these genes from all the MAGs independently. This Anvi’o program makes use of MUSCLE (version 3.8.1551; Edgar 2004) to return an alignment of the extracted sequences. Once all the amino acid sequence alignments were extracted, they were trimmed using Clipkit (version 1.1.5; Steenwyk et al. 2020). A super matrix was then built using the SCGs of the intersection. In cases where a MAG lacked a gene (due to sequencing limitations, assembly and/or binning related challenges and so on), alignment gaps were filled with dashes; the initial and the trimmed per gene alignments, the final super matrix alignment and a table showing exactly which SCGs were present in which MAG are available on the project’s GitHub repository under the SCG folder. Using IQ-TREE2 (Hoang et al. 2018, Minh et al. 2020) the phylogeny of the reconstructed MAGs was built using 1000 bootstrap replicates (-B 1000) and 1000 bootstrap replicates for Shimodaira–Hasegawa-like approximate likelihood ratio test (SH-aLRT; -alrt 1000). The best-fit model (LG+R10) was retrieved using ModelFinder (Kalyaanamoorthy et al. 2017). Using Barrnap (Seemann 2014), the 16S rRNA gene was extracted from the retrieved MAGs. The phylogeny of the MAGs and their relative abundances were integrated and visualized using GraPhlAn (Asnicar et al. 2015). All bioinformatics analyses were supported by the IMBBC High Performance Computing system (Zafeiropoulos et al. 2021).

Results

Taxonomic composition from 16S rRNA amplicon analysis

The DNA concentrations of the samples are shown in Table S1 (Supporting Information) and the results of the processing of the sequences are shown in Table S2 (Supporting Information). Sequencing of samples Elos08 and Elos11 was not successful and, therefore, they were not included in the following analyses. The final number of OTUs, after removal of the OTUs that were also found on the blank sample, was 2689 (see finalTable_noblanks.tsv). Overall, the most abundant phyla, as assessed by the relative abundance percentages of each replicate sample averaged per sampling station, were Bacteroidetes (∼17%), Euryarchaeota (∼16%), Proteobacteria (∼15%), and Halanaerobiaeota (∼10%, class Halanaerobiia; Fig. 2). Among the Bacteroidetes, the most abundant class was Bacteroidia (∼14%), followed by Rhodothermia (∼3%). Bacteroidetes had higher abundances in the aggregate samples (∼25% on average in Elos03, Elos06, and Elos07) and lowest in the mat samples (∼11% on average in Elos01 and Elos04). Euryarchaeota had very low abundances in samples Elos09, Elos10, and Elos13 and the higher abundances in the microbial mat samples (∼28% on average in Elos01 and Elos04); along with the Halanaerobiaeota (∼24% on average in Elos01 and Elos04) they were the dominant phyla of the mat samples. Halanaerobiaeota had low abundances in the orange and pink aggregates as well as in Elos09, Elos10, and Elos13. Among the Euryarchaeota, the most abundant class was Halobacteria (∼14%) followed by Thermoplasmata (∼2%).

Figure 2.

Balloon plot showing the abundances of the main microbial taxa, at the phylum level, at each sample, based on the 16S rRNA amplicon sequencing. Archaea and bacteria represented OTUs that could not reach higher taxonomic resolution.

Proteobacteria were almost equally distributed among all the three sample categories, i.e. mat, aggregates, and sediments (∼15%, ∼21%, and ∼13%, respectively), as well as among the classes Alphaproteobacteria (∼7%), Gammaproteobacteria (∼4%), and Deltaproteobacteria (∼4%; Figure S1, Supporting Information). Although Patescibacteria had a low average abundance (∼6%), they were dominant in Elos09 (∼46%). In addition, Chloroflexi were almost absent from the top layers and the aggregates, but they were found in the bottom sediment layer and in the combined sediment samples. Cyanobacteria were about ∼1% on average of all the samples.

The nMDS of the microbial OTUs (Figure S2, Supporting Information) showed that their spatial pattern differs both by their type and the year of sampling, which was also confirmed by the PERMANOVA results (Type: F.Model = 2.0396, P < .05; Year: F.Model = 2.3098, P < .05).

Coassembly, binning, and taxonomic composition from shotgun metagenomics analysis

Shotgun metagenomic sequencing of the chosen six samples resulted in 744 million reads totalling 112.3 Gbp, with each sample ranging between 16.77 and 21.78 Gbp. Coassembling of all the samples resulted in 1,5 million contigs totalling 5.04 Gbp (see assembly_report.html). The per-sample assemblies returned a total of 11.2 million contigs with a sum of 10.15 Gbp. Number of reads per sample, before and after the quality control, their length and the corresponding number of contigs are shown in Table S3 (Supporting Information).

Krona plots of the community profiles can be viewed through the kronagram.html. Based on the taxonomic profiles retrieved from Kraken2 (Figure S3, Supporting Information), after removing sequences belonging to Viruses and sequences that could not be classified (∼1%), Euryarchaeota (class Halobacteria) represent the majority of the total archaeal taxa (∼30% on average); however, they were almost absent from sample Elos10 (abundance ∼1%) while they were dominant in sample Elos01 (∼59%). As far as bacterial taxa are concerned, the most abundant ones were Alphaproteobacteria (∼19%), followed by Actinobacteria (∼13%), and Gammaproteobacteria (∼10%). Betaproteobacteria, delta/epsilon Proteobacteria subdivisions, and Bacteroidetes/Chlorobi group had similar abundances (∼5%, 4%, and 5%, respectively) with the latter having the highest abundance in sample Elos12 (∼10%). Cyanobacteria were limited in all samples (∼2% on average). Also, Kraken2 analysis did not identify any Nanoarchaeota; however, it identified in sample Elos01 the other archaeal taxa that are their hosts and namely (a) I. hospitalis, (b) Acidilobus sp. 7A, (c) Vulcanisaeta spp., (d) Pyrobaculum spp., (e) Metallosphaera spp., (f) Caldivirga sp., and (g) Sulfolobus sp.

Krona plots with the taxonomic profiles of each sample are available on GitHub. Prodigal predicted millions of genes per sample ranging from 2.1 (Elos03) to 4.3 million (Elos10). Their metabolic capacity/potential is further described in the ‘Biogeochemical cycles’ section. Based on the blobology results, among the 1 513 505 coassembled contigs a set of 102 250 were binned (see Figure S4, Supporting Information); according to megaBlast and the nt database of NCBI, among the binned contigs 53 536 were bacterial and 2230 archaeal while 60 contigs were assigned as viral and 739 as eukaryotic. The corresponding numbers for the case of the unbinned contigs were 1–2 orders of magnitude higher; thus, the number of unbinned contigs were 430 028 bacterial, 126 690 archaeal, 15 828 eukaryotic, and 1805 viral, respectively (Figure S5, Supporting Information).

MAGs phylogeny, functional annotation, and distribution across samples

In line with the quality definitions described in (Bowers et al. 2017), metagenome binning generated a total of 289 MAGs; details are shown in the Supplementary File 2. According to the CheckM software 194 MAGs were reconstructed with a completeness higher than 90% and a contamination lower than 5% and all the rest had a completeness > 70% and a contamination score < 6%. According to the anvio summary (using the coassembly as contigs database, the merged samples as profile and the reconstructed MAGs, i.e. the refined bins, as a collection) the redundancy of 10 (bin127, bin114, bin156, bin243, bin268, bin276, bin12, bin269, bin252, and bin226) of the reconstructed MAGs was >10% (see binning_results.png). After the manual refinement of these 10 MAGs, a total of: (i) 178 bacterial high quality (completeness > 90%, contamination <5%), (ii) 70 bacterial and 39 archaeal medium quality (50% < completeness < 90% and 5% < contamination < 10%), and (iii) two bacterial MAGs of low quality (bin263 and bin182 with a completeness score < 50%) were retrieved (see binning_reassembled.png). Combining the anvio summary results (see bin_by_bin folder in SECOND_SUMMARY) and the Barrnap outcome (see arc_rrnas and bac_rrnas on GitHub repo), the 16S rRNA gene was identified in 100 out of the total 250 bacterial MAGs. Likewise, from a total of 39 archaeal MAGs, 16S rRNA gene was found in 28 of them. Contigs included on those MAGs represented 1.03 Gbp of assembled reads. A set of 25 MAGs had a completeness of 100% and contamination less than 5% while 5 bacterial (MAG 143, MAG 66, MAG 129, MAG 189, and MAG 76) and 1 archaeal (MAG 232) MAGs among them had a contamination of 0%. Overall, bacterial MAGs had higher completeness scores.

Phylogenomic placement

The GTDB-Tk returned phylogenetic trees of the GTDB partition and the MAGs assigned to the corresponding domain for the cases of bacteria and archaea including the two low quality included (bin_182 : Proteobacteria and bin_263: Verrucomicrobiota). The phylogeny of the reconstructed MAGs (Fig. 3) was built based on single-copy genes present on both archaea and bacteria, using the total number of MAGs even if some of the MAGs did not have all the 25 single-copy genes. Although not all these 25 single-copy genes were found in every MAG, still, the mean number of occurrences of a gene among the 289 MAGs was 266.68, ranging from 211 to 278. In general, the archaeal MAGs had the fewer single-copy genes, most probably due to their lower completeness. Using the total number of MAGs, the phylogeny of the reconstructed MAGs highlights the robustness of the method as even for those MAGs the phylogenetic signal was enough to place them among the representatives of their phylum. Thorough investigation of the tree pointed out that the only two discrepancies were that the sole MAG of the RBG-13–61-14 phylum (bin_124) that was placed among the Myxococcota representatives and that a representative of the Patescibacteria phylum (bin_61) was not placed close to the rest of Patescibacteria but as the closest relative of the representatives of the Chloroflexota phylum.

Figure 3.

Concatenated marker gene phylogeny of the Karpathos’ marsh MAGs. Phylogeny of the 289 MAGs recovered from a hypersaline marsh in Karpathos island, based on 25 concatenated, single-copy genes present both in archaea and bacteria. From inside to outside, the concentric circles around the phylogeny indicate: the bin id, phylum level taxonomy, bin relative abundance in the (i) top layer of the sediment (i.e. the microbial mat), (ii) bottom layer of the sediment, (iii) orange aggregate, (iv) pink aggregate, (v) combined sediment layers (black upper layer), and (vi) combined sediment layers (orange upper layer). Stars indicate high quality bins. Branch colours stand for bootstrap values: (i) > 90 (green), (ii) 70–90 (black), and (iii) < 70 (red). The novel phylum (bin_202) is highlighted in grey.

The novel candidate phylum (bin_202) was placed within the same clade with Eisenbacteria (bin_31). In general, bootstrap values were > 90% with only exceptions a number of clades with representatives of the Nanoarchaeota phylum (of which the completeness and the number of single-copy genes present was relatively lower).

Distribution of MAGs across samples

Based on the classification of the retrieved MAGs [Figure S6, Supporting Information; MAG abundances per sample (metaWRAP): link], the most abundant phylum was Bacteroidota (∼28% on average), which almost dominated sample Elos03 (∼57%) and Elos07 (∼40%). The second most abundant phylum was Proteobacteria (∼13% on average), with abundances ranging from ∼2% in Elos02 to ∼23% in Elos01 and ∼22% in Elos07. Planctomycetota and Desulfobacterota were found at about ∼8% and ∼7%, respectively, with the latter being absent from Elos03 and very rare in Elos07 (∼2%). The only phylum that was present only in the microbial aggregates, i.e. in Elos03 and Elos07, and was absent from all the other samples was Myxococcota.

The most abundant archaeal phylum was Nanoarchaeota (∼5% on average), which was mostly found in Elos01 (∼16%) and in much lower abundances in the other samples. Thermoplasmatota and Asgardarchaeota were found in similar abundances (∼3% and ∼2% on average, respectively) and they were also absent from Elos03 and Elos07. Halobacteriota (∼2% on average) were not found in Elos10 and Elos12 and were mostly present in Elos01 (∼6%).

Elos10 was the sample with the highest number of bins (Figure S7, Supporting Information) even if it was the one with the lowest number of reads. In addition, it seems to be closer to Elos02, the bottom layer sediment sample from July, and to sample Elos12. The microbial aggregates (Elos03 and Elos 07) form another cluster, distinct from the other samples, but closer to Elos01, the microbial mat sample.

The MAG 202 that represents a novel candidate phylum was present in samples Elos12, Elos10, and Elos02.

Based on the results of the simper analysis, there were no MAGs that could significantly contribute to the differentiation of the sample categories (i.e. mat, aggregates, and sediment). The only comparison where certain MAGs were found to contribute, was the one between aggregates and combined sediment samples (Table S4, Supporting Information).

Functional annotation of MAGs

The reconstructed MAGs were annotated with KofamScan with a range of KEGG ORTHOLOGY terms ranging from 354 to 2879 terms (Figure S8, Supporting Information), leading to 1–87 complete KEGG modules (Figure S9, Supporting Information). The archaeal MAGs had, in general, lower completeness scores, lower number of KO terms assigned, and less complete modules.

As it is shown in Figure S10 (Supporting Information), MAGs form distinct clusters both based on the taxonomy, i.e. if they are bacterial or archaeal (PERMANOVA: F.Model = 35.767, P< .001), as well as based on their completeness (PERMANOVA: F.Model = 5.2156, P < .001). The modules that contribute most to this clustering, as identified by the simper analysis, as well as the significance of any given module’s contribution, are shown in Table S5 (Supporting Information). Examples of these modules were related to oxygenic photosynthesis and nitrogen, sulphur and carbon cycles. When examined separately, again they differ significantly by completeness (PERMANOVA: Bacteria: F.Model = 3.4053, P < .001; archaea: F.Model = 2.4452, P < .001).

Comparison of taxonomies between amplicon and metagenomic analysis

As expected, the taxonomic composition of the microbial communities in our samples depends on the analytical procedure that was followed in each case. However, when the similarity matrices of the samples were compared, the pattern deriving from the relative abundances of the microbial OTUs as derived from the amplicon survey, was highly correlated with the one deriving from the classification of the shotgun metagenomic bins (Mantel statistic: r = 0.83, P< .001). The pattern deriving from the Kraken2 analysis was also correlated both with the amplicon survey, as well as with the classification of the retrieved bins (Mantel statistic: r = 0.45, P < .05; r = 0.52, P< .05, respectively), but on a lesser degree.

Physicochemical analysis

The physicochemical variables are given in Table 2. According to the variation partitioning analysis (Table 3), the combination of oxygen and temperature explained 64% of the total variation in the Kraken2 community similarity matrix. For the other cases, i.e the amplicon data matrix and the classification of the retrieved MAGs, residuals were higher than 0.60 and, therefore, no significant models were retrieved.

Table 2.

The values of the physicochemical variables of the samples. O₂, oxygen concentration (mg/lt); Chl-a, chlorophyll-a concentration in the sediment (ug/g); CPE, chloroplastic pigment equivalents (ug/g); and TOC, total organic carbon (ug/g).

Sample	O2	Temperature (°C)	Salinity (psu)	TOC	Chl-a	Phaeopigments (ug/g)	CPE
Elos01	9.1*	32.1*	169.8*	59715.1	55.99	20.69	76.68
Elos02				72110.36	92.54	65	157.54
Elos04				54933.96	29.64	14.62	44.26
Elos05				67992.39	72.91	9.66	82.57
Elos08				44755.04	35.11	27.02	62.14
Elos09				58128.71	33.34	41.53	74.87
Elos10	14.74	24.2	63.9	61256.23	73.34	24.14	97.48
Elos11	10.7	24.3	63.9	66970.84	7.312	9.014	16.33
Elos12	9.72	24.3	166.5	45 000	16.04	56.74	72.78
Elos13	13.33	24.7	169.2	55528.85	48.63	30.03	78.66

*

: these variables were estimated on 8 July 2018, but for the variation partitioning analysis and db-RDA an assumption was made and they were used also for samples Elos02–Elos09.

Table 3.

The percentage of variation explained of each explanatory physicochemical variable, as well as their combinations. *: < .1, **: < .05.

	Adjusted R2	P
Oxygen + temperature	0.64	*
Oxygen with temperature as condition variable	0.46	**
Temperature with oxygen as condition variable	0.75	nonsign.

Functional profiles at the sample level

A set of 783 693 unique proteins were predicted from a total of 3 532 725 hits with NCBI COGs. The MEGAHIT assembler as implemented in the framework of DiTing, returned the assembly of each sample ranging from 1 323 538 (Elos03) to 2 773 933 (Elos10) contigs.

The relative abundances of metabolic and biogeochemical functional pathways in each sample are available in the DiTing folder on the GitHub repository. As shown in Fig. 4 and Figure S11 (Supporting Information), pathways belonging to the carbon cycle, central metabolism, and other metabolism were the most abundant (∼23%, ∼19%, and ∼18% of the total pathways on average, respectively). More specifically, the Reductive citrate cycle and the Dicarboxylate–hydroxybutyrate cycle dominate the carbon cycle (∼8% and ∼5%, ∼18%, respectively). In the central metabolism, pathways like the Embden–Meyerhof glycolysis pathway and tricarboxylic acid (TCA) cycle were most abundant (∼7% each). In contrast, the Entner–Doudoroff pathway, i.e. an alternative glycolytic pathway was found in abundances an order of magnitude lower than the Embden–Meyerhof glycolysis pathway (Figure S12, Supporting Information). Wood–Ljungdahl pathway that enables the use of hydrogen as an electron donor, was mostly found in bottom sediment layer sample (Elos02), but also in the combined sediments from the November 2019 sampling (Elos10 and Elos12) and to a lesser extent in the other samples. Bacterial chemotaxis, flagellar assembly, and dissimilatory arsenic reduction are also among the most abundant pathways (∼9% and∼7%, ∼3%, respectively).

Figure 4.

Balloon plot showing the abundances of the metabolic pathways for six biogeochemical cycles, at each sample, as retrieved from DiTing.

Regarding the methane cycle, methanogenesis pathways were found in the bottom sediment layer sample (Elos02), but also in the combined sediments from the November 2019 sampling (Elos10 and Elos12), as the Wood–Ljungdahl pathway; the most abundant methanogenesis pathway was the formation of methane from acetate. Interestingly, methane oxidation was almost absent in the samples.

Regarding the sulphur cycle, the most abundant pathways were the assimilatory and dissimilatory reduction of sulphate to sulphite (∼1% each). As shown in Fig. 5, thiosulphate oxidation as well as sulphite oxidation, but to a lesser extent, contribute to the sulphate pool. Sulphur disproportionation to sulphide and sulphite was absent in the aggregate samples, i.e. Elos03 and Elos07. In addition, although DMSO reduction was an abundant pathway in all the samples, DMS oxidation was very rare and mostly found in Elos03, i.e. the orange aggregate (Figure S13, Supporting Information).

Figure 5.

Relative abundances of the pathways involved in the sulphur cycle. The pie chart indicates the relative abundance of each pathway in each sample. The size of pie charts represents the total relative abundance of each pathway. ASR: assimilatory sulphate reduction; DSR: dissimilatory sulphate reduction.

As shown in Fig. 6, dissimilatory nitrate reduction to nitrite and nitrite to ammonia (DNRA) was prevalent in all samples, but it was mostly found in the combined sediment samples (Elos10 and Elos12). Denitrification (i.e. nitrite to nitric oxide and nitric oxide to nitrous oxide) was mostly found in sample Elos12. Nitrification, i.e. conversion of hydroxylamine to nitrite was almost absent from samples Elos01 (the microbial mat) and the microbial aggregates (Elos03 and Elos07); no amoABC genes were identified and hao genes were found only in bacterial MAGs. Nitrogen fixation was mostly abundant in the sample Elos07 (the pink aggregate).

Figure 6.

Relative abundances of the pathways involved in the nitrogen cycle. The pie chart indicates the relative abundance of each pathway in each sample. The size of pie charts represents the total relative abundance of each pathway. ANRA: assimilatory nitrate reduction to ammonium; DNRA: dissimilatory nitrate reduction to ammonium; and Anammox: anaerobic ammonium oxidation.

Anaplerotic genes were also very abundant in our samples (∼2%), and in particular the Pyruvate Carboxylase Pathway, which produces oxaloacetate from pyruvate and replenishes the intermediates of the TCA cycle.

Finally, KO terms that were not added to one of the DiTing pathways regarding, K04641 (bacteriorhodopsin) and K04642 (halorhodopsin) were more abundant in the microbial mat sample, as well as in the microbial aggregates, and absent from Elos10. Sensory rhodopsin (K04643) was found in all samples, but in much higher abundances in Elos01, Elos03, and Elos07.

Osmoadaptation and light harvesting in the reconstructed MAGs

A total of 215 KOs related to osmoadaptation were identified; out of those, 61 were found in the MAGs and 59 were found in the samples (see GitHub repo). Each MAG had at least one of those KOs, with some MAGs having several osmoadaptation KOs (Figure S14, Supporting Information). As shown in Figure S15 (Supporting Information), it seems that the preferred osmoadaptation strategy in all of the samples was related to the presence of alkali metal–cation transporters, while strategies related to uptake or production of compatible solutes were less abundant.

Regarding the light harvesting pigments/systems, we examined LH2 complex (‘peripheral light-harvesting complex’), chlorosomes, phycobilisome, and FMO protein (Figure S16, Supporting Information). The K08930 term, corresponding to light-harvesting protein B-800–850 alpha chain, was present in all the samples. It was highly abundant in the mat and the aggregates (Elos01, Elos03, and Elos07) and much less in the sediment samples. Proteins related to phycobilisomes (K02094, K02096, K02097, and K02290) were found in all of our samples and were most abundant in Elos 03 and Elos10.

KOs related to chlorosomes (various proteins of the chlorosome envelope), ranging from K08945 to K08954, were not detected in our samples. Similarly, the K08944 term, which corresponds to FMO protein, was absent from all the samples.

Discussion

Community composition and functional potential

Overall, there seems to be an agreement in all three ways that taxonomic composition has been derived for our samples. Even though 6 out of the 11 samples were sequenced using shotgun metagenomics, still we can compare them with the amplicon data and derive conclusions. In fact, the amplicon survey produced results that were highly correlated to the results of the classification of the retrieved bins; such congruence between the two methodologies has been reported in recent studies (Chan et al. 2015, Regalado et al. 2020) although there is also evidence for the opposite (Tessler et al. 2017). The results from Kraken2 were also correlated to the amplicon and the results of the classification of the retrieved MAGs, although the correlation was not as strong. However, Kraken2 is not restricted in the short amplicon length (Johnson et al. 2019) and, therefore, it was able to reach species level resolution (Lu and Salzberg 2020).

The top sediment layers, i.e. the microbial mat, was characterized by the presence of Halobacteria/Halobacteriota, Halanaerobiia, and Nanoarchaeota, which are all halophilic and were, therefore, able to withstand the salt crust that was on top of the mat (Norton et al. 1993, Casanueva et al. 2008, Çınar and Mutlu 2020, Akpolat et al. 2021). They were almost absent in the deeper sediment layers and they were absent in the winter, when there was no obvious microbial mat formed, despite the high salinity in sample Elos12. Therefore, the limiting factor for their presence could have been the lower temperature during the winter season since the optimum temperature for representatives of the Halobacteria is higher than 35°C (Grant 2015). Halobacteria, and their bacteriorhodopsin, were most likely responsible for the purple layer in the microbial mat, since they have been shown to cause striking red colours in salt flats (Stoeckenius et al. 1979). Another organism responsible for the red colouration could have been the algae Dunaliella salina (Oren 2016, Naghoni et al. 2017); however, Dunaliella sp. was very rare in our shotgun metagenomic data as analyzed by MGnify, therefore, the red colouration cannot be attributed to algae. Nanoarchaeota are obligate symbionts and they have found in association with Crenarchaeota/Thermoproteota (Huber et al. 2002, Podar et al. 2013, Munson-McGee et al. 2015, Wurch et al. 2016, Merkel et al. 2017). Out of the known host–symbiont pairs and the putative hosts that have been proposed, we identified two of the known hosts (Ignicoccus hospitalis and Acidilobus sp. 7A) and five of the putative hosts (Jarett et al. 2018). However, the Nanoarchaeotal representatives were assigned at taxonomies higher than the species level. Thus, we can only presume that the symbionts Nanoarchaeum equitans and Nanopusillus acidilobi of the aforementioned hosts were present in our data.

It was hypothesized that the microbial aggregates would be more closely related to the microbial mat samples and this was confirmed by our results. However, the microbial aggregates were mostly characterized by the presence of Bacteroidetes/Bacteroidota and Myxococcota and the absence, or very low abundance, of Desulfobacterota, Thermoplasmatota, and Asgardarchaeota.

In the photic zone of the microbial mats, phototrophic microorganisms transduce light into energy using either pigments (chlorophyll and/or bacteriochlorophylls) or retinal-based rhodopsins (Kurth et al. 2021). In hypersaline microbial mats, phototrophy is possible because light can penetrate salt crusts and, therefore, the only limiting factor is the capability of the microbial communities to actually perform phototrophy under salt-saturated conditions (Meier et al. 2021). Out of the phyla that have been reported to perform (bacterio)chlorophyll-based phototrophy, i.e. Cyanobacteria, Proteobacteria, Chlorobi/Chlorobia, Chloroflexi/Chloroflexota, Firmicutes, Acidobacteria/Acidobacteriota, Eremiobacterota and Gemmatimonadetes/Gemmatimonadota (Zeng and Koblížek 2017, Zheng et al. 2022), only Proteobacteria and Chloroflexi/Chloroflexota were found in high abundances. Cyanobacteria were very rare in our samples, which can be attributed to the high salinity of the marsh (DiLoreto et al. 2019), which can lead to osmotic stress and inhibition of photosynthesis (Sudhir and Murthy 2004). So, Cyanobacteria were not the foundation of the microbial mats in our study, as has been shown in other examples of hypersaline microbial mats (Bolhuis et al. 2014, Wong et al. 2016), which was expected as oxygenic photosynthesis is completely inhibited at saturation-level salinities (40%; Meier et al. 2021). However, cyanobacterial genera that are characterized by strategies and survival mechanisms that allow them to grow in such high salinities (Oren 2015), such as Euhalothece and Halothece, were present in our samples. Chloroflexi/Chloroflexota were most present in the bottom sediment layer and in the combined sediment samples. Representatives that were identified from our samples are Chloroflexus aurantiacus (Pierson and Castenholz 1974) and C. aggregans (Hanada et al. 1995); they can grow photoautotrophically and photoheterotrophically under anaerobic conditions but also chemotrophically under aerobic conditions, using sulphide or hydrogen as an electron donor (Tang et al. 2011, Kawai et al. 2019, 2021). Given that it is unclear if they can harvest light deeper in the sediment and considering the low abundance of Cyanobacteria, it is most probable that Chloroflexi/Chloroflexota grow chemotrophically in our samples. According to our data, and since it is shown that light harvesting proteins are mainly in the mat and aggregate samples and less in the deeper sediment samples, Proteobacteria seem to be taking over the photosynthetic pathways. More specifically, some of the species that have been found to perform photosynthesis and are present in our samples are (from the Alphaproteobacteria class) Citromicrobium sp. (Jiao et al. 2010), Rhodopseudomonas palustris, Rhodobacter sp., Rhodospirillum rubrum, Roseobacter denitrificans, Bradyrhizobium sp., Roseobacter sp. (Larimer et al. 2004, Bryant and Frigaard 2006) and (from the Gammaproteobacteria class) Marichromatium purpuratum (Shiung et al. 2018), Congregibacter litoralis (Fuchs et al. 2007, Spring et al. 2009), Allochromatium spp. (Kyndt and Meyer 2020). Overall, based on their taxonomic composition, the microbial mats under study seem to resemble microbial mats from an irregularly inundated tidal flat in Oman (Meier et al. 2021).

Regarding functional annotation of MAGs and their clustering according to taxonomy, it seems that it is driven by modules that are present in archaea and absent in bacteria, and vice versa. For example, regarding cysteine biosynthesis (M00021), this pathway is still unexplored in archaea and although it has been found in certain species, e.g. Methanosarcina barkeri, it is suggested that there might be a different cysteine biosynthesis pathway in archaea (Kitabatake et al. 2000). Likewise, there are no archaeal homologs for the bacterial pantothenate biosynthetic genes (Ronconi et al. 2008), therefore pantothenate biosynthesis (M00913) is one of the pathways that contributes to the differentiation between bacterial and archaeal MAGs. In addition, acetogen (M00618) is only found in bacteria.

Seasonality

Our second hypothesis was that samples from the deeper sediment layers collected in the summer would be more similar to the combined sediment samples collected in the winter. This was also confirmed by our results, as these samples were clustered together. However, despite the similarities between the samples, there were also certain differences; as expected, there is a high degree of spatial variation in the marsh under study (Dillon et al. 2009). Overall, in the deeper sediment layers or where microbial mat was not formed, acetogens such as Moorella thermoacetica, Clostridium aceticum, and Acetobacterium woodii (Schuchmann and Müller 2016), might have been utilizing the Wood–Ljungdahl pathway, i.e. using H₂ hydrogen as an electron donor and CO₂ as an electron acceptor. Acetogens are either competing directly with hydrogenotrophic methanogenic archaea or interacting syntrophically with acetotrophic methanogens (Ragsdale and Pierce 2008). Utilization of acetate for methanogenesis is present in the genera Methanosarcina and Methanothrix (Ferry 1992), which were both found in our samples. However, while methanogens were abundant in our samples, they might have been contributing little to anaerobic mineralization, since in salinities of 180‰ or less, they are inhibited by the increased activity of SRB (Sørensen et al. 2004). In addition, SRB might have also been using the Wood–Ljungdahl pathway in reverse (Ragsdale and Pierce 2008). It seems that other hypersaline microbial mats, such as Guerrero Negro and Shark Bay, are characterized by a high number of genes related to sulphur metabolism (Wong et al. 2016) and hypersalinity in coastal pans has been shown not to inhibit the activity of SRB (Porter et al. 2007). Therefore, the high salinity in our samples is not expected to inhibit sulphur cycling. Sulphate-reducing microorganisms (SRM) from the Euryarchaeota lineage (Muyzer and Stams 2008) were very abundant in the microbial mat and in the aggregates, while SRM belonging to Deltaproteobacteria (Muyzer and Stams 2008) were present in the deeper and the combined sediment samples. However, the occurrence of SRMs is not synonymous to the occurrence of sulphate reduction in the given habitat (Muyzer and Stams 2008).

It has been proposed that arsenic and sulphur cycling can sustain high microbial metabolic rates in permanently anoxic mats (Visscher et al. 2020). Bacteria capable of performing anoxygenic photosynthesis using arsenite (As(III)) as an electron donor, such as Ectothiorhodospira sp. and Halorhodospira halophila (Hoeft McCann et al. 2017), were present in our samples. It is suggested that might be performing oxidation of As(III) to arsenate (As(V)), which is afterwards reduced back to As(III) (Hoeft et al. 2004), which could explain the high occurrence of dissimilatory arsenic reduction in our samples. Arsenate can be an important electron acceptor in the biogeochemical cycling of carbon (Oremland et al. 2000); thus, arsenate reduction has a great potential to precipitate carbonates and it is energetically better than sulphate reduction (Visscher et al. 2020), although there is high potential for the latter to be also occurring our samples.

Thaumarchaeota, which are involved in nitrification in marine ecosystems (Veuger et al. 2013), were also present in our samples, according to the Kraken2 results, but in very low abundances, in contrast to other hypersaline microbial mats (Ruvindy et al. 2016). On the other hand, common nitrifying bacteria such as Nitrobacter where abundant in our samples; however, since the potential for nitrification was mostly present in the combined sediments and in the deeper layer, there seems to be a certain degree of hypersalinity limitation on the growth of nitrifying bacteria, as has been previously suggested (Jeffries et al. 2012). Interestingly, no amoABC genes were identified and, as expected, hao genes were found only in bacterial MAGs (Holmes et al. 2019); it has been suggested that ammonia can be oxidized to hydroxylamine by nirS (Liu et al. 2018), which could potentially explain the lack of amoABC genes in our samples. Potential for denitrification was also present in our samples, although mostly in the combined sediment, and not limited by the increased salinity (Laverman et al. 2007).

Osmoadaptation

It has been debated if hypersaline environments are thermodynamic limiting the occurrence of self-sustaining microbial communities (Oren 2011) or if they are biologically permissive (Lee et al. 2018). Cells need to implement strategies to counteract the osmotic stress (Galinski 1995, Gunde-Cimerman et al. 2018) but these strategies come with an energetic cost (Meier et al. 2021). It seems that the ‘salt in’ strategy prevails over the ‘compatible solute’ strategy in the samples of the present study; although they are both present, the former is much more abundant probably due to the high abundance of Halobacteria/Halobacteriota and Halanaerobiia (Mukhtar et al. 2020). This is in contrast to what has been observed in the hypersaline microbial mats of Shark Bay where choline and betaine uptake and betaine biosynthesis were the identified osmoadaptation strategies (Ruvindy et al. 2016).

It has been suggested that hypersaline microbial mats and, in particular, communities below salt crusts, cannot rely solely on primary production from anoxygenic phototrophy and mineralization from sulphate reduction (Meier et al. 2021). Instead, import of reduced substances and periods of reduced salinity are required, to allow the occurrence of oxygenic photosynthesis (Meier et al. 2021); in our study site this occurs in winter where evaporation is not as strong as in the summer, when combined with the increased precipitation, lowers the salinity of the marsh. During winter months, both the salt crust and the layering of the microbial mat disappears, as in Cardoso et al. (2019), and it seems that this temporal change and seasonal development of the microbial mats under study is the necessary element for the survival of the microorganisms. In addition, the potential for anaplerotic reactions, that were abundant in our samples, may play an important role in replenishing the intermediates of the TCA cycle, which is quite abundant in out samples, and thus allowing microbial growth with a carbohydrate as the sole carbon source (Tong 2013, Choi et al. 2016). Although it seems that hypersaline environments are ‘thermodynamically moderate’, DNA based studies can only identify the members of a community and not their metabolic activities. Therefore, future studies on hypersaline microbial mats should focus on the combination of metabolomics, metatranscriptomics, and metagenomics, in order to elucidate the functional repertoire of microbial communities, their metabolic potential and their metabolic and ecological interactions. Metabolic modelling of the microbial assemblages can shed further light on the effects of the environmental challenges on the mat construction as well as on which processes are taking place within each mat layer and among its different layers.

Conflict of interest statement

None declared.

Acknowledgments

We would like to thank the reviewers for their valuable comments and suggestions on our manuscript. We would also like to thank the Management Agency of the Dodecanese Protected Areas (MADPA) and especially Mr Dinos Protopapas and Mr Giorgos Prearis (captain of the R/V Saria) for providing assistance during our visit to the Tristomo bay. In addition, we would like to thank (a) Mr Antonis Potirakis and Mr Stelios Ninidakis, the best HPC system administrators, for their help and support during cluster maintenance and installation of third-party software, (b) Dr Christos A. Christakis (ORCID iD: 0000–0002-7075–0996) for his suggestion to choose a few samples for shotgun metagenomic sequencing so that higher coverage could be achieved, (c) Dr Paschalis Natsidis (ORCID iD: 0000–0002-2446–3208) for his insight on the creation of the phylogenomic trees, (d) Dr Jon Bent Kristoffersen (ORCID iD: 0000–0002-8804–1864) for performing the MiSeq sequencing at the IMBBC DNA Sequencing platform, and (e) Professor Jens Carlsson (ORCID iD: 0000–0002-9262–5627) for allowing us to use the bob server of Area 52 lab in University College Dublin. This research was supported in part through computational resources provided by IMBBC (Institute of Marine Biology, Biotechnology and Aquaculture) of the HCMR (Hellenic Centre for Marine Research). Funding for establishing the IMBBC HPC has been received by the MARBIGEN (EU Regpot) project, LifeWatchGreece RI and the CMBR (Centre for the study and sustainable exploitation of Marine Biological Resources) RI.

Funding

This work was supported by the project RECONNECT (MIS 5017160) financed by the Transnational Cooperation Programme Interreg V-B ‘Balkan-Mediterranean 2014–2020’ and cofunded by the European Union and national funds of the participating countries. There was no additional external funding received for this study. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Author notes

shared coauthorship