Oceanographic setting influences the prokaryotic community and metabolome in deep-sea sponges

Marine sponges (phylum Porifera) are leading organisms for the discovery of bioactive compounds from nature. Their often rich and species-specific microbiota is hypothesised to be producing many of these compounds. Yet, environmental influences on the sponge-associated microbiota and bioactive compound production remain elusive. Here, we investigated the changes of microbiota and metabolomes in sponges along a depth range of 1232 m. Using 16S rRNA gene amplicon sequencing and untargeted metabolomics, we assessed prokaryotic and chemical diversities in three deep-sea sponge species: Geodia barretti, Stryphnus fortis, and Weberella bursa. Both prokaryotic communities and metabolome varied significantly with depth, which we hypothesized to be the effect of different water masses. Up to 35.5% of microbial ASVs (amplicon sequence variants) showed significant changes with depth while phylum-level composition of host microbiome remained unchanged. The metabolome varied with depth, with relative quantities of known bioactive compounds increasing or decreasing strongly. Other metabolites varying with depth were compatible solutes regulating osmolarity of the cells. Correlations between prokaryotic community and the bioactive compounds in G. barretti suggested members of Acidobacteria, Proteobacteria, Chloroflexi, or an unclassified prokaryote as potential producers.


Overview
Here, we reported water masses as approximated by water depth to influence the prokaryotic community composition and holobiont metabolome in three deep-sea sponge species, Geodia barretti, Stryphnus fortis and Weberella bursa. By 16S rRNA gene sequencing and computational statistical analyses, we assessed and analysed the operational taxonomic units (ASVs) that increased or decreased with depth and found distinct groups of ASVs in samples associated with two water masses (Labrador Sea Water and Irminger Current) in the Davis Strait, North Atlantic. In parallel, by Ultra-high Performance Liquid Chromatography coupled with High-Resolution Mass Spectrometry (UPLC-HRMS) we acquired and analysed the metabolome of these sponge holobionts. Using computational analyses, we found that the metabolomes changed with depth. We extracted and identified metabolite signals driving the depth-dependent changes as well as signals of known bioactive compounds. Finally, we combined microbiota and metabolome data to generate a hypothesis about the producer of the antibiofounling compound barettin. This multi-omics data acquisition and analysis approach relied on the use of several instruments and methods that are detailed in the experimental section.
The supplementary methods described herein apply to the wet lab procedures. The code and comments for the subsequent computational analyses were deposited at the Paamiut GitHub repository.
The large amount of data produced in this study had several additional aspects and implications which were presented and discussed in the supplementary results section.
1 % formic acid (FA), and mobile phase B consisted of 40:60 MeCN/MQ water with 5 mM ammonium formate and 0.1 % FA. The gradient elution profile was as follows: mobile phase A was decreased non-linearly (slope factor 8, MassLynx) from 100 % A to 100 % B over 14 min, 100 % B was held for 2 min and then decreased back to 100 % A over 1 min. The column was re-equilibrated at 100 % A for 6 min for a total runtime of 23 min. Chromatographic separation in RP was performed on an Acquity UPLC BEH C18 column (1.7 m, 2.1 mm inner diameter × 50 mm, Waters Corp.). Mobile phase A consisted of MQ water with 0.1 % FA, and mobile phase B was MeCN with 0.1 % FA. The gradient elution profile started at 95 % A, was decreased linearly over 14 min to 5 % A, and 5 % A was held for 2 min before the column was re-equilibrated at 95 % A for 4 min. The flow rate was set to 0.4 mL/min, the column temperature was set to 40°C, the samples were kept at 8°C and the injection volume was 5 L in all experiments.

Data acquisition
Data acquisition was performed using MS E mode, and lock mass correction was applied using a solution of leucine enkephalin in both positive and negative mode. Ionization parameters were set as follows in positive/negative mode; the capillary voltage was 1kV/1.5 kV, the cone voltage was 30 V/25 V, the source offset was 50/60 and the source temperature was set to 120°C. Nitrogen was used as desolvation and cone gas with gas flows of 800 L/h and 50 L/h, respectively, and desolvation temperature was set to 500°C /450°C. For MS E acquisition a collision energy ramp from 20-45 eV was used with argon as collision gas.
The instrument was calibrated in the m/z range 50-1500 using sodium formate prior to each analysis. All study samples were analysed in both RP and HILIC, in positive and negative ionization mode, resulting in four metabolite datasets per sponge specimen. The column and sample cone were cleaned in between each analysis mode. Prior to each analysis 10 QC injections were made to condition the column, and to ensure stable retention times and signal intensities. The study samples were analysed in randomized order with QC injections interspaced every sixth injection.

