During the period of 2013 to 2017, we collected 101 biofilm samples developed on different material surfaces, including Petri dishes, natural rocks, and panels made of zinc, aluminum, titanium, stainless steel, polyether ether ketone, polytetrafluoroethylene, and polyvinyl chloride, and 24 seawater samples from eight locations in the Atlantic Ocean, the Pacific Ocean, the Indian Ocean, and the Red Sea (Zhang et al., 2019). The duration for the development of biofilm samples varied from 3 days to 1 month. Microbes in biofilms were scraped from panels with autoclaved cotton tips. Two liters of seawater was filtered from the adjacent water column where biofilms were developed and filtered through a 0.1 μm membrane filter to collect microbial cells. Microbes from biofilms and seawaters were stored in DNA extraction buffer as described by Zhang et al. (2019). The data sets deposited under Bioproject PRJNA438384 by Zhang et al. (2019) were used in the present study, and the detailed information for biofilm development and adjacent seawater collection is given in Supplementary Figure 1 and Supplementary Table 1.

DNA extraction and sequencing were described in our previous study (Zhang et al., 2019). For DNA sequencing, 350 bp libraries were constructed and sequenced on the Illumina HiSeq X Ten platform in Novogene Bioinformatics Institute (Novogene, Beijing, China) and the HiSeq 2500 System at Beijing Genomics Institute (BGI, Beijing, China). In total, 2.3 terabases of raw reads was generated for all subtidal biofilm metagenomes and adjacent seawater metagenomes. In addition, 67 surface seawater metagenomes were downloaded from Tara Oceans (Sunagawa et al., 2015) as a reference for further study.

Quality control of the biofilm and seawater metagenomes was performed using the next-generation sequencing quality control (QC) toolkit (version 2.0; Patel and Jain, 2012). Low-quality reads with a quality score <20 for >30% of the read length were removed. Subsequently, reads from biofilm and seawater samples were assembled into contigs by MEGAHIT (Li et al., 2015) with kmer values of 21–121, which increase in steps of 10.

To profile the type and abundance of signal transduction genes, metagenomic reads of each sample were mapped to NCBI Nr and Kyoto Encyclopedia of Genes and Genomes (KEGG; Kanehisa et al., 2012) using DIAMOND BLASTx (e-value < 1e−7; >60% sequence identity for >60% of the read length) (Buchfink et al., 2015). All the metagenomes were normalized into 10,000,000 reads per sample. The metagenomic reads were then classified on the basis of the SEED category (level 1 SEED functional category classification system was used here) in MEGAN (version 6.9.3) (Huson et al., 2016) and on KEGG PATHWAY Database.¹ The abundance of a given functional category is indicated by the number of metagenomic reads mapped to genes belonging to its given category. In SEED categories, the mapping results of reads classified into the category “Regulation and Cell Signaling” were selected for subsequent analyses, and the subcategories (i.e., level 2 SEED categories) were visualized. In the KEGG pathway, all the reads mapped to genes belonging to “09132 Signal Transduction” were collected. The significance of dissimilarities related to signal transduction systems between biofilm and seawater microbial communities was tested by two-tailed Student’s t-test using the SPSS software package (version 25.0) with relative abundance of the subcategories under “Regulation and Cell Signaling.” Principal coordinates analysis (PCoA) was performed on the basis of the Bray–Curtis distance matrix constructed in PAST (version 2.0) (Hammer et al., 2001).

To identify the source microbes of signal transduction genes, the reads mapped to the “Regulation and Cell Signaling” category (level 2, SEED subcategory) were then mapped to NCBI Nr database by using DIAMOND BLASTx (e-value < 1e−7; >60% sequence identity for >60% of the read length) (Buchfink et al., 2015) to identify the taxonomic affiliation. Each resulting read was subsequently classified into relevant phylum level (class level of proteobacteria) in MEGAN (version 6.9.3) (Huson et al., 2016). The contribution of each phylum to the taxonomic difference was assessed using permutational multivariate analysis of variance (PERMANOVA) and SIMPER analysis based on the source phylum composition of all biofilm and seawater microbiota by software PAST (version 2.0) (Hammer et al., 2001).

