Samples analyzed in this study were collected at noon (12.00 h) in May (spring, SP) and November (autumn, AU) during three consecutive years (2011–2013). Environmental temperatures in May and November ranged from 15–18°C and 13.5–15°C, respectively. The mats in all cases were flooded. The salinity of the water covering the mats was 58–62‰ in May and 55–65‰ in November. Mat samples were collected in cores (1 cm × 3 cm) and frozen in liquid nitrogen immediately. Then, cores were stored in the lab at -80°C until DNA extraction. We collected three cores as in previous works (Armitage et al., 2012; Harris et al., 2013), separated by 10 cm each, for each year and season. Our samples were taken each year at the same location. Dillon et al. (2009) sampled cores across 1 km. They observed that population structure diverged with increasing distance between sample sites, but positional replicates were highly similar among samples < 1 m distance. We pooled the extracted DNA for the three samples corresponding to year/season and layer to obtain a representative sample for each year/season/layer. We expected that if there were differences in microbial composition it could be due to seasonal environmental variables and not to location of sampling.

The frozen cores were sliced with a sterile blade in aseptic conditions horizontally in 2-mm increments (from the top to a depth of 6 mm): 0–2 mm (layer 1, oxic/photic layer), 2–4 mm (layer 2, oxic-anoxic transition layer), and 4–6 mm (layer 3, anoxic layer). A fresh blade was used at each interface. Then, pieces of microbial mat of approximately 3 mm³ were cut from each slice and suspended in 100 μl of TE buffer in 2.0-ml vials containing a capful of 0.1-mm glass beads. The mixture was homogenized for 1 min in a Mini-BeadBeater-8 (Biospec Products, Bartlesville, OK, United States), and centrifuged at high speed for 2 min. While avoiding transfer of the beads, ∼500 μl from each sample was pipetted into sterile 1.5-ml Eppendorf tubes. DNA was extracted using a phenol-chloroform mixture and precipitated in the cold using 95% ethanol. Three DNA extractions corresponding to every year, season and layer were performed. The DNAs obtained were mixed to correct for potential local heterogeneity effects to obtain a representative sample for each year, season and layer.

For years 1, 2, and 3, we performed amplicon sequencing of the bacterial 16S rDNA gene. The primers used for multiplex Roche 454 GS FLX pyrosequencing, contained a 25 nucleotide sequence adapter, 10-base-pair molecular barcode (multiplex identifier), and the universal bacterial sequence for the region V1–V2, 8F-338R (5′-AGAGTTTGATCCTGGCTCAG-3′ and 5′-TGCTGCCTCCCGTAGGAGT-3′) (Armitage et al., 2012; Harris et al., 2013; Yang et al., 2013). We used three different barcodes (each one for 0–1 mm, 2–4 mm, and 4–6 mm; ACGAGTGCGT, ACGCTCGACA, AGACGCACTC, respectively). Samples analyzed were: SP1-1, SP1-2, SP1-3; AU1-1, AU1-2, AU1-3; SP2-1, SP2-2, SP2-3; AU2-1, AU2-2, AU2-3; SP3-1, SP3-2, SP3-3; AU3-1, AU3-2, AU3-3 (SP and AU indicated the season analyzed; the first number, the year, and the second number, the layer). A PCR from each DNA was performed. The cycling conditions were 94°C for 3 min, followed by 30 cycles of 94°C for 30 s, 56°C for 40 s, 68°C for 40 s, and a final extension step at 68°C for 6 min. The resulting product was checked for size and purity on an agarose-SYBR Safe DNA gel that was subsequently stained (Invitrogen, San Diego, CA, United States). The amplicons were purified using a Pure Link kit (Invitrogen, San Diego, CA, United States) and quantified using Qubit and Bioanalyzer (Berlanga et al., 2016). A pool of amplicons was mixed in equimolar amounts (e.g., spring 1st year, amplicons obtained for 0–2, 2–4, and 4–6 mm), and then prepared for 454-pyrosequencing according to the manufacturer’s instructions. Pyrosequencing coverage (depth sequencing) resulted in 99,216 total raw reads that after quality control processing resulted in 44,787 reads (see bioinformatic analyses section) for the 18 samples.

Shotgun metagenomic analysis was performed on samples belonging to the third year (SP3-1, SP3-2, SP3-3; AU3-1, AU3-2; AU3-3). We repeated the DNA extraction several times to reach the approximate concentration of 500 ng to 1 μg of DNA for each sample. Random shotgun metagenomics was performed in the Unity of Genomics of Scientific and Technological Centers, University of Barcelona (CCiTUB). Number of sequences ranged from 61370 to 140208. Major scaffold distribution lengths were 390–470 bp.