MS data processing
Raw files were converted to netCDF files by Databridge (part of MassLynx, Waters Corporation, Milford, Massachusetts, USA). The netCDF files with the chromatographic spectra were sorted into folders according to species and processed with XCMS in R. Peak picking was performed using the centWave function with parameters ppm=8, peakwidth set to c(5,45) and the noise parameter set to 2000. Retention time alignment was performed with the obiwarp function and the response factor set to 10, grouping was performed with the "group" function and the "fillPeaks" function was used to impute a signal in cases where no matching 3 pseudospectra were detected. The data set was curated to remove features eluting in the void (retention time less than 45 s). A raw data set as well as two normalized data sets (Log10-transformed and median fold change normalized) were produced and filtered to only retain features with a coefficient of variation < 30% in the QC samples. After subsequent evaluation, raw (i.e. untransformed) data sets were used in subsequent statistics and modelling.

Peak picking with XCMS and annotation with CAMERA
We processed samples from all three sponnge species in random order with interspersed injection of a combined QC sample to monitor stability of the UPLC-HRMS run. The acquired signals/spectra were converted to netCDF format using the Program DataBridge, and thereafter sorted into four folders, three for the sponge species (Gb, Sf, Wb) and one for the QC samples (QC). Peak picking and combination of pseudospectra is performed with the R package xcms, the subsequent annotation of adducts and isotopes with the R package CAMERA. Code for XCMS and CAMERA data processing is found at the Paamiut GitHub repository s, 50°C for 20 s, 72°C for 20 s and a final extension of 7 min at 72°C. PCR products were visualised on a 1% (w/v) agarose gel. Five L of the PCR products were used as template in the second PCR reaction to incorporate eight nucleotide sample-specific barcodes, as previously described [3]. This second PCR was performed in triplicate for each sample. The PCR reactions contained 31 L nuclease-free water (Promega), 10 L of 5× HF buffer, 0.5 L of 2 U/ L Phusion hot start II high fidelity polymerase (Thermo Fisher Scientific), 5 L equimolar mixes of 10 M forward primer (barcode-linker-Unitag1) and reverse primer (barcode-linker-Unitag2), 1 L 10 mM dNTPs (Promega) and 2.5 L of the first PCR product as template for a total of 50 L.

