Biofilm sampling was conducted in May 2013 in the Thuwal seeps II (22°16N-38°53E) by the remotely operated underwater vehicle (ROV) Max Rover, DSSI, USA during the Red Sea exploration cruise. Biofilms were developed on six different types of materials in the brine pool for 3 and 6 days. The devices for biofilm development has been reported in Zhang et al. (2014) and the schematic diagrams are also shown in Figure S1 in the present study. The materials were deployed into the brine pool by the ROV. Environmental parameters, such as temperature, salinity, and dissolved oxygen (DO) were measured using conductivity-temperature-depth (CTD) and SBE43 DO sensors attached to the ROV.

The biofilms used in the present study were developed on aluminum and PVC in the brine pool for 6 days and were referred to as “biofilm_Al and biofilm_PVC.” Biofilms on these two materials were selected because biofilm_Al had the highest biomass and PVC represented a biofilm developed on plastic materials. Cells were recovered from the substrates using sterile cell scrapers and suspended in Tris-HCl buffer before extracting nucleic acid using the AllPrep DNA/RNA Mini Kit (Qiagen, Hilden, Germany). The extracted DNA was sequenced using an Illumina HiSeq 2000 platform. Metagenomic datasets of biofilm_Al and biofilm_PVC have been documented in Zhang et al. (2015), and the information of these two metagenomes are also listed in Table S2 in the present study. Assembly of the combined metagenomes of biofilm_Al and biofilm_PVC was performed using SPAdes Genome Assembler 3.6.1 (Bankevich et al., 2012). The specified K values for assembly included 21, 31, 41, 51, 61, 71, and 81, and the “—careful” and “—pe” options were used. The metagenome-based taxonomic structures were obtained by ribosomal database project (RDP) classification (Wang et al., 2007) of the extracted 16S rRNA sequences, and the coverage information.

Draft genome binning was conducted according to previously described methods (Albertsen et al., 2013; Tian et al., 2014). The genomes were separated mainly on the genome coverage, GC content, and tetranucleotide frequency. The Illumina reads were mapped to the contigs using Bowtie2 (version 2.0.0; Langmead and Salzberg, 2012). The genome coverage was calculated using SAMtools on a local server. The contigs were searched against a set of single-copy protein-encoding genes (Albertsen et al., 2013) using default cut-off values. The single-copy protein-encoding genes identified were searched against the National Center for Biotechnology Information non-redundant (NCBInr) database with BLASTP (e-value < 1e-07) with “xml” output format. Taxonomic information for the contigs was obtained by importing the BLAST results into MEGAN 5.0 (Huson et al., 2011). The metagenome pairs biofilm_Al and biofilm_PVC were selected to obtain the most effective separation and extraction of contigs. It must be noted that the diversity of the biofilm communities was not high and that many contigs belonging to different genomes formed clear, separate clusters based on coverage, facilitating genome binning and elimination of contamination (Figure S2).

Potential contig contamination was examined using CheckM (Parks et al., 2015), and PhylopythiaS, which is a sequence composition-based classifier that utilizes the hierarchical relationships between clades (Patil et al., 2012). Moreover, the Thaumarchaeota MAGs in the present study were compared with the eight genomes obtained using single-cell genomics by Ngugi et al. (2015), based on BLASTP (e-value < 1e-07) and BLASTN (e-value < 1e-07) with predicted open reading frames (ORFs) as input. The results were visualized using the Artemis Comparison Tool (ACT; Carver et al., 2005).

