Introduction

The advent of microbial life on earth predates that of animals by at least 3 billion years1. Even today, microorganisms account for most of the life on our planet, both in terms of diversity and biomass2. It is therefore not surprising that animals have evolved strategies for interacting with microbes1,3. Indeed, all animals engage in stable and highly-specific associations with microbial communities and these symbioses deeply impact animal ecology and evolution1,3. The recognition of microbes as evolutionary partners has changed the way we view animal systems and has opened new frontiers of research. A prominent example is the paradigm shift in our understanding of the immune system—from the classical view as conserved defence mechanism against pathogens to the emerging perspective of immunity as rudder that allows the host to navigate the microbial world, mediating both defence and tolerance4,5,6.

A common challenge for all animals is discriminating between microbes in order to maintain a specific microbiome, while also avoiding overgrowth, harmful infections, or energetically-expensive immune reaction to innocuous microbes. Upon microbial encounter, animals detect microbe-derived molecules (microbial-associated molecular patterns, MAMPs), such as lipopolysaccharide (LPS), peptidoglycan, or flagellin, which are absent in eukaryotic organisms7,8. Pattern-recognition receptors (PRRs) of the innate immune system recognise these MAMPs and transduce a signal that activates the corresponding immune response9,10. Detection of pathogen-derived MAMPs initiates pathogen destruction11,12,13, whereas detection of symbiont-derived MAMPs promotes tolerance4,7,14,15. Even in model animals, it is not yet fully understood how the identity of the microorganism shapes the down-stream interpretation of the microbial signal detected by the PRRs. It may be related to specific MAMP structures of certain microbes (e.g.16,17,) or to accompanying danger signals in pathogenic infections18. In any case, the appropriate response relies on specific recognition and fine-tuned down-stream regulation of the immune response. Due to the absence of an adaptive immune system, three mechanisms have been proposed as molecular basis for specific recognition in invertebrates19: (i) high genetic diversity of receptors or immune effectors, (ii) enhanced expression of relevant receptors upon microbial encounter, and (iii) synergistic interactions among immune components.

Several families of animal PRRs are characterized according to the distinct arrangement of conserved protein domains. The Toll-like receptors (TLRs) are membrane-bound receptors with an extracellular domain (leucine-rich repeats in canonical TLRs) that recognizes the MAMPs and an intracellular Toll/interleukin-1 receptor (TIR) domain that triggers a well-characterized signalling cascade. This signalling cascade is present and functional in early-diverging animals20. The nucleotide-binding domain and leucine-rich repeat containing receptors (NLRs) are mainly cytosolic receptors that detect signals from microbes, tissue damage, or cellular stress21. NLR-mediated activation of the mitogen-activated protein kinase (MAPK) signalling cascade (e.g., p38, JNKs) and caspases results in reactive oxygen species formation, inflammatory processes, production of antimicrobial peptides, as well as cell death22,23. Other receptor families, such as the scavenger receptor cysteine-rich (SRCR) and lectins, add to the diverse repertoire of immune receptors found in most animals10. Another abundant and diverse class of receptors is the G-protein coupled receptors (GPCRs)24. Although they are classically omitted from the PRR group, empirical evidence supports their role in the recognition of microbial signals in both invertebrates and vertebrates24,25.

Sponges (phylum Porifera) are among the earliest-diverging multicellular animals and thus considered key to understanding the origins of animal processes, including animal-microbe interactions26,27. Due to their sessile filter-feeder lifestyle, sponges constantly encounter microbes from the seawater, which serve as a food source, but at the same time maintain stable species-specific symbiotic communities28. The field of sponge symbiosis has consolidated in recent years29,30, but it remains largely focused on the microbial side, while host mechanisms for microbial recognition and control are still poorly explored. The genome of Amphimedon queenslandica showed, for the first time, the enormous complexity of the Poriferan genomic toolkit26. It comprised a high diversity of PRRs26,31,32, including expanded NLR and SRCR families10,32. Recent genomic and transcriptomic studies in other sponge species confirmed the complex repertoire of PRRs and the presence of key components of immune signalling cascades, such as the TLR-mediated signalling pathway33,34,35. The conserved domain architectures of PRRs, their similarity to vertebrate counterparts, and the striking expansion of PRR families in sponges collectively indicate conserved functions in MAMP recognition and signal transduction36. Still, empirical evidence of such functions remains scarce37,38,39.

