Introduction

Marine sponges (phylum Porifera) are diverse metazoans that arose in the Precambrian over 700 million years ago and are now globally distributed throughout marine and freshwater ecosystems in temperate, tropical and polar regions [1,2,3]. They possess a simple body plan divided into three layers of organisation: the outer pinacoderm, the inner choanoderm layer, and the mesohyl matrix in between [3]. As sessile filter-feeding organisms they are remarkably efficient in obtaining food from the surrounding water and thus draw down carbon, nitrogen and silicon from the pelagic environment into benthic biomass [4,5,6].

Marine sponges are well known for their wide range of interactions with microorganisms from all domains of life, from using them as food sources to housing them as symbionts [7]. Sponges are capable of differentiating between symbionts and food microorganisms [8], potentially through the recognition of specific eukaryote-like proteins (ELPs) presented by microbial symbionts [9,10,11,12] and/or immune receptors expressed by the sponge [13]. Sponges are often classified as either being of high microbial abundance (HMA), harbouring dense and often diverse microbial communities, or as being of low microbial abundance (LMA) [14, 15]. LMA sponges are typically dominated by only a few bacterial phyla, often being Proteobacteria [16]. Many bacteria found in sponges also fall into unique sponge-specific or sponge-enriched phylogenetic clusters [17, 18] as defined by being rarely found outside the sponge environment [19]. These unique clades have been proposed to have likely undergone specific evolutionary adaptations that allow them to only (or preferentially) thrive in association with their sponge host [7].

One of these unique bacterial groups, which was initially assigned to the class Betaproteobacteria (now order Burkholderiales), contains dominant and persistent members of the microbiota in various species of LMA sponges, in some instances comprising over 70% of the sponge’s microbial community [20,21,22,23,24]. First described in 2006 from the sponge Tethya aurantium [25], a few members of this bacterial group have been visualised via microscopy [20, 26, 27] and determined to be vertically transmitted from the adult sponge to the embryo [26, 28, 29]. It has also been proposed that this bacterial group co-evolved with their host [30] as each sponge species appears to harbour a distinct phylotype [20]. Despite the substantial interest in this group and its apparent broad involvement in sponge symbiosis, its phylogeny and taxonomy remain in disarray. Early studies placed this group in the family Nitrosomonadaceae (class Betaproteobacteria) based on the SILVA database 108 release [21] and assigned them to a sponge-specific cluster (SC112), which was related to other known betaproteobacterial clusters (SC110 and SC111) [22]. This classification was then superseded by placement in the betaproteobacterial order EC94 [31, 32], although in many studies these bacteria are simply described as unclassified betaproteobacteria [29, 33]. Furthermore, very little is known about this group’s functional features or role in sponge-microbe interactions. One metagenome-assembled-genome (MAG) of this group, called AqS2, has recently been described from the Great Barrier Reef sponge Amphimedon queenslandica and showed a likely heterotrophic metabolism [33].

The aim of this study is to characterise this widespread and often dominant group of sponge-associated bacteria and define evolutionary and functional properties of symbiosis. For this we employed metagenomics, in-depth phylogenetic analyses, metatranscriptomics, and fluorescence in situ hybridisation microscopy across several sponge species (Cymbastela concentrica, Crella incrustans, Crambe crambe, Scopalina sp. and Tethya stolonifera) that are known to harbour this bacterial group.

Methods

Sampling, metagenomic sequencing, genome reconstruction and annotation

Sponges were collected via SCUBA or snorkel from Bare Island, NSW, Australia; the Mediterranean Coast of Spain; and from the north-east coast of New Zealand (Supplementary Table S1). Sponge-associated microorganisms were enriched in each sample following the cell separation protocol described in Thomas et al. [9]. High-molecular-weight DNA was extracted from microbial cell pellets using the most appropriate method based on sample type and sequenced with an array of next generation sequencing technologies (Table 1). Reads from each sample were quality trimmed, de novo assembled and contigs over 1000 bp binned into MAGs (Table 2). CheckM [34] was employed to estimate genome completeness and contamination using the taxon-specific workflow. Further details are given in the Supplementary Information (File 2).

Table 1 Summary of DNA extraction, sequencing and analysis tools used to obtain MAGs and 16S rRNA gene sequences.

MAGs, including AqS2, were submitted to the Joint Genome Institute (JGI) Gold database [35] and annotated via the Integrated Microbial Genome (IMG) pipeline [36]. IMG annotations were manually refined by using the IMG “Find Genes” BLAST tool against National Centre for Biotechnology Information (NCBI) protein reference sequences. Hits close to 35% protein sequence similarity were further manually refined using the ExPASy BLAST tool (UniProtKB/Swiss-Prot databases only). An additional Pfam (v31.0) search was conducted with HMMER [37] to investigate the presence and diversity of ELPs in each MAG. Hits with bit scores less than 20 were disregarded and duplicate hits per predicted protein across multiple Pfam families were removed. Gene functions were also compared to publicly available isolate genomes, MAGs and single-cell amplified genomes in the IMG database (July 2020). Only genomes assigned by the IMG taxonomy to the class Gammaproteobacteria and betaproteobacteria (i.e. now Burkholderiales) that had ecosystem information containing the keywords “marine” as well as “oceanic” or “pelagic” were selected for comparison. This resulted in 141 genomes (Supplementary File 6) that represent relatives that are likely free living. Differences in the abundance of gene functions between the MAGs and free-living relatives were assessed by t tests with Benjamini–Hochberg correction for multiple testing.

Phylogenetic analysis