The MAGs were annotated by searching the Carbohydrate-Active EnZymes (CAZys) (Cantarel et al., 2009), Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa and Goto, 2000), and Clusters of Orthologous Groups (COGs) (Galperin et al., 2014) databases. The ORFs in the MAGs were predicted using Prodigal (Hyatt et al., 2010) using a local server. Genes were annotated by BLASTP searches against the above-listed databases using an e-value cut-off of 1e-07 on a local server. The online tool KAAS (KEGG Automatic Annotation Server) was also used to validate the KEGG annotation by local BLAST. Metabolic pathways were analyzed using online tools in KEGGMAPPER (http://www.genome.jp/kegg/mapper.html). Key enzymes revealed by KEGGMAPPER were rechecked against the NCBInr database based on the BLAST results.

In the phylogenetic analysis of the functional genes (e.g., polysaccharide deacetylase), the protein sequence alignment was determined by the molecular evolutionary genetics analysis using Muscle (MEGA 6.06; Tamura et al., 2013) with the following parameters: gap open -50, cluster method UPGMB and min diagonal length 24. Gblock software (Castresana, 2000) was applied to eliminate less informative sites in the alignments. The construction of maximum-likelihood trees was conducted using MEGA 6.06 with the Tamura-Nei model, the Nearest-Neighbor-Interchange (NNI) method and 500 bootstrap replicates.

Biofilms on different substrates were harvested (10 slides per sample) onboard immediately after recovery using sterile cotton tips and stored in RNAlater buffer. RNA extraction was performed using the AllPrep DNA/RNA Mini Kit (Qiagen, Hilden, Germany). Bacterial cells in RNAlater buffer were pelleted by centrifugation at 4000 g for 10 min and then lysed with lysozyme, proteinase K and RLT buffer provided in the kit. Total nucleic acid was extracted, and DNA was removed using the DNA column provided with the AllPrep DNA/RNA Mini Kit. The integrity and quantity of the RNA were checked by agarose gel electrophoresis and a Nanodrop device (ND-1000 spectrophotometer, DiaMed China Limited, Hong Kong), respectively. Prior to first-strand cDNA synthesis, genomic DNA contamination was further removed using the DNA-free™ DNA Removal Kit (Ambion, Austin, TX). The extracted RNA was then used as a template for whole-cDNA amplification using the Ovation Single-Cell RNA-Seq System, which enriches for non-rRNA during cDNA synthesis using proprietary whole transcriptome primers targeting non-rRNA sequences in the transcriptome. The synthesized cDNA was sequenced using an Illumina HiSeq 2000 platform.

The information for metatranscriptomic sequences are summarized in Table S2. The sequenced cDNA was subjected to quality control using the next-generation sequencing (NGS) QC toolkit (Patel and Jain, 2012). Assembly of the metatranscriptomic reads was performed using SPAdes Genome Assembler 3.6.1 (Bankevich et al., 2012) on a local server. The specified K values 21, 31, 41, 51, 61, 71, and 81 were used under the “—careful” and “—pe” options. The ORFs were predicted following the steps mentioned above and annotated by BLASTP against the COG -database. To depict the gene expression profiling for individual microorganisms, the reads were aligned to contigs of the MAGs using Bowtie2 (version 2.0.0; Langmead and Salzberg, 2012). Gene expression abundance is the number of transcripts mapped to the CAZys of each MAG that normalized to the length of the genes.

The MAGs of Cloacimonetes were deposited in the Whole Genome Shotgun (WGS) database under accession numbers of SAMN04423165, SAMN04423166, SAMN04423167, and SAMN04423168. The MAGs of Marinimicrobia were deposited in the WGS database under accession numbers of SAMN04423157, SAMN04423158, SAMN04423159, SAMN04423160, SAMN04423161, SAMN04423162, and SAMN04423164. The MAG of Bathyarchaeota was deposited in the WGS database under accession number of SAMN04423169.

Genome binning using two metagenomes, biofilm_Al and biofilm_PVC, generated 51 MAGs of the biofilm-inhabited microorganisms (TCS1-51; “TCS” served as the abbreviation for “Thuwal cold seep”), including members of Actinobacteria (two MAGs belonging to Candidatus Microthrix), Gammaproteobacteria, Deltaproteobacteria, Epsilonprtoteobacteria, Candidatus Marinimicrobia, Candidatus Cloacimonetes, Bathyarchaeota (previously referred to as Miscellaneous Crenarchaeotic Group, MCG) and Thaumarchaeota (Table 1). Based on the metagenome-derived 16S rRNA genes, these phyla were the dominant microbial groups in the biofilms (Figure S3). According to the results of CheckM, these MAGs had a maximum completeness value of 100% and a maximum potential contamination of 2.4%. Based on the results of PhylopythiaS, contigs of most of the MAGs were very consistent, as suggested by the percentage of contigs assigned to a defined taxonomy. For example, 97.2% of the Pseudomonas sp. nov (TCS13) contigs could be assigned to the genus Pseudomonas. Moreover, contigs of the newly defined bacterial phyla, were compared with their relatives with complete or single-cell derived genomes to assess the potential contamination. For example, 95.8% of the Thaumarchaeota sp. nov. 1 (TCS50) and 91.3% of the Thaumarchaeota sp. nov. 2 (TCS51) could be assigned to the eight single-cell Thaumarchaeota genomes, which were from another brine pool in the Red Sea (Ngugi et al., 2015). The assigned ORFs were distributed in all of the contigs of the TCS MAGs as exemplified by Thaumarchaeota sp. nov. 1 (Figure S4). These results indicated that the MAGs reconstructed in the present study were of high quality.

A large proportion of the genes identified in many of the MAGs were found to be involved in fermentative carbohydrate metabolism. The distribution of CAZys (Cantarel et al., 2009) responsible for the degradation of complex carbohydrates is listed in Table S1, such as β-glucosidase (GH3), α-amylase (GH13) and endo-β-1,4-galactanase (GH53). In total 17 unique CAZy categories corresponding to 1692 genes were identified in the 51 MAGs (on average 33 genes/genome). Genes encoding proteins playing a role in electron transport and energy production, including those encoding F₄₂₀H₂ dehydrogenase (fpo), cytochrome C (cytoC), the rhodobacter nitrogen fixation (rnf) complex, nitrate reductase and sulfate reductase, were identified in many of the MAGs. Annotation by COGs (Galperin et al., 2014) showed a high diversity in the gene arrangement of the F₄₂₀ biosynthesis gene clusters (Figure S5), which were present in eight of the MAGs belonging to different phyla. Moreover, putative bacterial microcompartment (BMC) clusters, types of protein-bound bacterial organelles that function in carbon catabolism and storage (Nobu et al., 2016), were present in five MAGs including those of Actinobacteria, Marinimicrobia and Desulfobacter (Figure S6). As reported by Nobu et al. (2016), genes related to aldehyde and sugar metabolism are present in the BMC clusters, including rnf-C like NADH dehydrogenase (COG4656), sugar isomerase (COG0363) and aldehyde dehydrogenase (COG1012). However, other genes were also identified in the present study, such as NTP pyrophosphohydrolases and oxidative damage repair enzymes (COG0494), making the BMC clusters different from those in previous study.

The prevalence and versatility of carbohydrate metabolism and energy generation was further supported by KEGG-based pathway reconstruction. The central metabolic strategies are summarized in Figure 1. Many of the MAGs encoded proteins associated with monosaccharide and polysaccharide transport, nitrate and sulfate reduction, electron transport systems (ETSs), glycolysis, and reoxidation of energy-conserving elements via acetyl-CoA reduction to ethanol or acetate. For example, monosaccharide and polysaccharide transport was indicated by the presence of a number of ABC transporters in all the bacterial genomes; dissimilatory nitrate reduction was indicated by both respiratory nitrate reductase (nar) and periplasmic nitrate reductases (nap); dissimilatory sulfate reduction was indicated by sulfur oxidation (SOX) gene cluster and adenylylsulfate reductase; and acetate generation was indicated by acetate kinase. Notably, the Marinimicrobia bacterium sp. nov. 2 (TCS39) and Bathyarchaeota sp. nov (TCS49) possessed both dissimilatory nitrogen and sulfur reduction pathways.

When the TCS MAGs were compared with reference genomes within the same phylum, unique features concerning carbohydrate metabolism and energy generation were observed. For example, genes encoding NADH dehydrogenase, such as nuoB, nuoL, ndufv2, NdhJ, and NoxF were only identified in four Cloacimonetes MAGs (TCS45-47) when compared with complete or single-cell Cloacimonetes genomes (Figure S7); similar results were observed for chitin deacetylase (EC3.5.1.41) and the α-D-glucose transformation pathway, which are important steps in glycolysis. Moreover, the complete pyruvate to acetate pathway was identified in Marinimicrobia bacterium sp. nov. 2 (TCS39) and Marinimicrobia bacterium sp. nov. 3 (TCS40), but the enzymes included in this pathway could not be found in the four reference genomes, the single-cell Marinimicrobia genomes from previous studies (Rinke et al., 2013; Nobu et al., 2015) (Figure S8). In addition, the Bathyarchaeota MAG (TCS49) harbored specific genes for chitin, cellulose, and β-D-glucose metabolism, such as a glycoside hydrolase family 3 gene (probably chitin deacetylase; EC3.5.1.41) and β-glucosidase (EC3.2.1.4) (Figure S9); the presence of aprAB and cysC genes, which are involved in dissimilatory sulfate reduction, differentiated TCS49 from the reference genomes.

Given that a number of unique genes were identified from these MAGs as compared with their relatives, we examined the phylogeny of these genes. The polysaccharide deacetylase (EC3.2.1.21) in Cloacimonetes bacterium sp. nov. 1 (TCS45) was phylogenetically close to the polysaccharide deacetylase from Lactobacillus ceti (Figure 2). Although, the bootstrap values in the tree were not high, BLAST analysis also revealed close relationship between polysaccharide deacetylase in Cloacimonetes bacterium sp. nov. 1 and those from firmicutes. The nitrate reductase α and β units in Marinimicrobia bacterium sp. nov. 3 (TCS40) were phylogenetically close to the nitrate reductases from Geothrix fermentans (Figures S10, S11). In addition, the glycoside hydrolase family 3 gene in Bathyarchaeota sp. nov (TCS49) was affiliated with Hungatella hathewayi and Clostridium hathewayi (Figure S12). The adjacent genes in the same contigs, which were affiliated with Cloacimonetes, Marinimicrobia or Bathyarchaeota, were displayed to confirm that the presence of these novel genes was probably not attributable to contig contamination.

To explore the major active functions, we obtained ~600 Mb of metatranscriptomic datasets for the biofilms developed on Al and PVC. The metatranscriptomic sequences were searched against the COG database, revealing an overall picture of the active functions. Transcripts of the genes involved in carbohydrate transport and metabolism were among the most abundant in the metatranscriptomes, particularly the glycosidases (COG0366), β-glucosidase-related glycosidases (COG1472) and β -galactosidase/beta-glucuronidase (COG3250) (Figure 3). In contrast, few peptidases that are essential for protein metabolism were identified, and rank analysis suggested relative lower expression levels for these genes (Figure S13).

The reads from these two metatranscriptomes were mapped to the 51 MAGs. Abundant active bacteria could be identified for nearly all of the phyla, whereas rare active bacteria were also included in each phylum. Key functional genes involved in nitrate and sulfate respiration were expressed, such as the napAB, narGH, nirK, norC, nosZ, and sulfur oxidation (SOX) gene cluster. The gene expression information was integrated in to pathways shown in Figure 1. Moreover, when the metatranscriptomic sequences were mapped to CAZys in each MAG, a heatmap showing expression pattern of CAZys was generated (Figure 4). The results confirmed the expression of diverse carbohydrate metabolism related genes, and indicated relative higher expression of the glycoside hydrolases (GHs), glycosyltransferase (GTs) than enzymes in other families.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.