We combined 16S rRNA amplicons, PICRUSt and shotgun metagenomics using the best of each method to obtain the maximal information to try to describe precisely the taxonomical structure and functionality of the samples. The advantages of using 16S rRNA amplicons sequencing had normally better taxonomic resolution than shotgun metagenomics (Tessler et al., 2017), and the availability of bioinformatic tools for prediction of functions (PICRUSt) is particularly attractive to microbial ecologists as it allows them to study the genes (functions) of complex microbial communities with reasonable accuracy at a high taxonomic resolution (Mukherjee et al., 2017). Random shotgun sequencing of environmental DNA provides a direct and potentially less biased view of the functional attributes of microbial communities (Klatt et al., 2013). 16S rRNA gene regions recovered from the shotgun metagenomic data can span the entire length of the genes; the PCR-based amplicon approach only targets the V1–V2 region. Therefore, the two approaches may not necessarily give identical results (Fierer et al., 2012).

For 16S rRNA amplicons, the raw data of each sample was preprocessed for demultiplex and quality control using a pipeline implemented in GPRO version 1.1 (Futami et al., 2011). Raw reads that contained < 150 nucleotides in size, ambiguities > 1, homopolimer > 8, as well as redundant sequences were removed from each metagenome dataset using screen.seqs and unique.seqs by Mothur1.31.2 (Schloss et al., 2009). Sequences were taxonomically classified using Silva database¹ (Quast et al., 2013). CD-HIT-EST from the CD-HIT 4.5.4 package (Fu et al., 2012) was used to define clusters of clones within each metagenome with a distance threshold of 0.03 (resulting in a cutoff at the species level). Alpha and beta diversity analyses of all samples were performed at 97% distance level of OTU. For diversity we rarified (normalized) samples to compare all the samples. Weighted UniFrac metrics was used to measure beta-diversity and to generate principal coordinates analysis plots, using the normalized OTU table. For the heatmap we used the OTU table at 0.10% genetic distance level (resulting in a cutoff at the family level) (Yarza et al., 2010), and make_otu_heatmap.py, and the script was modified by Stamp program. Hierarchical cluster analysis used for similarity measure was Pearson’s correlation, and for the clustering algorithms, Ward’s linkage.

Core microbiota was determined using compute_core_microbiome.py in qiime² (Caporaso et al., 2010). Core OTUs were defined as the OTUs that are present in at least 90% of the samples. From the set of OTUs that could be considered the core, we performed an ecological network of interactions. Ecological network was achieved by Molecular Ecological Network Approach (MENA)³ (Deng et al., 2012). Ecological network worked with RMT (random matrix). To visualize the network it was used Cytoscape 3.5.1.

Metagenomes were predicted from the 16S rRNA data using PICRUSt (Langille et al., 2013) for samples corresponding to years 1, 2, and 3. This was prepared by the predict_metagenomes.py script against functional database of KEGG Orthology. Functional contributions of various taxa to different KOs were computed with the script metagenome_contributions.py (Mukherjee et al., 2017). For the third-year samples, gene annotation of the shotgun method was analyzed by the United States Department of Energy Joint Genome Institute⁴ (Nordberg et al., 2014).

The DOE-JGI Metagenome Annotation Pipeline (MAP) supports the annotation of metagenomic sequences and it is organized in three stages: sequence data pre-processing, structural annotation, functional annotation and phylogenetic lineage prediction. Some of the processing methodology used by MAP was as follows: Unassembled 454 reads containing more than five occurrences of ‘N’s are removed. Sequences shorter than 150 bp after trimming are also removed. When two or more sequences are at least 95% identical, with their first 3 bp being identical as well, those sequences are considered to be replicates and only the longer copy is retained. For genomic assembler it is used the Velvet algorithm package. A good kmer size is just over half a read length, which prevents sequencing errors from forming bubbles. Ribosomal RNA genes (5S, 16S, 18S, 23S) are predicted using hmmsearch tool from the package HMMER 3.1b2. The pipeline runs against curated models, derived from full-length genes within IMG, while keeping the best scoring models. The identification of protein-coding genes is performed using a consensus of four different ab initio gene prediction tools: prokaryotic GeneMark.hmm (v.2.8), MetaGeneAnnotator (v. Aug 2008), Prodigal (v. 2.6.2) and FragGeneScan. Protein-coding genes with translations shorter than 32 amino acids are deleted. Assignment was made at 90% of the KO gene sequence that was covered by the alignment (Huntemann et al., 2016).