Annotation of SpongeEMP according to Dat et al.
The sequences of the abundant ASVs for each sponge species (relative abundance > 0.25%) were subjected to a BLAST search against a curated sponge Earth Microbiome Project (EMP) database (https://github.com/ amnona/SpongeEMP). The sponge microbiome project subASV sequences that were identical or had one nucleotide mismatch with sequences of the most abundant ASVs were uploaded to the spongeEMP online server (www.spongeemp.com) to identify ASVs that were significantly enriched in sponges [4].

Metabolomics: OPLS models
OPLS models are prone to overfitting when there are more variables than samples. This premise is met by our data sets that have vastly more variables (metabolomic features) than sponge samples. Therefore, in order to validate the model, a sevenfold cross validation is performed in the process of generating it, i.e. when calling the opls() function. The predictive power of a model is denoted Q2. For the actual data model it is "Q2 cum", and for the cross-validation models with permutated data it is "pQ2". In case of overfitting models (effectively modelling random noise), the predictive power of the cross-validated models should be similar or higher. In our case, in Tab. S3, we see that in two models ( S. fortis: RP neg: ion, and G. barretti: RP neg: ion) the predictive power of permutated data is higher than the original data. There are few additional cases where models on permutated data performed similar to original data ( S. fortis: HILIC pos: ion/RP pos: ion, W. bursa: RP pos: cleaned/pc group). However, there is no exact cut-off to define overfitting. We conclude thus that for most of the models presented in the table, the models seem to model data and not noise and that our conclusions that depth affects the metabolome are valid.

Metabolomics: Choline and serotonin/histamine
Biofouling inhibition is an activity tested and shown for several of the compounds from G. barretti and S. fortis belonging to different compound classes (diketopiperazines, peptides, bromotyrosin-derivatives [5], [6]. Two hypothetic modes of action have been proposed for the biofouling inhibition of barettin and 8,9dihydrobarettin [7]. One of these acts via 5-HT serotonin receptors as especially barettin has a high affinity to some of the receptor subtypes [8] which are also present in barnacles. The other mode involves inhibition of acetylcholinesterase (AChE), as barettin, 8,9-dihydrobarettin and stryphnusin inhibit the enzyme at concentrations close to other inhibitors of marine origin [7], [9]. AChE in turn was found in larval organs involved in obtaining settlement cues [10].
Given the implicated effect on choline; and the hypothesised connection with serotonin, we isolated their signals in the metabolomes and correlated them with depth. Serotonin was absent or present in very low relative abundance in G. barretti, present in low amounts in W. bursa and in low to moderate relative abundance in S. fortis, whereas choline was abundant in all three sponges (Fig. S9). Serotonin ( =0.58, p=0.04) and choline ( =0.72, p=0.0075) increased with depth in S. fortis. (Fig. S9). Serotonin [11], [12] and serotonin-derived compounds [13] have been previously identified with NMR in demosponges but genes coding for serotonin receptors and serotonin-synthesizing enzymes have so far not been found in sponge genomes [14]. The near absence of serotonin in G. barretti and the possible absence of serotonin receptors in sponges suggest that barettin targets are more likely to be macro-or microfoulers with serotonin-like receptors, and not the sponge itself, but its presence in the other sponges does not rule out the possibility entirely. As serotonin is a neurotransmitter, we further searched for histamine (Fig. S9). The signal intensity of histamine did not correlate with depth in any of the sponges.

Metabolome: uranidine
One VIP compound that was tentatively identified as Uranidine, a quinolone yellow pigment, originally isolated from the yellow shallow water sponge Aplysina aerophoba [15]. Uranidine was shown to have anti-HIV activity [16]. The compound has never been recorded outside of the verongiid sponges and, if it is confirmed, this is only the third record of this compound, this time in the deep-sea sponge G. barretti. We should stress that this compound was identified based on m/z 194.0448 [M+H]^+ and fragments (MSMS), without a reference. Uranidine is known to be an unstable compound that oxidises easily to a blue quinone that polymerizes and turns black [15]. On damaged areas of G. barretti, purple spots can be observed on the usually whitish surface and whitish to peach coloured choanosome (P. Cárdenas, pers. obs.). The reasons for the decrease of this compound with depth is currently unclear.

Metabolome: osmolytes
The role and function of the poorly understood arsenobetaine in marine environments has been investigated previously [17]. Arsenate might be taken up indiscriminately from the marine environment through phosphate membrane transport systems of microbes. Due to the structural similarity with betaine, another osmolytes, it has been hypothesised to have the same function. However, the signals of both compounds in the three different sponges as well as their depth responses did not warrant unequivocal conclusions (Fig.   S9). In G. barretti arsenobetaine decreased with depth (r=-0.705, p= 0.003) while betaine increased (r= 0.589, p= 0.021). In S. fortis, the signal remained stable and in W. bursa, arsenobetain increased with depth (r= 0.625, p= 0.01). The opposing responses to depth might be due to differences in the holobionts' depth response. However, if indiscriminately taken up by the environment, it would be expected that the signals of betaine and arsenobetaine would correlate. This was only the case in W. bursa (r= 0.631, p= 0.009).

Microbiome
Thaumarchaeota are one of the important constituents of sponge-associated prokaryotic communities, associated with nitrogen cycling, and particularly ammonia oxidation [18], [19]. ASV146 and ASV138 have only one bp difference, ASV138 is identical (140 bp query coverage) to a Thaumarchaeota clone which belongs to a sponge-enriched clade including Cenarchaeum symbiosum [20]; this clone was present in the sponges G.
barretti and Phakellia ventilabrum (Bergen area, Norway, 200-300 m) (e.g. JQ612474), as well as Stelletta normani KF597128 (Porcupine Bank, 1350 m). ASV146, showed generally higher relative abundances in specimens collected above 1000 m, but was absent in specimens collected below 1000 m, in both G. barretti (Fig. 4) and S. fortis (average relative abundance; shallow Gb 0.20 and Sf 0.14 versus deep Gb 0 and Sf 0.004). Conversely, ASV138 relative abundance is significantly increasing with depth, and is more abundant below 1000 m , but only in G. barretti (average relative abundance; shallow Gb 0.21 and Sf 0.16 versus deep Gb 0.55 and Sf 0.02) (Fig. 4). The inverse relationship of the relative abundance of the two crenarchaeal ASVs 138 and 146 with depth indicates that they are "sister-strains" and may represent two ecotypes from the two water masses: an IC-derived sponge-specific ecotype and a deeper LSW-derived sponge-specific ecotype.
In a similar way, the Thaumarchaeota ASV004 (genus Candidatus Nitrosopumulis) has several sister-strains with opposite trends. ASV004 abundance is decreasing with depth in G. barretti, while ASV097 (1 bp difference) increase with depth in G. barretti. The fact that ASV004 matches the sequence of the surface water mass Thaumarchaeota ecotype [20] (16S sequence courtesy of O. Müller) suggests that these sponge ASVs are this time seawater derived and may reflect the water mass Thaumarchaeota community [20]. Hence, both sponge-enriched and seawater-derived ASVs are related to the water masses, and lead to similar overall phylogenetic composition of prokaryotic communities at the species level.

Supplementary figures
All figures were produced with R v 3.5.   Barrettin: filled= positiv significant correlation, empty=negative correlation and/or not significant, absent=not enought data Figure S4: Annotated cladogram of all ASVs in G. barretti. Green arches indicate sister ASVs (n=19) as defined by t-test (pFDR<0.5) comparing relative abundance above and below 1000 m. The figure was produced with iTol and deposited online iToL (https://itol.embl.de/shared/karin_steffen).