Differential expression of immune receptors in two marine sponges upon exposure to microbial-associated molecular patterns

The innate immune system helps animals to navigate the microbial world. The response to microbes relies on the specific recognition of microbial-associated molecular patterns (MAMPs) by immune receptors. Sponges (phylum Porifera), as early-diverging animals, provide insights into conserved mechanisms for animal-microbe crosstalk. However, experimental data is limited. We adopted an experimental approach followed by RNA-Seq and differential gene expression analysis in order to characterise the sponge immune response. Two Mediterranean species, Aplysina aerophoba and Dysidea avara, were exposed to a “cocktail” of MAMPs (lipopolysaccharide and peptidoglycan) or to sterile artificial seawater (control) and sampled 1 h, 3 h, and 5 h post-treatment for RNA-Seq. The response involved, first and foremost, a higher number of differentially-expressed genes in A. aerophoba than D. avara. Secondly, while both species constitutively express a diverse repertoire of immune receptors, they differed in their expression profiles upon MAMP challenge. The response in D. avara was mediated by increased expression of two NLR genes, whereas the response in A. aerophoba involved SRCR and GPCR genes. From the set of annotated genes we infer that both species activated apoptosis in response to MAMPs while in A. aerophoba phagocytosis was additionally stimulated. Our study assessed for the first time the transcriptomic responses of sponges to MAMPs and revealed conserved and species-specific features of poriferan immunity as well as genes potentially relevant to animal-microbe interactions.


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 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 Invitrogen TM Qubit TM 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.35 46  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.5 53 , 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.
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 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 Eumetazoan 54 ) but with immunoglobulins instead of LRRs as extracellular domain 31 . 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 species 35,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 log 2 |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. 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 invertebrates 55 . Potential receptors according to sequence similarity, but without the corresponding conserved domains, are included in Tables 4, 5 and Supplementary Table S5. 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.  (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 pathway 56 . 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).

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 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 receptors 35 . Both A. aerophoba and D. avara constitutively express Immunoglobulin-TIR receptors, as found in other sponges 54 . 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 challenging 27 . Consequently, the set of Poriferan-unique and species-specific traits remain misrepresented 27,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 avoided 61 . 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 microbes 10,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 experiments 38 . 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 field 37 . 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 microbes 63,64 . Further studies have reported the up-regulation of these receptors upon bacterial exposure in other invertebrates 65 . 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 sponges 36 , 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 scarce 66 . 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 inflammasome 67 .
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 animals 25,69 , including sponges 68,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 TLRs 25,71 . In invertebrates, their role in defence has been suggested for Caenorhabditis elegans and Drosophila melanogaster 24 . 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 type 72 . Thus, our results provide additional support for the conserved role of GPCRs in animal-microbe interactions. In D. 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. 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 invertebrates 55,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 communication 75 , these domains are common in immune receptors 76 and are involved in microbial recognition in corals 77 , snails 78 , or nematodes 79 . Moreover, a ficolin-like gene was up-regulated in the sponge Cliona varians when "reinfected" with Symbiodinium compared to the aposymbiotic tissue 80 . 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 invertebrates 55,59,74,81 and potentially mediate a high diversity of cellular responses, such as cell death 81 , phagocytosis 82 , and metabolism regulation 55 . 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 species 83,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 microbes 38,85 . In Suberites domuncula, MyD88 expression was up-regulated 12 h after exposure to the same E.coli-derived LPS we used in our study 85 . However, before treatment, these sponges were kept in cultivation for a long period of time and their symbiotic bacterial load was reduced 85 , 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 . 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".  38 . 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 diversity 40,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 sponges 34 , while Germer et al. found NLRs to be absent from the transcriptome of the HMA sponge Vaceletia sp 35 . 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 considered 61,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 immunity 18 , 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 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". 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.

Data availability
Raw reads with the corresponding metadata and gene quantification matrices generated during the current study are available in the ArrayExpress database at EMBL-EBI archive (www.ebi.ac.uk/arrayexpress) under accession number E-MTAB-6757. De novo reference transcriptomes and their full annotation are available from the corresponding author upon request. Further processed data are included in this article and its Supplementary Information files.