The numbers of the analysis projects in the JGI were Ga0197827, Ga0197828, Ga0197830, Ga0197833, Ga0197836, Ga0197838. For the 16S rRNA amplicons, sequence data were deposited on the NCBI database by the Bioproject PRJNA416849.

Camargue microbial mats composition of Bacteria, Archaea and Eukarya were based on data obtained by shotgun metagenomics. Microbial communities were dominated by Bacteria (92.4–94%), while Archaea and Eukarya represented 4–5% and 1–1.6%, respectively. The distribution of Archaea phyla in spring and autumn was similar but there were several differences across the three depths sampled (0–2, 2–4, and 4–6 mm). Thus, Archaeal relative abundances at those depths were 4.2, 4.6, and 5.1%, respectively. The major phyla were Euryarchaeota (80.6%), followed by Crenarchaeota (8%), Candidatus Micrarchaeota (3.9%), and Thaumarchaeota (ammonia-oxidizing archaea, 2.8%).

The eukaryotic diversity of the Camargue microbial mats was sparse, in contrast to the vast bacterial diversity. This was probably due to the broad metabolic capabilities of Bacteria, which enable them to occupy a broad range of chemical niches, whereas the metabolic versatility of eukaryotes is more limited, despite their ability to survive under high sulfide, fermentative, and anoxic conditions. Eukarya represented 1% of the total relative abundance of microorganisms from the mat, with the most representative eukaryotes those related to algae (Chlorophyta), plants (Streptophyta), fungi (Ascomycota), and Arthropoda (insects, mainly Anopheles). This result contrasts with the findings in the Guerrero Negro microbial mats, where the dominant eukaryotic organisms are bacterivorous nematodes (Feazel et al., 2008).

More than 30 phyla of Bacteria were recovered from amplicon sequencing of the 16S rRNA gene and shotgun metagenomics isolated from the Camargue microbial mats. Among the distinguished phyla, there were six that dominated: Proteobacteria (40.2–75.2%), Bacteroidetes (2.6–16.6%), Firmicutes (1.2–21.5%), Actinobacteria (1.0–10.4%), Cyanobacteria (1.2–36.6%), and Spirochaetes (1.4–6.6%) (Figure 1). The use of shotgun metagenome analysis in the 3rd-year samples did not yield additional phyla compared to the 16S rRNA amplicons analyzed using Silva database. However, it did detect a relatively high abundance of Actinobacteria, Bacteroidetes, and Firmicutes, and a lower relative abundance of Proteobacteria. The distributions of several phyla and their families depended on the layer analyzed and the season (Figure 2 and Supplementary Figures S1A,B). Cyanobacteria were more abundant in the upper layer (0–2 mm) and in autumn. Alphaproteobacteria were the most abundant Proteobacteria, followed by Gammaproteobacteria and Deltaproteobacteria. Alphaproteobacteria, especially the family Rhodobacteraceae, were abundant in the upper layers (0–2 and 2–4 mm) and in the spring. Gammaproteobacteria were represented, in order descendent of relative abundance, by Chromatiaceae, Ectothiorhodospiraceae, and Pseudoalteromonadaceae. The Chromatiaceace family was present in all layers with a slightly increased abundance at 2–4 and 4–6 mm, but no difference between spring and autumn. Ectothiorhodospiraceae were more abundant in the upper layer (0–2 mm) than in the other, deeper layers. Pseudoalteromonadaceae were detected only in the autumn samples. Among the Deltaproteobacteria, Desulfobacteraceae was the most abundant family detected and their distribution in the different layers and in the two seasons was similar. Chloroflexales (phylum Chloroflexi) were more abundant in the autumn samples and in the upper layer. Several genera, such as Thioalkalivibrio, Desulfotigum, Roseovarius, etc., were detected through the years analyzed, but their distribution depended on the layer (Figure 2).

