The R/V Tansuoyihao TS-7 cruise was carried out in April 2018 in the South China Sea (SCS). A single sponge specimen (hereafter referred to as SQW35) was collected by dive No. 35 of the manned submersible Deep-Sea Warrior from a north SCS site (18° 41′ N, 113° 22′ E) at ~983 m depth. Upon arrival at the main deck of the R/V, the specimen was rinsed several times with 0.22-μm membrane-filtered seawater to remove loosely attached microbes and debris. The cleaned samples were placed into separate sterile plastic bags, transported to the laboratory, and stored at −80°C. Two pieces of tissue (~0.5 cm³ per piece) from an inner part of the sponge body as technological replicates (SQW35-1 and SQW35-2) were separately put into 1 ml of DNA extraction buffer (50 mM Tris–HCl, 40 mM EDTA, 500 mM NaCl, 0.75 M sucrose, pH = 8) and fully cut into tiny pieces with sterile scissors. Total suspensions of two replicates were separately subjected to DNA extraction with the PowerSoil^® DNA Isolation Kit (MoBio Laboratories, Carlsbad, CA, USA) following the default experimental procedure. The extracted DNA was quantified by the Qubit^® dsDNA HS Assay Kit with Qubit 2.0 fluorometer (Invitrogen, Carlsbad, CA, USA) and stored at −80°C for further processing.

One hundred nanograms of genomic DNA for each replicate were randomly fragmented to ~350 bp by Covaris M220 Focused-ultrasonicator (Covaris, Massachusetts, USA). Two metagenomic libraries were constructed using the TruSeq^® Nano DNA LT Kit (Illumina, San Diego, CA, USA) and sequenced on a HiSeq2500 platform (Illumina, San Diego, CA, USA) to produce 2 × 150 bp paired-end reads. Raw reads were trimmed and filtered by Fastp v0.20.0 (Chen et al., 2018) with default settings. The data were further processed with FastUniq v1.1 (Xu et al., 2012) to remove duplicated reads. 16S rRNA gene sequences were predicted from clean reads using the rRNA_HMM program for analysis (subsequently referred to as 16S miTags) (Huang et al., 2009) and subjected to QIIME1 pipelines for microbial community analyses (Caporaso et al., 2010). Briefly, the 16S rRNA reads with a shared similarity of 97% were clustered to operational taxonomic units (OTUs) using UCLUST. The longest read of each OTU was selected as the representative for further taxonomic classification with the SILVA 138 database as a reference (Quast et al., 2013). OTUs that were annotated as chloroplasts, mitochondria, and eukaryotes were excluded. False positive OTUs that were not assigned to any taxa were also removed. Finally, the taxonomic relative abundance of the microbial communities in two replicates was calculated at the phylum level.

