First Record of Microbiomes of Sponges Collected From the Persian Gulf, Using Tag Pyrosequencing

The Persian Gulf is a special habitat of marine sponges whose bacterial communities are under-investigated. Recently, next-generation sequencing technology has comprehensively improved the knowledge of marine sponge-associated bacteria. For the first time, this study aimed to evaluate the diversity of the Persian Gulf sponge-associated bacteria using tag pyrosequencing in Iran. In this study, 10 sponge samples from 6 different taxonomic orders were collected from the Persian Gulf using SCUBA diving. The diversity of the bacteria associated with the marine sponges was investigated using the 16S rRNA gene PCR-tagged pyrosequencing method. A total of 68,628 high-quality sequences were obtained and clustered at a 97% similarity into 724 unique operational taxonomic units (OTUs), representing 17 bacterial phyla. Cyanobacteria was the most abundant phylum in the sponges, followed by Proteobacteria, Chloroflexi, Acidobacteria, and Actinobacteria. Other phyla were detected as minor groups of bacteria. Bacterial community richness, Shannon, and Simpson indices revealed the highest diversity in sponge S11 (Dictyoceratida sp.) compared to other sponges. This study showed a diverse structure of bacterial communities associated with the Persian Gulf sponges. The dominance of Cyanobacteria may suggest an ecological importance of this phylum in the Persian Gulf sponges.