Non-random co-occurrence analysis was performed according to the methods described in a previous publication (Freilich et al., 2010). The source phyla of signal transduction genes were used to construct a correlation coefficient matrix in the R studio using the psych (Revelle, 2015) and reshape2 packages (Wickham, 2007). The network of co-occurring phyla was established by setting Spearman’s correlation coefficient at >0.6 and p-value at <0.01. Then, the matrix was constructed on the basis of the Benjamini-Hochberg standard false discovery rate correction (FDR-BH) to achieve multiple corrections prior to visualization by the interactive platform Gephi (0.9.2) (Bastian et al., 2009).

Open reading frames (ORFs) were predicted on the basis of contigs by MetaGeneMark (version 2.8) (Zhu et al., 2010) using the following parameter settings: gene prediction on both strands; gene overlaps were allowed; probability of initiation and termination non-coding state 0.5. Redundant protein sequences (minimum length of 100 bp) derived from the Red Sea biofilm and seawater metagenomes were removed by CD-HIT (>95% sequence identity for >60% of the length of the shorter sequences) (Fu et al., 2012).

Predicted proteins were annotated by COG database (Galperin et al., 2015) based on COG Hidden Markov Models (Dibrova et al., 2017) using HMMER hmmscan (e-value < 1e−7) (Finn et al., 2011), and the identified ORFs were assigned to COG categories in “T” function (“T” function in COG: signal transduction). Protein sequences in COG2204 (response regulator containing CheY-like receiver, AAA-type ATPase, and DNA-binding domains) were selected as the reference, and the genes mapped to this category were used to construct the phylogenetic tree with ClustW method and 1,000 bootstraps by MEGA 6 (Tamura et al., 2013).

Open reading frames from the Red Sea biofilm and seawater metagenomes were identified as the signal transduction genes in the COG database. The biofilm-derived proteins were then compared with seawater-derived ones by Blastp (e-value < 1e−7; >60% sequence identity for >60% of the read length) to identify top alignments (Johnson et al., 2008). The generated pairs over 70% identity and 100 bp alignment were selected as signal transduction genes of high similarity.

The profile of signal transduction genes in the marine environment was presented based on metagenomic tools (Figure 1). In the metagenomes of 101 biofilm samples and 91 seawater samples, the visualized signal transduction and regulation profiles at level 2 SEED functional classification were divided into 42 categories under “Regulation and Cell Signaling” (Figure 1). These categories included several quorum sensing (QS) gene families (e.g., QS in Vibrio, QS regulation in Pseudomonas, QS in Yersinia, AI-2 transporters, and AHL inducers), two-component regulator families (e.g., two-component regulatory systems in Campylobacter, oxygen and light sensor PpaA-PpsR, and mazE-mazF system), and unknown regulatory genes related to biofilm formation (e.g., biofilm formation in Staphylococcus) and virulence (e.g., Streptococcus pyogenes virulence regulators). The composition and relative abundance of these signal transduction categories displayed distinct patterns between the biofilm and seawater microbial communities (Figure 1), indicating the specific profile of signal transduction systems in marine biofilms.

Principal coordinates analysis based on the “Regulation and Cell Signaling” (subsystems and abundances) was conducted to confirm the differentiation of the signal transduction genes in the marine environment, which clearly separated the biofilms and seawater samples in the plot (Figure 2). Moreover, all the biofilm samples were clustered in a group, whereas the seawater samples were scattered (Figure 2), indicating smaller variations of signal transduction genes in biofilm microbial communities.