In the rarefaction curves for samples with respect to depth and season, at 0.03 similarity the samples did not reach an asymptote, suggesting their insufficient sequencing depths (Supplementary Figure S2). In our analyses we used the abundance estimator Chao1 and abundance-based coverage estimator (ACE), Shannon and Simpson diversity indexes, and Berger-Parker dominance index (Table 1). The highest diversity was found in the third layer (4–6 mm). The surface layer in all samples exhibited the lowest diversity, with a few strongly dominant OTUs, especially photosynthetic bacteria. A principal component analysis indicated that the population community structure of the upper layer (oxic/photic layer) differed from that of the transitional oxic-anoxic layer (2–4 mm) and the anoxic layer (4–6 mm), which were relatively close together (Figure 3A). These results were independent of the season and the sampling year. Samples for the oxic/photic layer (0–2 mm) had the most distant community distribution, probably due to the high relative abundances of Cyanobacteria and Chloroflexi. The phylogenetic P-test in Unifrac indicated that the microbial communities were not significant different (P > 0.05). But pairwise significance tests using the t-Student based on taxa detected showed significant differences between oxic/photic layer and the other layers (transition oxic-anoxic layer and anoxic layer). Significant differences between two samples were based on a 95.0% confidence level. P < 0.05 was considered to indicate statistical significance. No significant differences were observed between samples from the same layer. In addition, we could observe a clear difference in distribution respect of the functional annotation KO metagenomes to layers (data from shotgun metagenome for the third year samples) (Figure 3B). The t-Student from taxonomical data showed the same results that those just mentioned.

To understand the metabolic potential of the Camargue microbial mats and identify their many different functional features, we used PICRUSt (based on 16S rRNA gene amplicon) and random shotgun metagenomics methodologies. The predicted proteins were classified as KEGG orthologs (KOs). The nearest sequence taxon index (NSTI) values is a measure of how closely related the OTUs in each sample are to the reference genomes in the database. In our case, the “nearest sequence taxon index” (NSTI) values per sample ranged for 0.072–0.172. The taxonomical classification could be accurate at family level, in few cases to genera level, but it was difficult to achieve the species level. This result could explain the values observed at NSTI. Respect to the shotgun metagenomics, KEGG pathways via KO (percentage) ranged from 14.85 to 17.17%; and KO genes ranged from 24.97 to 29.30% respect to the number of sequences (total sequences analyzed ranged from 61,370 to 140,208), and assignment was made at 90% of the KO gene sequence that was covered by the alignment (Huntemann et al., 2016).

The biological processes identified were essential for sustaining prokaryotic life in the environment. They include transcription and translation functions (8.7–9.8% relative abundance genes, based on the total number of genes detected in the sample) and replication and repair functions (9–10.2%). Other functional processes were related to cellular processes such as cell motility (4.3–5.9%). Genes related to membrane transport (17.2–19.5%) and to metabolic functions (56.3–58.4%), which included the metabolism of carbohydrates, lipids, amino acid, cofactors and vitamins, xenobiotic biodegradation, and energy metabolism. PICRUSt and shotgun metagenomic analyses revealed similar functional biological processes in all samples analyzed, except of carbohydrate metabolism and energy metabolism, which contained more genes detected by shotgun than by PICRUSt analyses (Figure 4).

Gene content analysis provides a basis for inferring the possible metabolic functions of dominant populations present in the community. Cell motility, represented by chemotaxis genes, such as cherA, cheBR, motA, mcp, pixJ, etc., peaked at the oxic-anoxic transition zone, but they were also important in the oxic/photic zone. These genes were associated to phototrophic organisms (Cyanobacteria, Alpha-Gammaproteobacteria, and Chloroflexi), but also to heterotrophic members, such as Bacteroidetes and Spirochaetes. Ferredoxins have a negative redox potential act as electron distributors in various metabolic pathways. The genes that codify different ferredoxins were detected in all layers, but they were especially abundant in the oxic-anoxic transition zone. Ferredoxins likely reflected diversification of redox reactions required for respiration (Supplementary Table S1). Osmotic regulation is required for microbial survival in hypersaline environments. Accumulation of osmoprotective molecules, in particular glycine betaine, is an adaptive mechanism to pawn the high salinity conditions. We searched for genes implicated on the glycine-betaine biosynthesis, such as betA, betB, gbsA (Wong et al., 2015). These genes were distributed through the layers (especially on the oxic/photic layer), and they were associated with different taxa, showing that the microbial mat community could be adapted to salinity conditions (Supplementary Table S1).

Genes associated to oxygenic photosynthesis and bacteriochlorphylls were detected in the upper layer (photic zone) in autumn and spring samples (Table S1). Regarding the photosynthetic reaction center in the anoxygenic photho-breaksystem, pufL and pufM genes were detected and they belonged to Gammaproteobacteria (purple sulfur bacteria, Chromatiaceae). The possibility of alternative light energy usage by (bacterio)rhodopsin in different prokaryotic members of the mat cannot be confirmed because we could not detect related genes, even though retinal-based phototrophy could contribute as energy source in layers with low irradiance (Thiel et al., 2017).

