In search of the Aplysia immunome: an in silico study
BMC Genomics volume 23, Article number: 543 (2022)
The immune repertoires of mollusks beyond commercially important organisms such as the pacific oyster Crassostrea gigas or vectors for human pathogens like the bloodfluke planorb Biomphalaria glabrata are understudied. Despite being an important model for neural aging and the role of inflammation in neuropathic pain, the immune repertoire of Aplysia californica is poorly understood. Recent discovery of a neurotropic nidovirus in Aplysia has highlighted the need for a better understanding of the Aplysia immunome. To address this gap in the literature, the Aplysia reference genome was mined using InterProScan and OrthoFinder for putative immune genes. The Aplysia genome encodes orthologs of all critical components of the classical Toll-like receptor (TLR) signaling pathway. The presence of many more TLRs and TLR associated adapters than known from vertebrates suggest yet uncharacterized, novel TLR associated signaling pathways. Aplysia also retains many nucleotide receptors and antiviral effectors known to play a key role in viral defense in vertebrates. However, the absence of key antiviral signaling adapters MAVS and STING in the Aplysia genome suggests divergence from vertebrates and bivalves in these pathways. The resulting immune gene set of this in silico study provides a basis for interpretation of future immune studies in this important model organism.
Making inferences of immune function from genomics and transcriptomics data can be challenging with non-mammalian or non-ecdysozoan models due to a paucity of functional research and gene characterization outside these two groups. Immune research in mollusks has historically centered on commercially important species such as the Pacific oyster Crassostrea gigas [1,2,3], or on models of parasite transmission such as the bloodfluke planorb Biomphalaria glabrata [4,5,6]. Recent work has ventured to describe immune function in other mollusks, mainly in other commercially important bivalves [7,8,9], but also in a few important molluscan models such as the pond snail Lymnaea stagnalis  and the common octopus Octopus vulgaris . These organisms often lack reference genomes or transcriptomes, and researchers must rely on de novo assembly and homology-based transcript annotation. Many other mollusk models for which publicly available reference sequences are available have not been comprehensively surveyed for immune related genes, such as the marine model Aplysia californica.
Despite being a premier neural model for more than half a century and boasting a reference genome and transcriptome, the immune repertoire of Aplysia is still poorly understood [12,13,14]. Historically, the extent of interest in the Aplysia immune response was a model for the role of inflammation and immune response in neuropathic pain [15,16,17,18]. Recently, however, a novel virus in the order Nidovirales was discovered in the marine model Aplysia californica . This virus appears to have a broad tropism, with particularly high viral expression levels in the nervous system . This discovery presents an opportunity to further develop Aplysia as a new model for neurotropic viral infection, particularly in the context of neural aging. However, the dearth of functional annotation and lack of even a putative immune gene set makes future inference into the immune response to this novel virus difficult. While some putative immune associated genes have been identified in Aplysia previously in the context of the neuronal transcriptome  or aging associated neuro-inflammation , these studies did not perform a comprehensive survey of the Aplysia immunome. To address this gap in the Aplysia literature, an in silico survey of Aplysia reference genome was conducted to identify and describe a putative immunome to aid in future immunological studies in this well-established model.
A combination of InterProScan (IPS) protein domain annotation and OrthoFinder gene homology search hierarchical clustering was used to identify genes in the Aplysia genome with potential for immune function (see Methods). Candidate immune genes were first identified by conserved protein domains typical of immune genes in other taxa and then vetted for ortholoy to said immune genes with OrthoFinder. This was particularly necessary for genes which cannot be easily identified by diagnostic conserved domains alone, such as inhibitors of NFkB. Of the 26,658 predicted proteins in the AplCal3.0 gene feature format file (gff) available from the NCBI RefSeq FTP site, 26,209 (98.3%) were annotated with at least one analyses incorporated into IPS, and 22,511 (84.4%) proteins were annotated with at least one InterPro accession.
From the 5 species used for OrthoFinder analysis (Human, Drosophila melanogaster, Crassostrea gigas, Biomphalaria glabrata, and Aplysia californica), proteins were sorted into 19,353 orthogroups. Of these orthogroups, 12,594 (65.1%) contain genes from at least 2 species while 4606 (23.8%) contained genes from across all 5 species surveyed. Comparing among groups, 1,723 orthogroups were common to mollusks and humans but not the ecdysozoan Drosophila. A further 1501 orthogroups were common only to the mollusks surveyed and 1918 were common only to the two gastropods. A total of 849 orthogroups were unique to Aplysia.
By mining the IPS results for immune domain signatures (see Methods) and vetting the resultant set with OrthoFinder results, 2241 Aplysia proteins from 1669 genes were identified as having potential immune function. These putative immune genes are described in detail below.
IPS annotation of the Aplysia californica RefSeq predicted protein models can be found in supplemental file SFile1. Full OrthFinder results can be found in supplemental file SFile2. The fully compiled immunome can be found in supplemental file SFile3.
Pattern recognition receptors (PRRs)
Several PRRs are adapted to detect the glycans of various pathogens, including peptidogylcans, beta 1,3-glucan from fungi, and lipoglycans like bacterial liposaccharide (LPS) [22,23,24]. Many such PRRs were identified in Aplysia, including: five proteins in the CD36-like family, three gram-negative bacteria binding-proteins (LOC101850763, LOC101851233, LOC101861522), and two peptidoglycan recognition proteins (LOC101855568, LOC101858995). Orthologs of Beta 1,3-glucan recognition proteins (BGRP) which are antimicrobial PRRs known from insects, were not identified, however. In vertebrates, detection of bacterial liposaccharides (LPS) is facilitated by an extracellular cascade of the LPS binding proteins LBP, CD14, and MD2. LPB binds LPS and facilitates its binding to CD14, which then delivers LPS to the MD-TLR4 receptor complex . A single LBP protein, LOC101845376, as well as 15 proteins with MD-2-related lipid-recognition domains were identified, however potential orthologs of CD14 were not (Table 1).
Lectins are a diverse group of multifunctional proteins defined by their carbohydrate recognition domains (CRD). Many lectins have been implicated in immune sensing of pathogen associated molecular patterns (PAMPs) across taxa, particularly C-type lectins . The Aplysia genome encodes a diverse collection of lectins including: 15 Apextrins, 146 C-type lectin domain containing proteins (CTLD), 33 Chi-type lectins, 5 F-type lectins, 96 fibrinogen domain containing proteins (FBGDC), 5 galectins, 5 H-type lectins, 100 i-set immunoglobulin proteins which may contain i-type lectins, 2 L-type lectins, 6 M-type lectins, 13 P-type lectins, 4 pentraxins, 21 R-type lectins, and 3 SUEL lectins were identified (Table 2). No mannose-binding lectins or f-box lectins were identified.
A unique subclass of lectins that pairs an immunoglobulin superfamily (IgSf) domain with a carbohydrate recognition domain, called Variable Immunoglobulin and Lectin domain containing molecules (VIgL), have been shown to play an important role in innate immunity in gastropods . The first identified VIgL in Biomphalaria was dubbed a fibrinogen related protein (FREP). Later, genes with similar domain architecture but galectin or C-type lectin carbohydrate recognitions domains were dubbed GREP and CREP respectively.
Based on conserved domains identified in IPS, 2 possible CREP genes (LOC101860578 and LOC101862063) were identified. It was not possible to find any FREPs based on conserved domains, although two Aplysia FREPs have been reported previously. Because IgSf domains of mollusk VIgLs are difficult to detect for automated software such as those used by InterProScan , BLAST homology searches to known gastropod VIgLs from Biomphalaria were used to identify these subsets of FBGDC, CTLDC, and galectin domain containing genes. This approach recovered the two previously identified FREPs  and aforementioned CREPs. Multiple sequence alignment of the two putative CREPs with previously described Aplysia FREPs via Clustal Omega identified potentially distantly related IgSf regions (Fig. 1) .
To detect infection by pathogens, animals employ a diverse array of nucleotide sensing PRRs to alert intracellular immune cascades to the presence of pathogenic DNA or RNA in the cytosol. The quintessential cytosolic PRR sensor is RIG-I, which detects viral RNAs . Three proteins with RIG-I-like c-terminal domains (CTD) involved in viral RNA recognition were identified. Two of these, LOC101847784 and LOC101853474, contain DExD/H-box helicase domains in addition to RIG-I-like CTDs, whereas LOC101853805 lacked any identifiable domain other than the RIG-I-like CTD. As in Biomphalaria, all three proteins lacked CARD domains typical of vertebrate and bivalve RLRs (Fig. 2). Two additional genes were identified as RLRs by OrthoFinder homology search (LOC101858480 and LOC101846979). These genes contain RIG-I/MAVS-type CARD (IPR031964) domains but lack any DExD/H-box helicase domains or RLR CTD. One of these, labeled as an IFIH1/MDA5 ortholog in the REFSeq database (LOC101858480), contains two CARD domains typical of vertebrate/bivalve RIG-I and MDA. Although RIG-I acts to detect viral RNAs, RNA Polymerase III can allow RIG-I to detect DNA viruses as well by converting AT-rich viral DNAs into RNAs detectable by RIG-I . The RPC34-like subunit of the RNA Polymerase III was identified: LOC101860996.
In addition to RLRs, several other DExD/H-box helicases have been implicated in immune sensing . Fifty-nine other DExD/H-box helicase genes were identified in Aplysia. Among these are specific DExD/H-box helicases known to play a role in the immune response of vertebrates, including DDX1 (LOC101857175), DDX41(LOC101864496), DHX9 (LOC101848799), and DHX15 (LOC101860724). Furthermore, the Aplysia genome contains 12 DEAD-box OB fold, also called domain of unknown function 1605, containing DExD/H-box helicases. This group contains other known vertebrate immune-associated genes including DHX29/30/33/36 and the aforementioned DHX15. Similarly, the Aplysia genome also contains 26 Q-motif containing DExD/H-box helicases, a group that contains still other known vertebrate immune-associated genes including DDX 3/6/17/19/23/24/39/46/56 in addition to aforementioned DDX1 and DDX41.
In parallel to RIG-I and DExD/H-box helicases, zinc finger NFX1-type containing 1 (ZNFX1), has also been implicated as a dsRNA sensor and antiviral protein in metazoans . Using conserved domains, PANTHER annotation from IPS (PTHR10887:SF470), and homology search, three potential ZNFX1 genes (LOC101862822, LOC101857249, LOC101851452) were identified.
Several DNA-damage detecting proteins have also been implicated in cytosolic sensing of viral nucleotide [32, 33]. Putative orthologs of several of these, including all three subunits of the DNAPK complex  (DNAPKc LOC101854982, Ku70 LOC101860888, and Ku80 LOC101855590); both components of the MRN complex [35, 36] (MRE11 LOC101858081 and RAD50 LOC101847050); and several Mab-21 domain containing proteins that may be orthologs of vertebrate cGAS (LOC101848176, LOC101852616, LOC101863483) were identified. In addition to DNA-damage associated sensors, orthologs to several other potential cytosolic viral sensors were identified, including LRRFIP1 (LOC101862386) and its downstream signaling partner β-catenin (LOC100533383); and 32 High Mobility Group Box proteins (Table 3).
Notably, no putative orthologs of DNA-damage sensors DAI/ZBP, Pyrin and hematopoietic interferon-inducible nuclear (HIN) domain (PYHIN) proteins , or viral RNA sensors inducible oligoadenylate synthetase (OAS) and Protein Kinase R [37, 38] were identified (Fig. 3).
Broad Spectrum PRR families
While several families of PRRs are dedicated to a specific class of PAMPs, the Toll-like receptor family of PRRs can detect a diverse range of PAMPs. For example, in mammals TLRs can detect bacterial lipopeptides (TLR1,2,3), glycans (TLR4), nucleic acids (3,7,8,9), and bacterial proteins (TLR2,5) . In Aplysia 39 genes with Toll/interleukin-1 receptor homology (TIR) domain and Leucine Rich Repeat (LRR) domains characteristic of TLRs were identified. OrthoFinder divided these putative TLRs in 17 orthogroups. Of these, four orthogroups were unique to Aplysia (OG0001270, OG0009332, OG0016038, and OG0009441). Of the 17, only three were not unique to gastropods, with OG0011109 and OG0012371 containing a protein from C. gigas, and OG0003873 containing proteins from C. gigas and Drosophila.
Neighbor joining tree with 1000 boostraps of TIR domain region of these putative TLRs broadly supports OrthoFinder’s division of TLR groups (Fig. 4). Only TIR domains from LOC101851910 (XP_005112480.1) and LOC101863390 (XP_012942443.2) clustered away from their assigned orthogroups (OG0001276 and OG0007258 respectively). However, these divergent TIR assignments corresponded to weekly supported branches with low bootstrap values.
Interestingly, LOC101863390 contains 2 TIR domains. In addition, LOC101863390 also contains 3 transmembrane domains, and two sets of LRRs that give the appearance of an end-to-end fusion of two TLRs.
Phylogenetic clustering of TLR sequences from Aplysia californica, Biomphalaria glabrata, Crassosteraea gigas, Nematostella, Drosophila melanogaster, Strongylocentrotus purpuratus, and humans revealed lineage specific trends in TLR diversification (Fig. 5). Molluscan TLRs clustered together with other protostome TLRs and exhibited a unique molluscan radiation within this group. This protostome radiation clustered apart from a massive radiation in S. purpuratus. The number of publicly available reference TLR protein sequences for C. gigas was substantially less than the manually annotated number previously reported by Zhang et al. (18 vs. 83) , making direct comparison difficult. However, these results do suggest a unique molluscan TLR radiation. Furthermore, TLRs from then two gastropods primarily originate from this uniquely molluscan radiation or its sister branch, unlike available C. gigas TLRs.
While most of the Aplysia TLRs contain only TIR and LRR domains, several contain additional domains. Two of these putative TLRs (LOC101850809) and LOC101845154) contain detectable Cysteine-rich c-terminal flanking region in proximity of their LRR (Fig. 4). Another two putative TLRs contain the ROC and COR domains typical of ROCO GTPases. One (LOC101848126) contains an accessory Mbt repeat region while the other (LOC101845043) contains an EF-hand domain (Fig. 6). While these two genes are assigned to different orthogorups, tree analysis demonstrates they form a clade. Only two TLRs were identified as having c-terminal flanking regions (LOC101850809 and LOC101845154).
Like TLRs, Nucleotide-binding oligomerization domain-like receptors (NOD-like receptors, NLRs) are a diverse family of PRRs that detect a similarly broad range of PAMPs in vertebrates . In Aplysia, only 4 proteins with the characteristic NACHT domain of NLRs were identified, only one of which contained both a NACHT domain and LRR domains typical of NLRs. This NLR-like protein also contains a Death-like domain and a DZIP3-like HEPN domain.
Like TLRs and NLRs, Receptor Cysteine-Rich (SRCR) Domains have been implicated in innate immune function across taxa against a broad suite of pathogens . SRCR PRRs have been described in bivalves where they are found on the surface of haemocytes and respond to PAMPS . In Aplysia, 11 SRCR genes were identified (Table 4).
PRRs in several signaling cascades associate with adaptor proteins to facilitate transduction of pattern recognition into activation of kinases and ultimately transcriptional changes to facilitate an immune response, most notably the TLRs. Adaptor proteins of TLRs complex with their target receptor via their TIR domains. The most conserved TLR adaptor among animals is Myeloid differentiation primary response 88 (MyD88), which complexes with TLRs via its TIR domain and with Interleukin-1 receptor-associated kinases (IRAKs) via its Death-like domain (DLD) in a complex called a Myddosome . In Aplysia 5 MyD88-like proteins were identified. One of these contained only the TIR and DLD typical of vertebrate MyD88 (LOC101855709), while three contained an additional Armadillo-type fold (LOC101850562, LOC101850562, and LOC101858605), and the last contained the same ROCO domain as described earlier in some AcTLRs (LOC101861518). OrthoFinder homology search also identified LOC101847182 as a MyD88-like protein based on TIR domain homology, but this gene lacked any death-like domains. In addition to MyD88, 4 other TIR domain containing (TIRDC) adaptors have been identified: MyD88-adaptor-like (Mal), TICAM1 and 2 (also known as TRIF and TRAM), and sterile-alpha and Armadillo motif protein (SARM) . Although orthologs of MAL, TRIF, or TRAM could not be identified, three Aplysia SARMs (LOC101849069, LOC101850713, and LOC106011731) were identified. In addition, a diverse group of TIRDC proteins with as of yet unknown function were also identified. Of these, 10 contained Armadillo-type folds, two contained EGF-like domains, one contained COR domain, and the remaining 17 contained only a TIR domain. Potential orthologs of other, non-TIRDC TLR adaptor proteins Toll-Interacting Protein (TOLLIP, LOC101845996) and evolutionarily conserved signaling intermediate in Toll pathway (ECSIT, LOC101857763) were also identified.
One of the several downstream kinase complexes associated with PRR activation is the TAK1-TAB complex. The TAB proteins TAB1-3 facilitate the association of the kinase TAK1 with ubiquitin chain scaffolds to colocalize with downstream targets . A TAB1 ortholog was difficult to identify based on protein domain alone, as the only InterPro domain in vertebrate TAB1 is the PPM-type phosphatase domain (IPR001932) which is found in 14 Aplysia genes. However, OrthoFinder homology search identified a gene with two isoforms (LOC101854409) listed as TAB1 on the NCBI. The other TAK1 adaptors TAB2 and 3 were similarly difficult to identify based on IPS alone. It was not possible to identify vertebrate-like TAB2/3 proteins that contain both a CUE domain and RAnBP2 zinc finger. OrthoFinder homology search identified two isoforms of the same gene (LOC101854409), which is identified as TAB2 on the NCBI. Notably, this gene lacks any CUE domain but does contain a RAnBP2 zinc finger needed for ubiquitin chain binding and subsequent signaling . In fact, 19 RanBP2-type zinc finger domain containing genes were identified (Table 5).
Notably, an ortholog of the downstream adaptor of RLRs, mitochondrial antiviral signaling protein (MAVS) could not be identified (Fig. 7). Similarly, an ortholog of Stimulator of Interferon Genes (STING), could not be identified (Fig. 7). STING acts an adaptor downstream of many DNA-damage associated PRRs including cGAS, DAI, PYHIN proteins, DNAPK, and the MIRE11-RAD50 complex .
E3 ubiquitin ligases
Ubiquitylation of various adaptor proteins and kinase complexes is a critical component of PRR signaling cascades downstream of many PRRs, notably TLRs, NLRs, and RLRs .
Tumor necrosis factor receptor (TNFR)–associated factor (TRAF) proteins play a key role in activating kinases during PRR signaling events . In Aplysia, 6 TRAF-like proteins were identified: two resembling TRAF2 (LOC101856795 and LOC101844983), one resembling TRAF3 (LOC101846107), one resembling TRAF4 (LOC101859360, and two resembling TRAF6 (LOC101856862 and LOC101859878).
The LUBAC complex, comprised of proteins HOIP, HOIL1, and Sharpin, also plays a key role in RLR and TLR signaling, particularly in the activation of the IkB complex in NFkB signaling . In Aplysia, putative orthologs of LUBAC complex components HOIP (LOC101849805) and HOIL1 (LOC101863809) were identified. Although an individual gene representing the third LUBAC component, Sharpin, was not identified, one of the 4 AcHOIL1 isoforms (XP_035825204.1) resembled Sharpin more then HOIL1 in conserved domain composition (i.e. lacking a TRIAD domain) and thus may represent a functional analog of Sharpin.
Another immune associated E3 ligase complex is comprised of SKIP1, cullin-1, and f-box proteins, or SCF . In Aplysia, an ortholog of SKIP1 (LOC101854298), 6 putative cullins, and 60 f-box proteins, 14 of which contained WD40 domains (FBXW) and 29 of which contained leucine rich repeats (FBXL) were identified.
Other E3 ligases function to modulate or inhibit pro-inflammatory immune signaling. Pellino proteins modulates the activity of several immune associated kinases such as RIPKs and IRAKs . In Aplysia only one Pellino gene (LOC101863733) was identified. Similarly, the diverse family of tripartite-motif contain proteins (TRIMs) regulate immune and proinflammatory signaling cascade via Ubiquitylation of various kinases . TRIMs generally contain a b-box domain and a RING type zinc finger domain . Twenty-two such genes in Aplysia were identified.
While the baculoviral Inhibitor of Apoptosis (IAP) repeat (BIR) family of proteins are known to modulate apoptosis, BIRC2/3, also known as cIAP1/2, are also known to modulate immune associated NFkB signaling . The Aplysia genome encodes 22 BIR domain containing proteins, 11 of which also contained RING zinc finger domains that represent potential BIRC E3 ligases. Three of these also contained Death-like domains typical of cIAPs (LOC101851721, LOC101860034, and LOC101852760).
Signaling downstream of RLRs also depends on ubiquitylation via TRIMs and other E3 ligases such as MEX3C and RIPLET . In Aplysia, a single potential MEX3C ortholog (LOC101854000) and two RIPLET-like genes (LOC101858054 and LOC101862807) were identified (Table 6).
Interleukin-1 Receptor Associated Kinases (IRAKs) are serine/threonine kinases that complex with MyD88 via their Death-like domains and phosphorylate downstream adaptors in TIR signaling . Two IRAK genes were identified, one (LOC101845636) is listed as IRAK4 on the NCBI, and the other (LOC101845840) assigned to the Pelle-family of IRAKs by PANTHER (PTHR24419:SF33).
In vertebrates, IRAK activation leads to activation of TRAF6, which ubiquitylates and activates the TAK1-TAB complex which contains the mitogen activate protein kinase (MAPK) kinase kinase 7 (MAP3K7/ TAK1) . In Aplysia one ortholog of TAK1 which has 10 protein isoforms (LOC101852077) was identified. TAK1 phosphorylates two subsequent kinase chains: the IKK complex and mitogen activated protein kinases (MAPKs), that lead to translocation of transcription factors NFkB and AP-1 to the nucleus .
In canonical NFkB signaling, the IKK complex comprises IKKa/CHUK, IKKb, and IKKg/NEMO, and phosphorylates NFkB inhibitor proteins (NFkBIA, Cactus, etc.) to allow translocation of NFkB to the nucleus . The non-canonical IKKs, IKKe and TBK1, complex together and subsequently with adaptors like SINTBAD or TANK to activate Interferon regulatory factors (IRFs) . In Aplysia two canonical IKKS (LOC101856649, LOC106012158) and adaptor NEMO (LOC101858005) were identified, alongside two non-canonical IKKs (LOC101853683 and LOC101847449) that likely represent orthologs of TBK1 and IKKe. The only non-canonical IKK adaptor identified was a putative ortholog of NAP1 (LOC101852476) . Orthologs for adapter proteins of another key immune associate kinase, TBK1, TBKBP1/SINTBAD and TANK, were not identified.
Parallel to IKKs, TAK1 and other MAP3Ks activate MAP2Ks (MKKs or MEKs), which activate MAPKs which in turn activate transcription factors, notably those in the AP-1 family . In Aplysia several putative orthologs of several immune associated kinases in this three-tiered cascade were identified.
In addition to TAK1, one MAP3K in the Apoptosis signal-regulating kinase (ASK) family (LOC101862124) was identified, a MAP3K10 (LOC101864517), and a MAP3K12/13 (LOC101864115). Based on InterProScan annotation, it was not possible to identify orthologs of immune relevant MAP3Ks MAP3K14/NIK or MAP3K8/Tpl2. Annotation by PANTHER identified an additional 17 MAP3Ks.
MAP2Ks were difficult to identify purely using IPS domain annotation. Using the CDD analysis generated by IPS, a MAPK/ERK Kinase (MEK) with two isoforms (LOC101850764) was identified. Annotation by PANTHER further identified an additional 4 MAP2Ks (LOC101845548, LOC101846645, LOC101846695, LOC101860284).
Based on the MAPK-kinase conserved site (IPR003527), 6 MAPKs were identified, including orthologs of stress-activated kinases p38/MAPK14 (LOC100533304) and JNK, and an Extracellularly Regulated Kinase (ERK, LOC100533212). An additional p38 (LOC101857669) which lacks the conserved site but is identified as in the p38 family (IPR008352). PANTHER annotation identified 7 additional putative MAPKs that lacked the conserved site.
MAP kinase-activated protein kinase 2 (MK)-like gene (LOC101861689) and four MAPK phosphatases (MKP) that regulate MAPKs also were identified (Table 7).
Notably, it was not possible to identify any orthologs of receptor-interacting serine/threonine protein kinases (RIPKs) which also play an important role in immune and inflammatory signaling .
The ultimate action of PRR signal transduction cascades is the activation and translocation to the nucleus of pro-survival/pro-inflammatory transcription factors NFkB and AP-1 and interferon (IFN) regulatory factors (IRFs) .
Nuclear factor kappa beta (NFkB) transcription factors are the quintessential and evolutionarily conserved immune/inflammation transcription factors . Based on Rel homology domains (IPR032397 and IPR011539) a p105-like protein with two isoforms (LOC101855313) and a RelB/Dorsal-like (LOC101860399) protein were identified. A putative ortholog of Nuclear factor of activated T cells (NFAT, LOC101858233), which are a related class of immune transcription factors, was also identified. NFkB proteins are characteristically bound by inhibitor proteins (NFKBI/IkB) which are the target of IKKs. Aplysia contains putative orthologs of NFKBIA with 4 isoforms (LOC101852802) and of NFKBI Cactus (LOC101863680).
Just as NFkB is the target of the IKK branch of immune signaling, activator protein 1 (AP-1) is the downstream target of the MAPK cascade. AP-1 is a dimeric transcription factor made up of any two of a related group of proteins families known as JUN, FOS, MAF, and ATF . In Aplysia, 8 AP-1 type transcription factors were identified, including one identified as ATF2-like (LOC101856191). Homology search further identified the Aplysia c-Jun gene (LOC101863323).
Parallel to NFkB and AP-1, IRFs are activated by immune signaling such as TLRs or cytosolic nucleic acid sensors and facilitate the translation of interferon cytokines and interferon stimulated genes (ISG) to counteract infection . Seven IRFs were identified: 2 IRF1/2-like that lack SMAD domains (LOC101855936 and LOC101850966) and 5 SMAD domain containing IRFs. Translocation of these transcription factors is facilitated by major vault proteins (MVP) , 4 of which were identified in Aplysia (Table 8).
Successful activation of NFkB, AP-1, and IRFs transcription factors leads to the production of cytokines that recruit another family of critical immune transcription factor: Signal Transducer and Activator of Transcription (STAT) proteins. Cytokine receptors for Interleukins and Interferons complex with Janus Kinases (JAKs) which phosphorylate STAT transcription factors to promote transcription of interferon stimulated genes (ISGs) . This pathway is then autoinhibited by Suppressor of Cytokine Signaling (SOCS) proteins. In Aplysia, a single JAK (LOC101861981), a single STAT (LOC101861517), and 5 SOCS (LOC101845832) were identified.
Another major transcription factor that plays a role in immune response is the pro-proliferative Myelocytomatosis oncogene (MYC). While MYC is generally known for its role in modulation of the cell cycle and cellular proliferation, it has also been demonstrated to play a role in not only immune cell proliferation but also immunomodulation . Conserved domain search failed to identify an Aplysia MYC based on conserved leucine zipper (IPR003327) and n-terminal domains (IPR002418). Homology search also failed to identify an Aplysia ortholog of MYC. However, 52 genes were identified with a myc-type basic helix-loop-helix (bHLH) domain typical of transcription factors related to MYC, leaving the possibility a MYC-like transcription factor may yet be identified.
Successful activation of immune receptor transduction results in the secretion of pro-inflammatory cytokines . Activation of pro-inflammatory cytokine receptors and associated pathways then stimulates the proliferation of immune cells in an inflammatory immune response . Seven putative orthologs of the pro-inflammatory cytokine Macrophage migration inhibitory factor (MIF), a progranulin ortholog (LOC101846478), 4 putative tumor necrosis factors (TNF), and 4 interleukin-17 -like cytokines (IL-17) was identified. Three putative IL-17 receptors and 15 TNF receptor-like genes also were identified (Table 9).
TNF receptors associate with adaptor proteins via Death domains (DD) to transduce signaling downstream to the pro-inflammatory NFkB signaling pathway or to the pro-apoptotic caspase pathway. These DD containing signaling adaptors are TRADD, CRADD/RAIDD, and FADD. Putative CRADD (LOC101855728) and FADD (LOC101861195) orthologs were identified, but not TRADD.
Although not adapter proteins per se, A disintegrin and metalloproteinase (ADAM) proteins have been implicated in immune function, notably in cytokine signaling, particularly of TNF and Interleukins . Twenty-five ADAMs were identified based on conserved peptidase domain (IPR001590).
Apoptosis is the programmed self-destruction of cells and plays key roles in diverse biological processes such as development and the immune response. In mollusks, apoptosis is believed to play a critical role in the immune response by destroying pathogen infected cells and limiting the spread of infections without inducing inflammation . Indeed, apoptosis serves as a key mechanism for slowing the spread of viral infections across taxa . In addition to BIR containing apoptosis inhibitors mentioned in previous sections, several genes in the Bcl-2 family of proteins also act as important modulators of apoptosis . In Aplysia,6 Bcl-2 family proteins, including putative ortholog of anti-apoptotic Bcl-2 (LOC101848156), as well as pro-apoptotic proteins Bax (LOC101854270), Bak (LOC101851280), and Bok (LOC101854628) were identified. Two potential Fas Apoptosis Inhibitor Molecules (FAIM, LOC101863121 and LOC106012994), which act as apoptosis inhibitors and play a role in IFN signaling  were identified. Apoptosis signaling is most notably mediated through pro-inflammatory and pro-apoptotic proteases called caspases . Eight caspase genes were identified. Other apoptotic regulators Death factor and mitochondrial metabolism regulator AIF (LOC100533557)  and IAP inhibitor and pro-apoptotic factor HTRA2  was also identified (Table 10).
In mollusks, often the first line of defense against pathogens is the mucus, which contains a diverse array of antimicrobial compounds including lectins and mucins [75, 82]. In addition to previously discussed lectins, the Aplysia genome also encodes 13 mucins. Among the most ancient antimicrobial proteins are those that punch holes in the cellular membranes of pathogens, called pore forming proteins (PFP) [83, 84]. In Aplysia, 14 Aerolysin type, 5 ETX/MTX-2 type, and three MACPF containing MPEG-like pore forming proteins were identified. Lysozymes are similarly ubiquitous components of the innate immune system that act to lyse pathogen membranes . In Aplysia three lysozymes were identified (Table 11). Other antimicrobial proteins known from bivalves such as defensins, mytilins, or mytimicins were not identified.
In the case of virus infection specifically, cells mobilize a diversity of antiviral proteins that hamper virus entry, uncoating, replication, translation, assembly, and release . Two dynamin family Mx-like proteins (LOC101848065 and LOC101863954), one RNA Specific Adenosine Deaminase (ADAR) with three isoforms (LOC101855902), four viperins/RSAD, 8 caveolins, 5 Gilt-like proteins, three APOBEC3-like proteins, and 10 Interferon-induced transmembrane protein/Dispanin family (IPR007593) proteins were identified. OrthoFinder homology search identified one further potential ADAR gene (LOC101863647), which contains dsRNA binding domains (IPR014720) and adenosine deaminase/editase domains (IPR002466) but lacks a Z-binding domain (IPR042371). Orhtologs to other antiviral proteins including tetherins, SAMHD, Zinc finger antiviral protein (ZAP), IFI41, or RNAse-L orthologs were not identified.
Similarly to DNA damage response genes being implicated as PRRS, poly(ADP-ribose) polymerases (PARPs) which are associated with DNA damage response have also been demonstrated to have broad spectrum anti-viral and pro-inflammatory effects [86, 87]. In Aplysia, 17 PARP genes were identified (Table 12).
RNA interference (RNAi) is another potent antiviral mechanism of particular importance in arthropods . Key components of the RISC RNAi complex, including Dicer (LOC101848232), two argonaute proteins (LOC101846762 and LOC101863465), and two additional PIWI domain containing proteins that lack argonaute linker domains were identified (Table 13). An ortholog of RISC-loading complex subunit TARBP2 could not be identified via InterPro conserved domains alone. OrthoFinder homolog search, however, identified LOC101846211, named TARBP2 on the NCBI, as a potential TARBP2. LOC101846211 contains several dsRNA binding domains (IPR014720) typical of TARBP2, as well as a Staufen c-terminal (IPR032478).
ROS and Antioxidants
A prime function of immune cells is the capacity to recognize, phagocytose, and destroy pathogens during an immune response. Phagocytosed pathogens are destroyed by the generation of reactive oxygen species (ROS). ROS also play a role in potentiating pro-inflammatory signaling cascades during immune signaling . In Aplysia a diverse set of ROS generating enzymes were identified, including: 3 dual oxidases (DUOX) and 2 associated assembly factors (DUOXA), 7 nitric oxide synthases (NOS), 5 NADPH oxidases, and one peroxidasin (LOC100533347). Furthermore, 20 L-amino-acid oxidases (LAAO) were identified, among which was antimicrobial peptide Aplysianin-A (LOC100533295) .
Due to the cytotoxic nature of ROS, levels of reactive species are tightly regulated in the cell via antioxidant enzymes. Additionally, antioxidant enzymes have also been demonstrated to act as immunomodulators of proinflammatory cytokines . In Aplysia, a robust complement of antioxidant enzymes was identified, including: 8 Cu-Zn superoxide dismutases (SOD), one Fe-Mn superoxide dismuates (MnSod, LOC101852344), two catalases (LOC101854685 and LOC101860167), and 14 peroxiredoxins. In the glutathione antioxidant pathway, 3 glutathione peroxidases, 8 glutaredoxins, a glutathione reductase (LOC101847425), and 68 glutathione S-transferases were identified (Table 14).
Heat shock proteins
Stress signaling during an immune response is known to induce the transcription of heat shock protein (HSP), which also have been suggested to have immunomodulatory function . In Aplysia, many HSPs across several families were identified (Table 15), including: 20 HSP20s, 12 HSP70s, 4 HSP90s, and one heat shock factor (HSF, LOC101851429).
In arthropods, part of the immune response involves the production of toxic quinone compounds, wound healing factors, and pathogen encapsulating melanin in a process known as melanization. Although most well described in arthropods, melanization as an immune response has also been described in several bivalves  and in the gastropod B. glabrata . While 5 laccases and 6 tyrosinases that participate in melanization were identified in Aplysia, an ortholog of the critical Phenol Oxidase central to arthropod melanization could not be identified.
A major component of the vertebrate immune response is the complement system, which facilitates the formation of the pathogen destroying membrane attack complex (MAC) and the recruitment of phagocytic cells . The complement system is triggered after PAMP/DAMP binding to one of three lectins: Complement 1 (C1q), MBL, or fibrinogen and collagen domain containing lectins (ficolins) . While neither MBL nor ficolins were identified in Aplysia, 11 C1q domain containing genes (C1qDC) were identified. In addition, two orthologs for C1q receptor Multiple epidermal growth factor-like domains protein 10 (MEGF10, LOC101848617 and LOC106013429) were identified. Activation of complement PRRs recruits a group of thioester containing proteins called complement factors.
Thioester containing proteins (TEPs) play an important role in pathogen defense across the animal kingdom. TEPs are subdivided into complement factors in the complement system and alpha 2 Macroglobulins (a2Ms) which act as protease inhibitors . In Aplysia 13 TEPs based on conserved thioester domains (IPR011626) were identified. Two of these, LOC101848728 and LOC101862906, could be classified as complement factors C3 and C4-A respectively based on conserved domains. LOC101854858 and LOC101861335 were further classified as CD109-like proteins based on conserved domains. On the other hand, a factor B ortholog (Bf) based on conserved domains could not be identified. However, unlike vertebrate Bf which contains type A von Willebrand factor, Sushi, and Trypsin domains, the C. gigas gene identified as Bf by Zhang et al. 2015 contains von Willebrand type A (IPR002035), Sushi (IPR000436), EGF (IPR000742), and Pentraxin domains (IPR001759) . Based on these domains, an Aplysia Bf (LOC101854327) was identified from among the Aplysia genes identified as Pentraxins (Table 16).
In this in silico analysis, a robust immune complement in the Aplysia genome was proposed. This includes PRRs and adaptor proteins, signal transducers, transcription factors, and subsequent cytokine signaling components [97,98,99]. Among these are NFkB transcription factors, the IKK complex, a diverse compliment of MAPKs and associated kinases, the TAB-TAK1 complex, E3 Ubiquitin ligases like Pellino, TRAFs, and the LUBAC complex, and a diverse collection of TIR domain adaptors and receptors.
Functional studies in C. gigas have demonstrated that the oyster TLRs function as immune sensors and signal through MyD88 down to NfKB as in vertebrates . Similarly. Studies have shown that the oyster IKK complex does indeed form from two IKKs and adapter NEMO as in mammals. Furthermore, oyster IKK interacts with TRAF6 and MyD88 and induces the transcription of NF-kB and IRF target genes, likely through NFkB and IRF8 orthologs . Indeed, mollusks in general appear to retain a TLR signaling cascade like that of vertebrates . As in other mollusks, the Aplysia genome encodes essentially all the major components of the canonical TLR signaling cascade (Figs. 8). Given the broad conservation of genes from the TLR cascade in mollusks, and the conserved function of many of those genes as described in oyster, Aplysia possibly also retain a vertebrate-like TLR signaling cascade (Fig. 9). However, molluscan TLRs appear to have undergone massive expansion compared to vertebrates, suggesting a more diverse TLR signaling pool beyond the canonical pathway conserved with vertebrates.
Previously, a survey of the C. gigas genome demonstrated a massive expansion of the short protostome type TLRs and a smaller expansion of TLRs containing leucine-rich repeat c-terminals. This expansion led to a comparatively high TLR count, though was different in type compared to S. purpuratus which exhibited an expansion but in vertebrate-type TLRs. These TLRs were upregulated in response to boiotic challenge, suggesting a unique molluscan immune TLR cascade. Compared to C. gigas and S. purpuratus, Aplysia is comparatively poor in TLRs (39 in Aplysia vs 83 in C. gigas and 222 in S. purpuratus). Nevertheless, like other invertebrates, the Aplysia genome encodes many more TLRs than in vertebrates. While direct comparison to previously reported C. gigas numbers was not possible due to differences in custom vs reference protein number for C. gigas, the neighbor joining tree generated in this study corroborates earlier results that demonstrate that the massice expansion of TLRs in S. purpuratus and in mollusks represent distinct expansion of different TLR families . Interestingly, unlike the C. gigas TLRs, the gastropod TLRs predominantly clustered within the proteosome specific branch of the tree, with the majority clustering in the mollusk specific branch. These results agree with OrthoFinder results which suggest most of the gastropod TLRs form a unique group apart from arthropods, veterbrates, and even bivalves. This may suggest a unique gastropod TLR expansion that may hint at unique gastropod TLR signaling pathways.
Concomitantly, The Aplysia genome further contains a diversity of TIRDC that may represent novel or expanded TLR signaling adaptors. Among these are several MyD88 and SARM-like TIR adaptor proteins, both of which are known from diverse lineages and are considered the most ancient of the TLR adapters [44, 102]. However, many contain domains unknown from vertebrate TIRDC adapter. The presence of Armadillo-like folds in several Aplysia TLRs and TIRDCs may represent a novel TLR signaling cascade. Similarly, the presence of ROCO domains in two AcTLRs, two AcMyD88s, and an AcTIRDC may represent another novel GTPAse TLR signaling cascade. The inability to identify orthologs of TIR adapters TICAM1/TRIF and TICAM2/TRAM was expected as this MyD88-independent TLR signaling cascade arose only in the chordate lineage . This suggests that the TRIF/TRAM/TBK1 signaling axis of some vertebrate TLRs is absent in Aplysia, or that another set of TIRDC adaptors fill an analogous role of TICAMs in non-canonical NFkB signaling. Downstream of TIRDCs, Aplysia retains two IRAKs as in other invertebrates, like the arthropod Drosophila and basal chordate Amphioxus. Given that these two IRAK proteins transduce TIR activation in different ways between taxa, it is unclear how Aplysia IRAKs interact with Aplysia TRAFs or Cactus orthologs .
Unlike in arthropods, Aplysia also appears to retain several key components of the complement system, including C1q lectins, many thioester proteins including factors C3 and C4, and factor B. While the complement system has substantially changed in the vertebrate lineage since the emergence of jawless fish, evidence suggests that the basic components of this pathway have ancient roots in early metazoa . In this context, the Aplysia complement system likely serves to opsonize pathogens and induce inflammation . In oyster, the C1q domain containing family underwent an extensive expansion like the TLRs . In comparison, Aplysia is comparatively poor in C1qDC proteins, both compared to oyster and human.
Of particular interest for future studies of the Aplysia Abyssovirus, the Aplysia genome encodes many anti-viral genes. Among these are 5 Gamma-interferon-inducible lysosomal thiol reductases (Gilt/IFI30), which are key ISGs that inhibit viral entry in mammals . Crustacean orthologs of Gilt have been demonstrated to have conserved viral entry-inhibiting function . Although reported as an Mx ortholog, DQ821497 of Haliotis discus appears to be an ortholog of Aplysia Gilt based on homology and conserved domains. This abalone Gilt has been demonstrated to be induced by viral RNA mimic Poly I:C, suggesting conserved antiviral function of gastropod Gilt orthologs as well . Vertebrate Mx proteins are unrelated to Gilt proteins and are instead related to dynamins and serve to prevent entry and replication of a diverse array of viruses . Several dynamin-like genes were identified in Aplysia, but tree-based result from OrthoFinder suggests these Aplysia dynamins are more like vertebrate dynamins than vertebrate Mx proteins, and so may not exhibit any antiviral activity.
The Aplysia genome also encodes orthologs of antiviral deaminases ADAR and APOBEC3. These ISGs inhibit viral replication and function by introducing point mutations into viral genomes [109, 110]. RNA editing of viral genomes by molluscan ADAR genes has been demonstrated to play an important role in the immune response to malacoherpseviruses . The Aplysia genome also encodes a large complement of PARP proteins, which have been demonstrated to have strong anti-viral action against a diverse suite of viruses in mammals . If these defenses are insufficient to protect cells from active infection, the Aplysia genome also encodes several viperins to prevent viral spread. Viperins are ISGs that serve to inhibit viral budding and release from lipid bodies in diverse taxa . Conservation of function for molluscan viperins has been demonstrated in oysters . Although the Aplysia genome encodes an immune compliment with substantial conservation with known innate immune genes, several major immune pathways or components appear to be missing.
Despite containing two non-canonical IKKs, the Aplysia genome lacks many of the components associated with TBK1 signaling in vertebrates. The adaptor proteins STING and MAVS, neither of which have identifiable orthologs in Aplysia, are both known to act as scaffolds that facilitate TBK1 activation of IRF3 in response to antiviral PRR activation .
While present in bivalves, the absence of STING has been described previously in heterobranch gastropods like Aplysia . STING acts as an adaptor that facilitates phosphorylation of IRF3 by TBK1 in a diverse suit of nucleotide sensing pathways, many of which are either absent in Aplysia (DAI, PYHIN proteins), or non-productive in the heterobranch lineage (cGAS). The Aplysia genome does however encode for the MRN and DNAPK complexes, which have also been described as able to activate TBK1 in a STING dependent and independent fashions. Aplysia also lack the upstream components of TLR associated TBK1 signaling; notably TICAM proteins and TBK1 adaptors TANK and SINTBAD that form an analogous signaling cascade to the MyD88-TAB-TAK1 axis. TRADD, which is also known to complex with TRAF2/3 and activate the TANK-TBK1-IRF axis is also absent. Interestingly, Aplysia does retain an ortholog of RIOK3, which has been described as an essential adaptor of TBK1 and as an inhibiter of RLR signaling [62, 63, 116]. Given the diversity of yet undescribed TIRDC proteins in Aplysia, perhaps a convergent non-canonical IKK signaling pathway has emerged in the heterobranch lineage.
Given the absence of STING in heterobranchs, despite its presence in the bivalve lineage, it is perhaps not surprising that Aplysia also lack a MAVS ortholog [117, 118]. Nevertheless, the Aplysia genome does still include 3 putative RLRs that act upstream of MAVS in vertebrates and presumably bivalves. Notably, the Aplysia RLRs lack any CARD or DLD which act to form complexes between RLRs and MAVS in canonical RLR signaling . This condition of Aplysia RLRs is conserved with those in Biomphalaria, suggesting a functional divergence of heterobranch RLRs from those in bivalves. Perhaps heterobranch RLRs follow a similar pattern to the highly divergent immunoglobulin domains of VIgLs, where the functional domain is conserved but so highly divergent that it is difficult for automated systems to detect. Perhaps a similarly highly divergent MAVS-like protein has yet to be identified in heterobranchs, or a novel signaling axis for RLRs has emerged in the absence of MAVS. Future studies coprecipitating Aplysia RLRs with interacting proteins may be able to uncover such a highly divergent MAVS, if it indeed exists.
Another notably absent immune axis is that of OAS1 and its downstream effector RNAse L. OAS acts as a sensor for a diverse repertoire of viruses, including RNA viruses such as SARS-CoV2 .
In addition to missing some PRRs and associated signaling cascades, the absence of certain pro-inflammatory pathways may suggest Aplysia exhibits an inflammatory response different from vertebrates. A major component of the vertebrate inflammatory response is the formation of the NLRP3 inflammasome [121, 122]. While originally believed to have evolved in vertebrates, recent evidence suggests that the NLRP3 inflammasome may have antecedents in invertebrate innate immunity as orthologs of NLRP3 exhibit immune function in echinoderms and arthropods [123, 124]. The only Aplysia NLR identified does contain a Death-like domain, considered a possible predecessor to the CARD domain of NLRP3 in vertebrates, suggesting perhaps some conserved function . However, the dearth of NLRs, lack of PYHIN proteins, and absence of RIPKs suggest that Aplysia does not generate inflammasomes in response to pro-inflammatory signaling .
Studies in bivalves suggest an ancestral molluscan innate immune system with a high degree of conservation with that of deuterostomes [3, 8, 127,128,129]. However, the absence of several PRRs and their adaptors suggests that the immune response in Aplysia may be substantially different not only from vertebrates but also from bivalves. Antiviral immunity in arthropods acts through different mechanisms than in vertebrates, notably lacking RLR antiviral signaling. Arthropods detect and respond to viruses via the RNA interference, TLR, JAK/STAT, and IMD/PGRP pathways [88, 130]. Perhaps heterobranchs gastropods, exhibiting a reduced immune complement compared to their bivalve relatives, have developed alternative, lineage specific immune pathways as occurred in arthropods.
Together these results illustrate the difficulty of interpreting transcriptomic and genomics data in understudied lineages. Due to the overrepresentation of mammals and ecdysozoans in immune research, there is risk of erroneous interpretation of gene function based on similarity searches alone in mollusks. In these results, areas of broad conservation, notably in canonical toll-like receptor signaling were identified. However, Aplysia also exhibits areas of stark deviation from mammalian and even bivalve immunity, notably the absence of STING and MAVS. This illustrates the need for a healthy measure of caution when investigating immunity in Aplysia. By combining homology searches with annotation of conserved protein functional domains, an annotation set has been created that will provide a more nuanced basis for interpretation of immune-associated data in Aplysia. However, due to the in silico nature of this study, it is not possible to fully address potentially novel elements in Aplysia immunity. This dataset will serve as a springboard for future immune studies in Aplysia to further elucidate immune function in this venerable marine model.
In order to identify a suite of potential immune associated genes in Aplysia californica, conserved protein domain annotation via InterProScan was combined with homology search and tree based ortholog identification via OrthoFinder.
InterProScan (IPS) is a software suite that annotates queried proteins via several algorithms for in silico functional characterization [131, 132]. Annotations are given not only for the InterPro database, but also participating databases such as PFAM, SMART, PANTHER. For this analysis, the publicly available AplCal3.0 RefSeq reference genome assembly (GCF_000002075.1), generated from whole genome shotgun sequencing by the Broad Institute, was used . Aplysia RefSeq predicted protein models for this assembly (GCF_000002075.1_AplCal3.0_protein.faa) were annotated for conserved domains, motifs, and protein families via IPS (version 5.52–86.0) on the University of Miami Center for Computation Science PEGASUS supercomputer . IPS config file used can be found in supplemental file SFile4. The resulting database was then searched for proteins annotated with conserved domain suites that characterize described immune-associated genes in other organisms. The list of genes of interest and their requisite domains was constructed by combining previous in silico immune gene studies in Pacific oyster Crassostrea gigas  and pond snail Lymanaea stagnalis . This list was further expanded upon based on other reviews of immune-associated genes in vertebrates, mollusks, and arthropods [2, 4, 5, 9, 29,30,31, 33, 48, 55, 61, 63,64,65, 69, 84, 98, 110, 129, 134,135,136,137,138,139,140,141]. The full list of genes and requisite domains used can be found in supplemental file SFile5.
Conserved domain family screening of a proteome can recover putative members of known gene families with high confidence. Nevertheless, for proteins that are difficult to distinguish from a domain level alone (ie IkB with only Ankyrin repeat domains), this approach alone is insufficient. Furthermore, in several protein families of interest, molluscan proteins lack several domains or contain extra domains relative to their vertebrate and arthropod orthologs. In order to avoid missing these genes, previously identified molluscan representatives of these gene families were used in homology-based searches and validated with conserved domain data.
For more traditional homology searches for immune genes, reference predicted protein models for Aplysia, Crassostrea gigas (GCA_000297895.1_oyster_v9_protein.faa), and Biomphalaria glabrata (GCF_000457365.1_ASM45736v1_protein_Bglabrata.faa) from the RefSeq database [142, 143], and reference proteomes for Homo sapiens (UP000005640_9606_Homo_Sapies) and Drosophila melanogaster (UP000000803_7227_DroMel) from UNIPROT were analyzed with OrthoFinder. OrthoFinder combines reciprocal, all-against-all BLAST homology search with hierarchical clustering to identify probable orthologs and clusters of orthologs called orthogroups . These collections of orthologs were then queried for genes identified as immune associated in C. gigas , B glabrata , human, and Drosophila (SFile2).
The two sets of putative immune genes generated by these two analyses were then inspected and compared by hand to verify that hits contained requisite functional domains and were likely relatives of target genes. Data manipulation was done in the R statistical environment using RStudio and the tidyverse software suite [145, 146].
TIR domains of the longest protein isoform of each putative Aplysia TLR gene were extracted from RefSeq predicted protein fasta files using predicted TIR domains from IPS using tidyverse and the Biostrings R package . Extracted TIR domain sequences were aligned using the msa R package with default clustal parameters . A neighbor joining tree from Aplysia TLR TIR multi-sequence alignments was constructed with the Phangorn R package using the standard scripts for amino acid trees from the developer [149, 150]. Briefly, tree was built with a k of 4, a proportion of invariable sites of 2, and using the Jones-Taylor-Thornton model and bootstrapped with 1000 iterations. The resulting tree was visualized with the treeio and ggtree R packages [151,152,153,154].
TLR proteins for phylogenetic tree of TLRs were extracted from publicly available InterPro annotations of proteins from the UniprotKB data base for Biomphalaria glabrata (Bglab, UP000076420_IPS), Crassosteraea gigas (Cgig, UP000005408), Nematostella vectensis (Nvect, UP000001593), Drosophila melanogaster (Dromel, UP000000803), Strongylocentrotus purpuratus (Spurp, UP000007110), and human (Hs, UP000005640). These files were queried for TLRs based on TIR domain and LRR domains as was done with AplCal3.0 IPS annotation described above. Predicted TLR sequences were download from the UnirpotKB (or RefSeq for Aplysia and Nematostella) and used to generate a tree as described above for TIR domain sequences. The tree was not bootstrapped due to long iteration time. Tree was imported and visualized with treeio and ggtree.
Availability of data and materials
All data generated or analyzed during this study are included in this published article and its supplementary information files. All scripts used in this analysis are available on GitHub [https://github.com/Nicholas-Kron/Aplysia_Immunome].
Alpha-2-macroglobulin containing proteins
ATP-binding cassette sub-family F member-like
A disintegrin and metalloproteinase
Double-stranded RNA-specific adenosine deaminase
Mitochondrial apoptosis-inducing factor
Pyrin and hematopoietic interferon-inducible nuclear (HIN) domain (PYHIN) proteins
Activating protein 1 family
Apolipoprotein B Editing Complex 3 proteins
E3 ubiquitin-protein ligase TRIM23
Armadillo fold and TIR domain containing proteins
Apoptosis signal-regulating kinase (MAP3Ks)
Apoptosis regulator Bcl2- like
Tetherines (Bone marrow stromal antigen 2-like)
C1q domain containing proteins
Complement proteins 3/4/5
Caspase recruitment domain
Cluster of Differentiation 109
Monocyte differentiation antigen CD14
CD36/scavenger receptor class B member 1-like
Cyclic GMP-AMP synthase
C-terminal of Roc domain
Caspase and RIP adapter with death domain
Carbohydrate recognition domain
VIgL _C-type lectin-related protein
C-type lectin domain containing proteins
DNA damage-inducible transcript 3 protein (CHOP, c/EBPz/GADD153)
ATP-dependent RNA helicase with DEAD Q motif
ATP-dependent RNA helicase DDX1
ATP-dependent RNA helicase DDX21
ATP-dependent RNA helicase DDX41
DExD/H-box helicases - other
ATP-dependent RNA helicases with Domain of unknown function 1605
ATP-dependent RNA helicase DHX15
ATP-dependent RNA helicase DHX9
DNA-dependent protein kinase-like
Dual oxidase maturation factor
Evolutionarily conserved signaling intermediate in Toll pathway
Epidermal growth factor and TIR domain containing proteins
Extracellularly Regulated Kinase 1/2
Extracellularly Regulated Kinase 3/4
Estrogen receptor/oestrogen-related receptors
FAS-associated death domain protein
Fas apoptotic inhibitory molecule
Fibrinogen domain containing protein
Tyrosine-protein kinase Fyn
GTPase IMAP family member
Gram-negative bacteria binding-protein
VIgL _Galectin-related protein
Beta 1,3-glucan recognition proteins
Glutathion S transferase
- H-type lectin:
High mobility group box proteins
RanBP-type and C3HC4-type zinc finger-containing protein 1 RBCK1
E3 ubiquitin-protein ligase RNF31
Heat shock factor protein
Heat shock 20 kDa protein
Heat shock 70 kDa protein
Heat shock 90 kDa protein
Serine protease HTRA2
- I-type lectins:
Baculoviral IAP repeat-containing proteins
Interferon-induced protein 44
Interferon-induced transmembrane protein-like
Interferon alpha/beta receptor
Imunoglobulin and TIR domain containing proteins
Inhibitor of NFkB
Inhibitor of nuclear factor kappa-B kinase subunit alpha
Inhibitor of nuclear factor kappa-B kinase subunit epsilon
NF-kappa-B essential modulator
Interleukin 17 receptor
Interleukin-1 receptor-associated kinase 1
Interferon regulatory factor
Interferon stimulated gene
Tyrosine-protein kinase JAK
Mitogen-activated protein kinase 8
Activating protein 1 JUN
ATP-dependent DNA helicase II/X-ray repair cross-complementing protein 5/6
- L-type lectin:
- LAAO :
Lipopolysaccharide-binding protein/Bactericidal/permeability-increasing protein
Lipopolysaccharide-induced tumor necrosis factor-alpha factor
Leucine Rich Repeat
Leucine rich repeate and immunoglobin domain containing protein
Leucine rich repeat containg protein
Leucine-rich repeat flightless-interacting protein 1-like
Linear ubiquitin chain assembly complex
Membrane attack complex
MACPF domain containing protein (Perfoin-like)
Mitogen-activated protein kinase kinase 5 (MEK5)
Dual specificity mitogen-activated protein kinase kinase
Mitogen activated protein kinase kinase kinase 2/3
Mitogen activated protein kinase kinase kinase 10 (MST, MLK2)
Mitogen activated protein kinase kinase kinase 12/13 (DLK, MEKK12/13)
Mitogen activated protein kinase kinase kinase 4 (MEKK4, MTK1)
Dual specificity mitogen-activated protein kinase kinase kinase
Mitogen activated protein kinases
Mitochondrial Antiviral-Signaling Protein
Myeloid Differentiation factor 2/Lymphocyte antigen 96
Multiple epidermal growth factor-like domains protein 10
RNA-binding E3 ubiquitin-protein ligase MEX3C
Macrophage migration inhibitory factor
MAP kinase-activated protein kinase 2
Mitogen activated protein kinase phosphatase
Superoxide Dismutase (Mn/Fe)
Double-strand break repair protein MRE11-like
Major Vault Protein
Interferon-induced GTP-binding protein Mx
Myeloid differentiation primary response protein
NAIP, CIITA, HET-E and TP1 NTPase domain
Nuclear factor kB
- NXN :
- OAS :
- OrTIR :
Orphan TIR proteins
- P-type lectin :
Mannose 6-phosphate receptor homology domain containing
- p38 :
p38 MAPK/ MAPK14
- PAMP :
Pathogen associated molecular pattern
Pore forming protein
Peptidoglycan recognition proteins
Serine/threonine-protein phosphatase PP1
Pattern recognition receptor.
Pyrin and hematopoietic interferon-inducible nuclear (HIN) domain
DNA repair protein RAD50
RanBP2 domain containing proteins
RIG-I and MAVS type CARD domain containing proteins
RIO Kinase 3
Receptor-interacting serine/threonine-protein kinase
- RLR :
DNA-directed RNA polymerase III
E3 ubiquitin-protein ligase RNF135/RING finger protein leading to RIG-I activation
Ras of complex domain
ROC and COR domain
Viperins (Radical S-adenosyl methionine domain-containing proteins)
Retinoid X receptor
Deoxynucleoside triphosphate triphosphohydrolase SAMHD
NAD(+) hydrolase SARM1
SKIP1, Cullin, F-box protein complex
S-phase kinase-associated protein 1
- SMAD :
“Small” worm phenotype and “Mothers Against Decapentaplegic” homology domain
Suppressor of cytokine signaling
Superoxide Dismutase (Cu/Zn)
Scavenger receptor cysteine-rich domain containing protein
Receptor Cysteine-Rich (SRCR) Domain
Signal transducer and activator of transcription
Stimulator of interferon genes
TGF-beta-activated kinase 1 and MAP3K7-binding protein 1
TGF-beta-activated kinase 1 and MAP3K7-binding protein 2/3
MAP3K7/TGF-beta-activated kinase 1
TRAF family member-associated NF-kappa-B activator
TANK-binding kinase 1/ NFkB activated Kinase
TANK-binding kinase 1-binding protein 1
Thioester containing protein (thioester)
Toll/interleukin-1 receptor homology (TIR) domain
TIR and tetratrico peptide repeat containing proteins
Toll/interleukin-1 receptor domain-containing adapter protein
TIR domain containing
Tumor necrosis factor alpha-induced protein 3/A20-like
Tumor necrosis factor
Tumore necrosis factor receptor
Mitogen-activated protein kinase kinase kinase 8
Tumor necrosis factor receptor type 1-associated DEATH domain protein
TNF receptor-associated factors
TIR domain-containing adapter molecule 2
RISC-loading complex subunit TARBP
TIR domain-containing adapter molecule 1
E3 ubiquitin-protein ligase TRIM
E3 ubiquitin/ISG15 ligase TRIM25/65
E3 ubiquitin-protein ligase TRIM31
Tripartite motif-containing protein 5
Thioredoxin domain-containing protein
Variable Immunoglobulin and Lectin domain containing molecules
Zinc-finger antiviral proteins
Zinc finger NFX1-type containing 1.
Lafont M, Vergnes A, Vidal-Dupiol J, de Lorgeril J, Gueguen Y, Haffner P, et al. A Sustained Immune Response Supports Long-Term Antiviral Immune Priming in the Pacific Oyster, Crassostrea gigas. mBio. 2020;11(2):e02777–19.
Green TJ, Montagnani C. Poly I:C induces a protective antiviral immune response in the Pacific oyster (Crassostrea gigas) against subsequent challenge with Ostreid herpesvirus (OsHV-1 muvar). Fish Shellfish Immunol. 2013;35(2):382–8.
Zhang L, Li L, Guo X, Litman GW, Dishaw LJ, Zhang G. Massive expansion and functional divergence of innate immune genes in a protostome. Sci Rep. 2015;5:8693.
Dheilly NM, Duval D, Mouahid G, Emans R, Allienne JF, Galinier R, et al. A family of variable immunoglobulin and lectin domain containing molecules in the snail Biomphalaria glabrata. Dev Comp Immunol. 2015;48(1):234–43.
Li H, Gharamah AA, Hambrook JR, Wu X, Hanington PC. Single-cell RNA-seq profiling of individual Biomphalaria glabrata immune cells with a focus on immunologically relevant transcripts. Immunogenetics. 2022;74(1):77–98.
Mitta G, Galinier R, Tisseyre P, Allienne JF, Girerd-Chambaz Y, Guillou F, et al. Gene discovery and expression analysis of immune-relevant genes from Biomphalaria glabrata hemocytes. Dev Comp Immunol. 2005;29(5):393–407.
Pan B, Ren Y, Gao J, Gao H. De novo RNA-Seq analysis of the venus clam, Cyclina sinensis, and the identification of immune-related genes. PLoS One. 2015;10(4):e0123296.
McDowell IC, Nikapitiya C, Aguiar D, Lane CE, Istrail S, Gomez-Chiarri M. Transcriptome of American oysters, Crassostrea virginica, in response to bacterial challenge: insights into potential mechanisms of disease resistance. PLoS One. 2014;9(8):e105097.
Philipp EE, Kraemer L, Melzner F, Poustka AJ, Thieme S, Findeisen U, et al. Massively parallel RNA sequencing identifies a complex immune gene repertoire in the lophotrochozoan Mytilus edulis. PLoS One. 2012;7(3):e33091.
Seppälä O, Walser J-C, Cereghetti T, Seppälä K, Salo T, Adema CM. Transcriptome profiling of Lymnaea stagnalis (Gastropoda) for ecoimmunological research. BMC Genomics. 2021;22(1):144.
Castellanos-Martinez S, Diz AP, Alvarez-Chaver P, Gestal C. Proteomic characterization of the hemolymph of Octopus vulgaris infected by the protozoan parasite Aggregata octopiana. J Proteome. 2014;105:151–63.
Moroz LL. Aplysia. Curr Biol. 2011;21(2):R60–1.
Castellucci V, Pinsker H, Kupfermann I, Kandel ER. Neuronal mechanisms of habituation and dishabituation of the gill-withdrawal reflex in Aplysia. Science. 1970;167(3926):1745–8.
Moroz LL, Edwards JR, Puthanveettil SV, Kohn AB, Ha T, Heyland A, et al. Neuronal transcriptome of Aplysia: neuronal compartments and circuitry. Cell. 2006;127(7):1453–67.
Farr M, Mathews J, Zhu DF, Ambron RT. Inflammation causes a long-term hyperexcitability in the nociceptive sensory neurons of Aplysia. Learn Mem. 1999;6(3):331–40.
Clatworthy A, Castro G, Budelmann B, Walters E. Induction of a cellular defense reaction is accompanied by an increase in sensory neuron excitability in Aplysia. J Neurosci. 1994;14(5):3263–70.
Clatworthy AL. A simple systems approach to neural-immune communication. Comp Biochem Physiol A Physiol. 1996;115(1):1–10.
Clatworthy AL, Grose E. Immune-mediated alterations in nociceptive sensory function in Aplysia californica. J Exp Biol. 1999;202(Pt 5):623–30.
Bukhari K, Mulley G, Gulyaeva AA, Zhao L, Shu G, Jiang J, et al. Description and initial characterization of metatranscriptomic nidovirus-like genomes from the proposed new family Abyssoviridae, and from a sister group to the Coronavirinae, the proposed genus Alphaletovirus. Virology. 2018;524:160–71.
Debat HJ. Expanding the size limit of RNA viruses: Evidence of a novel divergent nidovirus in California sea hare, with a ~35.9 kb virus genome. bioRxiv. 2018:307678.
Kron NS, Fieber LA. Co-expression analysis identifies neuro-inflammation as a driver of sensory neuron aging in Aplysia californica. PLoS One. 2021;16(6):e0252647.
Jin P, Zhou L, Song X, Qian J, Chen L, Ma F. Particularity and universality of a putative gram-negative bacteria-binding protein (GNBP) gene from amphioxus (Branchiostoma belcheri): insights into the function and evolution of GNBP. Fish Shellfish Immunol. 2012;33(4):835–45.
Silverstein RL, Febbraio M. CD36, a scavenger receptor involved in immunity, metabolism, angiogenesis, and behavior. Sci Signal. 2009;2(72):re3.
Dziarski R, Gupta D. The peptidoglycan recognition proteins (PGRPs). Genome Biol. 2006;7(8):232.
Krasity BC, Troll JV, Weiss JP, McFall-Ngai MJ. LBP/BPI proteins and their relatives: conservation over evolution and roles in mutualism. Biochem Soc Trans. 2011;39(4):1039–44.
Liu Y, Liu J, Pang X, Liu T, Ning Z, Cheng G. The roles of direct recognition by animal lectins in antiviral immunity and viral pathogenesis. Molecules. 2015;20(2):2272–95.
Gorbushin AM, Panchin YV, Iakovleva NV. In search of the origin of FREPs: characterization of Aplysia californica fibrinogen-related proteins. Dev Comp Immunol. 2010;34(4):465–73.
Madeira F, Pearce M, Tivey ARN, Basutkar P, Lee J, Edbali O, et al. Search and sequence analysis tools services from EMBL-EBI in 2022. Nucleic Acids Res. 2022;50(W1):W276–9.
Goubau D, Deddouche S. Reis e Sousa C: cytosolic sensing of viruses. Immunity. 2013;38(5):855–69.
Su C, Tang YD, Zheng C. DExD/H-box helicases: multifunctional regulators in antiviral innate immunity. Cell Mol Life Sci. 2021;79(1):2.
Blasi G, Bortoletto E, Gasparotto M, Filippini F, Bai CM, Rosani U, et al. A glimpse on metazoan ZNFX1 helicases, ancient players of antiviral innate immunity. Fish Shellfish Immunol. 2022;121:456–66.
Nakad R, Schumacher B. DNA damage response and immune defense: links and mechanisms. Front Genet. 2016;7:147.
Paludan SR, Bowie AG. Immune sensing of DNA. Immunity. 2013;38(5):870–80.
Cavlar T, Ablasser A, Hornung V. Induction of type I IFNs by intracellular DNA-sensing pathways. Immunol Cell Biol. 2012;90(5):474–82.
Syed A, Tainer JA. The MRE11–RAD50–NBS1 complex conducts the orchestration of damage signaling and outcomes to stress in DNA replication and repair. Annu Rev Biochem. 2018;87(1):263–94.
Kondo T, Kobayashi J, Saitoh T, Maruyama K, Ishii KJ, Barber GN, et al. DNA damage sensor MRE11 recognizes cytosolic double-stranded DNA and induces type I interferon by regulating STING trafficking. Proc Natl Acad Sci U S A. 2013;110(8):2969–74.
Chakrabarti A, Jha BK, Silverman RH. New insights into the role of RNase L in innate immunity. J Interf Cytokine Res. 2011;31(1):49–57.
Pindel A, Sadler A. The role of protein kinase R in the interferon response. J Interf Cytokine Res. 2011;31(1):59–70.
Rauta PR, Samanta M, Dash HR, Nayak B, Das S. Toll-like receptors (TLRs) in aquatic animals: signaling pathways, expressions and immune responses. Immunol Lett. 2014;158(1–2):14–24.
Zhong Y, Kinio A, Saleh M. Functions of NOD-like receptors in human diseases. Front Immunol. 2013;4:333.
Sarrias MR, Gronlund J, Padilla O, Madsen J, Holmskov U, Lozano F. The scavenger receptor cysteine-rich (SRCR) domain: an ancient and highly conserved protein module of the innate immune system. Crit Rev Immunol. 2004;24(1):1–37.
Liu L, Yang J, Qiu L, Wang L, Zhang H, Wang M, et al. A novel scavenger receptor-cysteine-rich (SRCR) domain containing scavenger receptor identified from mollusk mediated PAMP recognition and binding. Dev Comp Immunol. 2011;35(2):227–39.
Deguine J, Barton GM. MyD88: a central player in innate immune signaling. F1000Prime Rep. 2014;6:97.
Belinda LW, Wei WX, Hanh BT, Lei LX, Bow H, Ling DJ. SARM: a novel toll-like receptor adaptor, is functionally conserved from arthropod to human. Mol Immunol. 2008;45(6):1732–42.
Xu YR, Lei CQ. TAK1-TABs complex: a central signalosome in inflammatory responses. Front Immunol. 2020;11:608976.
Zhang L, Ding X, Cui J, Xu H, Chen J, Gong YN, et al. Cysteine methylation disrupts ubiquitin-chain sensing in NF-kappaB activation. Nature. 2011;481(7380):204–8.
Jiang X, Chen ZJ. The role of ubiquitylation in immune defence and pathogen evasion. Nat Rev Immunol. 2011;12(1):35–48.
Xie P. TRAF molecules in cell signaling and in human diseases. J Mol Signal. 2013;8(1):7.
Tokunaga F, Iwai K. LUBAC, a novel ubiquitin ligase for linear ubiquitination, is crucial for inflammation and immune responses. Microbes Infect. 2012;14(7–8):563–72.
Lou Y, Han M, Song Y, Zhong J, Zhang W, Chen YH, et al. The SCF(beta-TrCP) E3 ubiquitin ligase regulates immune receptor signaling by targeting the negative regulatory protein TIPE2. J Immunol. 2020;204(8):2122–32.
Moynagh PN. The roles of Pellino E3 ubiquitin ligases in immunity. Nat Rev Immunol. 2014;14(2):122–31.
Ozato K, Shin DM, Chang TH, Morse HC 3rd. TRIM family proteins and their emerging roles in innate immunity. Nat Rev Immunol. 2008;8(11):849–60.
van Gent M, Sparrer KMJ, Gack MU. TRIM proteins and their roles in antiviral host defenses. Annu Rev Virol. 2018;5(1):385–405.
Bertrand MJ, Doiron K, Labbe K, Korneluk RG, Barker PA, Saleh M. Cellular inhibitors of apoptosis cIAP1 and cIAP2 are required for innate immunity signaling by the pattern recognition receptors NOD1 and NOD2. Immunity. 2009;30(6):789–801.
Manes NP, Nita-Lazar A. Molecular Mechanisms of the Toll-Like Receptor, STING, MAVS, Inflammasome, and Interferon Pathways. mSystems. 2021:e0033621.
Suzuki N, Suzuki S, Saito T. IRAKs: key regulatory kinases of innate immunity. Curr Med Chem Anti-Inflammatory Anti-Allergy Agents. 2005;4(1):13–20.
Li J, Chai QY, Liu CH. The ubiquitin system: a critical regulator of innate immunity and pathogen-host interactions. Cell Mol Immunol. 2016;13(5):560–76.
Israel A. The IKK complex, a central regulator of NF-kappaB activation. Cold Spring Harb Perspect Biol. 2010;2(3):a000158.
Shin CH, Choi DS. Essential roles for the non-canonical IkappaB kinases in linking inflammation to cancer, Obesity, and Diabetes. Cells. 2019;8(2):178.
Fujita F, Taniguchi Y, Kato T, Narita Y, Furuya A, Ogawa T, et al. Identification of NAP1, a regulatory subunit of IkappaB kinase-related kinases that potentiates NF-kappaB signaling. Mol Cell Biol. 2003;23(21):7780–93.
Huang G, Shi LZ, Chi H. Regulation of JNK and p38 MAPK in the immune system: signal integration, propagation and termination. Cytokine. 2009;48(3):161–9.
Shen Y, Tang K, Chen D, Hong M, Sun F, Wang S, et al. Riok3 inhibits the antiviral immune response by facilitating TRIM40-mediated RIG-I and MDA5 degradation. Cell Rep. 2021;35(12):109272.
Feng J, De Jesus PD, Su V, Han S, Gong D, Wu NC, et al. RIOK3 is an adaptor protein required for IRF3-mediated antiviral type I interferon production. J Virol. 2014;88(14):7987–97.
Newton K. RIPK1 and RIPK3: critical regulators of inflammation and cell death. Trends Cell Biol. 2015;25(6):347–53.
Wicherska-Pawlowska K, Wrobel T, Rybka J. Toll-like receptors (TLRs), NOD-like receptors (NLRs), and RIG-I-like receptors (RLRs) in innate immunity. TLRs, NLRs, and RLRs ligands as immunotherapeutic agents for hematopoietic diseases. Int J Mol Sci. 2021;22(24):13397.
Ghosh S, May MJ, Kopp EB. NF-κB AND REL PROTEINS: evolutionarily conserved mediators of immune responses. Annu Rev Immunol. 1998;16(1):225–60.
Shaulian E, Karin M. AP-1 as a regulator of cell life and death. Nat Cell Biol. 2002;4(5):E131–6.
Honda K, Takaoka A, Taniguchi T. Type I interferon gene induction by the interferon regulatory factor family of transcription factors. Immunity. 2006;25(3):349–60.
Wang W, Xiong L, Wang P, Wang F, Ma Q. Major vault protein plays important roles in viral infection. IUBMB Life. 2020;72(4):624–31.
Morris R, Kershaw NJ, Babon JJ. The molecular details of cytokine signaling via the JAK/STAT pathway. Protein Sci. 2018;27(12):1984–2009.
Gnanaprakasam JN, Wang R. MYC in regulating immunity: metabolism and beyond. Genes (Basel). 2017;8(3):88.
Holloway AF, Rao S, Shannon MF. Regulation of cytokine gene transcription in the immune system. Mol Immunol. 2002;38(8):567–80.
Gourbal B, Pinaud S, Beckers GJM, Van Der Meer JWM, Conrath U, Netea MG. Innate immune memory: an evolutionary perspective. Immunol Rev. 2018;283(1):21–40.
Lambrecht BN, Vanderkerken M, Hammad H. The emerging role of ADAM metalloproteinases in immunity. Nat Rev Immunol. 2018;18(12):745–58.
Sokolova IM. Apoptosis in molluscan immune defense. Isj-Invert Surviv J. 2009;6(1):49–58.
Everett H, McFadden G. Apoptosis: an innate immune response to virus infection. Trends Microbiol. 1999;7(4):160–5.
Banjara S, Suraweera CD, Hinds MG, Kvansakul M. The Bcl-2 family: ancient origins, Conserved Structures, and Divergent Mechanisms. Biomolecules. 2020;10(1):128.
Kaku H, Rothstein TL. Fas apoptosis inhibitory molecule expression in B cells is regulated through IRF4 in a feed-forward mechanism. J Immunol. 2009;183(9):5575–81.
Van Opdenbosch N, Lamkanfi M. Caspases in cell death, inflammation, and disease. Immunity. 2019;50(6):1352–64.
Bano D, Prehn JHM. Apoptosis-inducing factor (AIF) in physiology and disease: the tale of a repented natural born killer. EBioMedicine. 2018;30:29–37.
Vande Walle L, Lamkanfi M, Vandenabeele P. The mitochondrial serine protease HtrA2/Omi: an overview. Cell Death Differ. 2008;15(3):453–60.
Pales Espinosa E, Koller A, Allam B. Proteomic characterization of mucosal secretions in the eastern oyster, Crassostrea virginica. J Proteomics. 2016;132:63–76.
Dal Peraro M, van der Goot FG. Pore-forming toxins: ancient, but never really out of fashion. Nat Rev Microbiol. 2016;14(2):77–92.
Lacomel CJ, Dunstone MA, Spicer BA. Branching out the aerolysin, ETX/MTX-2 and Toxin_10 family of pore forming proteins. J Invertebr Pathol. 2021;186:107570.
Ragland SA, Criss AK. From bacterial killing to immune modulation: recent insights into the functions of lysozyme. PLoS Pathog. 2017;13(9):e1006512.
Zhu H, Zheng C. When PARPs meet antiviral innate immunity. Trends Microbiol. 2021;29(9):776–8.
Ke Y, Wang C, Zhang J, Zhong X, Wang R, Zeng X, et al. The role of PARPs in inflammation-and metabolic-related diseases: molecular mechanisms and beyond. Cells. 2019;8(9):1047.
Ruckert C, Bell-Sakyi L, Fazakerley JK, Fragkoudis R. Antiviral responses of arthropod vectors: an update on recent advances. Virusdisease. 2014;25(3):249–60.
Kohchi C, Inagawa H, Nishizawa T, Soma G. ROS and innate immunity. Anticancer Res. 2009;29(3):817–21.
Kamiya H, Muramoto K, Yamazaki M. Aplysianin-a, an antibacterial and antineoplastic glycoprotein in the albumen gland of a sea hare, Aplysia kurodai. Experientia. 1986;42(9):1065–7.
Coppo L, Ghezzi P. Thiol regulation of pro-inflammatory cytokines and innate immunity: protein S-thiolation as a novel molecular mechanism. Biochem Soc Trans. 2011;39(5):1268–72.
Zininga T, Ramatsui L, Shonhai A. Heat shock proteins as Immunomodulants. Molecules. 2018;23(11):2846.
Allam B, Raftos D. Immune responses to infectious diseases in bivalves. J Invertebr Pathol. 2015;131:121–36.
Le Clec'h W, Anderson TJ, Chevalier FD. Characterization of hemolymph phenoloxidase activity in two Biomphalaria snail species and impact of Schistosoma mansoni infection. Parasit Vectors. 2016;9:32.
Smith LC, Azumi K, Nonaka M. Complement systems in invertebrates. The ancient alternative and lectin pathways. Immunopharmacology. 1999;42(1–3):107–20.
Shokal U, Eleftherianos I. Evolution and function of thioester-containing proteins and the complement system in the innate immune response. Front Immunol. 2017;8:759.
Montagnani C, Kappler C, Reichhart JM, Escoubas JM. Cg-Rel, the first Rel/NF-κB homolog characterized in a mollusk, the Pacific oyster Crassostrea gigas. FEBS Lett. 2004;561(1):75–82.
Gerdol M. Immune-related genes in gastropods and bivalves: a comparative overview. Invertebr Surviv J. 2017;14(1).
Qiao X, Wang L, Song L. The primitive interferon-like system and its antiviral function in molluscs. Dev Comp Immunol. 2021;118:103997.
Zhang Y, He X, Yu F, Xiang Z, Li J, Thorpe KL, et al. Characteristic and functional analysis of toll-like receptors (TLRs) in the lophotrocozoan, Crassostrea gigas, reveals ancient origin of TLR-mediated innate immunity. PLoS One. 2013;8(10):e76464.
Huang B, Zhang L, Xu F, Tang X, Li L, Wang W, et al. Oyster versatile IKKalpha/betas are involved in toll-like receptor and RIG-I-like receptor signaling for innate immune response. Front Immunol. 1826;2019:10.
Yang M, Yuan S, Huang S, Li J, Xu L, Huang H, et al. Characterization of bbtTICAM from amphioxus suggests the emergence of a MyD88-independent pathway in basal chordates. Cell Res. 2011;21(10):1410–23.
Elvington M, Liszewski MK, Atkinson JP. Evolution of the complement system: from defense of the single cell to guardian of the intravascular space. Immunol Rev. 2016;274(1):9–15.
Nonaka M. Evolution of the complement system. Subcell Biochem. 2014;80:31–43.
Chen D, Hou Z, Jiang D, Zheng M, Li G, Zhang Y, et al. GILT restricts the cellular entry mediated by the envelope glycoproteins of SARS-CoV, Ebola virus and Lassa fever virus. Emerg Microbes Infect. 2019;8(1):1511–23.
Izumida M, Hayashi H, Smith C, Ishibashi F, Suga K, Kubo Y. Antivirus activity, but not thiolreductase activity, is conserved in interferon-gamma-inducible GILT protein in arthropod. Mol Immunol. 2021;140:240–9.
De Zoysa M, Kang HS, Song YB, Jee Y, Lee YD, Lee J. First report of invertebrate mx: cloning, characterization and expression analysis of mx cDNA in disk abalone (Haliotis discus discus). Fish Shellfish Immunol. 2007;23(1):86–96.
Verhelst J, Hulpiau P, Saelens X. Mx proteins: antiviral gatekeepers that restrain the uninvited. Microbiol Mol Biol Rev. 2013;77(4):551–66.
Stavrou S, Ross SR. APOBEC3 proteins in viral immunity. J Immunol. 2015;195(10):4565–70.
Lamers MM, van den Hoogen BG, Haagmans BL. ADAR1: “editor-in-chief” of cytoplasmic innate immunity. Front Immunol. 2019;10:1763.
Rosani U, Bai CM, Maso L, Shapiro M, Abbadi M, Domeneghetti S, et al. A-to-I editing of Malacoherpesviridae RNAs supports the antiviral role of ADAR1 in mollusks. BMC Evol Biol. 2019;19(1):149.
Shang JL, Smith MR, Anmangandla A, Lin HN. NAD(+)-consuming enzymes in immune defense against viral infection. Biochem J. 2021;478(23):4071–92.
Mattijssen S, Pruijn GJ. Viperin, a key player in the antiviral response. Microbes Infect. 2012;14(5):419–26.
Green TJ, Speck P, Geng L, Raftos D, Beard MR, Helbig KJ. Oyster viperin retains direct antiviral activity and its transcription occurs via a signalling pathway involving a heat-stable haemolymph protein. J Gen Virol. 2015;96(12):3587–97.
Liu S, Cai X, Wu J, Cong Q, Chen X, Li T, et al. Phosphorylation of innate immune adaptor proteins MAVS, STING, and TRIF induces IRF3 activation. Science. 2015;347(6227):aaa2630.
Takashima K, Oshiumi H, Takaki H, Matsumoto M, Seya T. RIOK3-mediated phosphorylation of MDA5 interferes with its assembly and attenuates the innate immune response. Cell Rep. 2015;11(2):192–200.
Huang B, Zhang L, Du Y, Xu F, Li L, Zhang G. Characterization of the Mollusc RIG-I/MAVS pathway reveals an archaic antiviral Signalling framework in invertebrates. Sci Rep. 2017;7(1):8217.
Liu Q, Li F, Liu W, Huang B, Li L, Wang X, et al. Transcriptional expression analysis reveals multiple effects of nonylphenol exposure on scallop immune system. Fish Shellfish Immunol. 2022;123:290–7.
Ren Z, Ding T, Zuo Z, Xu Z, Deng J, Wei Z. Regulation of MAVS expression and signaling function in the antiviral innate immune response. Front Immunol. 2020;11:1030.
Wickenhagen A, Sugrue E, Lytras S, Kuchi S, Noerenberg M, Turnbull ML, et al. A prenylated dsRNA sensor protects against severe COVID-19. Science. 2021;374(6567):eabj3624.
Zhao C, Zhao W. NLRP3 Inflammasome-a key player in antiviral responses. Front Immunol. 2020;11:211.
Kelley N, Jeltema D, Duan Y, He Y. The NLRP3 Inflammasome: an overview of mechanisms of activation and regulation. Int J Mol Sci. 2019;20(13):3328.
Zheng Z, Tang D, Zhao W, Wan Z, Yu M, Huang Z, et al. NLRP3-like protein negatively regulates the expression of antimicrobial peptides in Penaeus vannamei hemocyates. Fish Shellfish Immunol Rep. 2021;2:100039.
Lv Z, Wei Z, Zhang Z, Li C, Shao Y, Zhang W, et al. Characterization of NLRP3-like gene from Apostichopus japonicus provides new evidence on inflammation response in invertebrates. Fish Shellfish Immunol. 2017;68:114–23.
Rosenstiel P, Philipp EE, Schreiber S, Bosch TC. Evolution and function of innate immune receptors--insights from marine invertebrates. J Innate Immun. 2009;1(4):291–300.
Ozaki E, Campbell M, Doyle SL. Targeting the NLRP3 inflammasome in chronic inflammatory diseases: current perspectives. J Inflamm Res. 2015;8:15–27.
Qiao X, Zong Y, Liu Z, Wu Z, Li Y, Wang L, et al. The cGAS/STING-TBK1-IRF regulatory Axis orchestrates a primitive interferon-like antiviral mechanism in oyster. Front Immunol. 2021;12:689783.
Green TJ, Raftos D, Speck P, Montagnani C. Antiviral immunity in marine molluscs. J Gen Virol. 2015;96(9):2471–82.
Green TJ, Speck P. Antiviral defense and innate immune memory in the oyster. Viruses. 2018;10(3):133.
Xu J, Cherry S. Viruses and antiviral immunity in Drosophila. Dev Comp Immunol. 2014;42(1):67–84.
Blum M, Chang HY, Chuguransky S, Grego T, Kandasaamy S, Mitchell A, et al. The InterPro protein families and domains database: 20 years on. Nucleic Acids Res. 2021;49(D1):D344–54.
Jones P, Binns D, Chang HY, Fraser M, Li W, McAnulla C, et al. InterProScan 5: genome-scale protein function classification. Bioinformatics. 2014;30(9):1236–40.
Knudsen B, Kohn AB, Nahir B, McFadden CS, Moroz LL. Complete DNA sequence of the mitochondrial genome of the sea-slug, Aplysia californica: conservation of the gene order in Euthyneura. Mol Phylogenet Evol. 2006;38(2):459–69.
Buddenborg SK, Bu L, Zhang SM, Schilkey FD, Mkoji GM, Loker ES. Transcriptomic responses of Biomphalaria pfeifferi to Schistosoma mansoni: investigation of a neglected African snail that supports more S. mansoni transmission than any other snail species. PLoS Negl Trop Dis. 2017;11(10):e0005984.
Green TJ, Vergnes A, Montagnani C, de Lorgeril J. Distinct immune responses of juvenile and adult oysters (Crassostrea gigas) to viral and bacterial infections. Vet Res. 2016;47(1):72.
Meagher JL, Takata M, Goncalves-Carneiro D, Keane SC, Rebendenne A, Ong H, et al. Structure of the zinc-finger antiviral protein in complex with RNA reveals a mechanism for selective targeting of CG-rich viral sequences. Proc Natl Acad Sci U S A. 2019;116(48):24303–9.
Sparrer KMJ, Gableske S, Zurenski MA, Parker ZM, Full F, Baumgart GJ, et al. TRIM23 mediates virus-induced autophagy via activation of TBK1. Nat Microbiol. 2017;2(11):1543–57.
Hibino T, Loza-Coll M, Messier C, Majeske AJ, Cohen AH, Terwilliger DP, et al. The immune gene repertoire encoded in the purple sea urchin genome. Dev Biol. 2006;300(1):349–65.
Nauta AJ, Raaschou-Jensen N, Roos A, Daha MR, Madsen HO, Borrias-Essers MC, et al. Mannose-binding lectin engagement with late apoptotic and necrotic cells. Eur J Immunol. 2003;33(10):2853–63.
Huang J, Tilly D, Altman A, Sugie K, Grey HM. T-cell receptor antagonists induce Vav phosphorylation by selective activation of Fyn kinase. Proc Natl Acad Sci U S A. 2000;97(20):10923–9.
Park SH, Choi HJ, Yang H, Do KH, Kim J, Lee DW, et al. Endoplasmic reticulum stress-activated C/EBP homologous protein enhances nuclear factor-kappaB signals via repression of peroxisome proliferator-activated receptor gamma. J Biol Chem. 2010;285(46):35330–9.
DeJong RJ, Emery AM, Adema CM. The mitochondrial genome of Biomphalaria glabrata (Gastropoda: Basommatophora), intermediate host of Schistosoma mansoni. J Parasitol. 2004;90(5):991–7.
Zhang G, Fang X, Guo X, Li L, Luo R, Xu F, et al. The oyster genome reveals stress adaptation and complexity of shell formation. Nature. 2012;490(7418):49–54.
Emms DM, Kelly S. OrthoFinder: phylogenetic orthology inference for comparative genomics. Genome Biol. 2019;20(1):238.
Wickham H, Averick M, Bryan J, Chang W, McGowan L, François R, et al. Welcome to the Tidyverse. J Open Source Software. 2019;4(43):1686.
R Core Team: R. A language and environment for statistical computing. Vienna: R Foundation for Statistical Computing; 2013.
Pagès H, Aboyoun P, Gentleman R, DebRoy S. Biostrings: Efficient manipulation of biological strings. 2.64.0 ed; 2022.
Bodenhofer U, Bonatesta E, Horejš-Kainrath C, Hochreiter S. msa: an R package for multiple sequence alignment. Bioinformatics. 2015;31(24):3997–9.
Schliep K, Potts AJ, Morrison DA, Grimm GW, Fitzjohn R. Intertwining phylogenetic trees and networks. Methods Ecol Evol. 2017;8(10):1212–20.
Schliep KP. Phangorn: phylogenetic analysis in R. Bioinformatics. 2010;27(4):592–3.
Yu G, Smith DK, Zhu H, Guan Y, Lam TT. ggtree: an R package for visualization and annotation of phylogenetic trees with their covariates and other associated data. Methods Ecol Evol. 2017;8:28–36.
Yu G. Using ggtree to visualize data on tree-like structures. Curr Protoc Bioinformatics. 2020;69(1):e96.
Yu G, Lam TT-Y, Zhu H, Guan Y. Two methods for mapping and visualizing associated data on phylogeny using Ggtree. Mol Biol Evol. 2018;35(12):3041–3.
Wang L-G, Lam TT-Y, Xu S, Dai Z, Zhou L, Feng T, et al. Treeio: an R package for phylogenetic tree input and output with richly annotated and associated data. Mol Biol Evol. 2019;37(2):599–603.
The author gratefully acknowledges the helpful insight and support of Dr. Lynne Fieber and Dr. Michael Schmale.
This work was funded by the National Institutes of Health Grant (P40OD010952). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Ethics approval and consent to participate
Consent for publication
The author declares that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Full results of InterProScan annotation of the Aplysia reference proteome (AplCal3.0 GCF_000002075.1).
OrthoFinder Results. OrthoFinder Results. First page contains the overall statistics from OrthoFinder about orthogroups. The second page contains pairwise ortholog counts among species surveyed. Page three contains per species statistics. The remaining pages contain orthologs and ortholog families (orthogroups) identified by OrthoFinder using reference proteomes for Aplysia, Crassostrea gigas (GCA_000297895.1_oyster_v9_protein.faa), and Biomphalaria glabrata (GCF_000457365.1_ASM45736v1_protein_Bglabrata.faa) from the RefSeq database, and reference proteomes for Homo sapiens (UP000005640_9606_Homo_Sapies) and Drosophila melanogaster (UP000000803_7227_DroMel). Each page labeled is labeled with the accession number and organism name for each organism. The last page is a complete collection of all proteins assigned to all computed orthogroups. The complete output of OrthoFinder including gene trees are available from the author upon request.
Putative Aplysia Immunome gene set. A Complete list of RefSeq gene and protein identifiers of genes identified as putative immune associated gene based on InterProScan and OrthoFinder results.
InterProScan configuration file. This file contains the configuration parameters used to annotate the Aplysia AplCal3.0 reference proteome.
Target genes and requisite InterPro domains. The list of genes of interest and their requisite domains used to identify putative immune genes in Aplysia californica. The list combines data from previous mollusk immunome in silico studies in Pacific oyster Crassostrea gigas  and pond snail Lymanaea stagnalis  and reviews of immune-associated genes in vertebrates, mollusks, and arthropods [2, 4, 5, 9, 29,30,31, 33, 47, 54, 60, 62,63,64, 68, 88, 98, 117, 121,122,123,124,125,126,127,128,129]. Genes are categorized according to their broad immune function and more specific gene family. Required InterPro domains are listed as strict (as appears in human) or relaxed (required domains as seen in non-human or non-mammalian systems). Human UNIPROT identifiers of example genes are listed and were futher used in OrthFinder homology search comparisons. For some genes that exist only in ecdysozoans, the Drosophila UNIPROT identifier is also given. The DOI and link to the source literature is also provided.
About this article
Cite this article
Kron, N.S. In search of the Aplysia immunome: an in silico study. BMC Genomics 23, 543 (2022). https://doi.org/10.1186/s12864-022-08780-6