Ribosomal RNA genes were reconstructed from the metagenomic reads from each sponge species using the tools listed in Table 1. Briefly, rRNA reads were identified and then assembled to form partial- or full-length 16S rRNA genes or rRNA operons. Reconstructed genes were searched (BLASTn) against the MAGs to aid in identification. Reconstructed 16S rRNA gene sequences were also aligned with the SINA web aligner [38] and imported into the ARB software package [39] containing the SILVA 132 SSU Ref Nr99 database [40]. Alignments were further refined manually in ARB and sequences were trimmed before insertion by parsimony into the phylogenetic tree. This identified neighbouring sequences, which were then aligned with SINA to sequences over 850 bp from the sponge-specific clusters SC112, SC110 and SC111 [18] as well as sequences from type strain examples of various orders of the Gammaproteobacteria and Alphaproteobacteria (outgroup). A maximum-likelihood tree was constructed from this alignment with RAxML v. 8 [41] with 1000 bootstrap resampling and visualised using iTOL [42].

The MAGs from C. concentrica, C. incrustans, C. crambe, Scopalina sp., and T. stolonifera were analysed with the GTDB-Tk v0.2.1 Classify Workflow [43] containing the Genome Taxonomy Database (GTDB) [44] with 123,110 bacterial genomes. A genome tree and associated taxonomies were inferred using a set of 120 ubiquitous single-copy proteins and relative-evolutionary distance (RED) values were calculated. Relatedness of the MAGs was also assessed by average amino acid identity (AAI) [45] using the “AAI: Average Amino acid Identity calculator” [46].

Distribution analysis