In the studied metagenomes we identified the four known autotrophic carbon fixation pathways (the Calvin-Benson cycle, the reverse tricarboxylic acid cycle, the Wood–Ljundahl pathway, and the 3-hydroxypropionate bi-cycle) (Hügler et al., 2002, 2005; Ragsdale and Pierce, 2008; Berg, 2011; Thiel et al., 2017), which suggested the occurrence of a relatively diverse autotrophic community (Supplementary Table S1). Debris of predated bacteria by Bdellovibrionaceae and viruses may be another carbon source for other heterotrophs. Bdellovibrionaceae represented 1–9% of the relative abundance of Deltaproteobacteria. Bdellovibrionaceae are predatory bacteria upon a variety of Gram-negative bacteria. Viruses in the mat could be also involved in cell-lysis processes, based on CRISPR systems detected in the metagenomes.

Regarding nitrogen metabolism, we detected genes associated to nitrogen fixation. The oxic/photic zone contained the more diverse and abundant amount of nitrogen fixation genes (Supplementary Table S1). For the ammonium oxidation in nitrification, the main enzyme is the ammonia monooxygenase (amoA) that is present in both ammonia-oxidizing archaea and ammonium-oxidizing bacteria (Fan et al., 2015). However, amoA was not identified in the studied metagenomes. Nitritification provides the oxidant for anaerobic ammonium oxidation (anammox). We examined hzoA/hzoB genes because their ubiquity and high expression in anammox bacteria (Planctomycetes) (Hirsch et al., 2011), but no records of those genes were found, although the phylum was detected in the Camargue microbial mats.

Sulfate reduction genes were present in the metagenome dataset and distributed similarly through the different layers (Supplementary Table S1). They affiliated to Deltaproteobacteria and Gammaproteobacteria. Sulfur oxidation activity was also found in the Camargue microbial mats, based upon the presence of the enzyme sulfide:quinone oxidoreductase (sqr gene).

Finally, to identify potential biotic interactions within the dominant, prokaryotic communities in the Camargue microbial mats, we constructed a network based on the core OTUs (Figure 5). The core OTU were determined by the shared OTUs at 90% in all samples. Several minority populations were not included as the “core community,” and probably they could play important functions. Core community was performed by layer. We observed that there were no significant differences for one layer (e.g., 0–2 mm or 2–4 mm or 4–6 mm) among years and season. In the upper layer (oxic/photic zone), the core microbiota were represented Alphaproteobacteria (Rhodobacterales), Gammaproteobacteria (Marinicellales and Chromatiales—Ectothiorhodospiraceae) and Cyanobacteria, (Coleofasciculus [formerly, Microcoleus], Oscillatoriales). In the middle layer (transition oxic-anoxic zone), there were OTUs belonging to Alphaproteobacteria (Rhodobacterales), to Gammaproteobacteria (Marinicellales, Chromatiales, Thiotrichales), to Deltaproteobacteria (Desulfobacterales), to Bacteroidetes (Flavobacteriales and Cytophagales), to Spirochaetes, and to Cyanobacteria (Oscillatoriales, Coleofasciculus). In the bottom layer (anoxic zone), taxons were represented by Alphaproteobacteria (Rhodobacterales), OTUs to Gammaproteobacteria (Marinicellales and Chromatiales), Deltaproteobacteria (Desulfobacterales), Bacteroidetes (Flavobacteriales. Bacteroidales and Cytophagales), Spirochaetes, Planctomycetes, and Gemmatimonadetes. The network obtained showed that interactions among taxa could be done through different layers (Figure 5). The network was probably incomplete because there were not represented other populations, such as Firmicutes, Actinobacteria, Archaea, or other minor population with less than 0.1% of relative abundance (Planctomycetes, Nitrospinae, Saccharibacteria, etc.), which they could contribute and participate on metabolically interactions within the microbial mat. Cyanobacteria stablished the more diverse interactions with different population’s members of the microbial mat. Positive correlations based the thickness of the lines (in the figure marked by a purple line) and were observed between Cyanobacteria and Deltaproteobacteria; Cyanobacteria and Bacteroidetes; Cyanobacteria and Rhodobacterales; and Spirochaetes and Deltaproteobacteria.