Six adult sponges (Amphimedon queenslandica) were collected from Shark Bay, Heron Island, Great Barrier Reef, QLD (23° 27′S, 151° 55′ E), three in September 2012 (HSep1-3), and three in December 2012 (HDec1-3). Each sponge was removed from the substrate with a sterile scalpel and pieces were stored in tubes of 0.22 μm filter sterilized seawater (FSW) half submerged in iced water and quickly transported back to the Heron Island Research Station for immediate processing. Enrichment of sponge-associated bacteria from sponge cells followed the protocol of Thomas et al. (2010) with slight modifications, and total genomic DNA was extracted from these cell preps using a standard phenol-chloroform method. Details are provided in text of the Supplementary Material.

Each of the six samples was multiplexed and sequenced by Macrogen (http://www.macrogen.com) in one lane using the HiSeq Illumina platform. Pair-end reads were 100 bp long and had an insert size ranging between 150 and 250 bp. Reads were filtered with Trimmomatic (Bolger et al., 2014), removing the first five bases from the start of the read, performing a sliding window trimming with a window of four bases and an average quality threshold of 15, and eliminating resulting reads that were less than 50 bp long. We also did not retain ambiguous reads with Ns, nor reads that aligned to contigs that unambiguously belonged to the genome of the host sponge, Amphimedon queenslandica, using Bowtie2 (Langmead and Salzberg, 2012). Metabinning and genome construction were performed following the protocol of Albertsen et al. (2013) with slight modifications (Supplementary Material text; Figure S1).

Small subunit (i.e., 16S) rRNA gene-containing reads were identified from the filtered FASTQ files using Metaxa v. 2.0.1 (Bengtsson et al., 2011). Extracted reads were assembled with IDBA-UD (Peng et al., 2012) using a 99% overlap identity. We first performed assemblies on these reads keeping samples separate, and then another assembly pooling reads from all samples. We used CAP3 (Huang and Madan, 1999) on all assemblies to merge overlapping contigs, and checked manually the resulting singletons and collapsed contigs. These derived 16S rRNA sequences were broadly classified using BLASTN and the NR database, excluding uncultured and environmental samples. After confirming that we recovered most reads that had been initially extracted by Metaxa within our final list of assembled contigs, we then used the original FASTQ files and our derived 16S rRNA sequences as reference, and applied the closed-reference OTU picking strategy of the QIIME program (Caporaso et al., 2010b) to derive 16S rRNA gene abundance estimates for each sample.

We obtained gene annotation for the two bacterial genomes that we deemed had adequate coverage (AqS1 and AqS2). Contigs for each genome were uploaded on the RAST server (Aziz et al., 2008; Overbeek et al., 2014) and annotated using the Classic annotation scheme, the RAST gene caller, and FIGfams v. 70. We opted for errors and frame shifts to be fixed automatically, and for the pipeline to BLAST large gaps for missing genes. We compared the gene content of AqS1, AqS2, Thioalkalivibrio sp. HK1, three of AqS1 closest relatives and other known sulfur-oxidizing symbionts using the function based comparative tool available on the SEED viewer within the RAST platform. For each genome, we either used predicted protein-coding sequences that were publically available or obtained these ourselves using the RAST server. Although a gene annotation for T. sp. HK1 is available on NCBI, we derived a new one via the RAST server for ease of comparison with the AqS1 genome. This annotation can be made available upon request.

We searched for the presence of domains associated with proteins postulated to facilitate host-symbiont interactions in the predicted protein-coding sequences of AqS1, AqS2, other known sulfur-oxidizing symbionts and the three closest free-living relatives to AqS1. To do so, we scanned for specific PFAM domains using hmmsearch from the program hmmer 3.1 (hmmer.org), using the trusted cutoff bit score thresholding option (–cut_tc). We used Interpro (http://www.ebi.ac.uk/interpro/), NCBI conserved domains (http://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi) and TMHMM Server v. 2.0 (http://www.cbs.dtu.dk/services/TMHMM/) to predict domains and transmembrane regions of the high-affinity branched-chain amino acid ABC transporter proteins. An example of putative 2D topology for one of the transmembrane protein complex was derived using TOPO2 (http://www.sacs.ucsf.edu/TOPO2/topo2-help.html). Finally, we used the CRISPRfinder program (http://crispr.u-psud.fr/Server/) to deduce the organization of the CRISPR in AqS1. The 16S rRNA gene sequences of AqS1, AqS2, and AqS3 are available from the NCBI GenBank database under accession numbers KX772780-2. The draft genomes of AqS1 and AqS2 have been deposited as Whole Genome Shotgun projects at DDBJ/ENA/GenBank under the accession numbers MEID00000000 and MEIE00000000, respectively. The versions described in this paper are MEID01000000 and MEIE01000000, respectively. The genome browsers of AqS1 and AqS2 are hosted on the RAST server (http://rast.nmpdr.org/seedviewer.cgi) under genome IDs 6666666.207661 and 6666666.208689, respectively, and are accessible via username 'guest' with password 'guest').

Phylogenetic analyses were performed for the three most abundant 16S rRNA sequences in our metagenome data. Nearest related and described sequences were retrieved from the 2013 Greengene and NCBI NR databases. The partial sequence from T. sp. HK1 reported in Tian et al. (2014) was kindly provided by Dr Ren-Mao Tian. A sequence alignment was derived with PyNAST (Caporaso et al., 2010a) using the template file “core_set_aligned.fasta.imputed” from QIIME v. 1.8.0 (Caporaso et al., 2010b) and manually checked using AliView (Larsson, 2014). A maximum-likelihood tree was constructed with RAxML v. 8 (Stamatakis, 2014) with 1000 bootstrap resamples, and visualized using FigTree v. 1.4.0 (http://tree.bio.ed.ac.uk/software/figtree/). For trees with low bootstrap support at given nodes, we checked whether some taxa had ambiguous or insufficient phylogenetic signal by deriving their leaf stability indices using the program RogueNaRok (Aberer et al., 2013).

Since a large proportion of genes that endow the AqS1 symbiont with heterotrophic capabilities are most closely related to genes of several alphaproteobacteria, we wanted to test the hypothesis that some of these have been acquired horizontal gene transfers (HGT). To do so, we performed a phylogenetic analysis on the concatenated alignment of high-affinity branched-chain amino acid ABC transporter proteins. Nearest related and described sequences were retrieved from the NCBI NR database. We also included two species from the order Chromatiales as an outgroup. We derived alignments using AliView and checked for the presence of ambiguously aligned positions using Aliscore and Alicut (Kück et al., 2010). We used prottest3 (Darriba et al., 2011) to choose a model of evolution for each gene. We concatenated our sequence alignment into one supermatrix using FASconCAT (Kuck and Meusemann, 2010).

We sequenced six bacterial-enriched samples from the sponge Amphimedon queenslandica and recovered between 23,906,452 and 27,938,004 pairs of Illumina reads from each library. Filtering out low quality reads and contaminating sponge host DNA resulted in 14,563,169 to 16,046,247 read pairs contributing to bacterial genome assemblies (Table S1). To isolate the symbiont genomes of interest—AqS1, AqS2, and AqS3—we derived scaffold coverage plots with GC content and taxonomic information overlayed, using a slight modification of the protocol developed by Albertsen et al. (2013) (Supplementary Material; Figures S1, S2). The AqS1 genome was assembled into 127 contigs (median length 13 Kbp, total length 4, 2 Mbp, with 59% average GC content) and the AqS2 genome into 239 contigs (median length 3 Kbp, total length 1,6 Mbp, with 69% average GC content). We also recovered a third genome (AqS3) that we were able to assemble into 233 contigs (median length 5 Kbp, total length 3,17 Mbp, with 62% average GC content) (Table 1). To evaluate the quality of the assemblies, we checked for the presence of 105 conserved single-copy genes in each genome (Supplementary Material; Table S2). In general, we found that the total number of genes recovered in each of the three symbionts was on par with the unique essential gene list. Specifically, we estimate that the reconstructed AqS1 genome is 99% complete (we recovered 104 out of 105 genes) and that of AqS2 and AqS3 are 85% and 94% complete, respectively. However, since the coverage of AqS3 was low (<10X read coverage) and a principal component analysis on the tetranucleotide frequencies of its genome scaffolds indicated the presence of possible contamination that could not be removed, we elected to omit this genome from our downstream gene annotation, which is thus restricted to the two most abundant symbionts, AqS1 and AqS2.

Based on maximum-likelihood phylogenetic analysis of the 16S rRNA gene, the AqS1 symbiont clusters admist symbionts of various other sponge species (identity range 95–99%); sister to these is a small clade of marine environmental samples (Figure 1; Figure S3). These clades fall at the base of the family Ectothiorhodospiraceae, order Chromatiales (commonly known as sulfur-oxidizing bacteria). The closest cultured bacterial relatives (identity range 90–92%) are from Ectothiorhodospira and Thioalkalivibrio, both in the family Ectothiorhodospiraceae; the top BLAST NR hit was to the free-living Thioalkalivibrio paradoxus ARh 1. The AqS1 symbiont falls into the same sponge-specific clade as bacteria derived from the demosponges Callyspongia vaginalis and Gelliodes carnosa; also in this clade is a symbiont recently—but, we suggest, perhaps prematurely—named Thioalkalivibrio sp. HK1, whose genome was assembled from a metagenome, obtained from the sponge Haliclona cymaeformis (Tian et al., 2014). The other sponge-specific clade includes symbionts from a homoscleromorph (Plakortis halichondrioides) and several Demospongiae. For AqS1, we also performed phylogenetic analysis based on 30 concatenated essential marker genes, and including the genome of the close relative T. sp. HK1 plus a selection of 11 other publically-available genomes from free-living Chromatiales bacteria. In agreement with the 16S rRNA phylogenetic analysis, both the AqS1 and the HK1 sponge symbionts were recovered in a cluster of bacterial species belonging mostly to the family Ectothiorhodospiraceae (data not shown).

The 16S rRNA sequence of the AqS2 symbiont forms part of a large, well-supported cluster of demosponge and one coral-derived bacteria (identity range 84–99%), that is quite divergent from recognized families within the class Betaproteobacteria. Similar to AqS1, the AqS2 sequence clades with bacteria whose genome was assembled from metagenomic data derived from C. vaginalis and G. carnosa, as well as another Amphimedon sp. from the Red Sea (Figure 1; Figure S4). The closest cultivated relatives of the 16S rRNA of AqS2 are from the genera Azoarcus and Thauera (identity range 82–84%), both in the family Rhodocyclaceae; the top BLAST NR free-living hit was to the Azoarcus sp. CIB, but only with 84% identity.

The 16S rRNA gene sequence of symbiont AqS3 shows highest percent sequence similarity to bacterial clones originating from several demosponges (identity range 84–99%). The closest described relative is the Chromatiales Thioprofundum hispidum strain gps61, but the identity score is low (86%). In our phylogenetic analysis, AqS3 clustered separately with gammaproteobacteria derived from other demosponges, but this apparently novel clade has low bootstrap support, and thus its placement within the gammaproteobacteria is not yet resolved (Figure 1, Figure S3).

Since bacterial genomes can contain more than one copy of the 16S rRNA gene, the relative abundance of detected phylotypes from metagenomic DNA is potentially biased. We were interested in validating the binning strategy of our assembly by establishing whether the three assembled genomes would be also recovered as the most prevalent in the 16S rRNA pool, and whether any other prominent genome had not been accounted for. We also wanted to assess whether the overall composition of the microbial population found in the sponge A. queenslandica was broadly consistent across the six individual sponges that we sampled, and between the two different sampling seasons. It is also worth noting that we recovered the 16S rRNA gene in our single-copy gene list derived from our binned contigs for each of the three AqS genomes (Table S2).

Our 16S rRNA sequence reconstruction from the metagenomic Illumina data confirmed that the three symbionts for which we retrieved a genome assembly were also the most highly represented in the small subunit rRNA (ssrRNA) sequence pool derived from both separate and combined sponge metagenomic DNA samples (Table S3). In all of the individuals tested, and regardless of season, the microbiome was dominated either by the AqS1 or the AqS2 phylotype, which together represented as much as 90% of the overall microbial 16S rRNA pool. Reads matching to the ssrRNA sequence of the AqS3 symbiont were detected across all samples but recovered on average 20-fold less than those matching to the AqS1 phylotype. The remaining phylotypes that we detected all occurred at much lower abundance, with most displaying a relative abundance of <1% (Figure S5; Table S3), and could represent either rarer members of the A. queenslandica symbiont community or of the seawater, or contaminants introduced during sampling or processing. Because we could not discriminate between these scenarios, we do not discuss these further. Although we also recovered 16S rRNA sequences matching to members of the domain Archaea in three sponge individuals, these were from single pair-ended reads that could not be assembled into contigs.

We observed very little variation between the microbiome composition of individual sponges within a season; the two dominant phylotypes represented a similar proportion of the microbial community across sponge individuals. The AqS1 phylotype was more prevalent than the AqS2 phylotype in the samples collected in September, while the relative abundance of the two symbionts was more equivalent in those sponges collected in December, suggesting a slight seasonal variation (Figure S5; Table S3). Altogether, these data indicate that A. queenslandica individuals support a stable microbiome, the bulk of which comprises a coexisting gammaproteobacterium and betaproteobacterium. These observations are consistent with that reported in a companion paper, Fieth et al. (in review) in which OTU1 and OTU4 correspond to AqS1 and AqS2, respectively. The genome contents of both of these symbionts are described below.

A total of 4778 protein-encoding genes are predicted in the draft genome of AqS1, 54% of which are hypothetical at time of writing, while 29% are incorporated into 347 subsystems derived from the RAST server. The genome contains one copy of the 5S rRNA, the 16S rRNA, and the 23S rRNA genes – the latter bridging two scaffolds (scaffold 49 and 73), and 47 predicted transfer RNAs (tRNA). The gene repertoire of AqS1 is very similar to that of T. sp. HK1 (Figure 2; Table S4), and their contigs show high synteny (data not shown). Compared to AqS1, AqS2 displays a much smaller gene repertoire that reflects its reduced genome assembly size. A total of 1349 protein-encoding genes are predicted in the genome of AqS2, but only 22% of these are hypothetical, while 50% are incorporated into 244 subsystems. For AqS2, we recovered a single copy each of the 5S rRNA, the 16S rRNA, and the 23S rRNA genes, plus 42 putative tRNAs.

AqS1, along with its close relative T. sp. HK1 (Tian et al., 2014), show the highest phylogenetic similarity to Ectothiorhodospiraceae and like most bacteria belonging to this order, the two bacteria appear capable of sulfur oxidation. The genome of AqS1 encodes for a Sox enzyme complex (for periplasmic oxidation of various reduced sulfur compounds) coupled with a reverse dissimilatory sulfate reduction pathway (for sulfide oxidation in the cytoplasm) (Figure 3). The Sox pathway encoded by AqS1 is truncated (that is, it includes a soxBHXAYZ but no soxCD), as is the case for T. sp. HK1, although a point of difference is that no soxH gene was recovered in the latter. The Sox genes are spatially separated in the genome of AqS1, with soxB and soxH not occurring contiguously with soxAXYZ (Table S5). In the absence of the soxCD genes, obligatory intermediate sulfur globules are predicted to form as pathway byproducts, and these can be subsequently oxidized to sulfate by the dissimilatory sulfite reductase (dsr) pathway (Pott and Dahl, 1998). Homologues to the DsrABC proteins, as well as associated proteins needed to build the DsrMKOJP complex, are present in both the AqS1 and T. sp. HK1 genomes. In the current AqS1 assembly, the dsr genes are organized in a single polycistronic gene cassette (dsrABEFHCMKLJOP) that spans two contigs.

AqS1 and T. sp. HK1 also encode a gene encoding for ATP sulphurylase (sat) and for adenosine-5′-phosphosulphate (APS) reductase (aprAB), both of which enable the conversion of sulfide to sulfate via APS. We did not detect the membrane-binding protein or the QmoABC complex that participate in the electron transfer during this conversion step in other bacteria (Frigaard and Bryant, 2008). However, the genomes of both these symbionts possess an operon-like gene cluster, adjacent to the sat gene, that comprises a qmoA—and a qmoB-like gene followed by a heterodisulphide reductase (hdr) C and hdrB gene (Figure 3; Table S5). This putative QmoAB-HdrCB protein complex could provide an alternative means of cytoplasmic electron transfer during sulfide oxidation, as described in another sulfur-oxidizing bacterium Thiolapillus brandeum (Nunoura et al., 2014).

Another enzyme involved in sulfur metabolism and that catalyzes the formation of sulfur globules is sulfide:quinone oxidoreductase, but we could not identify this gene in either of the AqS1 or the T. sp. HK1 genomes. However, we did detect the presence of a thiosulphate sulphurtransferase (rhodanese) gene and two rhodanese-related sulphurtransferases in both symbionts, with one copy located near the soxH gene in AqS1. The proteins encoded by these genes could participate in thiosulphate import to the cytoplasm (Aussignargues et al., 2012). Finally, both genomes include genes whose products can process extracellular taurine (i.e. tauA, B and C), alkanesulfonate (ssuB) and dimethylsulfoniopropionates (dmdA and dmdB), all of which can contribute to sulfide metabolism (Figure 3; Table S5).

Our RAST-based gene annotation analyses indicate that components of the Calvin-Benson-Bassham (CBB) pathway are present (Table S5), but both the AqS1 and the T. sp. HK1 genomes are missing homologs to key enzymes: (a) ribulose-1,5-bisphosphate carboxylase/oxygenase (RubisCO) form I and II, and (b) others known to be essential in alternative carbonfixation routes (Hügler and Sievert, 2011).

To establish the heterotrophic capabilities of AqS1, we focussed on annotated genes found in the AqS1 symbiont but absent in its three closest free-living relatives, all of which are known to be obligate autotrophs (Figure 2B). From this list of 285 genes (86% of which are also present in the genome of T. sp. HK1), we recovered several genes capable of utilizing a wide range of sugars (monosaccharides, polysaccharides, organic acids and sugar alcohol subcategories) as well as amino acids and derivatives (e.g., polyamine, aromatic compound and creatine metabolism) as energy sources. We also found multiple genes encoding membrane transporters. These include four tripartite adenosine triphosphate (ATP)-independent periplasmic (TRAP) type transporters that are required for dicarboxylic acid import, as well as 29 tripartite adenosine triphosphate (ATP) binding cassette (ABC) transporter genes responsible for carbohydrate, amino acid and peptide uptake (Table S6A). Focusing on the high-affinity branched-chain amino acid ABC transporters, the AqS1 genome encodes multiple copies of the polycistronic LivFGHMK gene cassette, five of which are complete, and one of which is missing a livM homolog. Each cassette displays a different gene organization. Nonetheless, based on the transmembrane domain and protein domain predictions for each of the five complete cassettes, they all have the capacity to form functional transmembrane complexes (Figures 4A,B; Table S5). Finally, both D-octopine dehydrogenase and opine oxidases A, B and C are present in the genome of AqS1 (Table S5; Table S6A), which could enable AqS1 to use methyl-group donors like opines for carbon and nitrogen source. Therefore, as was observed for the symbiont T. sp. HK1 (Tian et al., 2014), the symbiont bacterium AqS1 harbors a large repertoire of genes that play a role in heterotrophic metabolism. This suggests that it might not need to rely on carbon fixation, as it is able to acquire and oxidize organic compounds.

As further evidence of the versatile metabolism of AqS1, we also found genes encoding small, medium and large subunits of the aerobic-type carbon monoxide (CO) dehydrogenase, a CO dehydrogenase F and a CO oxidation accessory protein (Table S5). These imply that the symbiont can generate reductive energy from the oxidation not only of sulfide, as described above, but also of CO. The RAST-based gene model predictions suggest that for catabolism, AqS1 is equipped with a complete Embden-Meyerhof-Parnas (glycolysis) pathway for the metabolism of hexose sugar, and a near-complete TCA cycle for the production of intermediates in cofactors and amino acid biosynthesis. A near complete pentose phosphate pathway is also present, providing a mechanism for production of NADPH and ribose sugars for nucleotide biosynthesis.

AqS1 is well equipped for phosphorus metabolism, including mechanisms to regulate the assimilation of inorganic phosphate (Pi). Specifically, the AqS1 genome encodes for components of the Phosphate (Pho) regulon, namely a Phosphate regulon sensor protein (PhoR) that detects levels of Pi in the environment, and downstream response genes (i.e., Phosphate regulon transcriptional regulatory protein [PhoB], Phosphate-starvation-inducible protein, predicted ATPase [PhoH], Phosphate transport system regulatory protein [PhoU], and two genes encoding response regulators in a two-component regulatory system) (Table S5). Also present is a gene encoding a polyphosphate kinase (EC 2.7.4.1) that catalyzes the formation of long-chain polyphosphate (polyP) (Table S5).

Finally, this symbiont has retained either partial or complete biosynthesis pathways for most essential amino acids (Figure 5A) and, like most bacteria, is capable of assimilating ammonia (genes encoding EC 1.7.1.4, EC 1.4.1.13, and EC 1.4.1.2 all are present) (Table S5). On the other hand, the absence of ureases and their associated proteins suggests that AqS1 is unable to use urea to cover its nitrogen needs. AqS1 also appears to lack a denitrification pathway; we could not detect genes encoding the intermediate enzymes of the anaerobic respiratory chain, nitrite reductase (NO forming) (EC 1.7.2.1) or nitric oxide reductase (EC 1.7.99.7).

The identification of 285 genes found in the AqS1 symbiont, but not in its closest free-living relatives, led us to examine the origin of these genes. Interestingly, 136 of the 285 genes showed highest BLASTP similarity to various species in the class Alphaproteobacteria, with the order Rhodobacterales being the most represented (please refer to numbers that are reported in brackets in Figure 2B; Tables S7A,B). This suggests that the versatile metabolic capacity of the AqS1 symbiont could be, at least in part, a result of horizontal gene transfer (HGT) from one or more alphaproteobacteria. To further explore this possibility, we performed a phylogenetic analysis focusing on the multiple high-affinity branched amino acid ABC transporter protein cassettes found in both AqS1 and T. sp. HK1 (see above). These all clade with different representatives of Alphaproteobacteria, and specifically within the orders Rhodobacterales and the Rhodospirilales. This same association was also discussed by Tian et al. (2014) for T. sp. HK1 (Figure 4C).

Further support for an HGT origin, albeit not alone compelling, is that the AqS1 genome encodes genes implicated in HGT mechanisms, including four copies of the ISSpo3 transposase, mobile element proteins and phage-associated proteins (i.e., adenine DNA methyltransferase, phage-associated; putative prophage protein ps3). Genetic systems for conjugation and transformation were not detected, but DNA recombination and repair enzymes are present (i.e., Rec F, N, O and R), as well as protein families involved in general DNA-modifying activities (i.e., HNH endonuclease), and restriction-modification (R-M) and toxin–antitoxin systems (Type I, Type II, and Type III R-M systems, proteins from the Doc/Phd family, and VapI of the HigAB system) (Table S5).

Interestingly, mechanisms to control excessive genetic exchange and to minimize the introduction of foreign DNA are also encoded by the AqS1 genome. Single clustered, regularly interspaced, short, palindromic repeats (CRISPRs) and their associated proteins form adaptive immunity systems that occur in most archaea and many bacteria and act against invading genetic elements, such as viruses and plasmids (Makarova et al., 2011). A single CRISPR is found in the genome of AqS1. It is composed of at least fifteen direct repeat regions that are each 36 nucleotides in length, separated by at least thirteen spacers that are each 30 nucleotides in length. A gap present in our assembly, which bridges the first recovered spacer at the 5′ end of the predicted CRISPR, precludes the recovery of its leader sequence. The AqS1 CRISPR system has a canonical genomic arrangement (Makarova et al., 2011), with CRISPR-associated (CAS) genes (Csn1 and Cas1) located upstream of and directly adjacent to the CRISPR array (Figure S6; Table S5).

Our RAST-based gene annotation analyses recovered complete or near-complete key functional pathways in AqS2, despite its smaller genome assembly size. These included glycolysis, the TCA cycle and the pentose phosphate pathway. This symbiont does not seem to have nitrogen reduction, carbon fixation or sulfur oxidation capacity and, in stark contrast to the AqS1 genome, there are no genes associated with phosphorus metabolism. However, AqS2 does possess a gene operon that includes the three subunits of carbon monoxide (CO) dehydrogenase form I, D and F, and CO oxidation accessory protein (CoxG), suggesting that, like AqS1, it can generate reductive energy from the oxidation of CO. AqS2 might also be capable of assimilating ammonia, as it encodes NADP-specific glutamate dehydrogenase (EC 1.4.1.4) (Table S5). It also (like AqS1) has retained partial or complete biosynthesis pathways for most essential amino acids (Figure 5A), and genes involved in carbohydrate metabolism and membrane transport are also present, indicating AqS2 has heterotrophic abilities (Table S7). On the other hand, the genome of this symbiont does not show any evidence of past horizontal gene transfer (including no transposable insertion elements, or CRISPR-Cas resistance systems). We did recover genes involved in protection against foreign DNA (e.g., R-M systems) and cellular stress response (e.g., chaperones, membrane proteases) (Table S5).

It has been postulated that co-occurring symbionts avoid competition by utilizing different carbon and energy sources (Duperron et al., 2006; Kleiner et al., 2012a). To establish preliminary evidence of resource partitioning between AqS1 and AqS2, we scanned for differences in their overall repertoire of predicted proteins involved in carbohydrate metabolism and transport. AqS1 encodes a wider range of membrane transporters than AqS2 (e.g., ABC transporter branched-chain amino acid, ABC transporter dipeptide, ABC transporter tungstate, copper transport system) (Table S7), and a genome-wide comparison clearly shows that each symbiont is equipped with a unique suite of proteins that would permit them to import or export different substrates (Figure 6). These include sulfur-related compounds, ammonium, vitamin B12, putrescine, spermidine, glutamate, aspartate, L-proline, glycine betaine, benzoate and chromate in AqS1; and biotin, methionine and manganese in AqS2.

Differences also exist in the means by which the two bacteria exploit carbohydrates (Table S7). AqS2 can utilize N-acetylglucosamine and chitin, both of which are important sources of carbon and nitrogen for marine organisms (Driscoll et al., 2007; Souza et al., 2011). It is also predicted to metabolize and catabolise a larger number of monosaccharides than AqS1, with subcategories, such as D-galactonate catabolism, D-gluconate and ketogluconates metabolism, deoxyribose and deoxynucleoside catabolism, mannose metabolism and xylose utilization being only recovered in AqS2. On the other hand, AqS1 is able to use the organic acid tricarballylate, which is a citrate analog and a potential source of carbon and energy (Gutnick et al., 1969) (Table S7).

While both the AqS1 and AqS2 genomes have the capacity to independently synthesize de novo most amino acids, they also lack key enzymes required for the production of others (Figure 5A). For instance, 3-dehydroquinate dehydratase (EC 4.2.1.10), an enzyme-encoding gene that participates in the production of aromatic acids (phenylalanine, tyrosine and tryptophan), is not present in the AqS2 genome. AqS2 also lacks a homolog for phosphoribosylanthranilate isomerase (EC 5.3.1.24) that catalyzes the third step of tryptophan synthesis. Interestingly, however, these missing genes are present in the genome of AqS1, suggesting that AqS2 could synthezise all three aromatic acids provided that AqS1 metabolism supplies the intermediate substrates, 3-dehydroshikimate and N-(5-phospho-beta-D-ribosyl)-anthranilate (Table S5; Figure 5B). The alternative—that AqS2 could could directly take up the synthesized amino acid provided by AqS1—is not supported, because we are unable to detect transporters for phenalanine, tryptophan or tyrosine amongst the predicted proteins of AqS2. These observations suggest that both essential and non-essential amino acid synthesis could depend on within-pathway complementation between the two symbionts.

Sponge-associated bacteria are predicted to use characteristic eukaryote-like domain-containing proteins to facilitate host-microbe interactions; these include ankyrin-repeats, tetratricopeptide-repeats, leucine-rich repeats, fibronectin domain IIIs, and cadherins (Thomas et al., 2010; Fan et al., 2012; Hentschel et al., 2012). Our comparison of AqS1 with its free-living relatives provided support for this prediction in that the genome of AqS1 is enriched in ankyrin-repeat containing proteins (Figure 7); a similar trend was reported for T. sp. HK1 (Tian et al., 2014). On the other hand, the AqS2 genome does not show enrichment for any of these domains (Figure 7). There is little evidence of virulence factors in either of the symbiont genomes; they share just a small number of predicted protein-coding genes conferring resistance to copper, cobalt-zinc-cadmium, and/or virulence. In AqS1, the predicted proteins belonging to the virulence resistance category are involved in the production of bacteriocin Colicin V, universal stress proteins, such as heat shock protein UspA, as well as PotD, PotA and PotH. Resistance to chromium is predicted only in AqS1 (Table S5).

Finally, AqS1 and AqS2 have the genomic capacity to produce products that could be beneficial to the host, although it must be noted that many non-symbiotic bacteria share the same capacity. For instance, the AqS1 genome encodes genes associated with the synthesis of coenzyme A, as well as the essential pyrophosphate (vitamin B1), which animals must obtain from their diet (Begley et al., 2008) (Table S5). The ability to synthezise vitamin B1 has also been observed in other sponge associated bacteria and archaea (Fan et al., 2012). The two Aq symbionts also contain genes involved in pyridoxine (a form of vitamin B6), menaquinone/phylloquinone and folate biosynthesis. AqS2 is additionally equipped with genes that are involved in biotin (vitamin B7) synthesis (Table S5).

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