In total, 60 healthy adult specimens of H. tubulosa (length 11 ± 0.98 cm and body weight 46 ± 7.5 g) were collected in the Gulf of Palermo (Sicily, Italy). The animals were acclimatized for 1 week in the laboratory tanks in a constantly aerated aquarium at a temperature of 15 ± 2°C and were fed with a commercial invertebrate feed until 24h before the sampling (Algamac 3000, Aquafauna BioMarine Inc., Hawthorne, CA, USA).

The coelomic fluid (CF) was collected by making a 3–5 cm incision on the antero-dorsal side of each individual using a scalpel, allowing the fluid, composed of both humoral and cellular components, to percolate naturally into 50 mL Falcon tubes pre-filled with Etylenedimine tetra-acetic acid (ISO-EDTA) anticoagulant (20 mM Tris, 0.5 mM NaCl, 70 mM EDTA; pH 7.5). This procedure ensured the prevention of clotting and the preservation of cellular integrity. The CF was collected, kept on ice, and immediately centrifuged at 1,000 g for 10 minutes at 4°C to separate the cellular component (cell pellet) and the humoral component (cell-free coelomic fluid). The cell pellet was pottered and sonicated with 1X RIPA buffer added with 1:200 antiprotease. The sample was centrifuged at 20,000 g for 10 min at 4°C. The supernatant (cell pellet lysate) and the cell-free coelomic fluid were frozen at -80°C and lyophilized. To obtain a sufficient amount of proteins for the analysis, all cell lysates were combined to form a single sample. The same procedure was followed for the cell-free coelomic fluids. The proteins presented in the humoral fraction correspond to extracellular proteins that have been secreted by the coelomocytes and the cellular component refers to the intracellular proteins that have been isolated from the coelomocytes.

The freeze-dried samples were dissolved in 0.25% RapiGest (Waters Co., Milford, MA, USA) and passed through Microcon 10 kDa filters (Merck KGaA, Darmstadt, Germany) to remove the RIPA buffer (the flowthrough was discarded). The nominal protein concentration was then measured for each sample using a NanoDrop 2000 UV-VIS spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). Moreover, 10 µg aliquots of each sample were taken in 30 µl of 5% methanol, reduced using RapiGest and dithiothreitol for 30 minutes at 60°C (Thermo Fisher Scientific, Waltham, MA, USA), and subsequently alkylated using iodoacetic acid (Thermo Fisher Scientific) for 30 minutes at room temperature in the dark in 25 mM ammonium bicarbonate buffer (Thermo Fisher Scientific). Samples were digested in-solution using endoproteinase LysC-trypsin (Mass Spec grade, Promega, Madison, WI, USA) at a 1:100 ratio for 1 h and then trypsin (Mass Spec grade, Promega) at a 1:25 ratio for 3 h at 37°C. Proteolysis was stopped by adding 1 µl of 100% formic acid (Thermo Fisher Scientific), and the samples were dried down and cleaned up using C18 spin columns (Thermo Fisher Scientific) according to the manufacturer’s protocol. The cleaned peptide extracts were dried down and stored at -20°C until further analysis. Aliquots of a nominal 1 µg of the tryptic digests were analyzed using a Dionex Ultimate 3000 nanoRSLC (Dionex, Sunnyvale, CA, USA) coupled to a Bruker Maxis II ETD mass spectrometer (Bruker Daltonics GmbH, Bremen, Germany) via a CaptiveSpray nanobooster ion source. The samples were first loaded on the trap column using 0.1% trifluoroacetic acid at a flow rate of 5 µl/min for 8 minutes using an Acclaim PepMap100 C-18 trap column (100 µm × 20 mm, Thermo Scientific). For the gradient elution, an ACQUITY UPLC M-Class Peptide BEH C18 column was used (130 Å, 1.7 µm, 75 µm × 250 mm, Waters) at a 300 nL/min flow rate and 48°C column temperature using a linear gradient from 4% B to 50% B in 120 minutes. Solvent A was 0.1% formic acid, and solvent B was acetonitrile with 0.1% formic acid. The cycle time for data-dependent acquisition was 2.5 s. MS spectra were acquired at 3 Hz, while MS/MS spectra were acquired at 4 or 16 Hz, depending on the intensity of the precursor ion. Singly charged ions were excluded from the analysis (36). The protein quantification was based on the mass-to-charge ratio, which represents the ratio of the mass number and the charge number, z.