Here we utilized an experimental approach in order to characterise the suite of PRRs and immune genes involved in the response of sponges to microbial elicitors. We aimed to induce an immune response that would reveal the gene toolkit that is relevant for sponge immunity in the context of microbial recognition. We challenged the sponges Aplysina aerophoba and Dysidea avara with MAMPs (LPS and peptidoglycan) under controlled conditions in aquaria and assessed their response by way of RNA-Seq analysis. These two Mediterranean sponge species illustrate a long-accepted dichotomy in sponge symbiosis40—sponges termed “high microbial abundance” (HMA), like A. aerophoba, harbour symbiotic communities in densities that are two to four orders of magnitude higher than in the “low microbial abundance” (LMA) sponges such as D. avara. The HMA-LMA dichotomy involves, in addition to differential symbiont densities, differences in microbial diversity and metabolic features of the sponges41,42. Moreover, a recent genomic analysis on HMA and LMA sponge representatives from the Red Sea suggested a more expanded repertoire of immune-related domains in the LMA than the HMA sponges34. Previous works reported that sponges can rapidly take up seawater bacteria but are unable to take up their own symbionts, which suggests that sponges are capable of differentiating microbes43,44. Therefore, we hypothesise that both sponges rely on differential expression of PRRs and signalling genes to recognize and respond to MAMPs. We also expect species-specific strategies according to their different immune repertoires and HMA-LMA status.

Methods

Specimen collection

Specimens of the Mediterranean sponge species Aplysina aerophoba and Dysidea avara were collected via SCUBA diving at the coast of Girona (Spain) in March 2015 (42.29408 N, 3.28944 E and 42.1145863 N, 3.168486 E; respectively). A. aerophoba was collected at a depth ca. 3 m and the water temperature at the time of collection was 11 °C. D. avara was collected at a depth ca. 15 m and the water temperature at the time of collection was 12 °C. Collection was performed in a way that a part of the sponge remained in the substrate, allowing the regeneration of the individual. Sponges were then transported to the Experimental Aquaria Zone (ZAE) located at the Institute of Marine Science (ICM-CSIC) in Barcelona (Spain). Sponges were placed in separated 6 L aquaria in a flow-through system with direct intake of seawater and a circadian cycle of 12 h light/12 h dark using artificial light sources. Sponges were acclimated under these conditions for one week prior to experimentation.

MAMP challenge

The same experimental design was applied to each sponge species and experiments were conducted consecutively. Before the experiments, sponges were kept overnight in 1µm-filtered seawater and an additional 0.1 µm-filter was applied for 3 h before the experiments. The flow-through was stopped during the experiment and small aquarium pumps were applied to ensure mixing of the water in the aquarium. Sponges were randomly assigned to each treatment (n = 5 individuals per treatment). In the MAMP treatment, sponges were injected with LPS (source: Escherichia coli O55:B5, Sigma L2880) and peptidoglycan (source: Staphylococcus aureus, Sigma 77140) (500 µL of a final concentration 1 mg/mL in sterile artificial seawater, 1:1), with the aim of triggering an acute immune response. Sponges in control treatment were injected with sterile artificial seawater (500 µL). Treatments were directly injected into the tissue at 3–5 different spots. Sponge pumping activity was assessed visually (i.e., open oscula). For each individual, one tissue sample from one of the injection sites was collected at 1 h, 3 h and 5 h post-injection. Samples were placed in RNAlater, maintained overnight at 4 °C, and stored at −80 °C until processed. For further analysis, 3 samples per time point and treatment were randomly selected.

Extraction and sequencing of eukaryotic mRNA

Eukaryotic mRNA was obtained following the protocol described by Moitinho-Silva et al.45. Briefly, cells were mechanically lysed and total RNA was extracted using the AllPrep DNA/RNA kit (Qiagen, Germany). Contaminating genomic DNA was removed using the RQ1 RNase-free DNase (Promega, USA). RNA quantity and integrity were analyzed using InvitrogenTM QubitTM fluorometer and Experion System (Bio-Rad, USA). Sponge mRNA was isolated from ca. 100 µg of total RNA (obtained from pooling 6–10 extractions from the same biological replicate) using a Poly(A) Purist MAG kit (Ambion, USA) with two round of poly(A) purification. Library preparation (including the reverse transcription of the mRNA into cDNA) and sequencing was performed at the IKMB Kiel (Germany). The cDNA libraries were prepared using the Illumina TruSeq stranded mRNA kit and paired-end sequenced on the HiSeq. 2500 platform using HiSeq v4 reagent kit (Illumina, Inc., USA).

Data filtering, de novo transcriptome assembly and functional annotation