Clean data of two technical replicates, SQW35-1 and SQW35-2, were separately assembled using both SPAdes v3.1.1 with the recommended settings for metagenomes (Bankevich et al., 2012) and MEGAHIT v1.2.6 with default settings (Li et al., 2015). Four assembled datasets (SQW35-1-spades, SQW35-2-spades, SQW35-1-megahit, and SQW35-2-megahit) were subjected to MetaWRAP v1.2.1 (Uritskiy et al., 2018) for genome binning with default settings, which generated 38 MAGs. These retrieved MAGs were then de-replicated together using dRep v.2.2.3 (Olm et al., 2017), which resulted in 14 high-quality MAGs. CheckM v1.0.11 was used to evaluate the completeness and contamination of the MAGs (Parks et al., 2015). Those with completeness > 50% and contamination < 10% were retained. Taxonomic annotation and relative evolutionary distance (RED) calculation of the MAGs were carried out by GTDB-tk v0.2.2 (Chaumeil et al., 2019). The nomenclature of the MAG was based on the List of Prokaryotic names with Standing in Nomenclature (LPSN) (https://lpsn.dsmz.de/).

The relative abundance of a MAG was estimated as the percentage of mapped metagenomic reads with the removal of eukaryotic reads. First, metagenomic contigs with a length of >1,000 bp were subjected to EukRep v.0.6.7 to identify eukaryotic contigs (West et al., 2018). Then, the metagenomic reads were mapped to the eukaryotic contigs by Bowtie2 (Langdon, 2015), and the mapped reads were removed from the metagenomes. Finally, the remaining reads were subjected to coverM v0.2.0 (https://github.com/wwood/CoverM) to calculate the relative abundance of sponge-associated MAGs. To investigate the global distribution of a sponge symbiont, the V4 region of the 16S rRNA gene extracted from the MAG was searched by BLASTN (e-value cutoff: 1e-05) against the datasets of 16S rRNA gene sequences in the Sponge Microbiome Project (SMP) (Moitinho-Silva et al., 2017b) and the Deep-sea Sponge Microbiome Project database (D-SMP) (Busch et al., 2022). Target sequences with >96% identity were retained to calculate their relative abundance in respective sponges.

Close relatives of symbiotic MAGs were obtained by the phylogenomic analysis of the GTDB-tk program, and their genome sequences were downloaded from the NCBI and Integrated Microbial Genomes (IMG) databases. Open reading frames (ORFs) of our MAGs and reference genomes were predicted by Prodigal v2.6.2 (Hyatt et al., 2012). The predicted genes were annotated by KofamScan v1.1.0 (Aramaki et al., 2020) against KEGG databases. Eukaryotic-like protein (ELP) domains were annotated by PfamScan script against the PFAM database (El-Gebali et al., 2019). ORFs were searched by BLASTP v2.5.0 against the COG_2019_v11.0 database to find genes encoding a transposase. Clustered regularly interspaced short palindromic repeats (CRISPRs) and Cas proteins were predicted using the online CRISPRminer2, a toolkit comprising CRISPRCasFinder, CRT, and PILER-CR programs (Zhang et al., 2018). Pairwise ANI values of MAGs were calculated by PyAni v0.2.10 (Pritchard et al., 2016) with the ANIb model. AAI values of MAGs were calculated by CompareM v0.0.23 (https://github.com/dparks1134/CompareM).

The coxI (cytochrome c oxidase subunit I) gene of the deep-sea sponge sample was annotated from metagenome contigs by Prokka v1.13.7 (Seemann, 2014). 16S rRNA genes of symbiotic MAGs were predicted using the RNA_HMM program (Huang et al., 2009). For phylogenetic analysis of both coxI and 16S RNA genes, the targeted gene sequences were searched by the BLASTN program against the NCBI GenBank database to identify close relatives. The collected gene sequences were aligned using MAFFT v7.427 (Katoh and Toh, 2010) and trimmed using trimAl v1.4 (Capella-Gutiérrez et al., 2009). Phylogenetic trees were built using IQ-TREE v1.6.10 (Nguyen et al., 2015) with the “TIM2+F+I+G4” model. For phylogenomic analysis, the alignment of 43 concatenated conserved proteins deduced from our MAGs and reference genomes was produced using the CheckM program with default settings and further treated with trimAl v1.4 (Capella-Gutiérrez et al., 2009) to remove poorly aligned regions. The maximum likelihood (ML) tree was built using IQ-TREE v1.6.10 (Nguyen et al., 2015) with the “LG+F+R7” models. Bootstrap values of the trees were calculated based on 1,000 replicates.

Sponge metagenomic contigs in size of >5,000 bp were imported into DeepVirFinder (Ren et al., 2020), VIBRANT (Kieft et al., 2020), VirFinder (Ren et al., 2017), VirSorter (Roux et al., 2015b), and VirSorter2 (Guo et al., 2021) to identify viral genomes. CheckV v0.8.1 (Nayfach et al., 2021) was used to assess the viral genome quality, to identify and remove potential host contamination in integrated proviruses, and to match closely relative viral genomes from public databases. Viral genome completeness was estimated by searching against a database that comprises 76,262 complete viral genomes from publicly available metagenomes, metatranscriptomes, and metaviromes by IMG/M (Chen et al., 2019), MGnify (Mitchell et al., 2020), and the study of the human microbiome (Nayfach et al., 2019) and ocean virome (Gregory et al., 2019). Taxonomic assignment of positive viral contigs was performed using vConTACT2 (Bin Jang et al., 2019), which was designed to cluster the protein sequences with a RefSeq database based upon shared protein clusters. The relative abundance of each viral contig in a deep-sea sponge sample was calculated by mapping clean reads to assembled contigs using BWA-MEN with default settings (Li and Durbin, 2010).

A thick-walled tubular deep-sea sponge was collected from the South China Sea at a depth of 983 m. The top of the sponge body is a moderately sized, trumpet-shaped osculum. The body color is pale beige. Prostalia protrudes over the body by several centimeters (Figure 1A). Phylogenetic inference based on the coxI gene (1,575 bp in size) suggested that the sponge is closely related to Bathydorus laniger from the coast of California (Kahn et al., 2013) and Bathydorus spinosus in the Weddell Sea (Dohrmann et al., 2011), two known species of the genus Bathydorus in the glass sponge class Hexactinellida (Figure 1B). The coxI gene of our sponge sample shared an identity of 94.83 and 93.27% with its two Bathydorus relatives. The coxI genes had respective matching thresholds at order, family, and genus levels (45, 73, and 91%, respectively) based on a sponge identification protocol proposed by a previous study based on 37 sponge species belonging to 10 orders from South Australia (Yang et al., 2017). Thus, our sample should represent a novel sponge species in the genus Bathydorus and is preliminarily named Bathydorus sp. SQW35.

Taxonomic classification of 16S miTags demonstrated highly consistent prokaryotic communities in two technical replicates of the sponge Bathydorus sp. SQW35. The bacterial symbionts were mainly affiliated with the phyla Pseudomonadota (synonyms Proteobacteria) (Oren and Garrity, 2021), Planctomycetota, Nitrospirota, and Bacteroidota, and archaeal symbionts were dominated by the phylum Nitrososphaerota (synonyms Thermoproteota) (Supplementary Figure 1; Supplementary Table 1). Among them, the highest proportion belonged to Nitrososphaerota, which accounted for 81.81 ± 2.09% of the prokaryotic communities. Further taxonomic analysis revealed the dominance of the Nitrosopumilaceae family in Nitrososphaerota, known as AOA, which can obtain energy from ammonia oxidation and use CO₂ as a carbon source for chemoautotrophy (Wang et al., 2019). Gammaproteobacteria that mainly belonged to SOBs were the second dominant group in the sponge Bathydorus sp. SQW35 and accounted for 10.60 ± 0.85% of the prokaryotic communities. Alternatively, Nitrospirota that is involved in nitrite oxidation carbon fixation (Lücker et al., 2013) occupied 0.85 ± 0.11% of the prokaryotic communities.

In total, 14 metagenome-assembled genomes (MAGs) with >50% estimated completeness and <10% contamination were obtained from the sponge metagenomes of the same sponge individual, and five of them were more than 90% complete (Table 1). The MAGs have a genome size ranging from 0.54 to 3.63 Mbp and GC content of 31.93–52.87%. Classification by the GTDB-tk program revealed that these MAGs were affiliated with the phyla Nitrososphaerota, Pseudomonadota, Nitrospirota, Bdellovibrionota, SAR324, Bacteroidota, and Patescibacteria (Supplementary Table 2). Relative evolutionary divergence (RED) that infers a phylogenetic distance between the last common ancestor (set to RED = 0) and all extant taxa (RED = 1) can be used to establish taxonomic ranks (Parks et al., 2018). RED values indicated that only MAG B05 fell into known species and the others likely represent novel species. Among these novel species, five MAGs have RED values of 0.71–0.89, indicative of novel genera. Three MAGs with RED values of 0.64–0.69 are suggested to present novel families (Parks et al., 2018).

Phylogenomic inference also revealed that the 14 symbiotic MAGs fell into seven phyla, and the far distance of these MAGs with known species further indicated their taxonomic novelty (Figure 2A). Five MAGs (B01–B05) fell into the clade composed of AOAs in the phylum Nitrososphaerota, whereas only B03 was closely related to sponge-associated AOAs. ANIs between B03 and its close relatives Ca. Nitrosopumilus sp. ESC from the sponge Hymedesmia (Stylopus) methanophila (Haber et al., 2020) and Ca. Nitrosopumilus sp. LS from the glass sponge Lophophysema eversa (Tian et al., 2016) were 90.40 and 90.86%, respectively (Supplementary Table 3). The affiliation of MAG B14 with sulfur-oxidizing Gammaproteobacteria Gsub from the cold seep sponge Suberites sp. (Tian et al., 2017) in the phylogenomic tree was in line with their ANI of 70.80%. Nitrite-oxidizing bacterium (NOB) MAG B08 in the phylum Nitrospirota was affiliated with the genus SPGG5 and was placed together with Nitrospinae bacterium UBA9942 (Parks et al., 2018), which shared an ANI of 67.37%. MAGs B12 and B11 were classified as the phylum Bdellovibrionota and shared the ANIs of 68.07 and 67.91%, respectively, with its relatives Bdellovibrionales RBG_16_40_8 (Kauffman et al., 2018) and Bdellovibrionales bacterium SXSP01 (Xing et al., 2020).

The relative abundances of the sponge symbionts in metagenomes of the sponge Bathydorus sp. SQW35 are summarized in Figure 2B. The AOA MAG B01 recruited 70.31 ± 1.15% of the prokaryotic metagenome reads and represented the sole dominant symbiont (Figure 2B; Supplementary Table 4). The SOB MAG 14 as the second most abundant symbiont was mapped by 2.83 ± 0.40% of the prokaryotic metagenomic reads. The NOB MAG B08 could also be mapped with less abundant metagenomic reads. Relative abundances of the representative MAGs including AOA B01, SOB B14, and NOB B08 were consistent with the result of the 16S miTags analysis. In contrast, the Bdellovibrionota MAGs B11 and B12 that accounted for 0.29 ± 0.29% and 0.07 ± 0.07% metagenomic reads were not detected in 16S miTags analysis.

Our detailed phylogenomic inference revealed that the dominant symbiont AOA B01 was closely related to the Cenarchaeum clade that is now composed of four lineages including Cenarchaeum symbiosum A from the demosponge Axinella mexicana, Cenarchaeum sp. bin74s, Cenarchaeum sp. bin90o, and Cenarchaeum sp. bin143o (phylum: Nitrososphaerota) from Vazella pourtalesii, Nitrososphaerota archaeon TS (phylum: Nitrososphaerota) from Theonella swinhoei, and Candidatus Cenporiarchaeum stylissum S13–S15 (phylum: Nitrososphaerota) from Stylissa flabelliformis (Figure 3). ANIs between B01 and these reference genomes range from 68.96 to 69.97% (Supplementary Table 3). With its RED value of 0.87 (Table 1), MAG B01 should represent a new genus. Global distribution analysis against the SMP database (Moitinho-Silva et al., 2017b) revealed that the V4 region of the 16S rRNA gene in MAG B01 showed 100% identity with amplicons from three glass sponge individuals in the family Dactylocalycidae (Figure 4), with relative abundances of 23.25, 8.78, and 8.60% in the respective communities (Supplementary Table 5). B01 also shared 97–99% identities with sponge symbionts in three other sponge species, Aphrocallistes beatrix, Axinella sp., and an unnamed species, but was only 96% similar to its closest relatives from seawater samples (Supplementary Table 5). By a query against the D-SMP database (Busch et al., 2022), B01 shared more than 99% identity with symbionts in Pheronema carpenteri, Bathydorus sp., Aphrocallistes beatrix, Saccocalyx tetractinus, Rossellinae indet, and Lophocalyx sp. The relative abundances of B01 in these sponges span between 0.007 and 0.657% (Supplementary Table 5). Because the D-SMP database employed bacterial primers for amplification, the relative abundance of archaeal relatives of MAG B01 in this database is likely underestimated. Sponge symbionts can be classified into generalists (found in a wide range of sponge species) and specialists (living in a small number of sponge species) according to their distribution patterns (Haber et al., 2020). In brief, the distribution patterns against both the SMP and D-SMP databases indicate that AOA B01 is a specialist.

Phylogenomic inference showed that both MAGs B11 and B12 fell in the Bdello-group2 in the phylum Bdellovibrionota but form two separate deep branches (Figure 5). Bdello-group2 comprised of the representative Bdellovibrio predators including Bdellovibrio exovorus JSS, Bdellovibrio bacteriovorus HD100, Bdellovibrio bacteriovorus Tiberius, and Bdellovibrio bacteriovorus 109J (Li et al., 2021), suggesting that B11 and B12 may also represent bacterial predators. The Bdellovibrionota MAGs B11 and B12 have remarkably smaller genome sizes and lower GC contents than those free-living Bdellovibrionota species (Supplementary Figure 2). The comparative genomic analysis further revealed that the free-living Bdellovibrionota bacteria harbored genes encoding chemotaxis proteins and mobility systems (flagellum and type IV pili), yet MAGs B11 and B12 lacked these genes (Supplementary Table 6). Predatory Bdellovibrio bacteria invade the periplasm of bacterial prey cells by penetrating the peptidoglycan layer to form transient structures there. Penicillin-binding proteins (PBPs) are important for Bdellovibrionota to lyze the cell wall of a prey bacterium (Lerner et al., 2012). Annotation analysis showed that symbiotic Bdellovibrionota MAGs lack the PBP1C coding gene compared to free-living relatives but retain PBP1A, PBP1B, PBP2, PBP3, and PBP4 coding genes (Supplementary Table 6), which suggests that MAGs B11 and B12 may still retain the predatory potential.

We searched the three main prokaryotic autotrophic carbon fixation pathways: 3-hydroxypropionate/4-hydroxybutyrate (3HP/4HB) cycle, reductive citric acid (rTCA) cycle, and the Calvin–Benson–Bassham (CBB) cycle. The 3HP/4HB cycle was identified in sponge-associated AOA (family Nitrosopumilaceae) with the presence of key genes encoding 3-hydroxypropionyl-CoA synthetase (K18594) and 4-hydroxybutyrate-CoA ligase (K18593) (Supplementary Table 7). Almost all genes of the rTCA cycle were identified in Nitrospirota B08 (family UBA8639, genus SPGG5), including those coding for citrate synthase (K01647), fumarate hydratase (K01676), and succinate dehydrogenase (K00240), which were considered as evidence for the rTCA cycle. Furthermore, more than 90% of the genes for carrying out the CBB cycle including those encoding subunits of ribulose-bisphosphate carboxylase (K01602 and K01601) were found in the SOB MAGs (class Gammaproteobacteria). These findings indicate the existence of various carbon fixation pathways in the microbial consortium of the sponge Bathydorus sp. SQW35.

The amoABC genes (K10944, K10945, K10946) encoding ammonia monooxygenase subunits were identified in three Nitrososphaerota MAGs (B01, B04, and B05) but not in any bacterial MAGs (Supplementary Table 7). This suggests that AOA is the only symbiotic group in charge of ammonia oxidation in the sponge Bathydorus sp. SQW35. Nitrospirota B08 had the gene encoding nitrate/nitrite transporter (K02575) but lacked genes encoding nitrite oxidoreductase. In addition, the B08 has a gene set involving urea decomposition that produces ammonia for nitrification and biomass production, where urea would be imported from the extracellular milieu by the urea ABC transporter urtABCDE (K11959, K11960, K11961, K11962, and K11963) and be hydrolyzed by urease ureABC (K01428, K01429, and K01430) (Supplementary Table 7). These results further imply the potential role of NOB B08 in urea utilization.

Microbial sulfur oxidation is frequently reported in sponges to remove toxic sulfide produced by other bacterial symbionts that use sulfate as an electron acceptor (Tian et al., 2014). Gammaproteobacteria SOB B13 and B14 encode adenylyl-sulfate reductase A (K00394), adenylyl-sulfate reductase B (K00395), and sulfate adenylyltransferase (K00958) to oxidize sulfite to sulfate (Supplementary Table 7). SOB B14 also contains the soxAXYZ gene clusters (K17222, K17223, K17226, and K17227) that code for the SOX complex, confirming potential capacity for sulfur oxidation (Supplementary Table 7). Previous studies showed that taurine is a natural product of sponges, and sponge symbionts can import and utilize taurine as suggested by the presence of ABC transporter genes (tauABC) (Emura et al., 2006; Karimi et al., 2018). The SOB MAGs B13 and B14 also have genes encoding taurine ABC transporter (tauAC; K15551 and K15552), taurine-pyruvate aminotransferase (K03851), and sulfoacetaldehyde acetyltransferase (K03852), which can catalyze taurine to sulfite for further oxidation by the SOX complex.

We analyzed the distribution of ELPs including ankyrin repeats (Ank), tetratricopeptide repeats (TPRs), NCL-1, HT2A and Lin-41 repeats (NHL), fibronectin type III (Fn3) and cadherin (CAD) CUB, bacterial Ig-like domain (Big), WD40, and pyrroloquinoline quinone repeat (PQQ) in the deep-sea glass sponge microbiome (Supplementary Table 8). We found that five of the nine bacterial MAGs encoded Ank domains, including SOB B14 and NOB B08. Notably, Bdellovibrionota B11 has 12 Ank domains, which was the highest number among these bacterial MAGs. However, none of the AOA MAGs had Anks. The bacterial and archaeal MAGs, except B10, had TPR domains, especially in B08 and B09. As the most prevalent symbiont, AOA B01 encoded the maximum number of NHL domains. NOB B08 also showed the richness of NHL domains (n = 13).

One major strategy of prokaryotic defense against phages is dependence on the Clustered Regularly Interspaced Short Palindromic Repeats (CRISPR)–Cas system (Makarova et al., 2013). Here, all the MAGs were analyzed for their CRISPR–Cas systems. Almost all the symbiotic MAGs in the sponge Bathydorus sp. SQW35 has genes encoding Cas proteins (Supplementary Figure 3A). CRISPRminer2 analysis further indicated that CRISPR arrays could be found in six MAGs (B01, B08, B09, B11, B13, and B14), which represented almost all the dominant microbial groups (Supplementary Figure 3B; Supplementary Table 9). Interestingly, the dominant symbiont AOA B01 has a CRISPR array with the highest number of spacers (n = 85). The second abundant symbiont SOB B14 also has a CRISPR array with 81 spacers, followed by a CRISPR array in Bdellovibrionota B11 consisting of 51 spacers.

To understand the high prevalence of CRISPR–Cas systems in the dominant symbiotic inhabitants, we analyzed the phages in the deep-sea glass sponge Bathydorus sp. SQW35. In total, 125 contigs in SQW35-1 and 148 contigs in SQW35-2 were predicted as potential phage genomes (>5,000 bp) (Supplementary Table 10). The potential viral contigs of SQW35-1 and SQW35-2 were predicted to harbor 237 and 236 viral genes by CheckV (Nayfach et al., 2021), respectively. There are 232 positive virus sequences after redundancy removal of the contigs, which indicated a large phage community in the sponge Bathydorus sp. SQW35. Among them, the most abundant phage fragment accounted for 1.48 and 2.19% of metagenome reads, respectively, in the two replicate samples (Supplementary Table 11).