Extensive studies have been conducted to investigate the bacterial communities associated with different sponge species using both culture-dependent and culture-independent techniques (Giles et al., 2013;Jeong et al., 2014). As is common with most environments, <1% of bacteria present in sponge tissues can be successfully cultivated (Nam et al., 2011;Jackson et al., 2012). In the past two decades, culture-independent methods (mainly based on 16S rRNA gene) have led to a deeper understanding of the microbial diversity in sponges (Alex and Antunes, 2015). Numerous sponge-associated bacteria have been identified using culture-independent techniques such as denaturing gradient gel electrophoresis (DGGE) (Li et al., 2006), fluorescent in situ hybridization (FISH) , terminal restriction fragment length polymorphism (TRFLP) analyses (Zhang et al., 2006;Lee et al., 2011), PCR cloning and sequencing (Kennedy et al., 2008). The 454 tag pyrosequencing is a next-generation sequencing (NGS) technology that provides a faster and simpler way to analyze the microbial communities associated with marine sponges (Webster et al., 2010;Lee et al., 2011;Schmitt et al., 2012;Gao et al., 2014a;Gaikwad et al., 2016). This new method enables hundreds of thousands of nucleotide sequences from multiple samples to be examined in a single 10 h reaction (Lee et al., 2011;Nam et al., 2011;Jeong et al., 2013Jeong et al., , 2014. The next-generation sequencing techniques have revealed the presence of more than 25 different bacterial phyla and 2 archaeal lineages in marine sponges around the world (Moitinho-Silva et al., 2014;Rodríguez-Marconi et al., 2015;Verhoeven et al., 2017). Members of the phyla Actinobacteria, Acidobacteria, Cyanobacteria, Chloroflexi, Proteobacteria, Bacteroidetes, and Firmicutes have been described in association with different marine sponges (Jackson et al., 2012;Giles et al., 2013;Bayer et al., 2014). However, the marine bacterial communities can vary in different sponges with respect to both microbial richness and diversity.
The Persian Gulf, a small, shallow, semi-enclosed body of water bordered by the Arabian Peninsula and Iran, is a unique and greatly underexplored marine ecosystem. There are about 55 sponge genera recorded in the Persian Gulf (Najafi, 2012). However, next-generation sequencing technology has not been applied to identify the sponge-associated bacterial communities in Iran and the Middle East.
The present study was aimed at characterizing the bacterial community associated with the marine sponge species collected from the Persian Gulf, Iran using 454 pyrosequencing.

Sponge Sampling
Sponge sampling was performed in May to September 2016 at a depth of 2-3 m offshore Bushehr, Persian Gulf, Iran by SCUBA diving. In this study, the sponges were living in an area where it was exposed to light. Sponge samples were placed in sterile plastic Ziploc bags containing seawater and immediately transported to the Persian Gulf of Marine Biotechnology Research Center. Sponge tissues were rinsed with 0.22-µm-membrane-filtered seawater (FSW) to remove exogenous materials and loosely attached microbes (Jackson et al., 2012). The samples were stored at −80 • C until further processing.

Ethics Statement
In this study, the sponges collected did not involve endangered or protected sponge species. No specific scientific research permission was required to collect the sponges from the Persian Gulf.

Sponge Identification
Sponge taxonomic identifications were confirmed with a combination of multilocus DNA markers. The cytochrome oxidase subunit 1 (COI) and partial 28S rDNA fragments (ITS) were amplified using specific primers as previously reported . PCR amplifications were carried out on a thermal cycler PeQlab, peqSTAR 96X Universal Gradient, Germany under the following conditions: 94 • C for 30 s; followed by 35 cycles of 94 • C for 5 s; 50 • C for 5 s; 72 • C for 12 s; followed by 72 • C for 1 min . PCR products were purified and sequenced by Macrogen Inc. (Seoul, Korea).
Also, morphological and spicule examination was carried out by Dr. Yusheng M. Huang (National Penghu University of Science and Technology, Taiwan). List of sponge species collected from different locations of the Persian Gulf is shown in Table 1.

Metagenomic DNA Extraction From Sponges
Frozen sponge tissues were defrosted and washed with sterilized and filtered seawater. Then, they were cut into small pieces (about 1 cm 3 ) and ground to fine powder under liquid N 2 using a sterile pestle and mortar (Jackson et al., 2012;Bayer et al., 2014). DNA was extracted using a hexadecyltrimethylammonium bromide (CTAB) method. Briefly, a subsample of approximately 100 mg of each sponge tissue was suspended in lysis buffer [100 mM Tris, 100 mM EDTA, 1.5 M NaCl (w/v), 1% CTAB (w/v), 2% SDS (w/v)]. Then they were disrupted in the presence of factors such as proteinase K and extraction buffer containing chloroform: isoamyl alcohol (24:1). Samples were washed with a solution of phenol-chloroform in a few steps for DNA purification. Finally, DNA was precipitated with sodium acetate (3 M, pH 5.2) and isopropanol, then washed in 70% ethanol, dried and redissolved in TE buffer (Jackson et al., 2012;Schmitt et al., 2012). Metagenomic DNA was qualified by agarose gel electrophoresis [1% w/v agarose in Tris-acetate-EDTA (TAE) buffer]. The quantitative assessment of the isolated DNA was carried out using a NanoDrop 1000 spectrophotometer (Thermo Scientific, Wilmington, DE, USA). The high-quality DNA was stored at −20 • C until use (Schmitt et al., 2012;Jasmin et al., 2015).

Pyrosequencing of Barcoded 16S rRNA Gene Amplicons
In this study, the universal primers 27F and 518R were used to amplify a ∼ 400 bp fragment of the bacterial 16S rRNA gene targeting the V1 to V3 hyper-variable regions. These regions were amplified using primer sets (V1-27F: 5 ′ -X-MID-GAGTTTGATCMTGGCTCAG-3 ′ and V3-518R: 5 ′ -X-MID-WTTACCGCGGCTGCTGG-3 ′ ), in which X indicates the  Jeong et al., 2013Jeong et al., , 2014. This approach allowed for the mixing of multiple samples in parallel and re-sorting the sequences into order (Nam et al., 2011;Jackson et al., 2012). PCR amplification was performed in a volume of 50 µL containing 3 mM of MgCl 2 , 0.2 mM of dNTPs, 2.5 U of Pfu Turbo DNA polymerase (Stratagene, La Jolla, CA, USA), 1X Pfu reaction buffer, 0.1 µM of each pair of barcoded primers, and 20 ng of metagenomic DNA (Lee et al., 2011;Gao et al., 2014a). PCR was conducted in a Thermal Cycler (Applied Biosystems ABI Perkin Elmer 9600 GeneAmp) using the following conditions: an initial denaturation at 95 • C for 5 min, followed by 35 cycles of denaturation at 95 • C for 60 s, annealing at 55 • C for 60 s, extension at 72 • C for 60 s, and a final extension at 72 • C for 5 min (White et al., 2012). The PCR amplicon libraries were purified using the NucleoSpin R Gel and PCR Clean-up (Macherey-Nagel, Germany) and quantified using NanoDrop 1000 spectrophotometer (Thermo Scientific, Wilmington, DE, USA). Pyrosequencing was performed through a GS FLX Titanium system (454 Life Sciences) according to the manufacturer's instructions (Roche, Germany) by a commercial sequencing provider (Macrogen, Seoul, Korea).

Processing of 454 Tag Sequences Data
In the present study, the sequences generated from pyrosequencing were analyzed as previously reported (Jeong et al., 2014). The low-quality sequences were filtered from the raw reads using Trimmomatic v0.30 (Gaikwad et al., 2016). Briefly, sequences with a read length of less than 172 bp or with mismatches on primer or barcode, a quality score of less than 25 with ambiguous bases N and homopolymers longer than 6 nucleotides were removed from further analysis (Gao et al., 2014a;Moitinho-Silva et al., 2014;Rodríguez-Marconi et al., 2015;Gaikwad et al., 2016). Final sequences that passed the quality checks were then assigned to respective samples based on the barcodes and subjected to Denoiser to increase the accuracy of the sequence processing (Gao et al., 2014a). The sequences were merged into one file and clustered into operational taxonomic units (OTUs) using the Quantitative Insights Into Microbial Ecology Toolkit (QIIME), version 1.8.0. (Caporaso et al., 2010). Chloroplast and mitochondria sequences and chimeric reads were excluded from downstream analyses. In this study, the chimeric reads were removed through the CD-HIT-OTU program (http://cd-hit.org). Taxonomic assignment of representative OTU sequences was performed using the UCLUST (version 1.2.22) taxonomy assigner method (Edgar, 2010) against the SILVA release 119 database as a reference.
To have more information about the strains, each OTU was compared to the most closely related 16S rRNA gene sequences from the NCBI nucleotide databases using BLAST search. Phylogenetic analysis was inferred by using the Maximum Likelihood method based on the Tamura-Nei model (Tamura and Nei, 1993). Also, the evolutionary analysis was conducted using MEGA version 7.0 (Kumar et al., 2016).

Analysis of Microbial Richness and Diversity
Taxonomic abundance was presented in the phylum, class, order, family, and genus. Alpha diversity metrics were computed (observed species, Chao1, Shannon and inverse Simpson) among sponge samples using the QIIME package with a step size of 100 and 100 repetitions per step. These indices were presented to evaluate the richness and evenness of the associated bacteria within each sponge sample (Caporaso et al., 2010). To show whether the number of reads used in the analysis was sufficient in identifying species/OTU, rarefaction curves were calculated using the QIIME script alpha_rarefaction.py. (Caporaso et al., 2010). Good's coverage index was calculated as C = 1-(s/n), where "s" is the number of unique OTUs and "n" is the number of individuals in the sample (Naim et al., 2014). The beta diversity among the microbial communities in different sponges was evaluated using UniFrac analysis and the QIIME package (Caporaso et al., 2010). The phylogenetic tree was constructed with the FastTree program to cluster the samples by an unweighted-pair group method with arithmetic mean (UPGMA) using average linkages (Nam et al., 2011;Giles et al., 2013). Also, principal coordinate analysis (PCoA) plots were provided using the QIIME to visualize the effect of the microbial community on structuring the diversity in different sponges (Caporaso et al., 2010).

Sponge Taxonomic Identification
In this study, 10 sponge samples were collected from three different sites in the Persian Gulf to evaluate their bacterial diversity ( Table 1). Taxonomic identification of the sponges using a combination of multilocus DNA markers (28S rDNA and COI mtDNA) showed that the sponges belong to families Suberitidae, Pseudoceratinidae, Chondrillidae, Chalinidae, Halichondriidae, unclassified Dictyoceratida, and Irciniidae (Table 1).

Bacterial Richness and Diversity Analyses
A combined total of 134,495 raw pyrosequencing reads of the bacterial 16S rRNA gene fragment comprising 52,835,851 bases were obtained from the sponge samples. Trimming and quality filtering of the raw reads derived 68,628 high-quality sequences. These denoised sequences were clustered at a 97% similarity into 724 unique OTUs. The highest number of OTUs was obtained from the sponge Dictyoceratida sp. (S11), representing 165 OTUs and minimum of 29 OTUs from the sponge Suberites diversicolor (S02). In this study, 7.43% of 16S rRNA gene fragments were unassigned OTUs at the phylum level. To have more detailed confirmation of the OTU representative sequences corresponding to unassigned, we extracted the sequences from each OTU and performed a BLAST search to check the sequences corresponding to the conserved region of the target region (16S rDNA). But the closest known sequences had less than 89% similarity rate (Supplementary Table 1) and the phylogenetic tree was not informative. Therefore, these OTUs were not included in subsequent analyses. The total number of reads retrieved and OTUs from each sponge are shown in Table 3.
Rarefaction plots were constructed based on OTUs at a 97% sequence similarity cut-off value by the QIIME package. Alpha rarefaction curve showed that a reasonable number of reads have been used in analysis and identifying species/OTU. The sequencing depth of the sponge samples indicated the microbial communities were very well sampled (Figure 1). However, additional reads may be required to discover more OTUs for the samples such as S02 to show their total bacterial diversity. OTU based alpha diversity measures and Chao1 estimation of species richness revealed the highest richness of bacterial species in the sponge S11 and lowest in the sponge S02 (Table 3). Microbial community diversity, Shannon and Simpson indices displayed the highest community diversity and obviously distinguished the sponge S11 when compared to other sponge species collected from the same location (Table 3). In this study, mean and s.d. Expected richness (Chao1) was 80.242 ± 37.647, non-parametric Shannon (H ′ ) was 3.28 ± 1.070 and Simpson (D) was 0.774 ± 0.141. Also, the obtained average of the Good's coverage index was 99.8 % ± 0.001 for all the sponge species.
Diversity sponge-associated bacteria at the lower taxonomic levels showed that 45 classes and 87 orders were recovered from all datasets. Cyanophyceae (in Cyanobacteria), Gammaproteobacteria, Alphaproteobacteria were the major classes in the sponges, making up 43.78, 15.82, and 7.79%, respectively. The bacterial communities in sponges at the order level were heavily loaded with Synechococcales (43.78%) in Cyanobacteria and Rhodobacterales (16.32%) and Rhizobiales (9.42%) in Proteobacteria. At all the three levels of taxonomic classification, the bacterial community in the sponge S11 was more diverse than those associated with other sponges. Interestingly, sponge S11 has contained more rare bacteria compared with other 9 sponges. In addition, the bacterial compositions in the same sponges (e.g., S05 and S06) from similar sites did vary substantially (Figure 2).

BLASTN and Phylogenetic Relationships of Highly Abundant OTUs
To have more information about the sponge-associated bacteria strains, the representative sequences of each OTU were compared against the nucleotide database in GenBank. In this study, the first 27 OTUs with average proportions of more than 0.5% among sponge samples are shown in Figure 3. The most predominant OTU was OTU_denovo0 in the phylum Cyanobacteria, which accounted for proportions of 16.60% in the sponge-associated bacteria and shared the highest similarity (99%) to Synechococcus rubescens ( Table 4). The second most dominant OTU was OTU_denovo2 in the phylum Proteobacteria, which accounted for proportions of 5.30% and shared 93% identity with the closest relative Proteobacteria bacterium Sinobacterium caligoides ( Table 4). The phylum Proteobacteria was divided into more than 100 OTUs mainly belonging to Gammaproteobacteria, Alphaproteobacteria, and Betaproteobacteria. Following these symbionts in terms of relative abundance were OTU_ denovo1 in the phylum Cyanobacteria, OTU_denovo4 in Acidobacteria (read count: 1614), and OTU_denovo9 in Chloroflexi (read count: 1472) (Figure 3). Another highly abundant OUT in the phylum Cyanobacteria was related to OTU_denovo11 that affiliated with Prochlorococcus marinus with an identity of 97%. Details of all cyanobacterial and proteobacterial sequences included in the study are summarized in Table 4. Also, the evolutionary analysis of these two dominant phyla involved 127 nucleotide sequences and are shown in Figure 4.

Clustering of Sponges According to Bacterial Diversity
Weighted UPGMA tree from the Unifrac analysis was constructed to show the relationships between different sponges according to their bacterial communities ( Figure 5). In this study, sponge S11 was always distinct from the others in all trees. Also, sponges S4 and S9 exhibited the same topology and the closest distance to each other. The association among other sponges was not noticeable in the Unifrac UPGMA tree topologies. Because the bacterial profiles in each sponge were highly different from each other. Further, the difference in the bacterial community structure of the sponge samples was evaluated using Principal coordinate analysis (PCoA) plot based on the unweighted unifrac distance metrics (Figure 6). Sponge  Table 1 for sponge abbreviation. S11 formed a distinct clade and showed notable differences from other sponges. Three pairs of sponges belonged to the same species (Table 1). But each pair was clustered in different clades except for the species Chondrila sp. (S04 and S09). Also, only the sponges S02, S10, and S13 contained Cyanobacteria as a major phylum.

DISCUSSION
Recently, marine sponges have been a major target of different studies due to their abundant and diverse microbial communities, ecological roles, production of novel bioactive natural compounds, and biotechnological significance (Menezes et al., 2010;Alex and Antunes, 2015). However, our understanding of bacteria-sponge interactions, nature, and diversity of bacteria associated with marine sponges is still incomplete. There are significant gaps in research on the bacterial composition, function, and maintenance of the symbiotic relationships (Menezes et al., 2010;Verhoeven et al., 2017). Deep sequencing approaches such as 454 tag pyrosequencing can be used to explore the microbial diversity of sponges with high efficiency rather than was possible by clone library construction and Sanger sequencing methods (Gao et al., 2014a;Moitinho-Silva et al., 2014).

Taxonomic Richness and Bacterial Community Diversity
In this study, the number of unique OTUs (724) is in accordance with a study on 12 different marine sponge species from the Atlantic coast, where only 686 OTUs at 97% sequence similarity reported (Alex and Antunes, 2015). Here, the low diversity of OTUs in bacterial communities might be related to the bias of selected primers (V1-V3), in comparison to V3-V4 and V6 regions in other studies (Lee et al., 2011;Gaikwad et al., 2016;Souza et al., 2017). It will be led to inefficient amplification of the bacterial 16S rRNA genes (Gao et al., 2017). On the other hand, differences in sampling depth and the water temperature may affect the number of OTUs in our sponge-associated bacterial communities than other studies (Jackson et al., 2012;Alex and Antunes, 2015;Souza et al., 2017). Another possibility may be because of the difference in using various bioinformatics patterns to analyze the sponge data. In CD-HIT-OTU program, when defining the OUT, there is a step to remove a cluster with fewer reads, such as a singleton or doubleton, as noise. This is not the same OTU picking method in MOTHUR or QIIME. It looks like the total number of OTU appears to be low when removing low size clusters. In the present study, the cutoff X value defined as the low size cluster was applied as 7, and the cluster consisting of less than 7 reads were removed without being picked by OTU. At this time, the total number of reads removed is 46,318. In this study, alpha rarefaction curve showed that a reasonable number of reads have been used in analysis and identifying species/OTU. The bacterial richness estimation reported in the Persian Gulf sponges ranged from 29 to 165 OTUs. It is concordant with the range observed (29-370 OTUs) in a study targeting 12 marine sponge species sampled from the Atlantic coast (Alex and Antunes, 2015). In contrast, higher bacterial richness (570-3,013 OTUs) was found in the sponge-specific bacterial communities in Irish waters (Jackson et al., 2012) and also in the sponges from the Red Sea (251-444 OTUs) (Moitinho-Silva et al., 2014). Also,  showed contributing marine sponges to the total bacterial diversity of the world's oceans, with a bacterial richness accounted for 50-3,820 OTUs in each sponge . In the present study, diversity index and number of OTUs were higher in sponge S11 compared to other sponges. Sponge S11 was obviously distinguished from the other sponges, in terms of containing taxa from 15 different bacterial phyla and candidate phyla. At all the three levels of taxonomic classification, the bacterial community in the sponge S11 was more diverse than those associated with other studied sponges. Interestingly, sponge S11 has contained more rare bacteria compared with other sponges.
Taxonomic identification of sponge S11 using a combination of multilocus DNA markers showed that the sponge belongs to the genus of Dictyoceratida sp. Many studies have considered the order Dictyoceratida as the largest producer of new marine natural products, contributing more than 20% of all spongederived novel compounds (Mehbub et al., 2014(Mehbub et al., , 2016. Different bacterial communities associated with the order Dictyoceratida are reported to produce a wide range of natural compounds with a variety of biological activities (Thakur et al., 2005;Mehbub et al., 2014Mehbub et al., , 2016. Also, Sponge S11 formed a distinct clade in the weighted UPGMA tree from Unifrac analysis and showed notable differences from other sponges. Interestingly, in this study, the bacterial compositions in the same sponges (S05 and S06) from similar sites did vary substantially. There was a 100 m distance between the sampling sites of these two sponges. This finding confirms the previous studies, indicating that two cohabiting sponges may have different bacterial signatures (Jasmin et al., 2015;Jeong et al., 2015).
In the present study, a small number of OTUs (7.43%) were recorded as unassigned at the phylum level, after quality filtering and removal of chimera. When alignment is performed with UCLUST in the Reference DB, it means that there are no more than 90% references to match (OTU representative sequence). Unassigned OTUs may be a sequencing error element (Chimera, etc.) that could not be removed in the previous step (OTU picking step). Because there is no match result due to high cut off value, 90%. The frequency of unassigned OTUs was also maintained, after re-analysis of the OTU representative sequences corresponding to unassigned in a BLAST search for the conserved region of the target region (16S rDNA). These OTUs were not included in further analyses, but their presence is noteworthy. However, this number was much lower than those reported from Florida (White et al., 2012) and Indonesia (Cleary et al., 2013), where 36 and 34% of OTUs could not be assigned to any bacterial phylum, respectively.

Community Composition of Sponge-Associated Bacteria
One of the highlights of our study was the high frequency of Cyanobacteria in the Persian Gulf sponges, contributing 44.22% of the total phylum-level diversity. Also, the presence of Cyanobacteria was confirmed in all the sponges studied with various abundance. Cyanobacteria were also the dominant phylum in the sponge samples reported from other tropical and subtropical regions (Alex et al., 2012;Gao et al., 2014aGao et al., , 2017Regueiras et al., 2017).
Within this phylum, our sponge samples contained a high proportion (97.37%) of free-living Synechococcus. BLAST search further revealed the dominant OTU in sponge samples and showed a high similarity with previously reported freshwaterspecific species "Synechococcus rubescens" isolated from the deep subalpine lakes (central Europe), Lake Biwa (Japan), Lake Balaton (Hungary), and the Baltic Sea (Ernst et al., 2003). The Cyanobacteria Synechococcus has considered as an autotrophic plankton community and a substantial fraction of marine primary production (Flombaum et al., 2013), as sponge feeding on Cyanobacteria has been extensively confirmed (Pile et al., 1996;Hadas et al., 2009). This genus of marine bacteria has widely distributed in many ocean regions, covering both polar and high-nutrient waters (Flombaum et al., 2013). Furthermore, the cyanobacterial Synechococcus lineage is believed to have several ecotypes that are adapted to different environmental conditions including light, temperature, nutrients, and chlorophyll a concentration (Flombaum et al., 2013). It is possible that the core OTUs illustrate bacterial sponge ecotypes that are matched to the niche sponges and are probably environmentally transmitted (Schmitt et al., 2012). As Cyanobacteria are the center of carbon fixation and provide necessary nutrients to photosynthetic sponge hosts (Taylor et al., 2007), the high abundance of Cyanobacteria further indicated the specific roles of photosynthetic bacteria and their profitability in sponge biology (Lemloh et al., 2009;Alex et al., 2012). The Persian Gulf sponges were recurrently exposed to light in low depth (<3 m) in this study. Therefore, the prominence of photosymbionts was predictable in these sponges. In spite of the important role of Cyanobacteria, they are typically absent in Antarctic sponges (Rodríguez-Marconi et al., 2015). It seems that environmental factors such as temperature, salinity or nutrient levels might impact the composition of bacterial community structures in different sponges (Giles et al., 2013;Cuvelier et al., 2014).
The co-evolution and functional aspects of spongecyanobacteria associations have not been revealed in details. Sponges may acquire their cyanobacterial symbionts by vertical, horizontal or combined transmission routes (Thacker and Freeman, 2012). The genome-level research on cyanobacterial symbionts of sponges showed that they have general and specific adaptations to life within the sponge host in comparison with free-living cyanobacteria (Gao et al., 2014b;Burgsdorf et al., 2015). Sponge symbionts have adapted mechanisms to actively seek out by their sponge hosts (Webster and Thomas, 2016). For example, cyanobacterial symbionts contains large number eukaryotic-like domains, such as ARs. These domains may be involved in avoiding digestion by the sponge host (Gao et al., 2017).
Lifestyle evolutionary and functional studies on other functions enriched and depleted in cyanobacterial symbionts of sponges compared to members of the closely related freeliving strains reveled the precise and smart adaptation of cyanobacteria to live in full of challenge sponge's intercellular environment. Sponge amoebocytes may not actively distinguish between food bacteria and their cyanobacterial symbionts (Webster and Thomas, 2016). Hence, the depleted genes involved in biosynthesis of LPS O antigen in cyanobacterial symbionts of sponges produces a defense mechanism against sponge predation and phage attack (Burgsdorf et al., 2015;Webster and Thomas, 2016). In functional profile studies of bacterial symbionts with sponges, common functions in similar niches were found, indicating functionally convergence of symbionts in the divergent hosts (Fan et al., 2012;Rua et al., 2015). The biological and ecological roles of these functional equivalences may be of general importance for the adaptation of cyanobacterial symbiont to the sponge host environment and other symbiotic interactions.
In the present study, Proteobacteria was the second most abundant phylum in the microbiome of sponge-associated bacteria. Our study is in accordance with other studies showing  proteobacteria as one of the most diverse phyla of spongeisolated bacterial communities, irrespective of the habitat (Menezes et al., 2010;Schmitt et al., 2012;Jeong et al., 2014Jeong et al., , 2015Gao et al., 2017). In Proteobacteria phylum, OTUs were mainly affiliated with Gammaproteobacteria, Alphaproteobacteria, Betaproteobacteria classes, respectively. Multiple studies have highlighted the presence of Gammaproteobacteria in marine invertebrates, such as sponges (Menezes et al., 2010;Giles et al., 2013;Graça et al., 2015;Rodríguez-Marconi et al., 2015), corals (Sun et al., 2014), and oysters (Garnier et al., 2007). In this study, Gammaproteobacteria accounted for 15.82% of the total bacterial community, including mainly isolates belonging to genus Nitrosococcus. Different species of the genus Nitrosococcus are well known as ammonia-oxidizing bacteria (AOB), inducing the process of nitrification in different sponges and removing the ammonia excreted by the sponge host (Hentschel et al., 2006;Gao et al., 2017). The combined action between AOB and nitrite-oxidizing bacteria (NOB), such as members of the phylum Nitrospirae might then be responsible for the conversion of ammonia to nitrate in sponges (Bayer et al., 2007;Han et al., 2013). In the present study, nitrite-oxidizing phylum Nitrospirae (0.19%) was identified as one of the minor phyla in the Persian Gulf sponges. Our study is in accordance to other pyrosequencing studies in which Nitrospirae constituted a small portion of reads, accounted for 0.01-3% among several sponge species (Webster et al., 2010;Lee et al., 2011;Bayer et al., 2014;Gaikwad et al., 2016).
Chloroflexi was ranked as the third most abundant group (8.67%) using our sequencing approach. In this study, sponge 04 (Chondrilla sp.) harbored the largest proportion of Chloroflexi in comparison to other sponges. The presence of this phylum was previously reported in the Mediterranean sponge Chondrilla nucula using a clone library of 16S rRNA gene sequences (Thiel et al., 2007). The Chloroflexi is one of the most common and diverse bacterial phyla associated with a wide range of sponges, with many sponge-specific lineages detected (Schmitt et al., 2011;Hardoim et al., 2012;Jeong et al., 2013;Gao et al., 2017). Different studies have revealed the important role of Chloroflexi in nutrition, defense (Hardoim et al., 2012) and carbon fixation through photosynthesis in marine sponges (Brück et al., 2010), indicating its autotrophic lifestyle. In the present study, the Chloroflexi OTUs were closer to the autotrophic lineages.
Actinobacteria was another phylum inhabiting the Persian Gulf sponges. This phylum has been widely reported in marine sponges (Schmitt et al., 2012;Giles et al., 2013;Bayer et al., 2014;Cuvelier et al., 2014;Naim et al., 2014). Different studies have been considered Actinobacteria as an important source of bioactive natural products (Izumikawa et al., 2010;Pimentel-Elardo et al., 2010;Abdelmohsen et al., 2012), protecting the sponge hosts against pathogens (O'Connor-Sánchez et al., 2014). There is a possibility that this phylum may provide new opportunities for novel marine drug discovery.
The candidate division TM7, also known as phylum candidatus Saccharibacteria, is a highly ubiquitous and uncultured phylum of bacteria, described through environmental 16S rRNA gene sequence and genome data only (Ferrari et al., 2014). The existence of this phylum has widely been reported from different sponges (Webster et al., 2010;Lee et al., 2011;Schmitt et al., 2012;Gao et al., 2014a;Montalvo et al., 2014;Alex and Antunes, 2015;Gaikwad et al., 2016). In the present study, TM7 was found at very low abundance in some sponge species. It highlights the importance of deep sequencing technology for detection of the sponge-associated uncultivated bacteria and rare microbial groups in sponges. Otherwise, they might have been not discovered by other routine molecular methods (Lee et al., 2011). However, because of the lack of cultivated representatives and minimal genomic sampling knowledge on the metabolism and biological activities of this enigmatic group has been remained unclear (Ferrari et al., 2014).
The candidate phylum Poribacteria is a sponge-specific phylum that has been widely detected and described in a variety of sponge species (Lafi et al., 2009;Cleary et al., 2013). This candidate phylum has typically been reported in high microbial abundance (HMA) sponge microbiomes (Hochmuth et al., 2010) and considered as "indicator species" for these group of sponges . However, unexpectedly, no Poribacteria was found in sponges of the present study.
One reason for the lack of Poribacteria in the Persian Gulf sponges could be related to the use of different primers and bioinformatics pipelines (Souza et al., 2017). Also, some studies have shown that Poribacteria has not been associated with the orders such as Halichondrida, Dictyoceratida, and Haplosclerida (Lafi et al., 2009;Jeong et al., 2013). Noteworthy, these orders have been among the six orders found in this study. In addition, in contrast to other studies conducted on sponges (Kennedy et al., 2008;Montalvo et al., 2014;Alex and Antunes, 2015), the Verrucomicrobia and Spirochaetes phyla were not observed in the bacterial community of the Persian Gulf sponges.
Our study showed that more than 90% of the microbial groups observed in the Persian Gulf sponges were also represented in other studies conducted on the seawater samples (Lee et al., 2011;Schmitt et al., 2011;Gao et al., 2014a;Alex and Antunes, 2015). Nitrospirae, Tenericutes, Armatimonadetes, and the two candidate phyla BD1-5 and TM6 were the only bacterial communities reported exclusively in association with the marine sponges.
Our results showed that the representative OTUs sequences were related to sequences mainly from marine environments such as seawater and marine sediments. Interestingly, sponges had a very low proportion of the isolation source in the closest relative bacteria. Also, some sequences were related to other marine invertebrates such as abalone, octocoral, sea cucumber, oyster, ascidian, starfish, etc. This finding supports the hypothesis of possible environmental acquisition and/or horizontal transmission of bacteria (Moitinho-Silva et al., 2014;Alex and Antunes, 2015). It seems that innate immune system in sponges is responsible for differentiation between symbionts and food microbes (Müller and Müller, 2003). As a result, some overlap between the bacteria in the surrounding seawater and marine sponges would occur. A further deep-sequencing approach is needed to improve our knowledge about the nature of the bacterial specificity among the marine sponges (Alex and Antunes, 2015).

CONCLUSION
In this study, evaluation of 16S rRNA gene amplicon tag pyrosequencing data showed a complex structure of the previously uncharacterized bacterial communities associated with the Persian Gulf sponges (class Demospongiae). OTU-based description of the bacterial communities exhibited altogether 17 different bacterial phyla, containing 14 formally described phyla and three candidate phyla. The dominance of Cyanobacteria may suggest an ecological importance of this phylum in the Persian Gulf sponges. More specifically, most of the bacterial symbionts were previously described as significant participants in the carbon, nitrogen cycle, and chemical defense of the studied sponges. Sponge S11 was highly diverse in comparison to other studied sponges and contained more rare bacteria. More research is needed to fully understand the composition of the Persian Gulf sponge-associated prokaryotic communities, specifically the functionality of these specific microbiomes as an important part of the marine ecosystem.

AUTHOR CONTRIBUTIONS
IN and AN designed the research. AN and MM performed the experiments. AN analyzed the data. AN and IN wrote the manuscript with contributions from all authors.