Given the lack of reference genomes for these sponges, a reference transcriptome was assembled de novo for each species. Raw Illumina reads were filtered to remove adapters and low-quality reads in Trimmommatic-version 0.3546 (filtering parameters- LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:75). Read quality was visualised in FastQC. Additional filtering of prokaryotic and microbial eukaryotic reads was performed in the classifier Kaiju47, in greedy-5 mode (version and database accessed in October 2016). The remaining reads of samples belonging to the same species were combined to create de novo reference assemblies in Trinity-version 2.2.048, following the general pipeline for stranded libraries. Statistics from the assemblies were obtained in Trinity and TransRate-version 1.0.249. Completeness was assessed by comparing the assemblies against the Metazoa reference data in BUSCO-version 1.2250, trans mode). Assemblies were annotated in Trinotate-version 3.0.1 (e-values < 1 e−5), a comprehensive suite that includes homology search to publicly available data (BLAST+/SwissProt), protein domain identification (HMMER/Pfam), protein signal peptide and transmembrane domain prediction (signalP/tmHMM), as well as eggNOG, GO and KEGG annotation. Those contigs with blast matches to Bacteria, Archaea, or Virus were further removed from the reference assembly. The annotation report was manually screened for the presence of the most common PRR families based on the PFAM annotation. Specifically, non-canonical TLR were identified by the presence of the TIR domain (PF01582), in combination with Ig-like domains (PF00047), NLRs by the presence of the NACHT domain (PF05729), in combination with leucine-rich repeat (LRR) domains (PF13516), and SRCRs by the presence of the SRCR domain (PF00530 or PF15494).

Transcript quantification and differential gene expression analysis

Following the Trinity pipeline, gene abundance was estimated separately for each sample by RSEM bowtie2-based quantification (version 1.2.19). Trinity outputs include the estimates for genes (Trinity components) and isoforms (Trinity transcripts). Distinguishing true isoforms from chimeras or fragmented genes remains a challenge; thus, the analysis presented here is based on gene (Trinity component) abundances. Differential gene expression analysis within each time point (i.e. 1 h, 3 h, and 5 h) was performed in edgeR (exact test mode) as implemented in the Trinity pipeline (default parameters). Differentially expressed genes (DEGs) in the MAMP compared to control treatment were defined by False Discovery Rate –corrected (FDR) p-value < 0.005 and log2|fold change| ≥ 2. For comparison, DeSeq 2 tool (as implemented in Trinity pipeline) was also tested for identification of DEGs in order to check for consistency with edgeR results. DESeq 2 found a higher number of DEGs than edgeR for the same significance threshold (Supplementary Fig. S3). Importantly, 91% and 100% of edgeR-DEGs (FDR p-value < 0.005) were consistently retrieved by DESeq 2 in A. aerophoba and D. avara, respectively. Therefore, we further explored and report here the edgeR-based results.