To further explore features of the signal transduction in marine biofilms, statistical analyses were performed to identify the abundances of signal transduction gene categories between the global biofilm and seawater samples. Most genes in the “Regulation and Cell Signaling” categories of SEED were more abundant in biofilms than those in seawater (two-tailed Student’s t-test, p < 0.01). A total of 35 subsystems were significantly (p < 0.01) enriched in the biofilm samples, whereas only two were significantly enriched in the seawater samples (Figure 3). Notably, almost all the genes related to QS were concentrated in biofilms. For example, the genes that belonged to the AHL autoinducer category were abundant (with up to 4,000 reads per biofilm metagenome), but almost no reads could be mapped by the seawater metagenomes (Figure 3). By contrast, the gene categories of G-protein-coupled receptor systems in Actinobacteria and of WhiB and WhiB-type regulatory proteins were enriched in the seawater samples (Figure 3). To confirm these results, we counted the number of reads that were mapped to “09132 Signal transduction” categories according to KEGG pathway classification. We found that only the HIF-1 signaling pathway was enriched in the seawater samples, whereas almost all signaling pathways were enriched in the biofilm samples (Supplementary Figure 2). Furthermore, the abundance of these genes varied slightly in biofilms but dramatically in seawater microbial communities (Figure 3), which was consistent with the results of PCoA (Figure 2).

To understand the source of signal transduction genes, we determined the taxonomic affiliation of these genes in biofilm and seawater. The results demonstrated that dominant source phyla (here and after, Proteobacteria fallen into the class level) in biofilm and seawater were different, which could be confirmed by PERMANOVA test (p < 0.001). In the biofilm microbial communities, Gammaproteobacteria contributed 46.38% signal transduction genes, followed by Alphaproteobacteria (23.58%) and Cyanobacteria (10.71%), whereas among microbes in seawater, Alphaproteobacteria had the most contribution (39.66%), followed by Cyanobacteria (26.49%) and Gammaproteobacteria (17.71%; Figure 4). SIMPER analysis unveiled the top five dominant source phyla that contributed to the taxonomic differentiation of signal transduction genes, accounting for 85.73% taxonomic difference between biofilm and seawater samples (Supplementary Figure 3).

Considering that the microbial community of biofilm seeded by surrounding seawater is not an insular and enclosed ecosystem in the marine environment, the signal transduction to connect the interspecies interactions of microbes from biofilm and neighboring seawater is important. Therefore, 12 biofilm samples and 12 seawater samples collected in pairs from the Red Sea were selected for the following analyses. The co-occurrence network was constructed based on all source phyla of signal transduction genes to explore the interspecies interactions between biofilms and adjacent seawater. Among the 331 potential correlations identified, 137 inter-phyla interactions between the Red Sea biofilm communities and their adjacent seawater communities were identified (Figure 5). In Red Sea biofilms, 30 phyla were involved in the inter-phyla relationships, of which 25 phyla could achieve inter-phyla correlations with the 31 phyla derived from adjacent seawater samples (Supplementary Figure 4). Acidithiobacillia contributed the most connections in biofilms, followed by Alphaproteobacteria and Deltaproteobacteria/Epsilon subdivisions, whereas Fusobacteria, Firmicutes, and Bacteroidetes/Chlorobi were found in seawater (Figure 5).

We hypothesized that if biofilm and adjacent seawater communities shared high sequence similarity of substantial signal transduction genes, the related signal transduction systems could be employed and responded by microbes both from marine biofilms and neighboring planktonic counterparts. To test the hypothesis, the sequence similarity of signal transduction genes between Red Sea biofilm and the nearby seawater samples were identified. There were 1,74,273 protein sequences in Red Sea biofilm samples (n = 12) annotated as signal transduction genes by COG database, and a total of 15,780 (9.05%) protein sequences shared over 70% identity and 100 bp alignment with those from the seawater samples (n = 12), which could be considered as highly similar signal transduction genes (Supplementary Figure 5). To confirm the biological basis of interspecies signaling interactions between surface-associated and free-living microbial communities, a phylogenetic tree of a protein ortholog was established. The proteins assigned to a subunit COG2204 (response regulator containing CheY-like receiver, AAA-type ATPase, and DNA-binding domains) of a two-component transduction system were selected and the phylogenetic tree showed a relatively independent branch for the response regulator proteins from the Red Sea biofilm and seawater (Figure 6), indicating the high protein similarity of signal transduction genes between marine biofilms and adjacent seawater.