As the group of bacteria investigated here has been reported to occur globally in a number of sponge species, we used the extensive dataset recently produced for the Sponge Earth Microbiome Project (SEMP) [16, 47] to explore its geographic distribution. For this, we first trimmed the V4 region of the 16S rRNA gene from each MAG to match the 100 bp regions of the SEMP dataset (https://github.com/amnona/SpongeEMP). These trimmed sequences were then BLASTn searched against the sequences used to construct the 16S rRNA gene phylogenetic tree shown in Fig. 1. From this, a 93% similarity threshold was found to encompass only sequences of the target group and exclude those from other orders. The trimmed V4 regions were then BLASTn searched against the operational taxonomic units (OTUs) of the SEMP database and hits with over 93% identity were matched with the relative abundance and metadata (host taxonomy and location) from Moitinho-Silva et al. [16]. As we were examining the group as a whole, relative abundances (RA%) of matched sequences were summed within samples and reported as the mean, standard deviation and range (min–max). Samples were retained for further analysis if the RA of the group was greater than 0.05%. Data were then subjected to statistical analysis (F test and one-way ANOVA).

Gene expression of MAG

Tissue of C. concentrica, C. incrustans, C. crambe, and Scopalina sp. individuals were collected via SCUBA from their respective sampling sites and brought to the surface within 20 min after being cut from the sponge individuals (Supplementary Table S2) and as part of previous studies involving Scopalina sp. [48] and C. crambe [49]. After surfacing, the sponge samples (~5–10 cm3) were immediately washed three times (<1 min) in 50 ml of 0.22 μm filtered seawater to remove planktonic microorganisms before being snap frozen in liquid nitrogen within 5 min after surfacing of the SCUBA diver. Frozen samples were transported to the laboratory. Total RNA extraction and bacterial mRNA enrichment followed the protocol established by Diez et al. [48] and Öztürk et al. [49] (see Supplementary Information File 1 for details). Expression of the coding sequences for each MAG was estimated with RSEM v. 1.3.0 [50] using default parameters and implementing the Bowtie aligner [51]. Transcripts per million (TPM) measures were used as estimates of expression [52, 53] and only genes with an average TPM values over 30 across replicates were considered. Expression of genes is reported as TPMav for those with multiple replicates. TPMsum is the combined expression of all genes in a pathway, while X − Y TPM signifies a range of expression values when multiple copies of a gene are present.

Bacterial visualisation with catalyzed reporter deposition fluorescence in situ hybridisation (CARD-FISH)

Sponges were collected from the original sampling sites (Supplementary Table S1). T. stolonifera was collected via snorkel at low tide in October 2018 and C. concentrica, C. incrustans and Scopalina sp. were collected via SCUBA in November 2018. Sponges were dissected and soaked in 4% paraformaldehyde for 4 h then washed in PBS and stored in 1:1 PBS and ethanol at −20 °C. Processing of sponge tissue was undertaken as described by Webster and Hill [54] and is described in more detail in the Supplementary Information File 1. CARD-FISH probes labelled with horseradish peroxidase (HRP) were designed for each target bacterium using the probe design tool in the ARB software package [39]. Probe specificity was assessed using the ARB probe match tool, the SILVA probe test tool and BLASTn searches against the NCBI database. Probe sequences can be found in the Supplementary Information (Table S3). Photo bleaching of C. incrustans’ autofluorescence was conducted by applying light from a mercury lamp onto tissue sections for 1 min prior to the CARD-FISH procedure [54]. CARD-FISH was carried out on sponge tissue sections as described by Croue et al. [20] with several modification (see Supplementary Information File 1 for details) and slides were visualised on a Nikon A1 Spectral confocal microscope. Conditions for probe specificity and stringency were assessed experimentally in an instance of two mismatches between a probe and non-target bacteria (Scopalina sp. probe on C. concentrica tissue). Each CARD-FISH experiment was conducted with a negative control slide that underwent the exact same procedure, but without adding the probe.

Results and discussion

For each sponge species analysed here, we obtained a single MAG each that was affiliated initially with the EC94 cluster (see below). The MAGs range in size from 1.1 to 2.2 Mbp, have completenesses of 66.66–95.69%, contamination rates of <1%, and an estimated genome size of 1.65–2.41 Mbp (based on the bacteria-specific markers of CheckM [34]). The observed GC contents of the genomes vary over a wide range from 50 to 68% (Table 2). Of the protein-coding genes, 22.2 ± 8.6% are hypothetical, i.e. without a predicted function (Table 2). With the exception of Cram cc1, the relative proportion of hypothetical proteins is in the lower range of reported values for MAGs from sponge symbionts [33, 55, 56], possibly because of the relatively small genome sizes and because the phylum Proteobacteria is comparatively well studied and annotated [57].

Table 2 Characteristics for MAGs affiliated with the EC94.

Phylogenetic analyses revealed a novel order within the Gammaproteobacteria

We recovered near full-length 16S rRNA genes from the metagenomes of C. concentrica (1494 bp), C. incrustans (1494 bp), C. crambe (1527 bp), Scopalina sp. (1499 bp), and a partial sequence from T. stolonifera (941 bp) (see Supplementary File 1, Table S4 for sequences). The AqS2 16S rRNA gene sequence (1519 bp) was already present in the SILVA 132 release (MEIE01000018). The 16S rRNA gene sequences from the symbionts of C. concentrica and C. crambe were entirely assembled in their respective MAGs, while a 240 bp region flanking the reconstructed 16S rRNA gene from Scopalina sp. overlapped with 97% identity with the end of a MAG contig and was thus retrospectively added to it. The 16S rRNA gene from the symbionts of C. incrustans and T. stolonifera did not overlap with any MAG contig, but we still manually included them based on comparative 16S rRNA and genome phylogenetic analysis (see below).

Insertion of these reconstructed 16S rRNA gene sequences by parsimony into the SILVA 132 SSU Ref Nr99 tree places them into the Betaproteobacteriales group EC94 [58], among sequences from other uncultured bacteria from marine samples. The 16S rRNA phylogeny (Fig. 1) provides a taxonomic assignment for previously described sponge-specific cluster SC112 to within the group EC94. Contrary to the recent SILVA 132 release, EC94 appears to form a sister clade sharing a common ancestor with the recently re-ranked Betaproteobacteriales (now Burkholderiales), rather than a subgroup [44]. The reconstructed 16S rRNA gene sequences appear to form two families within this order based on previously proposed percentage similarity thresholds [59]. Specifically, the 16S rRNA gene from AqS2 and the symbiont of C. concentrica share 90% similarity (i.e. same family) and sequences from the symbionts of C. incrustans, C. crambe, Scopalina sp., and T. stolonifera share between 89 and 93% similarity; however, the two families only share 83–85% sequence similarity. The order contains mostly 16S rRNA gene sequences derived from sponge samples, with the only exception being one sequence from a seawater sample and four sequences from two coral species (Fig. 1).

Fig. 1: 16S rRNA phylogenetic tree of gammaproteobacterial orders with Alphaproteobacteria as an outgroup.
figure 1

The displayed tree is a maximum-likelihood tree constructed with sequences longer than 850 bp. Circles represent bootstrap support of 70–100% (1000 bootstraps). Bar represents 10% sequence divergence. Sequences are named as: host species or source, length bp, accession number, order/sponge-specific cluster. Names in red are reconstructed rRNA gene sequences from this study, names in blue represent EC94 members of non-sponge origin.

The genome tree based on the GTDB reflects a similar phylogenetic arrangement (Fig. 2). The MAGs found here are clustered within the AqS2 order, separate from the Burkholderiales. AqS2 is the only genome representing EC94 in the GTDB. The accompanying classification places the genome Cym b2r in the same family as the AqS2 genome (RED value of 0.77), and the genomes Crel c29, Cram cc1, Sco b3, and Teth b1 (RED values of 0.58–0.59) in a separate family. The AAI analysis of the genomes supports this taxonomic discrimination. The genomes AqS2 and Cym b2r belong to the same family (AAI 53%), whilst the genomes Crel c29, Cram cc1, Sco b3, and Teth b1 belong to a second family (AAI 48.7–52.6%). Together both families belong to the same order (AAI 43–44.8%) using the thresholds proposed by Konstantinidis et al. [60].

Fig. 2: Phylogenomic tree of the Proteobacteria based on GTDB.
figure 2

Names in red are draft genomes of this study. As AqS2 was already in the database (GCA 001750625.1) it was used in tree construction. Circles represent bootstrap support of 70–100%. Bar represents 10% sequence divergence. Genomes are named as species name (genome bin ID).

Given these results and the refined placement of the order EC94, we propose taxonomic names following guidelines developed by the Genomics Standards Consortium [61], Konstantinidis et al. [60] and Chuvochina et al. [62]. As type material for this order we designate the MAG Cym b2r as it is highly complete (estimated to be 95%), contains a near full-length 16S rRNA gene (as well as 5S and 23S rRNA genes and 43 tRNA genes) and has the least number of contigs (Table 2) [62]. As all known members of EC94 are of marine origin, we chose to name them after the Oceanids (also known as sea nymphs), which are the 3000 daughters of the Titan Oceanus and Titaness Tethys in Greek mythology [63, 64]. In replacement of EC94 and the GTDB order-level placeholder AqS2, we propose the order name Ca. Tethybacterales, hailing from both Tethys and the sponge Tethya stolonifera studied here. The family Ca. Tethybacteraceae encompasses Ca. Tethybacter castellensis (represented by Cym b2r) and Ca. Amphirhobacter heronislandensis (represented by AqS2). The second family we propose is Ca. Persebacteraceae, encompassing Ca. Persebacter sydneyensis (represented by Crel c29), Ca. Beroebacter blanensis (represented by Cram cc1), Ca. Calypsobacter congwongensis (represented by Sco b3) and Ca. Telestobacter tawharanui (represented by Teth b1). A full etymological description can be found in the Supplementary Information (Table S5).

Ca. T. castellensis (Cym b2r), Ca. A. heronislandensis (AqS2), Ca. P. sydneyensis (Crel c29) and Ca. B. blanensis (Cram cc1) meet the minimum (>80% complete, <5% contamination, species genetic discreteness and ecological data; see below) and some additional (near full-length 16S rRNA gene and/or microscopy picture; see below) standards for descriptions of uncultured species as outlined by Konstantinidis et al. [60]. Ca. C. congwongensis (Sco b3) and Ca. T. tawharanui (Teth b1) only fall short on the genome completeness estimation; however, considering their symbiotic lifestyle (see below) and the fact that symbionts often have reduced genomes [62], it is possible that the draft genomes are more complete than CheckM predicts [34] Indeed, the Ca. Tethybacterales MAGs are all missing one CheckM bacterial single-copy marker gene (RecR). Konstantinidis et al. [60] noted that, for highly divergent species such as members of the proposed Ca. Tethybacterales, reliable matches may not be able to be retrieved against the universal protein-coding genes analysed by CheckM, which would result in an underestimated completeness. Given this, as well as the microscopic identification shown below, we feel confident in prescribing nomenclature.

Distribution analysis of Ca. Tethybacterales reveals global presence in both LMA and HMA sponges

In addition to world-wide reports of sponge-associated Ca. Tethybacterales, from Antarctic waters to tropical coral reefs [65, 66], our distribution analysis using the global data from the SEMP found Ca. Tethybacterales species belonging to the two described families to be present across a vast geographic range (Fig. 3). The BLASTn search found 157 OTUs in the SEMP belonging to the Ca. Tethybacterales, in 1077 samples of 159 identified sponge species (see Supplementary Information File 5).

Fig. 3: Global distribution of Ca. Tethybacterales-containing samples collected for the Sponge Earth Microbiome Project.
figure 3

Bubbles indicate collection sites for (a) marine sponges, (b) seawater, and (c) marine sediment samples. Bubble sizes are proportional to number of samples as indicated.

Of these 159 sponge species, 77 are classified or predicted to be LMA sponges and had a RA of 13.29 ± 23.23% (range: 0.05–93.40%), while 36 belonged to the HMA sponge group with RA of 1.66 ± 4.00% (range: 0.05–27.89%) [15], which were significantly different in RA (one-way ANOVA: p = 5.48797E−14). Ca. Tethybacterales dominance (e.g. >20% average relative sequence abundance) seems to be a feature predominately associated with some LMA sponges, as previously noted [20]. However, our analysis now clearly demonstrates that the Ca. Tethybacterales are also found in HMA sponges, contrasting with previous suggestions [20].

The global distribution patterns indicate a possible horizontal transmission of symbionts. Although there is also some strong evidence for vertical transmission of Ca. Tethybacterales symbionts in T. rubra [28], Tedania sp. [29], and of AqS2 in A. queenslandica (Ca. A. heronislandensis) [26], these two modes of transmission are not mutually exclusive. For example, Sipkema et al. [67] analysed the microbial communities of Corticium candelabrum, C. crambe and Petrosia ficiformis and identified vertical–horizontal transmission phylogenetic clusters (VHT clusters), where similar sponge-associated bacteria can be acquired either via vertical or horizontal transmission. It was proposed that successful horizontal transmission of bacteria could be facilitated through ELPs [67] (see further information below). Further, a recent study revealed that vertical transmission of symbiotic bacteria to sponge larva is fickle and the symbionts involved are often not as specific to a particular host species as previously assumed [68].

Members of the Ca. Tethybacterales are also found in marine invertebrates other than sponges (e.g. Fig. 1). Ca. Tethybacterales was found to make up 16–32% of the bacterial community of the coral Erythropodium caribaeorum [58, 69] and to be a dominant OTU in the coral Scleronephthya gracillimum [70]. Finding sponge-associated bacteria in corals is not uncommon and this has been defined by the concept of sponge/coral-specific clusters (SCCs), though SC112 was not originally considered one [18]. Shared and similar niches (e.g. biofilm-associated growth) in coral and sponge holobionts could facilitate the exchange of related bacteria (i.e. secondary symbiont transmission). Members of the Ca. Tethybacterales have also been recently reported in the healthy tissues of crown-of-thorns starfish Acanthaster cf. solaris at an abundance of 0.1–1% [71]. Since these starfish eat corals and sponges, it is possible that the presence of Ca. Tethybacterales is due to dietary uptake.

Our analyses also found Ca. Tethybacterales in 18% of seawater samples and 27% of marine sediment samples analysed in the SEMP (Fig. 3), although generally in very low abundances (RA%: 0.91 ± 2.36, range: 0.06–15.58 and RA%: 0.46 ± 0.71, range: 0.07–2.51, respectively). The discovery of sponge-associated bacteria in environmental samples is not surprising as many are thought to be widespread, but very rare. For example, Taylor et al. [19] investigated the prevalence of sponge-specific bacteria in the environmental sample analysed by the International Census of Marine Microbes (ICoMM). They found a total of 23 reads assigned to SC112 (Ca. Tethybacterales) in 412 million 16S rRNA gene pyrotags (V6 region) generated by ICoMM in four sample types (Amazon-Guianas Water (AGW) sediment, Baltic Sea Proper (BPS) seawater, Salt Marsh (LSM) and New Zealand (NZS) marine sediment), which is congruent with our findings. Ca. Tethybacterales has also been detected in seawater samples by Gantt et al. [23] (<0.01% of all reads) and Matcher et al. [30] (<0.05% of all reads). Given the predicted metabolic repertoire (below), it is possible that these bacteria could persist in the seawater and sediments and this would in turn form reservoirs for horizontal symbiont transmission.

Ca. Tethybacterales exhibit diverse localisation patterns within the sponge tissue

Visualisation by CARD-FISH confirmed the existence of each proposed species within the tissue of its respective host sponge and revealed remarkable diversity in term of cell morphology and localisation (Fig. 4).

Fig. 4: Detection of Ca. Tethybacterales species by CARD-FISH in sponge hosts.
figure 4

a Tethya stolonifera, b Crella incrustans, c Cymbastela concentrica and d Scopalina sp.. Each sponge’s target bacterium (Ca. Telestobacter tawharanui, Ca. Persebacter sydneyensis and Ca. Tethybacter castellensis Ca. Calypsobacter congwongensis, respectively) is labelled in green (Alexa Fluor™ 488) and counter-stained with blue (DAPI).

Ca. T. tawharanui cells appear large (2–2.5 µm) and clustered in structures that resemble bacteriocytes within the tissue of T. stolonifera (Fig. 4a). Individual spherical bacterial cells can be clearly seen surrounding an apparent sponge nuclei and probable vacuoles. T. stolonifera is considered a LMA sponge [72], though in this individual, bulging bacteriocytes are abundant throughout the mesohyl and in between megasclere tracts, forming a band from the inner pinacoderm to the outer cortex [73]. The intracellular lifestyle indicated here may also explain why the recovered Ca. T. tawharanui MAG had a relatively low read coverage compared to the other MAGs (Table 1), given that Ca. T. tawharanui comprises over 70% of T. stolonifera’s bacterial community based on 16S rRNA gene sequencing [22].

Ca. P. sydneyensis was difficult to localise due to the high autofluorescence of C. incrustans, even after photo bleaching the tissue. Figure 4b, however, shows a clear probe signal, where cells bear resemblance to Ca. T. tawharanui in shape, size and clustering, indicating a possible lifestyle associated with a bacteriocyte.

Ca. T. castellensis displays very large cell size (>3 µm) with a peculiar square shape (Fig. 4c). The cells are sparsely distributed on the pinacoderm side of the mesohyl in C. concentrica. The unusual shape and signal of these cells could be a result of active cell division and may be related to the minimum number of ribosomes required to emit a signal in both halves of a dividing cell [74].

Ca. C. congwongensis cells are comparatively small (1–1.75 µm) (Fig. 4d) and are more rod-like in shape. Ca. C. congwongensis is distributed abundantly and evenly throughout the mesohyl of Scopalina sp.. There is an apparent increase in density of Ca. C. congwongensis cells at the pinacoderm edge of the tissue section (see Supplementary Information, Fig. S1).

Localisation experiments are useful to support functional predictions, thus helping to provide deeper insights into sponge-microbe interactions. Two studies have also employed FISH to gather morphological and localisation data on Ca. Tethybacterales from C. crambe [20] (i.e. Ca. B. blanensis) and T. anhelans [27]. Both studies describe a homogenous distribution of Ca. Tethybacterales cells throughout the sponge mesohyl, and in the case of T. anhelans, particularly near choanocyte chambers, similar to what we observed in Scopalina sp. Bacterial cells in close proximity to choanocyte chambers will be exposed to oxygenated seawater pumped in by the sponge, which could support aerobic respiration and an opportunity to promptly sequester seawater-derived nutrients. In both cases the cells appear oval shaped and vary in size from 0.5 to 1 µm in C. crambe and 1.01 to 1.6 μm in T. anhelans.

Our analysis reveals considerably more diversity in Ca. Tethybacterales localisation patterns than previously observed, likely involving an intracellular lifestyle in bacteriocytes. Bacteriocytes are specialised host cells that encase bacteria, which likely evolved from a mutualistic infection [75]. Bacteriocytes have also been reported in a number of marine sponge species, mostly via transmission electron microscopy, although the taxonomic or phylogenetic nature of the intracellular bacteria is often undetermined [76,77,78]. The physiological reasons why sponges house bacteria within bacteriocytes is unclear, although various hypotheses have been proposed, such as that this type of cell structure can control the growth or metabolite production of bacteria that have been taken up from the environment [77]. Lastly, it could be the bacteria themselves require these specialised structures to survive and thus facilitate the specific metabolic or molecular interactions, such as those described in the following sections.

Metabolic features that unify the Ca. Tethybacterales

The five MAGs described in this study and Ca. A. heronislandensis share a number of predicted metabolic features, broadly summarised as being heterotrophic with the ability to use various carbon, nitrogen and sulfur sources, including taurine, spermidine and dimethylsulfoniopropionate (DMSP).

The genomes contain mostly complete and expressed tricarboxylic acid cycles (TCA) and glycolysis pathways (Supplementary Information File 3). Of the six genomes analysed, Ca. P. sydneyensis and Ca. C. congwongensis also encode both proteins required for the glyoxylate bypass, isocitrate lyase (Ca. P. sydneyensis 329 TPMav) and malate synthase (Ca. P. sydneyensis 788 TPMav, Ca. C. congwongensis 9069 TPMav), allowing acetate to be used for growth. The Ca. B. blanensis genome contains only the malate synthase (148 TPM). It important to note though that the absence of specific proteins should be interpreted with caution due to potential genome incompleteness.

The genomes also display a mostly complete and expressed electron transport chain for respiration. Two types of heme- and copper-containing terminal oxidase (HCO) superfamily cytochrome c oxidases are found in the genomes: a low-O2-affinity cytochrome aa3-type oxidase is encoded in Ca. A. heronislandensis, Ca. T. castellensis (682 TPMsum), Ca. B. blanensis (1360 TPMsum), Ca. C. congwongensis (CtaC: 3231 TPMav) and Ca. T. tawharanui; and a high-affinity cytochrome cbb3-type oxidase is found in Ca. B. blanensis (6210 TPMsum), Ca. C. congwongensis and Ca. T. tawharanui. The symbiont Ca. P. sydneyensis genome lacks one catalytic subunit (CtaD) but encodes and expresses subunit two (CtaC: 637 TPMav) of the aa3-type cytochrome complex IV, suggesting that the absence of a complete Complex IV is due to incompleteness of the genome. Ca. T. castellensis further encodes a non-HCO cytochrome bd oxidase (CydAB: 2188 TPMsum), which has a high affinity for O2 and is typically induced under O2-limited conditions [79]. The presence and expression of both high- and low-affinity cytochrome c oxidases is likely a reflection for the differing oxygen concentrations experienced in the sponge microenvironment [80] and has been observed in other sponge symbionts [56].

Ca. T. castellensis displays the potential for nitrate respiration via the NarGHJI system (1898 TPMsum) and also encodes a nitrate/nitrite transporter gene, NarK (Table 3); both were expressed. As with the presence of high-affinity cytochrome c oxidases, nitrate respiration would allow these bacteria to thrive under microaerobic conditions or when oxygen levels rapidly change, such as at the interface of choanocyte chambers and mesohyl tissues [81] where these bacteria have been localised (see Fig. 4). Nitrate respiration is also documented in an alphaproteobacterial symbiont of C. concentrica [56]. Ca. P. sydneyensis (2169 TPMsum), Ca. B. blanensis (NapC: 188 TPM), Ca. C. congwongensis (NapA: 1420 TPMav) and Ca. T. tawharanui encode a periplasmic NapABC nitrate reductase, which is possibly involved in dissimilatory nitrate reduction and redox balancing [82]. The Ca. B. blanensis (149 TPM) and Ca. T. tawharanui genomes also carry a copper-containing nitrite reductase gene (NirK) that can reduce nitrite to nitric oxide (NO). The enzymes required for subsequent reduction of NO are absent in all of the six genomes. In the absence of a mechanism for NO detoxification, NO accumulation would reach toxic levels [83], meaning that nitrite is unlikely to be used for respiration, and may solely be dissimilatory. The consistent lack of genes involved in NO reduction, despite the presence of nitrate and nitrite reductases, has been previously reported in the bacterial metagenomes from the marine sponges Rhopaloeides odorabile and C. concentrica [21].

Table 3 Transport systems (TS) present in Ca. Tethybacterales.

Transporter systems and metabolic pathways indicate differential utilisation of environmentally- or sponge-derived nutrients in the Ca. Tethybacterales

The majority of the transporters present in the Ca. Tethybacterales MAGs belong to the ATP-binding cassette (ABC) transporter superfamily, which can target compounds at low concentrations [84], of which several are significantly enriched when compared to free-living, marine Gammaproteobacteria (Supplementary File 6). Other transporter families present are the tripartite ATP-independent periplasmic transporters (TRAP) [85], the solute:sodium symporter family of the amino acid–polyamine–organo–cation superfamily [86] the major facilitator superfamily [87], the betaine–choline–carnitine transporter family (BCCT) [88], and the ammonium transport proteins family [89].

The MAGs display a diverse, but variable, range of transporters reflecting the heterotrophic metabolism mentioned above (Table 3). Sugar transporters are one of the most notable differences between families. MAGs of the family Ca. Tethybacteraceae (Ca. A. heronislandensis and Ca. T. castellensis) encode more sugar transporters, which are almost entirely absent in the genomes of the Ca. Persebacteraceae (Ca. P. sydneyensis, Ca. B. blanensis, Ca. C. congwongensis and Ca. T. tawharanui). Sugar utilisation has been previously noted as an adaptation by members of the phylum Poribacteria to the sponge environment, possibly by benefitting from the recycling or scavenging of the carbohydrate matrix of the host [90].

Taurine is often found in crustaceans and sponges [91] and cultured alphaproteobacterial symbionts from sponges have been recently shown to have the genomic potential to degrade this organosulfonate [92]. Ca. Tethybacterales show a differential ability to utilise taurine as a potential source of carbon, nitrogen, and sulfur. ABC transporters for taurine import are encoded in the four members of the Ca. Persebacteraceae (Table 3), which are rarely found in free-living, marine Gammaproteobacteria (Supplementary File 6). These four members also encode and express enzymes to degrade taurine to acetyl-CoA (via sulfoacetaldehyde and acetyl phosphate: Xsc: Ca. P. sydneyensis 1322 TPMav, Ca. B. blanensis 519 TPM, Ca. C. congwongensis 2048 TPMav, Pta: Ca. P. sydneyensis 109 TPMav and Ca. B. blanensis 352 TPM), while these genes and taurine transporters were not detected in the genomes of Ca. T. castellensis and Ca. A. heronislandensis. Since Ca. C. congwongensis, Ca. P. sydneyensis and potentially Ca. B. blanensis encode a glyoxylate bypass (see above), taurine could potentially be used as a source of carbon for growth. Taurine degradation also yields sulfite, which can be reduced to sulfide for biosynthesis using an assimilatory sulfite reductase (CysJI), which is present in Ca. B. blanensis (521 TPMsum), Ca. P. sydneyensis (1238 TPMsum), Ca. C. congwongensis and Ca. T. tawharanui. In contrast, genes are present in all six symbionts for the synthesis of cysteine from serine (CysE: Ca. T. castellensis 93 TPMav, Ca. P. sydneyensis 643 TPMav, Ca. B. blanensis 169 TPM, CysM: VT. castellensis 90 TPMav, Ca. P. sydneyensis 284 TPMav). However, none of the six MAGs encode any identifiable genes for the reductive assimilation of sulfate, suggesting that Ca. Persebacteraceae are dependent on sulfur derived from taurine degradation, while Ca. Tethybacteraceae may require other organic sulfur sources to meet their cellular sulfur demands. In Ca. B. blanensis, sulfite generated by taurine desulfonation may also be detoxified by a Soe-type sulfite dehydrogenase, and oxidised to sulfate (2569 TPM) [93].

The Ca. C. congwongensis genome uniquely encodes sulfur oxidation proteins, including those involved in the oxidation of sulfide (SoxF) and thiosulfate (SoxA, SoxB) [94, 95], similar to some sponge-associated Chromatiales [96, 97]. The function of these proteins in Ca. C. congwongensis may be to use sulfur compounds as energy sources, to supplement an organotrophic metabolism. Certain small organic substrates might also be oxidised as sources of reductant, including formate (FdhA: Ca. P. sydneyensis 463 TPMav, Ca. B. blanensis 613 TPM Ca. C. congwongensis, Ca. T. tawharanui), methanol (MoxF: Ca. T. tawharanui), and di/trimethylamine (Tmd: Ca. P. sydneyensis 86 TPMav). However, no symbiont exhibits an apparent genomic potential for methylotrophic growth. An aerobic carbon monoxide dehydrogenase (CoxLMS) is encoded in Ca. A. heronislandensis and Ca. T. tawharanui, suggesting CO oxidation can be used as an energy source. In these cases, it appears that oxidation is coupled to oxygen as the terminal electron acceptor (see above).

Marine phytoplankton often accumulate compatible solutes such as glycine betaine and DMSP in order to maintain favourable osmotic tensions and positive turgor [98]. These compounds are released into the marine environment upon cell lysis and are therefore also likely introduced into filter-feeding sponges. The symbionts Ca. B. blanensis, Ca. P. sydneyensis, Ca. C. congwongensis and Ca. T. tawharanui encode ABC transporters for the uptake of glycine betaine (Table 3), whereas Ca. A. heronislandensis and Ca. B. blanensis have BCCT family transporters for glycine betaine uptake [88]. In addition to being used as a potential osmolyte, glycine betaine may serve as a carbon and nitrogen source in some of these symbionts [99]. In support of this, Ca. P. sydneyensis, Ca. B. blanensis, Ca. C. congwongensis and Ca. T. tawharanui all encode enzymes for the degradation of glycine betaine to glycine through demethylation (Bhmt) via dimethylglycine (Dmgdh) and sarcosine (SoxBDAG) [100]. A glycine oxidase (ThiO) encoded in Ca. P. sydneyensis (240 TPMav), Ca. B. blanensis and Ca. C. congwongensis could then convert glycine to glyoxyl imine and H2O2. Glyoxyl imine can then be used for the biosynthesis of the thiazole ring of thiamine or spontaneously hydrolyses to glyoxylate and ammonia [101]. Ca. A. heronislandensis, Ca. T. castellensis, Ca. P. sydneyensis, Ca. B. blanensis, and Ca. T. tawharanui also have genes for the synthesis of glycine betaine from choline (BetA: Ca. T. castellensis 40 TPMav, Ca. P. sydneyensis 103–156 TPMav Ca. B. blanensis 308–465 TPM, BetB: Ca. P. sydneyensis 195 TPMav, Ca. B. blanensis 267 TPM), with Ca. B. blanensis additionally having choline sulfatase (310 TPM) for the conversion of choline sulfate to choline.

The sulfur analog of glycine betaine, DMSP [98], may also be imported by BCCT transporters [88]. DMSP can be demethylated by DMSP demethylase (DmdA: Ca. P. sydneyensis 192–17461 TPMav, Ca. B. blanensis 543–3244 TPM, Ca. C. congwongensis and Ca. T. tawharanui) or cleaved by DMSP lyase (DddP: Ca. P. sydneyensis 310 TPMav, Ca. B. blanensis 186 TPM and Ca. T. tawharanui). In the latter case, the gas dimethylsulfide is released as a by-product. DMSP degradation seems to be widespread in sponge symbionts given recent evidence for the presence of relevant degradation genes in genomes of some alphaproteobacterial sponge symbionts [102] and the metagenomes of two Antarctic sponges [103].

Sulfopropanediol (2,3-dihydroxypropane-1-sulfonate; DHPS) is another common metabolite in the marine environment, which is produced and released by diatoms [104]. The members of the family Ca. Persebacteraceae exhibit the genetic capacity to utilise DHPS as a carbon source, while we could find no support for this in the Ca. Tethybacteraceae. While the method of uptake in unclear, the four genomes of the Ca. Persebacteraceae encode sulfopropanediol 3-dehydrogenases (HspN: Ca. P. sydneyensis 430 TPMav, Ca. B. blanensis 342 TPM) and (2R)-sulfolactate sulfo-lyase (SuyAB: Ca. P. sydneyensis 997 TPMsum, SuyB: Ca. B. blanensis 345 TPM).

The Ca. Tethybacterales MAGs also contain genes for the import and catabolism of spermidine and putrescine, which are polyamines commonly found in animals, including sponges [105, 106]. The primary transport system PotABCD (Table 3) can import both spermidine and putrescine, but has preference for the former [107, 108]. Again, the genomes of the Ca. Persebacteraceae apparently distinguish themselves from those of the Ca. Tethybacteraceae by encoding enzymes to degrade putrescine to succinate via the transamine pathway (SpuC: Ca. P. sydneyensis 450 TPMav, Ca. B. blanensis 586 TPMGabD: Ca. P. sydneyensis 411 TPMav, Ca. B. blanensis 215 TPM) and the glutamated putrescine pathway (PuuA/PuuB/PuuC: Ca. P. sydneyensis 594 TPMsum, Ca. B. blanensis 683 TPMsum, PuuC: Sco 4782 TPMav) [109] yielding ammonia, alanine and glutamate. Succinate would then be available to enter the TCA cycle for energy conservation and biosynthesis [110]. Even though PotABCD prefers spermidine import over putrescine, genes for spermidine-specific degradation (e.g. SpdH) are absent [111]. This somewhat contradictory instance has also been observed in the marine alphaproteobacterium Ruegeria pomeroyi [110], where it was concluded that given the expression of homologues for putrescine degradation, exogenous spermidine is imported and converted to putrescine (then ultimately succinate) through a mechanism yet to be determined. Indeed, spermidine to putrescine conversion has been documented in other systems [112] and therefore could also be occurring here.

All six Ca. Tethybacterales also encode proteins for the synthesis of polyhydroxyalkanoate, a carbon and energy storage polymer produced from acetyl-CoA (PhaA/PhaB/PhaC: Ca. T. castellensis 3650 TPMsum, Ca. P. sydneyensis 2870 TPMsum, Ca. B. blanensis 2744 TPMsum and Ca. C. congwongensis 13418 TPMsum). Symbionts Ca. C. congwongensis and Ca. T. tawharanui also encode polyphosphate kinase, allowing for the accumulation of inorganic polyphosphate, for energy and/or phosphate storage [113]. Cyanobacterial symbionts of sponges have been previously found to sequester phosphate in this fashion, thought to be driven, in part, by oxic–anoxic perturbations within the sponge tissue, which ultimately impacts phosphate cycling within the benthos [114]. Polyphosphate accumulation has also been associated with bacterial environmental durability and virulence [115, 116] and while many species of endosymbiotic bacteria appear to have lost this ability [115], retaining this function in sponge symbionts may be related to the horizontal transmission (see above) and the requirement to persist outside of their host. Overall, these data indicate that members of the Ca. Tethybacterales have strategies to cope with carbon, nitrogen (see Supplementary Information File 1) and, for some species, phosphate imbalances.

Symbiosis effectors in the Ca. Tethybacterales

Given their intimate interaction with and within sponge cells (see Fig. 4), the Ca. Tethybacterales must contain and express effector molecules to mediate and control their symbiosis, especially as sponges feed on bacteria and nutrient particles through phagocytosis [8]. ELPs are predicted to be involved in mediating sponge-bacteria interactions and have been shown to inhibit eukaryotic phagocytosis, which would be an important feature for a symbiont to persist inside the sponge [10,11,12]. The Ca. Tethybacterales genomes contain ELPs belonging to seven classes: ankyrin repeats (ANK), tetratricopeptide repeats (TPR), pyrrolo-quinoline quinone repeats (PQQ), Sel1-like repeats (SEL1) fibronectin domain III repeats (FN3), and cadherin (a calcium-dependent cell adhesion molecule) (Table 4). Of these, ANK-type ELPs are more abundant in the Ca. Tethybacterales than in free-living, marine Gammaproteobacteria, while the other ELP types do not appear to be enriched in the Ca. Tethybacterales. The distribution and abundance of these proteins is generally consistent across all genomes (Table 4), with the exception of FN3 only being present in Ca. T. castellensis and cadherin only being found in Ca. B. blanensis. Expression (>40 TPM) of all ELPs is observed in Ca. T. castellensis and Ca. B. blanensis, whereas expression of only TPR and PQQ is observed in Ca. P. sydneyensis and only TPR has detectable expression in Ca. C. congwongensis (Supplementary Information File 3). Such subtle differences in the presence and/or expression of certain ELPs might determine whether a bacterium can reside extracellularly (e.g. by avoiding an initial uptake through a phagocytic cells) or end up in intracellular localisation (e.g. vacuole) by blocking the maturation of a phagosome [10, 11]. Totipotent archaeocytes are the main phagocytes in the sponge mesohyl [117] and an ELP-mediated accumulation of intracellular bacteria could contribute to the formation of bacteriocytes (e.g. structure in Fig. 4).

Table 4 Number of genes encoding ELP classes in MAGs.

Often considered a feature of pathogenicity, by definition another form of symbiosis, all Ca. Tethybacterales genomes also contain genes for the release of lipoproteins via the ABC-2 transporter LolABCDE complex (Table 3). The release of lipoproteins is reported to play a direct role in bacterial colonisation/invasion, evasion of host defences and immunomodulation [118] and strategies advantageous to a symbiont in the host environment. The genomes Ca. P. sydneyensis, Ca. B. blanensis and Ca. C. congwongensis also contain lipo-oligosaccharide exporting nodulation genes, NodIJ (Table 3). Better known from legume-rhizobia symbiosis, the NodIJ transport system secrets host-specific Nod factors important for the recognition of rhizobia during legume nodulation [119]. The possible role of nodulation genes in bacterial symbionts of sponges is completely unexplored but could represent another mechanism for mediating symbiosis.

Conclusion

Our analysis reveals that Ca. Tethybacterales are a large, widespread, and divergent clade of (almost exclusively) sponge-associated symbionts, which suggests that they have adapted to the sponge environment reflecting a degree of co-evolution, as also proposed by Matcher et al. [30]. We propose that the Ca. Tethybacterales entered a symbiosis with sponges relatively early in evolutionary time, as indicated by their deep branching within the Gammaproteobacteria. This early symbiosis is reflected in the common features that the Ca. Tethybacterales share, such as their heterotrophic metabolism and respiratory pathways that have allowed them to occupy specific niches in the sponge and derive energy and carbon either from their host or compounds in seawater. However, we also find that members of the Ca. Tethybacterales have subsequently functionally radiated, reflected in their quite distinct lifestyles, most notably their diverse morphologies, predicted substrate preferences and localisation within the sponge tissue. These differences may be the result of a more specialised co-evolution with their present-day sponge hosts, similar to the diverse characteristics observed in the Rickettsiales, which have evolved to have mutualistic, commensal, and parasitic relationships with an extensive array of hosts [120]. To expand on these hypotheses, further acquisition and analysis of additional Ca. Tethybacterales genomes are required, complemented with the inference of a molecular clock analysis [121].