The proteins were first identified by searching against the Uniprot Aechinodermata (downloaded: 07/03/2020) database using the Byonic software search engine (v3.8.13, Protein Metrics Inc, San Carlos, CA, USA) (RRID: SCR_016735) with the following parameters: 1% false discovery rate (FDR), 20 ppm peptide mass tolerance, 30 ppm fragment mass tolerance, two missed cleavages, trypsin as the enzyme, carbamidomethylation of cysteines as the fixed modification; and the following variable modifications: Oxidation/+15.994915 @ M, Deamidated/+0.984016 @ N, Deamidated/+0.984016 @ Q, Gln- >pyro-Glu/-17.026549 @ NTerm Q, Glu->pyro-Glu/-18.010565 @ NTerm E. Protein hits were filtered using Scaffold (version 4.11, Proteome Software, Inc., USA) (RRID: SCR_014345) using the same parameters stated above in addition to the following parameters: protein grouping strategy, Experiment-wide grouping with protein cluster analysis; peptide thresholds, 95.0% minimum; protein thresholds, 1% FDR and 2 peptides minimum. Subsequently, protein identification was also performed using BlastP comparison to non-redundant protein sequence and model organism databases (available at https://blast.ncbi.nlm.nih.gov/Blast.cgi?PAGE=Proteins; accessed in November 2024) (RRID: SCR_004870). An expected value of < 1 was set as cutoff. Figure 1 shows the workflow of proteomic data, starting from sample preparation to downstream data analysis ( Figure 1 ).

The Uniprot web tool (https://www.uniprot.org/) (RRID: SCR_002380) was used for protein annotation. It is a web resource for protein sequence and functional information through ID mapping (https://www.uniprot.org/id-mapping).

The Search Tool for the Retrieval of Interacting Genes/Protein (STRING) database (RRID: SCR_005223) (https://string-db.org/), which allows the visualization of complex networks (through clustering analysis), was used to retrieve the predicted interactions for the identified proteins. This web tool returns different information: first, it can reconstruct potential protein-protein interactions, aiming to evidence specific sub-networks from input data. Second, the web tool can perform clustering analysis. Finally, STRING can provide information on functional enrichment and pathway analysis for each cluster. We used the following settings for STRING analysis: full STRING network as the network type, in which edges indicate both functional and physical protein associations; evidence as the meaning of the network edges; different types of interaction sources, such as text-mining, neighborhood, experiments, co-occurrence, and more; and a medium confidence cut-off of 0.400 as the minimum required interaction score to highlight the more significant interactions. Proteins were clustered according to the k-means algorithm, an unsupervised clustering algorithm based on an adjacency matrix, and values with a cut-off score of FDR <0.05 were considered statistically significant and evaluated for further analysis. Each cluster produced by the algorithm was analyzed using Gene Ontology (GO). The STRING database is also linked to the SMART database (http://smart.embl-heidelberg.de/) (RRID: SCR_005026), which allowed us to compare different protein domains of proteome data of H. tubulosa and compare them with domains present in the human species.

Besides the STRING “clustering” function, the “homology search” function (string- db.org/cgi/proteinhomology) was used to calculate the percentage of protein sequence similarity between different species. H. tubulosa extracellular protein fractions were compared with those of the Homo sapiens species, and protein sequences were blasted and evaluated by using a STRING similarity bit score value. These are Smith–Waterman scores calculated for alignments between two random or unrelated sequences and are internally used by STRING as a proxy for protein homology. The bit score measures sequence similarity independent of query sequence length and database size and is normalized based on the raw-pairwise alignment score. The original scores were computed by the similarity matrix of proteins (SIMAP) project (SIMAP, http://mips.gsf.de/simap/) (RRID: SCR_007927). The bitscore is defined by “S” and it is determined by the following formula:

Where:

λ is the Gumble distribution constant,

S is the raw alignment score,

K is a constant associated with the scoring matrix used.

Clearly, the bit score (S) is linearly related to the raw alignment score (S). Thus, the higher the bit score, the more highly significant the match is. The bit score provides a constant statistical indicator for searching different databases of different sizes or for searching the same database at different times as the database grows (www.biostars.org/p/187230/). The “homology search” mode is the most popular strategy for inferring functional similarity between the same proteins of two different species. It is known that homologous sequences have similar structures, and frequently, they have similar functions as well. Scores less than 50 were not significant, and proteins with scores less than 50 were discarded from further analyses.

Proteomic analysis of extracellular fractions produced 174 unique annotated proteins. A workflow of proteomic data analysis is shown in Figure 1 where detailed steps of sample preparation and data analysis are summarized (see Figure 1 in the Materials and Methods section). The Uniprot web tool was used for protein annotation analysis (https://www.uniprot.org/).

Of the 174 proteins identified by mass spectrometry, 6 were common contaminants and were discarded from subsequent analysis, and 57 proteins were annotated as uncharacterized proteins. This means that, although they have been identified experimentally, they still do not have a known functional or structural domain. Moreover, these proteins do not have a homolog in humans, and they were also discarded from the analysis. Furthermore, 27 proteins did not have information about coverage percentage or log probability score from MS analysis, and thus 111 annotated proteins remained that were suitable for the subsequent analysis. They were analyzed using the STRING web tool to elucidate the potential interactions between the proteins found by the humoral fraction. When analyzing the remaining 111 proteins in the STRING database, the tool returned as the output of the analysis just 59 proteins, thus the remaining portion was lost during this last analysis.

Three main clusters ( Figure 2 , clusters A–C) were found, and they were composed as follows: the first cluster was composed of a few proteins linked to fibrinogen (cluster A, three proteins), the second cluster (cluster B, two proteins) contained Vitellogenin domain-containing protein (VTG domain-containing protein) and a putative endoplasmic reticulum resident protein 44 (ERp44), and the third cluster (cluster C, three proteins) consisted of Histone proteins H1-β, H2A, and H2B. Besides the STRING “clustering” function, we used the “homology search” (string- db.org/cgi/proteinhomology) function that calculates the percentage of protein sequence similarity between different species. Thus, the sequences of H. tubulosa extracellular protein fractions were compared with H. sapiens, and protein sequences were blasted and evaluated using a STRING similarity bit score value.

Table 1 shows the annotated proteins of the humoral fraction of coelomic fluid and their relative STRING similarity sequence analyses. Sequence homology analysis of extracellular protein fraction showed a high sequence conservation between H. tubulosa and humans. Indeed, the STRING similarity scores were between 639.4 (DMBT1), the highest homology score, and 73.2 (FREP-A). Most immune molecules showed a similarity score of approximately 150, indicating a high protein sequence similarity. Functional analysis of proteins in the humoral fraction of coelomic fluid allowed us to classify them into four main categories: immune receptor factors (IRFs), non-canonical immune proteins (NCIPs), signaling molecules, and effector proteins.

As shown in Table 1 , the first class of immune-related proteins can be divided into different sub-classes such as scavenger receptor cysteine-rich (SRCR) proteins [Putative deleted in malignant brain tumors 1 protein (DMBT1)], lectins (putative techylectin-5B-like), pattern recognition receptors (PPRs) (putative ficolin-2-like protein), fibrinogen-like family (FREP) proteins, and molecules with antibacterial activity or involved in neutrophil extracellular trap (NET) formation. The second class identified in the humoral component of coelomic fluid is represented by non-canonical proteins interacting with pathogens (NCIPs). To this second class belong proteins such as beta-actin and collagen. The third class of immunity-related proteins is involved in the signaling cascade. It includes vesicular trafficking molecules, calcium-binding proteins, cell polarity and cytoskeleton proteins, extracellular matrix (ECM) proteins, and finally polyubiquitination proteins. The fourth class of immune-related proteins found in the humoral compartment is represented by effector proteins. They belong only to the transferrin family: melanotransferrin-4 (MELTF-4) and Major York proteins 1 and 2.

Proteomic analysis of intracellular fractions produced 215 annotated proteins. STRING analysis of the intracellular fractions of H. tubulosa coelomic fluid proteins evidenced three different clusters ( Figure 3 ). Cluster A was mainly represented by proteins linked to cell signaling, cytoskeleton, and actin remodeling. Cluster B included proteins of different metabolic pathways. Cluster C was linked to ribosomal activity, ubiquitin proteasome system (UPS) machinery, and primary cilium. A graphical representation of some of the proteins present in the three clusters is also shown in Figure 4 .

Cluster A ( Figure 3 ) refers to proteins linked to different cellular mechanisms such as cell motility, calcium binding, cytoskeleton remodeling, or signaling networks. STRING analysis found GO terms linked to biological processes (BP) such as synaptic vesicle priming and actin filament network. Furthermore, molecular function (MF) and cellular component (CC) terms confirmed BP ontology, with MF terms such as actin filament binding, structural constituent of the cytoskeleton, and phosphatidylinositol-4,5-biphosphate binding, and CC terms such as actin and cytoskeleton ( Figure 5 ). This last cluster is mainly represented by proteins belonging to pathways such as phagosome, Salmonella infection, MAPK signaling, and endocytosis, showing the putative role of these proteins in processes linked to vesicle trafficking and cytoskeleton remodeling after pathogen stimuli ( Figure 5 ). A detailed list of proteins linked to GO and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway STRING analysis is in the Supplementary Material ( Supplementary Material Table , sheet GO and KEGG cluster A, sheet proteins A). A STRING homology search was also applied for all the proteins belonging to cluster A to evidence the sequence homology between H.tubulosa and H. sapiens. The results of STRING homology search are shown in Table 2 . The STRING similarity scores were between 944.1 (FLNA), which is the highest homology score, and 84.3 (PLS3). Most molecules showed a similarity score of approximately 300, indicating a high sequence similarity.

Cluster B ( Figure 3 ) is linked to BP terms such as pentose-phosphate shunt, nucleosome assembly, response to oxidative stress, and carbohydrate metabolism. MF terms were linked to different functions of the pentose-phosphate pathway, such as transketolase activity; cellular oxidant detoxification, such as glutathione peroxidase activity; and the carbohydrate metabolic pathway, such as alpha-amylase activity ( Figure 5 ). The KEGG analysis was in accordance with the GO findings previously shown ( Figure 5 ). A detailed list of proteins linked to the STRING analysis is in the Supplementary Material and a graphical representation of the main pathways in which proteins of cluster B are involved is shown in Figure 4 . ( Supplementary Material, Table , sheet GO and KEGG cluster B, sheet proteins B). A STRING homology search was also applied for all the proteins belonging to cluster B to evidence the sequence homology between H. tubulosa and H. sapiens. The results of STRING homology search are shown in Table 3 . The STRING similarity scores were between 750.0 (HAL), the highest homology score, and 95.5 (GPX3). Most of the molecules showed a similarity score of approximately 200.

The proteins of cluster C ( Figures 3 , 4 ) can be grouped into three sub-clusters: proteins linked to primary cilium and BBsome, proteins linked to UPS machinery, and proteins related to ribosomal machinery. The BP terms of this cluster were as follows: cytoplasmic translational elongation, axonemal dynein complex assembly, positive regulation of the proteasomal ubiquitin-dependent protein complex, sperm axoneme assembly, and non-motile cilium assembly. MF terms associated with cluster C are large subunit rRNA binding and structural constituents of the ribosome. Finally, the CC terms of cluster C are BBsome, cytosolic large ribosomal subunit, ciliary membrane, non-motile cilium, and centrosome ( Figure 5 ). Furthermore, the KEGG analysis was in accordance with the GO results previously shown ( Figure 5 ). The Supplementary Material ( Supplementary Material Table , sheet GO and KEGG cluster C, sheet proteins C) provides a detailed list of the proteins linked to the STRING analysis. A STRING homology search was also applied for all the proteins belonging to cluster C to evidence the sequence homology between H.tubulosa and H. sapiens. The results of the STRING homology search are shown in Table 4 . Most molecules showed a similarity score higher than 200. Finally, some of the proteins belonging to all three clusters of intracellular fractions, such as Tubulin-A, Rab, F-actin capping protein, Cdc42, and Polyubiquitin were linked to both the Wnt and Hedgehog pathways. In particular, RhoA, and Dynein were present exclusively in the Wnt pathway, and Calreticulin, Tubulin-beta, BBS, and Elav were present in the Hedgehog pathway.