For a DEG annotated as a GPCR in A. aerophoba, we confirmed its presence in other sponge species by performing a blast search (at protein level, 1e−5 threshold) against a custom local database constructed from publically available transcriptomic information for 17 sponge species (Amphimedon queenslandica, Ephydatia muelleri, Haliclona amboinensis, H. tubifera, Leucosolenia complicata, Oscarella carmela, Oscarella sp., Stylissa carteri, Sycon ciliatum, Xestospongia testudinaria, Chondrilla nucula, Corticium candelabrum, Ircinia fasciculata, Petrosia ficiformis, Pseudospongosorites suberitoides, Aphrocallistes vastus, and Sycon coactum). We also searched for similar genes (blast search at protein level, e-value < 1e−5) against other marine invertebrates available in the Ensembl Metazoa database (i.e., Mnemiopsis leidyi, Nematostella vectensis, Strongylocentrotus purpuratus) and against vertebrate species available in the Ensembl database (i.e., Homo sapiens, Danio rerio, and Xenopus tropicalis). The protein alignment was built in MAFFT version 7.402 as implemented in CIPRES Science Gateway, with E-INS-i strategy and default parameters, and further visualized in Jalview Desktop51. The resulting alignment was used for phylogenetic tree construction in RAxML version 8.2.1052 within CIPRES Science Gateway, with 500 rapid bootstrap inferences and maximum likelihood search under GAMMA and WAG substitution model. The phylogenetic tree was annotated in FigTree v1.4.3 (http://tree.bio.ed.ac.uk/software/figtree/).

The set of DEGs when applying a more relaxed significance threshold, FDR p-value < 0.05, was explored via interaction network analysis in STRING-version 10.553, accessed in October 2017. We used the protein name of the top blast hit (HUGO nomenclature) of Trinotate annotation as input for STRING. STRING searches for the corresponding COG annotations and depicts a network of COG-COG interactions based on multiple types of evidences (e.g. known interactions from curated databases and experiments or predicted interactions based on gene co-occurrence and gene neighbourhood)53. We applied a minimum interaction score of 0.700 (high confidence). For A. aerophoba, two networks were created: one for the set of up-regulated genes, the other for the down-regulated genes. For D. avara, the number of annotated genes was relatively low, and therefore, a single network combining both up-regulated and down-regulated genes was created.

Results

Sequencing and de novo transcriptome assemblies

The number of paired-end Illumina reads generated in this study is summarised in Table 1. They originated from a total of 18 samples from A. aerophoba and 17 samples from D. avara, corresponding to three biological replicates per treatment within each of the three time points (except for D. avara 1 h post-MAMP treatment, for which the library construction of one replicate failed). The surviving paired reads post-filtering (Table 1) were used for generating a de novo reference assembly for each species. The statistics of the resulting reference transcriptomes are summarised in Table 2. Those contigs with similarity (blast hits) with Bacteria, Archaea, or Virus-derived sequences were removed from the reference assembly (Table 2, filtering after annotation). BUSCO assessments revealed that 69% and 70% of the 843 core Metazoan genes were detected in A. aerophoba and D.avara reference assemblies, respectively, with 21% of the genes found as fragments.

Table 1 Number of read pairs (million reads).
Full size table
Table 2 Statistics of the de novo transcriptomic assemblies.
Full size table

Diverse repertoire of putative PRRs in reference transcriptomic assemblies

Based on the presence of conserved domains (Pfam annotation), we identified putative PRRs within the families of non-canonical TLRs, NLRs, and SRCRs in the reference transcriptomes of A. aerophoba and D. avara. Bona fide NLRs are characterised by the presence of NACHT and leucine-rich repeat (LRR) domains (as in Yuen et al.32). In the A. aerophoba reference transcriptome, only one gene (Trinity component TR172818_c2_g1) showed this architecture across a complete open reading frame (Supplementary Table S1). However, 75 additional genes contained a NACHT domain and could potentially belong to the NLR family but lacked the LRR domain (Supplementary Table S1). In D. avara, 80 bona fide NLR genes were detected of which 39 were complete (Supplementary Table S2). The number of additional NACHT domain-containing genes in the D. avara reference transcriptome extended to 390 (Supplementary Table S2). The reference transcriptomes of A. aerophoba and D. avara also included >250 genes containing single or multiple SRCR domains, sometimes in combination with other conserved domains such as fibronectin III, protein kinases, Sushi repeats, or epidermal growth factor-like domains (Supplementary Tables S1, S2). While sponges lack bona fide TLR, they do contain Immunoglobulin-TIR receptors characterised by an intracellular TIR domain (which is homologous to the TIR domain in TLR in Eumetazoan54) but with immunoglobulins instead of LRRs as extracellular domain31. We detected a single gene in A. aerophoba (TR170373_c0_g1, Supplementary Table S1) and two genes in D. avara (TR163581_c0_g2 and TR169736_c5_g2, Supplementary Table S2) presenting this architecture. In addition, KEGG annotation identified components of the TLR signalling pathway (Supplementary Figs S1, S2), as reported in other sponge species35,54.

Transcriptomic profiles in response to MAMPs

Overall, 83.35 ± 0.21% and 82.17 ± 0.26% of the reads in the samples aligned to the corresponding transcriptome reference in A. aerophoba and D. avara, respectively (average ± standard error). Next, gene expression levels in MAMP challenge treatment were compared to those in the control treatment at each time point (1 h, 3 h, and 5 h). DEGs were defined by log2|FC| ≥2 (4-fold change) and FDR p-value < 0.005. The DEGs were classified as up-regulated or down-regulated in the MAMP treatment when compared to expression levels in the control treatment. Overall, a higher number of DEGs was detected in A. aerophoba than in D. avara (Fig. 1). A total of 235 and 249 genes were identified as up-regulated and down-regulated, respectively, in A. aerophoba. In D. avara, the total number of DEGs was 29 up-regulated and 20 down-regulated. Most DEGs detected within a sponge species were time-specific (Fig. 1). In A. aerophoba, the highest number of DEGs was detected 3 h after MAMP challenge. In D. avara, the highest differential expression occurred 1 h after treatment; but only 2 replicates from the MAMP treatment were available for this time point, which could have influenced the observed trend. Heatmaps illustrate the consistency of DEG-expression profiles among biological replicates in each treatment and time point (Fig. 2). The full results from the differential expression analysis in edgeR are reported in Supplementary Tables S3 and S4 and the full annotation report for DEGs is available as Supplementary Tables S6 and S7.

Figure 1
figure 1

Numbers of DEGs those were either common or specific for each time point (1 h, 3 h, 5 h) in each sponge species upon MAMP treatment. Within each time point, DEGs were identified by comparing gene expression levels in MAMP relative to control treatment and according to the defined threshold FDR p-value < 0.005 and log2|FC| ≥ 2 expression, as calculated in edgeR. “Up-regulated” and “down-regulated” refers to genes with higher and lower expression in MAMP than in control treatment, respectively.

Full size image
Figure 2
figure 2

Heatmaps of differentially expressed genes upon MAMP treatment. Heatmaps show relative expression levels of each DEG (rows, hierarchically clustered) in each sample (columns) from A. aerophoba (left) and D. avara (right). DEGs are defined by FDR p-value < 0.005 and log2|FC| ≥ 2 expression (4-fold change), as calculated in edgeR. Expression values are log2-transformed median-centred TMM-normalised values.

Full size image

PRR expression and signalling in response to MAMPs (FDR p-value < 0.005)

Based on Pfam domain architectures, several putative PRRs were identified as differentially expressed in response to the MAMP challenge (Table 3). In A. aerophoba, the repertoire of receptors that were differentially expressed included one gene with a SRCR domain (TR13528_c0_g1, partial gene). We also include in this category a gene identified as a GPCR by the presence of a GPS motif (PF01825: GPCR proteolytic site). Further phylogenetic analysis of this gene suggests that it belongs to the group of adhesion GPCRs, with similarity to the vertebrate group I (ADGRL2 genes, also known as latrophilin-2) (Fig. 3). In D. avara, bona fide NLRs were significantly up-regulated upon MAMP challenge (Table 3). Within them, the TR172577_c0_g1 gene was among the 10 highest differentially expressed genes at each time point (in terms of fold change and FDR p-value) and contained a predicted transmembrane domain (Supplementary Table S6). Also, a leucine-rich repeat-containing gene and several genes containing fibrinogen-related domains were differentially expressed and included as putative PRRs (Table 3). The fibrinogen domain containing genes showed similarity to vertebrate ficolins and angiopoietin-related genes (blastp, e-value < 1e−5). Fibrinogen-like proteins have been proposed as potential immune receptors in molluscs and other invertebrates55. Potential receptors according to sequence similarity, but without the corresponding conserved domains, are included in Tables 4, 5 and Supplementary Table S5.

Table 3 Differential expressed genes identified as immune receptors in A. aerophoba and D. avara, according to the presence of conserved domains.
Full size table
Figure 3
figure 3

Phylogenetic analysis of the A. aerophoba GPCR up-regulated gene TR165761_c4_g1. A part of the alignment is reported. The star (*) denotes the beginning of the GPS motif. Conserved residues (identical in all sequences) are shown in dark blue, and those identical in at least 50% of the sequences are in light blue. A schematic representation of the domain architecture of each gene is provided. As TR165761_c4_g1 gene is incomplete, we removed the 7tm domain from the other protein sequences included in the alignment prior to tree construction by maximum likelihood analysis. Node labels represent bootstrap support greater than 50% of 500 pseudoreplicates.

Full size image
Table 4 Differentially expressed genes (FDR p-value < 0.005) in A. aerophoba.
Full size table
Table 5 Annotated DEGs (FDR p-value < 0.005) in D. avara.
Full size table

Genes involved in signal transduction (e.g., kinases), chaperones (i.e., hsp70), and genes related to adhesion and extracellular matrix were differentially expressed upon MAMP challenge in both species (Tables 4 and 5). We also detected differential expression of genes related to ubiquitination (i.e., ubiquitin ligases) and apoptosis (Tables 4 and 5). In A. aerophoba (Table 4), the set of DEGs included genes with conserved domains such as ankyrin repeats, immunoglobulin domains, Sushi and fibronectin III domains or tetrapeptide repeats that could be involved in recognition, adhesion, and cell-cell interactions. The A. aerophoba gene TR175974_c14_g10, which was identified as a GPCR by sequence similarity but not by Pfam domain architecture, was therefore excluded from Table 3 and included in Table 4. According to blast results, several genes potentially involved in GPCR signalling were also significantly differentially expressed upon treatment in this sponge (Table 4). Signalling transduction in A. aerophoba was further mediated by a DEATH -domain containing gene as well as by several mitogen-activated protein kinase kinase kinases (MAPKKK), which were all down-regulated (Table 4). In D. avara, the genes involved in recognition, adhesion and cell-cell communication were all up-regulated (Table 5). Signalling transduction was mediated by protein kinases and serine/threonine protein kinases, which were up-regulated too (Table 5). DEGs related with apoptosis were up-regulated 1 h post-treatment in D. avara. And this sponge up-regulated a gene annotated as phospholipase D, which may be involved in lipid metabolism and in the phosphatidylinositol signalling pathway.

DEGs in A. aerophoba included genes with functions in metabolic processes (Table S5), such as lipid metabolism (e.g., long-chain-fatty-acid-CoA ligases). Other functions under regulation in this species were chromatin remodelling and transcription (e.g., differential expression of DNA-binding proteins and transcription factors) (Table S5). Also, a gene with similarity to Dictyostelium discoideum DD3-3 gene (DDB_G0283095) was up-regulated 3 h after treatment (Table 4). Homologs of this gene are present in other invertebrates, including cnidarians and echinoderms, but are absent in Vertebrata. In D. discoideum, a DD3-3 knockout yields faster cell aggregation than in the wild type and compromised cAMP signalling pathway56. Another DEG in A. aerophoba contained a Reeler domain (PF02014), similar to insect defence proteins (Table 3), which may have antimicrobial activity. Several genes remained unidentified due to a lack of similarity with genes in public databases or conserved domains. For example, in A. aerophoba, the gene TR170260_c3_g2 was within the top DEGs at all time points (in terms of fold change and FDR p-value) and was identified as a non-transmembrane signalling peptide but no further annotation was available for this gene. Several DEGs within D. avara which lack annotation were identified as signalling peptides (Supplementary Table S7).

COG network analysis (FDR p-value < 0.05)

We also explored the set of DEGs when a more relaxed significance threshold was applied (FDR p-value < 0.05; annotation in Supplementary Tables S6 and S7) to probe for further support of the biological processes activated upon MAMP treatment. In both species, the complex network represented a signalling cascade mediated by kinases (Figs. 4 and 5). In A. aerophoba, the groups of serine-threonine protein kinases (COG0515) and the ankyrin repeat-containing genes (COG0666) occurred in multiple interactions in both the up-regulated and the down-regulated networks (Fig. 4). In the network of up-regulated genes (Fig. 4, left side), the central nodes (in terms of number of interactions) were leucine-rich repeat proteins (COG4886) and transcription factors involved in chromatin remodelling (COG5076). In the network of down-regulated genes (Fig. 4, right side), the category of phosphatidynositol-3 (PI-3) kinases (COG5032) was also a central node and it connected with other kinases as well as with a network of genes related with lipid metabolism (COG1022; COG1024; COG1562; COG4281). In D. avara, up-regulated and down-regulated genes were analysed in a single network (Fig. 5). Serine-threonine kinases (COG0515), as well as the category of leucine-rich repeat proteins (COG4886) were the COGs with the highest number of connections (Fig. 5). They interact with each other and with other protein groups, including GTPases (COG1100), and to COGs related to extracellular matrix (Fig. 5).

Figure 4
figure 4

COG association network analysis from annotated differentially expressed genes in A. aerophoba upon MAMP treatment. Networks of up-regulated and down-regulated DEGs (FDR p-value < 0.05; log2|FC| ≥ 2), as obtained in STRING. The reference protein names identified in Trinotate were used as input. STRING searches for COG annotations and calculates and depicts the association network. Edges represent protein-protein associations coded by colour according to the type of evidence for the shown interaction (see legend). Minimum required interaction score: 0.700 (high confidence). NOG means “non-categorised orthologous group”.

Full size image
Figure 5
figure 5

COG association network analysis from annotated differentially expressed genes in D. avara upon MAMP treatment. Network of annotated DEGs (FDR p-value < 0.05; log2|FC| ≥ 2), as obtained in STRING. The reference protein names identified in Trinotate were used as input. STRING searches for COG annotations and calculates and depicts the association network. Edges represent protein-protein associations coded by colour according to the type of evidence for the shown interaction (see legend). Minimum required interaction score: 0.700 (high confidence). The network includes both down-regulated and up-regulated genes. NOG means “non-categorised orthologous group”.

Full size image

Discussion

We investigated the transcriptomic profiles of two Mediterranean sponge species upon MAMP exposure (LPS and peptidoglycan). Previous genomic information for A. aerophoba and D.avara was lacking; thus, this study provides a valuable resource with the generation of a de novo-assembled reference transcriptome for these species. The reference transcriptomes of A. aerophoba and D. avara contain a complex inventory of PRRs. Both species harbour hundreds of genes containing single or multiple SRCR domains, sometimes in combination with other conserved domains such as fibronectin III or immunoglobulin domains. In D. avara, 80 bona fide NLRs are found in the reference transcriptome. In the A. aerophoba reference transcriptome, only one gene could be identified as a bona fide NLR and it was constitutively expressed in all samples. However, several incomplete transcripts contained NACHT domains and could potentially add to the repertoire of expressed NLRs in this species. The NLRs represent a PRR family that is highly expanded in the A. queenslandica genome (comprising 135 genes, in contrast to 20 genes in humans)32; however, the reference transcriptome of the sponge Vaceletia sp. lacks these receptors35. Both A. aerophoba and D. avara constitutively express Immunoglobulin-TIR receptors, as found in other sponges54. In organisms with limited amenability to genetic manipulation, such as sponges, gene function is typically inferred from data from distantly-related organisms as validation of functions is challenging27. Consequently, the set of Poriferan-unique and species-specific traits remain misrepresented27,57. Nevertheless, by adopting an experimental approach, we have identified receptors and other genes that are potentially relevant to the sponge response to microbes and have narrowed the list of target genes for future research.

MAMPs (mainly LPS, but also peptidoglycan or flagellin) have been broadly used as immune activators in multiple organisms (including plants, invertebrates, and vertebrates)8,55,58,59. The MAMP-triggered immune pathways are considered, besides physical barriers, as the first line of the response to microbes. As filter-feeders, sponges constantly encounter diverse microbes carrying different MAMPs. To increase the chances of inducing an immune response, we chose here commercially-available MAMPs (LPS and peptidoglycan) derived from non-marine organisms. We applied them simultaneously to increase the array of transcriptionally inducible PRRs and pathways in the same treatment. For example, Zhang et al.55 showed a stronger transcriptomic response (more number of DEGs) to LPS than to peptidoglycan and fucoidan in the snail Biomphalaria glabrata. Similarly, Weiss et al.60 reported little overlap in the transcriptomic response of the coral Acropora millepora to muramyl dipeptide and poly I:C as MAMPs. The MAMP challenge is preferable over challenge with live cells when the aim is to induce the transcriptionally inducible PRRs and their activated downstream response because interference with microbial-derived effector molecules is avoided61. We thus consider the MAMP challenge approach meaningful for unveiling animal-microbe molecular talk, although future studies addressing other microbial challenges would help to further identify the underlying molecular mechanisms.

In invertebrates, a high diversity of PRRs and their tuned expression upon microbial stimuli has been proposed as a mechanism for specific recognition of microbes10,19,36,62. Here, we detected sponge species-specific signatures in the expression profiles of these PRRs upon MAMP challenge (Table 3). A SRCR domain-containing gene was up-regulated in A. aerophoba in response to MAMPs (Table 3). In A. queenslandica juveniles, more than 30 SRCR domain-containing genes with diverse architectures were differentially expressed upon exposure to microbes in aquaria experiments38. The implication of SRCR on microbial recognition in sponges was first evidenced by the upregulation of a SRCR-domain containing gene in symbiotic vs aposymbiotic (i.e., cyanobacteria-free) Petrosia ficiformis in the field37. SRCR-domain containing genes are also expanded in echinoderm genomes as well as being highly expressed in their immune cells and activated in response to microbes63,64. Further studies have reported the up-regulation of these receptors upon bacterial exposure in other invertebrates65. In D. avara, two NLRs were differentially expressed upon MAMP treatment. The complex repertoire of NLRs in A. queenslandica already hinted towards their role in microbial recognition in sponges36, but our findings provide the first experimental evidence of enhanced expression of poriferan NLRs in response to microbial cues. Evidence of the role of NLRs in invertebrates is scarce66. However, in vitro studies in the cnidarian Hydra showed that a non-conventional NLR genes (lacking the LRR domain) are differentially-expressed in response to LPS and flagellin stimulation and yield the activation of caspases in a manner that may be analogous to the mammalian inflammasome67.

Our study also revealed other putative immune receptors. GPCRs were differentially expressed in both A. aerophoba (up-regulated; Table 3, Supplementary Table S6) and in D. avara (down-regulated; Supplementary Table S7). The phylogenetic analysis of the A. aerophoba gene TR165761_c4_g1 showed that it belongs to the adhesion GPCR family (Fig. 3), which is involved in adhesion and signalling. Krishnan et al.68 also classified a group of A. queenslandica adhesion GPCRs as basal of human Group I and Group II adhesion GPCRs, whereas the rest of A. queenslandica adhesion GPCRs were either sponge specific or more similar to other vertebrate GPCR families. GPCRs constitute a highly diverse receptor family in animals25,69, including sponges68,70. In vertebrates, they take part in crosstalk with microbes, by detecting microbial-derived metabolites (e.g., short-chain fatty acids) and interacting with other PRRs such as TLRs25,71. In invertebrates, their role in defence has been suggested for Caenorhabditis elegans and Drosophila melanogaster24. In addition, RNA-Seq analysis revealed that GPCR signalling played a role in the response of the sea anemone Aiptasia to symbiotic states and Symbiodinium type72. Thus, our results provide additional support for the conserved role of GPCRs in animal-microbe interactions. In D. avara, there is furthermore a noteworthy differential expression of several fibrinogen-domain containing genes. This domain is commonly found in the DEGs responding to microbial cues in invertebrates55,73,74. In addition, both species differentially expressed several genes containing immunoglobulin domains, LRR domains, DEATH domains and genes with sequence similarity to lectins (e.g. galectin). Besides their roles in cell-cell communication75, these domains are common in immune receptors76 and are involved in microbial recognition in corals77, snails78, or nematodes79. Moreover, a ficolin-like gene was up-regulated in the sponge Cliona varians when “reinfected” with Symbiodinium compared to the aposymbiotic tissue80. Therefore, GPCRs, fibrinogen-containing and lectin-like genes could add to the repertoire of genes key for immune recognition in sponges.

The response of both sponges to MAMPs involved the up-regulation of ankyrin repeat-containing genes, immunoglobulin-domain containing genes, DEATH-domain containing genes, CARD-domain containing genes and chaperones (hsp70), as well as regulation of collagen. Signalling transduction was also mediated by serine-threonine protein kinases, which were significantly down-regulated in A. aerophoba but up-regulated in D. avara. The network analyses in STRING (Figs 4, 5) show that the information available from other organisms supports the co-expression patterns reported in our study, but further studies on co-localization analysis and protein-protein interactions would be necessary to confirm these networks. These MAMP-triggered transcriptomic profiles resemble those found in other invertebrates55,59,74,81 and potentially mediate a high diversity of cellular responses, such as cell death81, phagocytosis82, and metabolism regulation55. Here, the activation of apoptosis in both species is indicated at the earliest time point (1 h). Moreover, the enhanced expression of a folate receptor (Table 4), SRCR and GPCR (Table 3) in A. aerophoba together with the differential expression of Ras family gene, dynamin and genes involved in cytoskeleton rearrangement (Table 4; Fig. 4) hints to the activation of a phagocytic response in this sponge species83,84.

We did not detect differential expression of genes encoding Immunoglobulin-TIR receptors or its adaptor protein MyD88 (myeloid differentiation primary response 88), even though both sponge species investigated here constitutively expressed Immunoglobulin-TIR domain receptors (Supplementary Tables S6, S7) and the MyD88-dependant downstream pathway (Supplementary Figs S1 and S2). In contrast, other sponge species activated MyD88 gene in response to LPS or microbes38,85. In Suberites domuncula, MyD88 expression was up-regulated 12 h after exposure to the same E.coli-derived LPS we used in our study85. However, before treatment, these sponges were kept in cultivation for a long period of time and their symbiotic bacterial load was reduced85, which could affect the immune reaction. Also, the combination of LPS and peptidoglycan may be a reason for the different responses reported in our study. In A. queenslandica juveniles, the up-regulation of Immunoglobulin-TIR receptors and components of the signalling pathway (including MyD88) was induced 2 h after exposure to bacteria38. The different results may be due to species-specific strategies, time-dependent responses, or the different experimental design (e.g., challenge with different microbial elicitors, different sampling points, or the use of adults vs juveniles).

The two species investigated here exemplify the HMA-LMA dichotomy in sponges, defined by differences in symbiont density and diversity40,42. Previously, Ryu et al. observed that SRCR, NLRs, Immunoglobulin-like, and fibronectin-3 containing genes were more abundant in the genomes of LMA than HMA sponges34, while Germer et al. found NLRs to be absent from the transcriptome of the HMA sponge Vaceletia sp35. Similarly, we observed a more abundant repertoire of NLRs in D. avara (LMA) than A. aerophoba (HMA) transcriptomes. However, comparative genome analysis would be necessary for further confirmation of this pattern between HMA and LMA sponges. In our study, both species showed certain similarities in the response to MAMPs; for example they activated apoptotic processes in the immediate response (1 h after treatment). However, the repertoire of PRR genes involved differed between species and the magnitude of the transcriptionally-regulated response (in terms of the number of DEGs) was more complex in A. aerophoba than in D. avara. In particular, further regulation of genes related with transcription and phagocytosis account for the greater transcriptomic response in A. aerophoba. These differences may point to species-specific features. For example, coral immune responses to LPS challenge and to thermal stress differ significantly depending on the species considered61,81. However, they may also reflect different immune strategies according to their differing HMA-LMA status. We propose that HMA sponges require a more fine-tuned regulated response to deal with potential conflicts between the signals from the MAMP stimulation and the symbiotic feedbacks from their highly dense microbial community. In line with the Danger model of immunity18, we further hypothesize that the host danger signals released upon apoptosis subsequently trigger an enhanced immune response and phagocytic activity. This hypothesis is supported in A. aerophoba by an increased expression of apoptosis genes after 1 h and of phagocytosis-related signalling pathways after 3 h of MAMP challenge. Further studies including more HMA-LMA species are on-going to elucidate whether the HMA-LMA status contributes to the variation in immune responses to microorganisms among sponge species.

Conclusions

The characterization of the innate immune response through experiments and functional studies remains limited to few animal groups and was previously lacking in the phylum Porifera. We exposed two Mediterranean sponge species to MAMPs (LPS and peptidoglycan) and described, to our knowledge for the first time, the response of the sponges to immune stimuli by RNA-Seq. The sponges responded by increased expression of a subset of relevant receptors (i.e., NLRs in D. avara, SRCR and GPCRs in A. aerophoba) and the transduction of signals by kinase cascades that likely yield apoptosis and regulation of metabolic processes. In addition, the magnitude of the transcriptomic response was higher in A. aerophoba and this was related to the regulation of additional processes such as phagocytosis. The differences between species in the subset of regulated receptors and pathways when exposed to MAMPs may relate to their different symbiont load (HMA/LMA status). We propose that the presence of a highly dense symbiotic community in A. aerophoba influences the signalling feedbacks and determines the more complex transcriptomic response upon MAMP challenge in this species. Our findings address a prominent gap in marine sponge research by providing novel information on the repertoire of genes involved in immune recognition and signalling in this ancient animal phylum.