A proteo-transcriptomic investigation of toxin evolution in planarians and their role in flatworm terrestrialization

The transition from aquatic to terrestrial environments represents a major evolutionary transition in animals, requiring significant adaptations in physiology and defense mechanisms to the challenges presented by the harsh terrestrial environment. Platyhelminthes, which include both aquatic and terrestrial species with a single terrestrialization event in the family Geoplaniidae, serve as excellent model organisms for studying the evolutionary adaptations required for terrestrialization. This study investigates the evolutionary dynamics of toxin orthologous groups (as a proxy to gene families) in aquatic and terrestrial flatworms, together with mucus composition, focusing on their role in terrestrialization from a molecular ecology perspective. Using a proteo-transcriptomic approach, we predicted and identified a broader toxin gene repertoire in terrestrial flatworms compared to freshwater ones. Although most toxins in flatworms arose before terrestrial planarian diversification—gaining a novel evolutionary origin .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 14, 2025. ; https://doi.org/10.1101/2025.01.09.632098doi: bioRxiv preprint 2 at the Tricladida and Continenticola nodes—the mucus protein repertoire appears to have a far older evolutionary origin in both species. Moreover, distinct orthologous groups underpin the toxin gene repertoire and mucus composition in each lineage, highlighting the contrasting evolutionary trajectories of these two functional components. While toxin families in both aquatic and terrestrial flatworms revealed overall common functions, including cytokine modulation and ion channel regulation, terrestrial flatworms exhibited specific expansions of lectin-like proteins and pro-inflammatory responses, highlighting their potential key role to respond to land-based threats. This study provides new insights into the differential evolutionary trajectories of toxin and mucus proteins in planarians, offering a deeper understanding of the genetic innovations that facilitated flatworm terrestrialization.

The transition from aquatic to terrestrial environments marks one of the most significant events in evolutionary history, occurring independently across multiple animal phyla. This transition demanded profound physiological adaptations, including changes in locomotion, respiration, reproduction, predation, and defense against novel pathogens (Little, 1983). Genomic changes by gene gain, duplication and loss have accompanied these evolutionary shifts, driving the diversification of terrestrial lineages across various animal groups (Aristide & Fernández, 2023; Meyer et al., 2021). Some orthologous groups (as a proxy to gene families) that expanded or contracted in terrestrial species are likely directly related to overcoming the challenges posed by this transition, as observed for instance in aquaporins (Martínez-Redondo et al., 2023), directly involved in osmoregulation, and more expanded overall in freshwater and terrestrial lineages compared to marine ones. A key challenge related to terrestrialization is the exposure to new pathogens, new predators and new prey. T oxins, typically peptides or proteins with harmful effects on other organisms, are key molecular players in predation and defense. These compounds have independently evolved from non-toxin proteins in at least eight different phyla (Casewell et al., 2013; Schendel et al., 2019). Despite their diversity, toxin molecular functions and domains often exhibit convergent evolution across lineages (Casewell et al., 2013; Fry et al., 2009; Zancolli et al., 2022). While venomous vertebrates and certain arthropod groups have traditionally received more research .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 14, 2025. ; https://doi.org/10.1101/2025.01.09.632098doi: bioRxiv preprint 3 focus, recent advances in venomics are increasingly uncovering toxins in under-studied invertebrates (Jenner et al., 2019; Verdes et al., 2018, 2022). Notably, invertebrates, which lack adaptive immune responses (Kloc et al., 2024; Yacoub et al., 2020), also use toxins as a defense against pathogens, making them promising sources of antimicrobial peptides (Yacoub et al., 2020). Despite this growing interest, the mechanisms driving toxin evolution across invertebrate lineages, particularly in response to terrestrialization, remain poorly understood. One promising group for studying these processes is the Platyhelminthes (flatworms). T errestrialization within this phylum is confined to the monophyletic clade Geoplaniidae, part of the Tricladida order, which also includes freshwater and marine species (Riutort et al., 2012). The transition to land in planarians was accompanied by notable gene gains at the Tricladida node, which facilitated adaptation to terrestrial habitats (Benítez-Álvarez et al. 2025a). These gene repertoire changes likely played critical roles in addressing the new environmental challenges, including potentially toxin evolution. Free-living flatworms are predators, and although all species utilize a muscular protrusible pharynx for feeding (Giribet & Edgecombe, 2020), hunting strategies vary between taxa. For example, while freshwater species like those in the Planariidae and Dugesiidae families often feed on weakened or immobile prey (Vila-Farré & C Rink, 2018), terrestrial species, such as those in the Geoplanidae, are active predators. These terrestrial flatworms use physical force and adhesive mucus to capture prey (Cuevas-Caballé et al., 2019; Sluys, 1999). In these species, mucus plays a key role in immobilizing prey, and it has been suggested to have neurotoxic properties (Thielicke and Sluys 2019; Sluys 2019; Cardoso et al. 2023), though the mechanisms remain unclear. Planarians lack specialized venom glands, complicating the study of their toxins. However, toxins have been identified across the Platyhelminthes phylum, with most studies focusing on parasitic species. T oxins have also been reported in free-living flatworms, including tetrodotoxin (TTX) in some terrestrial species of the genus Bipalium (Stokes et al., 2014). In these species, TTX is concentrated in the head and eggs, similar to other marine flatworms (Miyazawa et al., 1987), suggesting region-specific toxin expression, particularly in body parts involved in hunting, such as the pharynx or secreted mucus, and offspring protection. .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 14, 2025. ; https://doi.org/10.1101/2025.01.09.632098doi: bioRxiv preprint 4 Although toxins may be localized in specific regions, expression patterns likely vary across species, and the absence of clear venom delivery systems or specialized venom-producing tissues makes toxin characterization challenging. Traditional methods like BLAST are insufficient for predicting protein toxicity (Sharma et al., 2022). However, recent advances in machine learning and deep learning provide more reliable predictions (Cole & Brewer, 2019; Sharma et al., 2022; Vishnoi et al., 2020), offering new avenues for toxin discovery. In this study, we aim to characterize the toxin repertoire of a freshwater (Schmidtea mediterranea) and a terrestrial planarian (Obama nungara) using a proteo-transcriptomic approach, with the goal of further understanding how toxin diversity has evolved in response to different ecological pressures, particularly the transition from aquatic to terrestrial habitats. We employed two complementary methods: (i) in silico prediction of putative toxins and (ii) characterization of mucus secretions, focusing on identifying proteins with potential toxin activity. Additionally, we explored the evolution of putative toxins and mucus proteins across the Platyhelminthes phylum to understand their evolutionary origins and their roles in terrestrial adaptation, investigating how these proteins may have contributed to novel predatory strategies, defense mechanisms, and resistance to pathogens in land-dwelling species. This comprehensive approach sheds light on the molecular innovations that have accompanied the conquest of land within this animal lineage.

RNA and protein extraction Specimens of Obama nungara were collected in Sopelana, Vizcaya, while specimens of the freshwater planarian Schmidtea mediterranea were obtained from a lab culture established from specimens from Montjuic (Barcelona). All animals were maintained in starvation before sample processing (2-3 days for the O. nungara specimens; 7 days for the S. mediterranea specimens). Each specimen was divided into three body parts (assigned as head, body and pharynx) by fixing first the specimens with RNAlater, dissecting with a scalpel, and freezing immediately with liquid nitrogen. Samples were stored at -80ºC until the extraction. For O. nungara, a fragment from each body part was used per replicate for the RNA and protein extraction, while for S. mediterranea we combined four fragments from four different individuals per replicate. Both RNA and proteins were extracted by using TRIzol ® reagent (Invitrogen, USA). The frozen .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 14, 2025. ; https://doi.org/10.1101/2025.01.09.632098doi: bioRxiv preprint 5 samples were homogenized by adding 400μl of TRIzol ® and using an electrical pestle until homogenous. The RNA was extracted by following the manufacturer’s instructions, while for proteins the TRIzol ® protocol was modified (following (Simões et al., 2013)), and the final protein pellets were resuspended in 25μl of resuspension buffer (1%SDS and 8M Urea in Tris-HCl 1M, pH 8). The mucus proteins were collected by submerging one specimen of O. nungara or S. mediterranea in 5% N-Acetyl Cysteine (NAC) in PBS buffer and incubating for 8 minutes, followed by the collection of the supernatant and overnight precipitation with 4 volumes of acetone at -20C, centrifugation at 20,000g for 20 minutes and resuspension in 50μl for O. nungara or 30μl for S. mediterranea of resuspension Buffer. For each species, we sequenced 3 replicates of each body region, including mucus. The protein profiles of all the protein extracts were examined in a SDS page before preparing the sample for proteomic sequencing to make sure that there wasn’t any overrepresented protein masking other proteins. 4–20% Mini-PROTEAN® TGX™ Precast Protein Gels (Bio-Rad, Hercules, CA, US) were used with 10x Tris/Glycine/SDS as a running buffer. PageRuler™ Unstained Protein Ladder (Thermo Scientific, San Jose, CA, USA) was also loaded in the SDS page to serve as reference. De novo transcriptome assembly and differential expression analysis The concentration of all samples was determined using the Qubit RNA BR Assay kit (Thermo Fisher Scientific). RNA libraries underwent polyA enrichment for library preparation and were sequenced on a NovaSeq 6000 platform (Illumina, 2 × 150 bp) to achieve a coverage of 6Gb. Transcriptomes were assembled following the pipeline detailed in MATEdb (Fernández et al. 2022). In brief, adapters were removed with Trimmomatic 0.38 (Bolger et al., 2014), and the de novo assembly for each different species was performed using Trinity 2.11 (Haas et al., 2013). We used Transdecoder 5.5 (https://github.com/TransDecoder/TransDecoder) to perform the prediction of the coding regions, and we removed all non-metazoan sequences identified through Diamond 2.0.8 (Buchfink et al., 2015) and filtered with Blobtools (Challis et al., 2020). From these filtered predicted proteomes we also extracted the longest isoforms. The differential gene expression analyses across different body parts were carried out following the Trinity pipeline using edgeR (Challis et al., 2020; Haas et al., 2013). Detailed parameters and scripts are provided in the Github repository. .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 14, 2025. ; https://doi.org/10.1101/2025.01.09.632098doi: bioRxiv preprint 6 Sample digestion and proteomics sequencing Samples (10 µg) were reduced with dithiothreitol (100 mM, 37 ºC, 60 min) and alkylated in the dark with iodoacetamide (5 µmol, 25 ºC, 20 min). The resulting protein extract was washed with 2M urea with 100 mM Tris-HCl and with 50 mM ammonium bicarbonate for subsequent digestion with endoproteinase LysC (1:10 w:w, 37ºC, o/n, Wako, cat # 129-02541) and trypsin (1:10 w:w, 37ºC, 8h, Promega cat # V5113) following Wiśniewski et al. FASP procedure. After digestion, the peptide mix was acidified with formic acid and desalted with a MicroSpin C18 column (The Nest Group, Inc) prior to LC-MS/MS analysis. Samples were analyzed using a Orbitrap Fusion Lumos mass spectrometer (Thermo Fisher Scientific, San Jose, CA, USA) coupled to an EASY-nLC 1200 (Thermo Fisher Scientific (Proxeon), Odense, Denmark). Peptides were loaded directly onto the analytical column and separated by reversed-phase chromatography using a 50-cm column with an inner diameter of 75 μm, packed with 2 μm C18 particles (Thermo Fisher Scientific, cat # ES903). Chromatographic gradients started at 95% buffer A and 5% buffer B with a flow rate of 300 nl/min and gradually increasing to 25% buffer B in 79 min and then to 40% buffer B in 11 min. After each analysis, the column was washed for 10 min with 100% buffer B. Buffer A: 0.1% formic acid in water. Buffer B: 0.1% formic acid in 80% acetonitrile. The mass spectrometer was operated in positive ionization mode with nanospray voltage set at 2.4 kV and source temperature at 305°C. The acquisition was performed in data-dependent acquisition (DDA) mode and full MS scans with 1 microscans at resolution of 120,000 were used over a mass range of m/z 350-1400 with detection in the Orbitrap mass analyzer. Auto gain control (AGC) was set to ‘standard’ and injection time to ‘auto’. In each cycle of data-dependent acquisition analysis, following each survey scan, the most intense ions above a threshold ion count of 10000 were selected for fragmentation. The number of selected precursor ions for fragmentation was determined by the “T op Speed” acquisition algorithm and a dynamic exclusion of 60 seconds. Fragment ion spectra were produced via high-energy collision dissociation (HCD) at normalized collision energy of 28% and they were acquired in the ion trap mass analyzer. AGC was set to 2E4, and an isolation window of 0.7 m/z and a maximum injection time of 12 ms were used. Digested bovine serum albumin (New england biolabs cat # P8108S) was analyzed between each sample to avoid sample carryover and to assure stability of the instrument and QCloud (Chiva et al., 2018; Olivella et al., 2021) has been used to control instrument longitudinal performance during the project. .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 14, 2025. ; https://doi.org/10.1101/2025.01.09.632098doi: bioRxiv preprint 7 Proteomics data analysis Acquired spectra were analyzed using the Proteome Discoverer software suite (v2.5, Thermo Fisher Scientific) and the Mascot search engine (v2.6, Matrix Science). The data were searched against the newly-generated de novo assembled transcriptomes of O. nungara and S. mediterranea (see above) plus a list of common contaminants and all the corresponding decoy entries. For peptide identification a precursor ion mass tolerance of 7 ppm was used for MS1 level, trypsin was chosen as enzyme, and up to three missed cleavages were allowed. The fragment ion mass tolerance was set to 0.5 Da for MS2 spectra. Oxidation of methionine and N-terminal protein acetylation were used as variable modifications whereas carbamidomethylation on cysteines was set as a fixed modification. False discovery rate (FDR) in peptide identification was set to a maximum of 1%. Peptide quantification data were retrieved from the “Precursor Ions Quantifier” node in Proteome Discoverer (v2.5) using a mass tolerance of 2 ppm for the peptide extracted ion current (XIC). The obtained values were then used to calculate the T op3 value for each protein, which is the average of the three most abundant peptides for each protein. The final list of identified proteins of each species was constructed by merging the results of the 12 samples (3 replicates for each body part, including the mucus fraction), but only proteins identified with high confidence levels and designated as the lead protein within a protein group (master protein) were considered. A protein group was defined as a collection of proteins containing at least one unique peptide. We finally explored the relationships among sample replicates through a Principal Component Analysis (PCA) using as input the normalized abundances of all the master proteins. All statistical analyses were performed using Proteome Discoverer software suite. For specifically characterizing the mucus of each species, we selected a subset of proteins that were differentially abundant in the mucus, in comparison to the body samples, that were used as the background signal. In the protein subset we included (i) proteins that were differentially expressed between the mucus and the body background signal, exhibiting a significant adjusted p-value and higher abundance in the mucus, and (ii) proteins that were identified in at least 2 out of 3 samples of the mucus, and were not present in non of the remaining samples (i.e., proteins that were detected as unique in the mucus fraction). .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 14, 2025. ; https://doi.org/10.1101/2025.01.09.632098doi: bioRxiv preprint 8 Prediction of toxins and functional annotation The predicted proteomes of each species were used to predict candidate toxins. The candidate toxins were predicted by using two different softwares: T oxify, which follows a deep learning approach (Cole & Brewer, 2019) and T oxinpred2, which follows an hybrid method that combines BLAST-based similarity, motif search and prediction models (Sharma et al., 2022), with a threshold of 0.6 and 0.95 respectively. Candidate toxins predicted by both softwares were merged to generate a list, which was used to detect the secreted toxins by using SignalP 6.0. T o reduce the amount of false positives during toxin prediction, we compared the list of putative secreted toxins with the proteins identified in the proteomic analyses, and built a final list including only those toxins that were also detected in the proteomic analyses. The putative toxins and the subset of mucus proteins were annotated through BLAST , using as

database the non-redundant protein sequences (nr). Since many putative toxins yielded no hits, we further explored putative protein function with FANTASIA (Martínez-Redondo et al. 2024), which leverages natural language models to infer Gene Ontology (GO) terms. Finally, we explored which GO terms were enriched in the mucus and toxins subsets against the proteome of each species inferred through proteomics using T opGO (Alexa & Rahnenfuhrer 2024) and REVIGO 1.8.1 (Supek et al., 2011), with the Gene Ontology database dated from january 2024. Graphic representations were done using RawGraph (Mauri et al. 2017) and (R version 4.2.2) Toxins and mucus evolution across the phylum Platyhelminthes We used the dataset and methods described in Benitez-Álvarez et al. 2025a and Benítez-Álvarez 2025b. Briefly, we combined the longest isoforms from the de novo transcriptomes of O. nungara and S. mediterranea (this study), with data from 31 other platyhelminthes species (Benítez-Álvarez et al. 2025b). We also included 4 outgroups belonging to the phyla Mollusca, Nemertea, Annelida, and Gastrotricha (Benítez-Álvarez et al. 2025b), obtaining a final dataset of 37 species. We inferred the Hierarchical Orthologous Groups (HOGs) with OMA v2.6 (Altenhoff et al., 2013), identifying the origin of HOGs at each diversification event in the phylogeny through pyHam (Train et al., 2019). For all the HOGs containing genes identified as toxins only for one species, but also sequences from the other species we explored potential shifts in the selection pressure using Pelican .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 14, 2025. ; https://doi.org/10.1101/2025.01.09.632098doi: bioRxiv preprint 9 (Duchemin et al., 2023). We coded the aquatic lineages as background and terrestrial ones as foreground characters. A p-value at gene level was estimated using using the Gene-wise Truncated Fisher’s method considering the best k=10 p-values in each alignment (https://gitlab.in2p3.fr/phoogle/pelican/-/wikis/Gene-level-predictions). Additionally, we applied the adjustment method Benjamini & Hochberg (BH) (Benjamini & Hochberg, 1995) based on the false discovery rate to obtain an adjusted p-value. Finally, we used the IHam tool (Train et al., 2019) to visually inspect the clustering of gene members across all HOGs that included toxins and mucus sequences. From this inspection, we identified the HOGs that showed expansions or contractions in the terrestrial lineage. Detailed

and data from the species included in these analyses can be found in the Github repository and Benitez-Álvarez et al. 2025a, Benítez-Álvarez et al. 2025b.

Toxin prediction and proteomic identification For each species, we constructed a high-quality reference transcriptome and predicted its proteome by merging the 9 RNA libraries (see Github repository). These libraries corresponded to different body regions, including the head, pharynx, and the rest of the body, for both O. nungara and S. mediterranea. Based on the predicted proteome, we identified putative toxins using T oxify (Cole & Brewer, 2019) and T oxinPred2 (Sharma et al., 2022). T oxinPred2 predicted a greater number of toxins than T oxify. Regardless of the software used, we detected a higher number of toxins in O. nungara than in S. mediterranea (Figure 1A, 1B, Github repository). For both species, we performed a protein extraction from the previously mentioned body regions as well as the mucus, and sequenced them through a bottom-up proteomics approach. We unambiguously identified 3584 proteins in S. mediterranea and 5302 proteins in O. nungara (Github repository). T o reduce false positives in our toxins predictions, we compared our final toxin candidates with the proteins identified in the proteomic analyses. Again, the number of toxins found was higher in O. nungara than in S. mediterranea (Figure 1C). Most of the toxins .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 14, 2025. ; https://doi.org/10.1101/2025.01.09.632098doi: bioRxiv preprint 10 that were identified in the proteomic data corresponded not only to different isoforms, but also to different genes. Figure 1. T oxin predictions and identifications. The Venn diagram for (1A) Obama nungara and (1B) Schmidtea mediterranea, showing the number of toxins predicted by T oxify (purple) and T oxinPred2 (blue). Green circles represent proteins predicted to be secreted from the shared list of toxins identified by both software tools. 1C. Barplot showing the proteins predicted and identified for Obama nungara (brown) and Schmidtea mediterranea (blue). “Identified toxins (i)” refers to the total amount of toxins identified in the proteomic analyses, while "Identified toxins (g)" indicates the number of distinct genes corresponding to these toxins. Toxin characterization BLAST results revealed that, although functional annotations were similar in both species, the most frequent toxins differed between the aquatic and terrestrial planarians (Figure 2A, T able S1). In O. nungara, 12 out of 57 toxins were annotated as lectins or scavenger receptors, whereas this function was not found in any of the toxins predicted in S. mediterranea. Two functions were also exclusively found in S. mediterranea, which included one conotoxin and one follistatin. On the other hand, both species contained toxins annotated as phospholipases, .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 14, 2025. ; https://doi.org/10.1101/2025.01.09.632098doi: bioRxiv preprint 11 insulin-like proteins, granulins, and venom allergens. Finally, for a considerable proportion of toxins, we obtained no match or found the function to be unknown. In addition, some annotations were either too general or unclear. Enrichments performed with FANTASIA annotations showed that the amount of enriched categories was substantially higher in O. nungara than in S. mediterranea, likely reflecting the greater amount of toxins predicted in the terrestrial planarian. In general, most of the terms that were enriched in S. mediterranea also were found enriched in O. nungara. Although in the biological process category three out of five categories were common between both species, O. nungara exhibited a notable number of categories related to the immune system, defense mechanisms, and hunting, such as “killing cells of another organism”, “antibacterial humoral response”, or “melanization defense response”. We also identified categories related to respiration and oxidative stress (indicated by pale green squares). In S. mediterranea we found an additional category related with the nervous system (“sleep”), and also “chemotaxis” (Figure 2B,C). Regarding the molecular function (Figure S1), both species shared two common enriched categories: "cytokine activity" and "ion channel regulation activity", in concordance with the findings from the biological process enrichments. S. mediterranea had only one additional enriched GO term, "toxin activity". In contrast, in O. nungara we detected up to 10 more enriched GO terms, including among others several binding-related categories, such as “galactose binding” or “antigen binding”, as well as immune system-related categories. T o fully characterize the function and expression of these putative toxins, we examined whether there was an increased expression across the various body parts analyzed in both transcriptomic (head and pharynx) and proteomic data (head, pharynx and mucus) when compared against the body background signal. We found only one toxin that was significantly differentially expressed in the pharynx of S. mediterranea in the transcriptome (identified as “SMED3_DN4834_c0_g1”). This same protein showed higher abundance also at protein level, although the difference was not significant (T able S2). BLAST annotations did not retrieve any hits for this protein, but FANTASIA annotations indicated a potential relationship with the neuropeptide signaling pathway and some kind of calcium binding activity (Github repository). .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 14, 2025. ; https://doi.org/10.1101/2025.01.09.632098doi: bioRxiv preprint 12 .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 14, 2025. ; https://doi.org/10.1101/2025.01.09.632098doi: bioRxiv preprint 13 Figure 2. Functional annotation of all the putative toxins identified in the proteomic results. 2A. Alluvial plot showing the summarized BLASTp annotation results, with O. nungara in brown and S. mediterranea in blue. Numbers indicate the amount of toxins included in a specific category. 2B, 2C. Results of the enrichment analyses for the Biological process category in O. nungara (B) and S. mediterranea (C). At proteomic level, we did not find any putative toxin which was significatively more abundant in the regions analyzed compared with the body background signal. However, for both species, we found some toxins showing relatively higher concentrations in the head, mucus or pharynx, reflected in the abundance ratio of each region compared to the body background signal (T able S2, T able S3). Nonetheless, the normalized abundance values across replicates exhibited a considerable variability, and in some cases there was also a lack of quantification in samples from specific body regions. Despite the observed trends, with our results we cannot confirm nor discard differential expression of the putative toxins across the body of S. mediterranea and O. nungara, and further studies will be needed to determine if true differential expression exists. Mucus protein characterization SDS page gel showed that the mucus fraction presented a different protein composition pattern in comparison with the other body fractions (Figure S2), and PCAs revealed a similar result (Figure S3). T o go further in this, we created mucus protein subsets for each species containing (i) proteins that were significantly differentially expressed in mucus, and (ii) proteins that were only found in mucus samples. While our strict criteria aim to exclude proteins that are also abundant in the body of these animals - and therefore potential contaminants produced during the extraction of the mucus fraction - our subsets may not capture the full mucus composition. Instead, our subsets primarily include the most abundant and specific proteins in O. nungara and S. mediterranea mucus. Enrichment analyses of both species retrieved terms related to adhesion processes, as well as terms related with development, growth, regeneration and immune system (Figure 3A, B). In both species, we annotated through blast proteins that were related with cell structure (such as actin, myosin and other related proteins), adhesion (mucins, von Willebrand factor), immune system and defense (lectins, proteases) and proteins related with metabolism and nucleic acids (such as ribosomal proteins). Full results are detailed in the Github repository and Figure S4. .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 14, 2025. ; https://doi.org/10.1101/2025.01.09.632098doi: bioRxiv preprint 14 .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 14, 2025. ; https://doi.org/10.1101/2025.01.09.632098doi: bioRxiv preprint 15 Figure 3. Results of the enrichment analyses of mucus for O. nungara in the Biological process (A) and Molecular function (B) and S. mediterranea in the Biological process (C) and Molecular function (D). In addition to these conserved general patterns, in O. nungara the enrichments revealed GO terms that coincided with the most common toxin annotations, including 'scavenger receptor activity,' 'carbohydrate binding,' and 'cytokine activity', as well as a higher amount of GOterms related with immune system and defense, such as “defense response to other organism”, “immune effector process” or “protein destabilization”, which could also be related with the lectins identified in the mucus fraction. On the other hand, enrichment analyses of S. mediterranea retrieved a higher number of GO terms that were related with morphogenesis and development (Figure 3). Evolutionary origin of toxins and mucus proteins T o better understand the role of mucus and toxins during the transition from aquatic to terrestrial environments, we explored the evolutionary origin of toxin- and mucus-protein orthologous groups across the platyhelminth phylogeny. The general results of gene assignment to hierarchical orthologous groups (HOGs) are outlined in T able 1. Table 1: Number of genes of each species (O. nungara and S. mediterranea) and dataset (toxins and mucus) and its assignment to hierarchical orthologous groups (HOGs). Genes that were species-specific were not assigned to any HOG. Finally, the number of HOGs that were gained after the diversification of terrestrial planarians (O. nungara results) and freshwater planarians (S. mediterranea results) are detailed. A high proportion of HOGs including predicted toxins and mucus proteins originated prior to the diversification of terrestrial planarians, although mucus and toxin orthologous groups showed different patterns (Figure 4). A large proportion of toxin orthologous groups was gained in the branch leading to Tricladida (which includes marine, freshwater and terrestrial species) and .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 14, 2025. ; https://doi.org/10.1101/2025.01.09.632098doi: bioRxiv preprint 16 Continenticola (Github repository). On the other hand, most mucus orthologous groups of O. nungara and S. mediterranea had an ancient origin, with approximately half of the total HOGs emerging before the diversification of platyhelminthes and therefore being shared with other animal phyla. Interestingly, despite the majority of the orthologous groups originated before the terrestrial planarians diversification, most of the HOGs were not shared between species, indicating that O. nungara and S. mediterranea may rely on different orthologous groups in their toxin and mucus composition. However, while mucus composition was assessed by direct identification in the selected subset of proteins, and therefore reflects the most abundant and specific proteins in the mucus of each species, the toxin identification was based on in silico predictions. In order to further clarify this, we conducted a deeper investigation into the sequences within the HOGs containing toxin genes to determine whether the HOGs with predicted toxins in one species (O. nungara or S. mediterranea) also included orthologous sequences from the other species. Five out of 33 HOGs contained toxins that were both predicted and identified in both species. However, of the remaining 28 HOGs, we identified 15 that included toxins predicted and identified in only one species (O. nungara or S. mediterranea), yet also contained sequences from the other species, which in most cases were not predicted as toxins by T oxify (see T ableS4 for details). Results obtained with Pelican showed that 10 out of these 15 HOGs exhibited a significant shift in the selection at the amino acid level between aquatic and terrestrial planarians (T able S5). We also detected two HOGs containing mucus genes that specifically originated in the ancestor of terrestrial planarians. The first HOG contained up to 10 O. nungara genes, including at least one identified in the mucus subset (“ONUN2_DN347_c0_g1”). GO terms inferred with FANTASIA were again associated with immune system response, retrieving some interesting terms as “defense response to bacterium” or “innate immune response” (Github repository). The other HOG gained in the ancestor of terrestrial planarians contained only one O. nungara gene. In this case, annotations included functions as “lipid binding” and also functions related with entry to the cell (Github repository). .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 14, 2025. ; https://doi.org/10.1101/2025.01.09.632098doi: bioRxiv preprint 17 Figure 4. Evolutionary origins of toxins and mucus proteins across Platyhelminthes. Outgroups are highlighted in red, and the habitats of the different species are indicated by squares: dark blue for marine, light blue for freshwater, and brown for terrestrial. Gained HOGs are represented by circles at the nodes, with the colors and numbers inside indicating HOGs that contain sequences from S. mediterranea (blue), O. nungara (brown), or both species (grey). Finally, we used pyHam to visualize the toxin and mucus HOGs and to identify potential expansions and contractions that could be associated with the adaptation of the planarians to the terrestrial environment (Github repository). We did not identify any mucus orthologous groups that presented expansions or contractions in terrestrial planarians. However, among the HOGs containing toxin genes, one showed an expansion and another a contraction in terrestrial planarians. The first orthologous group originated in the Continenticola node and was expanded in terrestrial planarians, with 19 gene copies in O. nungara and up to 26 copies in the case of the terrestrial flatworm Bipalium kewense (Figure 4A). At least 6 sequences of O. nungara were predicted as secreted toxins by the two software tools and identified in the proteomic results. Interestingly, and in concordance with previous results, all the BLAST hits of these sequences retrieved hits with lectins and scavenger receptors (T able S1). .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 14, 2025. ; https://doi.org/10.1101/2025.01.09.632098doi: bioRxiv preprint 18 Figure 5. HOGs expanded or contracted in terrestrial flatworms. The habitats of the different species are represented by colored squares: dark blue for marine, light blue for freshwater, and brown for terrestrial. The number of gray squares aligned with each species name corresponds to the number of gene copies that the species has in the specified HOG. 5A. HOG (HOG00036620_64) is expanded in terrestrial flatworms. 5B. HOG (HOG00044669_65) is contracted in terrestrial flatworms. The second HOG was gained in the branch leading to Tricladida and contracted in the terrestrial planarians (Figure 4B). Two out of the 6 sequences of S. mediterranea included in this HOG were predicted as secreted toxins, with one being also found in the proteomic results (T able S4). The BLAST annotation indicated that this toxin was classified as a venom allergen. Summarized

for all the proteins classified as toxins and mucus components, including the HOGs to which they belong and a general category for the BLAST hits can be found in T able S1 and S6.

Our study provides new insights into the composition and evolutionary trajectories of toxins and mucus proteins in flatworms and their potential role in flatworm terrestrialization. By comparing the terrestrial flatworm O. nungara and the freshwater species S. mediterranea, we aimed to elucidate how the adaptive gene repertoires differ between aquatic and terrestrial species, particularly in relation to toxins and mucus proteins. .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 14, 2025. ; https://doi.org/10.1101/2025.01.09.632098doi: bioRxiv preprint 19 Evolutionary origin of the toxin repertoire in aquatic and terrestrial flatworms We found that a substantial proportion of the predicted toxin families in both species arose before the diversification of the flatworm families examined here (Dugesiidae and Geoplanidae). This observation is consistent with recent findings indicating that many genes differentially expressed under abiotic stress conditions also originated prior to their diversification, specifically in the branches leading to Tricladida and Continenticola (Benítez-Álvarez et al. 2025a). T aken together, these results suggest that the genetic framework for toxin production was largely established earlier in flatworm evolution than previously thought, providing a reservoir of ancestral genes that could be co-opted or fine-tuned in subsequent lineages. This ancient gene repertoire may have played a pivotal role in enabling flatworms to adapt to diverse ecological niches, including both freshwater and terrestrial habitats, by conferring pre-existing molecular mechanisms for defense, prey capture, or other adaptive functions. Despite a common evolutionary origin mostly in the branches leading to Tricladida and Continenticola, most of the toxin-containing HOGs were not shared between species. However, a substantial proportion of HOGs that included toxins predicted for only one species also contained sequences from the other. These results can be explained by two potential hypotheses: (i) although our protocol aimed to reduce the proportion of false positives by combining the results of two different prediction software tools, this approach may have increased the proportion of false negatives. Consequently, some sequences not predicted as toxins could indeed exhibit toxin activity; (ii) species may rely on distinct HOGs for their toxin repertoire, and genes belonging to the same HOG may have different functions depending on species-specific or environmental factors. In line with this, selection analyses revealed that approximately 66% of these HOGs exhibited a significant shift in selection between terrestrial and aquatic species, supporting the possibility that, at least in some cases, the latter hypothesis could be true. Hence, these findings suggest that differences in selection pressures may shape the specialized toxin repertoires of O. nungara and S. mediterranea, underscoring the complexity of molecular adaptations that emerge as flatworms occupy different environments. However, our data do not allow us to clearly distinguish between these two explanations, and further experimental studies will be required to definitively determine whether these proteins have indeed toxin activity. .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 14, 2025. ; https://doi.org/10.1101/2025.01.09.632098doi: bioRxiv preprint 20 A larger toxin gene repertoire in terrestrial flatworms and its potential role in land colonization We detected more putative toxin-related genes in O. nungara than in S. mediterranea, also reflected in the greater enrichment of significant categories in the terrestrial species. Specifically, the higher number of enriched GO terms related to immune response and inflammation in O. nungara suggests that these functions could have played a significant role in the adaptation to terrestrial life, where exposure to new pathogens, predators, and environmental stressors may have produced a more complex interaction between toxins and terrestrial ecosystems. Although many toxin protein families were detected in both species, we observed some characteristic patterns which differed between the terrestrial and aquatic flatworms. For example, in O. nungara, at least 12 predicted toxins were identified as lectin-like proteins. Lectins are proteins that can bind to carbohydrates, and previous research suggests that these proteins may have a key role in the pathogen recognition and the immune system of invertebrates (Hanington et al., 2010; Pees et al., 2016) and have been also identified in venomous animals from other phyla including cnidarians, mollusks and nemertea (Delgado et al., 2022; Verdes et al., 2022; von Reumont et al., 2020). The 12 lectins identified in O. nungara belonged to two different HOGs: the first originated within the terrestrial planarian lineage, while the second was gained in the branch leading to Continenticola and underwent a significant expansion in terrestrial species, (including O. nungara). This expansion aligns with previous observations in mollusks, where several C-type lectin families are expanded in terrestrial species and have been suggested to serve defensive roles in terrestrial environments (Aristide & Fernández, 2023) In addition, in planarians, several novel lectins have been also reported in freshwater flatworms (Gao et al., 2017; Shagin et al., 2002). One of these novel lectins, described in the freshwater planarian Dugesia japonica, has been associated with innate immune response against pathogens and with the first steps of regeneration and wound healing (Gao et al., 2017), further supporting its potential role in the adaptation to new environments and pathogens and its involvement in regeneration. We also identified HOGs containing group 1 (secreted) Venom Allergen-like proteins (VAL), which are known to be widely distributed across metazoans, including flatworms (Chalmers et al., 2008). One of these families has an ancient origin and is shared across multiple phyla. The other two emerged at the branch leading to Tricladida. One of these HOGs experienced a .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 14, 2025. ; https://doi.org/10.1101/2025.01.09.632098doi: bioRxiv preprint 21 contraction in terrestrial species, while the aquatic planarians, including S. mediterranea, showed a larger copy number. While the role of VAL proteins has been linked to parasitism and their function has been mainly studied in parasitic species (Chalmers & Hoffmann, 2012; Yoshino et al., 2014), the most expanded VAL orthologous group of this phylum has been precisely described in freshwater free-living flatworms (Chalmers & Hoffmann, 2012,Sipley 2019), in concordance with our results. Altogether, it seems that VAL proteins may have a relevant role in the biology of aquatic free-living planarians, with a potential shift in functional roles as flatworms adapted to terrestrial life. Functional convergence of toxins with different evolutionary origin in aquatic and terrestrial planarians Beyond the functional differences between the toxin gene repertoire in aquatic and terrestrial flatworms discussed above, our analyses revealed that several functions identified in O. nungara and S. mediterranea were shared between both species, despite being encoded by different HOGs. Common enriched categories included cytokine activity, ion channel regulation and categories related with the neuromuscular system, suggesting that, as observed in other species and phyla, some flatworm toxins may trigger inflammatory responses, target ion channels and/or disrupt the nervous system (Brodie, 2009; Zhang, 2015). In both species, we identified several toxins belonging to the CAP superfamily (Gibbs et al., 2008), including one conotoxin-like protein in S. mediterranea and Venom Allergen-like (VAL) proteins in both species, the latter of which have previously been described in different flatworm lineages (Chalmers & Hoffmann, 2012). CAP-domain proteins have been convergently recruited in venoms of many taxa, ranging from mammals to annelids or arthropods (Jenner et al., 2019; Verdes et al., 2018; von Reumont et al., 2020; Zhang et al., 2022). These proteins exhibit several biological functions, including disruption of ion channels, proinflammatory activity and the induction of allergic responses (Fry et al., 2009; Zhang et al., 2022) We also identified phospholipases, which are frequently recruited in venoms of other taxa (Lynch, 2007; Undheim et al., 2014; Verdes et al., 2018) and induce multiple effects on the targets, including myotoxicity, neurotoxicity, inflammation and changes in the reactive oxygen species production (El-Benna et al., 2021; Fry et al., 2009). Finally, insulin-like proteins were also predicted and detected in both species. This protein group is also found in venoms across .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 14, 2025. ; https://doi.org/10.1101/2025.01.09.632098doi: bioRxiv preprint 22 various taxa, including nemerteans (Verdes et al., 2022; von Reumont et al., 2020) and cone-snails (Guo et al., 2024), where it is used to facilitate prey capture. Ancient origins of mucus proteins encoded by distinct HOGs in aquatic and terrestrial flatworms In contrast to toxin genes, the mucus protein repertoire appears to have a more ancient evolutionary origin, with many mucus-related HOGs also present across other animal phyla. Despite this shared ancestry, O. nungara and S. mediterranea seem to rely on distinct sets of HOGs for mucus composition, suggesting that each species draws upon different genetic resources to produce its most abundant mucus proteins. For example, although we found proteins and GO terms related to immune defense in both species, the proportion of those in O. nungara mucus was higher, which could be due to an enhanced protective function of the mucus in terrestrial habitats. In both species, we found proteases and protease inhibitors (such as serpin proteins), including metalloproteases in the S. mediterranea mucus. Both proteases and proteases inhibitors have been repeatedly reported in the mucus of different species (Cerullo et al., 2023; Pales Espinosa et al., 2016), including S. mediterranea (Bocchinfuso et al., 2012), and it has been suggested that these proteins contribute to the immunity. We also detected lectins in the mucus of both species, aligning with previous studies in planarians (Shagin et al., 2002; Zayas et al., 2010) and other invertebrates, where lectins have been suggested to contribute to the mucus adhesion and microbial protection (Liegertová et al., 2022; Smith et al., 2021). Additionally, proteins involved in adhesion, such as mucins, spondin and fibrillin-like proteins, were present in the mucus of both species, underscoring their functional significance in aquatic and terrestrial habitats (Lang et al., 2007; Liegertová & Malý, 2023). Finally, in both species, we also detected proteins of intracellular origin, including those related with cell structure and cytoskeleton, such as actin or myosin, and ribosomal proteins. These results are highly consistent with previous mucus characterizations in species from other phyla (Lopes et al., 2024; Pales Espinosa et al., 2016; Schwaner et al., 2024) and have been also previously described in the S. mediterranea mucus (Bocchinfuso et al., 2012). .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 14, 2025. ; https://doi.org/10.1101/2025.01.09.632098doi: bioRxiv preprint 23

Overall, we successfully predicted and identified putative translated toxins in O. nungara and S. mediterranea, and characterized the most abundant proteins in their mucus secretions. Our findings revealed both similarities and contrasting dynamics in the toxins and mucus proteins of these two flatworm species, which we discuss in the context of flatworm terrestrialization. We predicted a larger number of toxins in O. nungara, also reflected in a higher number of enriched GO terms related to immune response and inflammation. In both species, most orthologous groups encoding its mucus and toxin repertoire originated before the major split of terrestrial flatworms, showcasing how genetic elements that originated prior to the colonization of terrestrial environments were co-opted to deal with threats related to life on land. However, while most toxin-related HOGs were gained in the branches leading to Tricladida and Continenticola, a significant proportion of mucus-related HOGs seem to have a more ancient evolutionary origin. Although the toxin families in O. nungara and S. mediterranea appear to share similar functions, terrestrial flatworms exhibit an expansion of a lectin orthologous group, suggesting that some lectin families may have a role in facilitating adaptation to land. These findings underscore the pivotal role of ancient genetic frameworks, supplemented by lineage-specific expansions like those of lectin families, in shaping the evolutionary success of flatworms on land. Altogether, these results highlight the importance of examining both ancient genomic elements and lineage-specific expansions to better understand how flatworms—and, by extension, other taxa—have successfully adapted to life on land. Supplementary information Figure S1: Results of the enrichment analyses for the Molecular function category in O. nungara and S. mediterranea. Figure S2: SDS Page results for the different body fractions sampled in O. nungara and S. mediterranea. 10 micrograms of protein per lane. A. SDS Page coomassie stained. Lane 1, PageRuler™ Unstained Protein Ladder (ThermoScientific), Lanes 2 to 5 proteins from O. nungara. Head shown in lane 2, body in lane 3, pharynx in lane 4, and mucus in lane 5. Lanes 6 to 10 proteins from S. mediterranea. Head in lane 6, body in lane 7, pharynx in lane 8, lane 9 and 10 shows mucus extracted by 4 minutes or 8 minutes incubation in NAC respectively, lane 10. Note that the mucus fraction in lane 5 is from a failed attempt of collecting mucus directly by scrapping the individuals, so we repeated it with NAC as for S. mediterranea and it is presented in B. B. SDS Page of O. nungara proteins from mucus and body. Lane 1 shows the protein ladder. Lanes 2 and 3 show the mucus extracted by 4 minutes or 8 minutes incubation in NAC respectively. Lane 4 shows protein from the body. Figure S3: PCAs showing the distinct proteomic compositions of the body fractions (body, head, mucus, and pharynx) for O. nungara (left) and S. mediterranea (right). Notably, the mucus fraction is clearly separated from the other samples in both species. .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 14, 2025. ; https://doi.org/10.1101/2025.01.09.632098doi: bioRxiv preprint 24 Figure S4: Barplot representing the summarized Blastp annotation results for the mucus fraction of O. nungara and S. mediterranea. Table S1: Integrated results of all the putative toxins identified in O. nungara and S. mediterranea. Includes HOG information and blast summarized annotations for each putative toxin. Table S2: Proteomics results for the toxins detected in S. mediterranea, including the normalized abundance of each sample. Table S3: Proteomics results for the toxins detected in O. nungara, including the normalized abundance of each sample. Table S4: Predicted toxins by T oxify and T oxinPred2 for the sequences included in HOGs that already contained at least one predicted and detected toxin in S. mediterranea or O. nungara. Table S5: Pelican results and statistics. Table S6: Integrated results of all the proteins identified in the mucus subset of O. nungara and S. mediterranea. Includes HOG information and blast summarized annotations for each mucus protein. Supplementary Materials include all intermediate files, final files, scripts and tables, that can be found at the Github repository. Data availability Full link to the Github repository: https://github.com/MetazoaPhylogenomicsLab/Garcia_Vernet_et_al_2024_Flatworm_toxin_evol ution

R.G-V was funded through a Margarita Salas grant, funded by the Spanish Ministry of Universities and the European Union Next Generation EU/PRTR. RF acknowledges support from the following sources of funding: the European Research Council (this project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement no. 948281), the OSCARS project (funding from the European Commission’s Horizon Europe Research and Innovation programme under grant agreement no. 101129751) and the Secretaria d’Universitats i Recerca del Departament d’Economia i Coneixement de la Generalitat de Catalunya (AGAUR 2021-SGR00420 and 2021-SGR01225). We acknowledge support of the Spanish Ministry of Science and Innovation through the Centro de Excelencia Severo Ochoa (CEX2020-001049-S grant funded by MCIN/AEI/10.13039/501100011033). The CRG/UPF Proteomics Unit is part of the Spanish Infrastructure for Omics T echnologies (ICTS OmicsT ech). We also thank Centro de Supercomputación de Galicia (CESGA) and the HPC Drago from the Centro Superior de Investigaciones Científicas for access to computer resources. .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 14, 2025. ; https://doi.org/10.1101/2025.01.09.632098doi: bioRxiv preprint 25

Alexa, A., & Rahnenführer, J. (2009). Gene set enrichment analysis with topGO. Bioconductor Improv, 27, 1-26. Altenhoff, A. M., Gil, M., Gonnet, G. H., & Dessimoz, C. (2013). Inferring hierarchical orthologous groups from orthologous gene pairs. PloS One, 8(1), e53786. Aristide, L., & Fernández, R. (2023). Genomic Insights into Mollusk T errestrialization: Parallel and Convergent Gene Family Expansions as Key Facilitators in Out-of-the-Sea Transitions. Genome Biology and Evolution, 15(10). https://doi.org/10.1093/gbe/evad176 Benítez-Álvarez, L., García-Vernet, R., Eleftheriadi, K., Vargas-Chávez, C., Escudero, N., Salces-Ortiz, J., Rojo, I., Fernández-Álvarez, F .A., Chiva, C., Sabidó, E., Mateos, E., Carbayo, F . & Fernández, R. (2025a). Genomic exaptation and regulatory landscape shifts as key mechanisms enabling flatworm terrestrialization. bioRxiv. Benítez-Álvarez, L., Escudero, N., Salces-Ortiz, J., Rojo, I., Fernández-Álvarez, F .A., Mateos, E., Carbayo, F ., & Fernández, R. (2025b). Disentangling the evolutionary history of terrestrial planarians through phylogenomics. bioRxiv. https://doi.org/10.1101/2025.01.07.631701 Benjamini, Y ., & Hochberg, Y . (1995). Controlling the false discovery rate: a practical and powerful approach to multiple testing. Journal of the Royal Statistical Society Series B-Methodological, 57(1), 289–300. Bocchinfuso, D. G., T aylor, P ., Ross, E., Ignatchenko, A., Ignatchenko, V., Kislinger, T ., Pearson, B. J., & Moran, M. F . (2012). Proteomic profiling of the planarian Schmidtea mediterranea and its mucous reveals similarities with human secretions and those predicted for parasitic flatworms. Molecular & Cellular Proteomics : MCP, 11(9), 681–691. Bolger, A. M., Lohse, M., & Usadel, B. (2014). Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics , 30(15), 2114–2120. Brodie, E. D., 3rd. (2009). T oxins and venoms. Current Biology : CB, 19(20), R931–R935. Buchfink, B., Xie, C., & Huson, D. H. (2015). Fast and sensitive protein alignment using DIAMOND. Nature Methods, 12(1), 59–60. Cardoso, J. C., Carbayo, F ., & Gonzaga, M. O. (2023). Land planarians (Platyhelminthes) also prey on web-building spiders. Neotropical Biology and Conservation, 18(3), 157-162. Casewell, N. R., Wüster, W., Vonk, F . J., Harrison, R. A., & Fry, B. G. (2013). Complex cocktails: the evolutionary novelty of venoms. Trends in Ecology & Evolution, 28(4), 219–229. Cerullo, A. R., McDermott, M. B., Pepi, L. E., Liu, Z.-L., Barry, D., Zhang, S., Yang, X., Chen, X., Azadi, P ., Holford, M., & Braunschweig, A. B. (2023). Comparative mucomic analysis of three functionally distinct Cornu aspersum Secretions. Nature Communications, 14(1), 5361. Challis, R., Richards, E., Rajan, J., Cochrane, G., & Blaxter, M. (2020). BlobT oolKit - Interactive Quality Assessment of Genome Assemblies. G3 , 10(4), 1361–1374. Chalmers, I. W., & Hoffmann, K. F . (2012). Platyhelminth Venom Allergen-Like (VAL) proteins: revealing structural diversity, class-specific features and biological associations across the phylum. Parasitology, 139(10), 1231–1245. Chalmers, I. W., McArdle, A. J., Coulson, R. M., Wagner, M. A., Schmid, R., Hirai, H., & Hoffmann, K. F . (2008). Developmentally regulated expression, alternative splicing and distinct sub-groupings in members of the Schistosoma mansoni venom allergen-like (SmVAL) gene family. BMC Genomics, 9, 89. Chiva, C., Olivella, R., Borràs, E., Espadas, G., Pastor, O., Solé, A., & Sabidó, E. (2018). QCloud: A cloud-based quality control system for mass spectrometry-based proteomics laboratories. PloS One, 13(1), e0189209. Cole, T . J., & Brewer, M. S. (2019). TOXIFY: a deep learning approach to classify animal venom .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 14, 2025. ; https://doi.org/10.1101/2025.01.09.632098doi: bioRxiv preprint 26 proteins. PeerJ, 7, e7200. Cuevas-Caballé, C., Riutort, M., & Álvarez-Presas, M. (2019). Diet assessment of two land planarian species using high-throughput sequencing data. Scientific Reports, 9(1), 8679. Delgado, A., Benedict, C., Macrander, J., & Daly, M. (2022). Never, Ever Make an Enemy… Out of an Anemone: Transcriptomic Comparison of Clownfish Hosting Sea Anemone Venoms. Marine Drugs, 20(12). https://doi.org/10.3390/md20120730 Duchemin, L., Lanore, V., Veber, P ., & Boussau, B. (2023). Evaluation of methods to detect shifts in directional selection at the genome scale. Molecular Biology and Evolution, 40(2). https://doi.org/10.1093/molbev/msac247 El-Benna, J., Hurtado-Nedelec, M., Gougerot-Pocidalo, M.-A., & Dang, P . M.-C. (2021). Effects of venoms on neutrophil respiratory burst: a major inflammatory function. The Journal of Venomous Animals and Toxins Including Tropical Diseases, 27, e20200179. Fry, B. G., Roelants, K., Champagne, D. E., Scheib, H., Tyndall, J. D. A., King, G. F ., Nevalainen, T . J., Norman, J. A., Lewis, R. J., Norton, R. S., Renjifo, C., & de la Vega, R. C. R. (2009). The toxicogenomic multiverse: convergent recruitment of proteins into animal venoms. Annual Review of Genomics and Human Genetics, 10, 483–511. Gao, L., Han, Y ., Deng, H., Hu, W., Zhen, H., Li, N., Qin, N., Yan, M., Wu, W., Liu, B., Zhao, B., & Pang, Q. (2017). The role of a novel C-type lectin-like protein from planarian in innate immunity and regeneration. Developmental and Comparative Immunology, 67, 413–426. Gibbs, G. M., Roelants, K., & O’Bryan, M. K. (2008). The CAP superfamily: cysteine-rich secretory proteins, antigen 5, and pathogenesis-related 1 proteins--roles in reproduction, cancer, and immune defense. Endocrine Reviews, 29(7), 865–897. Giribet, G., & Edgecombe, G. D. (2020). The Invertebrate Tree of Life. Princeton University Press. Guo, Q., Huang, M., Li, M., Chen, J., Cheng, S., Ma, L., & Gao, B. (2024). Diversity and Evolutionary Analysis of Venom Insulin Derived from Cone Snails. Toxins, 16(1). https://doi.org/10.3390/toxins16010034 Haas, B. J., Papanicolaou, A., Yassour, M., Grabherr, M., Blood, P . D., Bowden, J., Couger, M. B., Eccles, D., Li, B., Lieber, M., MacManes, M. D., Ott, M., Orvis, J., Pochet, N., Strozzi, F ., Weeks, N., Westerman, R., William, T ., Dewey, C. N., … Regev, A. (2013). De novo transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis. Nature Protocols, 8(8), 1494–1512. Hanington, P . C., Forys, M. A., Dragoo, J. W., Zhang, S.-M., Adema, C. M., & Loker, E. S. (2010). Role for a somatically diversified lectin in resistance of an invertebrate to parasite infection. Proceedings of the National Academy of Sciences of the United States of America, 107(49), 21087–21092. Jenner, R. A., von Reumont, B. M., Campbell, L. I., & Undheim, E. A. B. (2019). Parallel Evolution of Complex Centipede Venoms Revealed by Comparative Proteotranscriptomic Analyses. Molecular Biology and Evolution, 36(12), 2748–2763. Kloc, M., Halasa, M., Kubiak, J. Z., & Ghobrial, R. M. (2024). Invertebrate Immunity, Natural Transplantation Immunity, Somatic and Germ Cell Parasitism, and Transposon Defense. International Journal of Molecular Sciences, 25(2). https://doi.org/10.3390/ijms25021072 Lang, T ., Hansson, G. C., & Samuelsson, T . (2007). Gel-forming mucins appeared early in metazoan evolution. Proceedings of the National Academy of Sciences of the United States of America, 104(41), 16209–16214. Liegertová, M., & Malý, J. (2023). Gastropod Mucus: Interdisciplinary Perspectives on Biological Activities, Applications, and Strategic Priorities. ACS Biomaterials Science & Engineering, 9(10), 5567–5579. Liegertová, M., Semerádtová, A., Kocholatá, M., Průšová, M., Němcová, L., Štofik, M., Kříženecká, S., Malý, J., & Janoušková, O. (2022). Mucus-derived exosome-like vesicles from the Spanish slug (Arion vulgaris): taking advantage of invasive pest species in .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 14, 2025. ; https://doi.org/10.1101/2025.01.09.632098doi: bioRxiv preprint 27 biotechnology. Scientific Reports, 12(1), 21768. Little, C. (1983). The Colonisation of Land: Origins and Adaptations of Terrestrial Animals. Cambridge University Press. Lopes, D., Aveiro, S. S., Cruz, S., Cartaxana, P ., & Domingues, P . (2024). Proteomic analysis of the mucus of the photosynthetic sea slug Elysia crispata. Journal of Proteomics, 294, 105087. Lynch, V. J. (2007). Inventing an arsenal: adaptive evolution and neofunctionalization of snake venom phospholipase A2 genes. BMC Evolutionary Biology, 7, 2. Martínez-Redondo, G. I., Simón Guerrero, C., Aristide, L., Balart-García, P ., T onzo, V., & Fernández, R. (2023). Parallel duplication and loss of aquaporin-coding genes during the “out of the sea” transition as potential key drivers of animal terrestrialization. Molecular Ecology, 32(8), 2022–2040. Martínez-Redondo, G.I., Barrios-Núñez, I., Vázquez-Valls, M., Rojas, A. M., & Fernández, R. (2024). Illuminating the functional landscape of the dark proteome across the Animal Tree of Life through natural language processing models. bioRxiv. doi: https://doi.org/10.1101/2024.02.28.582465 Mauri, M., Elli, T ., Caviglia, G., Uboldi, G., & Azzi, M. (2017). RAWGraphs: A Visualisation Platform to Create Open Outputs. In Proceedings of the 12th Biannual Conference on Italian SIGCHI Chapter (p. 28:1–28:5). New York, NY , USA: ACM. https://doi.org/10.1145/3125571.3125585 Meyer, A., Schloissnig, S., Franchini, P ., Du, K., Woltering, J. M., Irisarri, I., Wong, W. Y ., Nowoshilow, S., Kneitz, S., Kawaguchi, A., Fabrizius, A., Xiong, P ., Dechaud, C., Spaink, H. P ., Volff, J.-N., Simakov, O., Burmester, T ., T anaka, E. M., & Schartl, M. (2021). Giant lungfish genome elucidates the conquest of land by vertebrates. Nature, 590(7845), 284–289. Miyazawa, K., Jeon, J. K., Noguchi, T ., Ito, K., & Hashimoto, K. (1987). Distribution of tetrodotoxin in the tissues of the flatworm Planocera multitentaculata (Platyhelminthes). Toxicon: Official Journal of the International Society on Toxinology, 25(9), 975–980. Olivella, R., Chiva, C., Serret, M., Mancera, D., Cozzuto, L., Hermoso, A., Borràs, E., Espadas, G., Morales, J., Pastor, O., Solé, A., Ponomarenko, J., & Sabidó, E. (2021). QCloud2: An Improved Cloud-based Quality-Control System for Mass-Spectrometry-based Proteomics Laboratories. Journal of Proteome Research, 20(4), 2010–2013. Pales Espinosa, E., Koller, A., & Allam, B. (2016). Proteomic characterization of mucosal secretions in the eastern oyster, Crassostrea virginica. Journal of Proteomics, 132, 63–76. Pees, B., Yang, W., Zárate-Potes, A., Schulenburg, H., & Dierking, K. (2016). High Innate Immune Specificity through Diversified C-Type Lectin-Like Domain Proteins in Invertebrates. Journal of Innate Immunity, 8(2), 129–142. Riutort, M., Álvarez-Presas, M., Lázaro, E., Solà, E., & Paps, J. (2012). Evolutionary history of the Tricladida and the Platyhelminthes: an up-to-date phylogenetic and systematic account. The International Journal of Developmental Biology, 56(1-3), 5–17. Schendel, V., Rash, L. D., Jenner, R. A., & Undheim, E. A. B. (2019). The Diversity of Venom: The Importance of Behavior and Venom System Morphology in Understanding Its Ecology and Evolution. Toxins, 11(11). https://doi.org/10.3390/toxins11110666 Schwaner, C., Barbosa, M., Haley, J., Pales Espinosa, E., & Allam, B. (2024). Transcriptomics, proteomics, and physiological assays reveal immunosuppression in the eastern oyster Crassostrea virginica exposed to acidification stress. Fish & Shellfish Immunology, 146, 109366. Shagin, D. A., Barsova, E. V., Bogdanova, E., Britanova, O. V., Gurskaya, N., Lukyanov, K. A., Matz, M. V., Punkova, N. I., Usman, N. Y ., Kopantzev, E. P ., Salo, E., & Lukyanov, S. A. (2002). Identification and characterization of a new family of C-type lectin-like genes from planaria Girardia tigrina. Glycobiology, 12(8), 463–472. .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 14, 2025. ; https://doi.org/10.1101/2025.01.09.632098doi: bioRxiv preprint 28 Sharma, N., Naorem, L. D., Jain, S., & Raghava, G. P . S. (2022). T oxinPred2: an improved

for predicting toxicity of proteins. Briefings in Bioinformatics, 23(5). https://doi.org/10.1093/bib/bbac174 Simões, A. E. S., Pereira, D. M., Amaral, J. D., Nunes, A. F ., Gomes, S. E., Rodrigues, P . M., Lo, A. C., D’Hooge, R., Steer, C. J., Thibodeau, S. N., Borralho, P . M., & Rodrigues, C. M. P . (2013). Efficient recovery of proteins from multiple source samples after TRIzol(®) or TRIzol(®)LS RNA extraction and long-term storage. BMC Genomics, 14, 181. Sipley, B. N. (2019). Venom allergen-like protein diversification in flatworms (Master's thesis, Auburn University). Sluys, R. (1999). Biodiversity and Conservation, 8(12), 1663–1681. Sluys, R. (2019). The evolutionary terrestrialization of planarian flatworms (Platyhelminthes, Tricladida, Geoplanidae): a review and research programme. Zoosystematics and Evolution, 95(2), 543-556. Smith, A. M., Huynh, P ., Griffin, S., Baughn, M., & Monka, P . (2021). Strong, Non-specific Adhesion Using C-Lectin Heterotrimers in a Molluscan Defensive Secretion. Integrative and Comparative Biology, 61(4), 1440–1449. Stokes, A. N., Ducey, P . K., Neuman-Lee, L., Hanifin, C. T ., French, S. S., Pfrender, M. E., Brodie, E. D., 3rd, & Brodie, E. D., Jr. (2014). Confirmation and distribution of tetrodotoxin for the first time in terrestrial invertebrates: two terrestrial flatworm species (Bipalium adventitium and Bipalium kewense). PloS One, 9(6), e100718. Supek, F ., Bošnjak, M., Škunca, N., & Šmuc, T . (2011). REVIGO summarizes and visualizes long lists of gene ontology terms. PloS One, 6(7), e21800. Thielicke, M., & Sluys, R. (2019). Prey capture and feeding behaviour in an endemic land flatworm from São T omé Island. Journal of Natural History, 53(23-24), 1385-1393. Train, C.-M., Pignatelli, M., Altenhoff, A., & Dessimoz, C. (2019). iHam and pyHam: visualizing and processing hierarchical orthologous groups. Bioinformatics , 35(14), 2504–2506. Undheim, E. A. B., Jones, A., Clauser, K. R., Holland, J. W., Pineda, S. S., King, G. F ., & Fry, B. G. (2014). Clawing through evolution: toxin diversification and convergence in the ancient lineage Chilopoda (centipedes). Molecular Biology and Evolution, 31(8), 2124–2148. Verdes, A., Simpson, D., & Holford, M. (2018). Are Fireworms Venomous? Evidence for the Convergent Evolution of T oxin Homologs in Three Species of Fireworms (Annelida, Amphinomidae). Genome Biology and Evolution, 10(1), 249–268. Verdes, A., T aboada, S., Hamilton, B. R., Undheim, E. A. B., Sonoda, G. G., Andrade, S. C. S., Morato, E., Marina, A. I., Cárdenas, C. A., & Riesgo, A. (2022). Evolution, Expression Patterns, and Distribution of Novel Ribbon Worm Predatory and Defensive T oxins. Molecular Biology and Evolution, 39(5). https://doi.org/10.1093/molbev/msac096 Vila-Farré, M., & C Rink, J. (2018). The Ecology of Freshwater Planarians. Methods in Molecular Biology , 1774, 173–205. Vishnoi, S., Matre, H., Garg, P ., & Pandey, S. K. (2020). Artificial intelligence and machine learning for protein toxicity prediction using proteomics data. Chemical Biology & Drug Design, 96(3), 902–920. von Reumont, B. M., Lüddecke, T ., Timm, T ., Lochnit, G., Vilcinskas, A., von Döhren, J., & Nilsson, M. A. (2020). Proteo-Transcriptomic Analysis Identifies Potential Novel T oxins Secreted by the Predatory, Prey-Piercing Ribbon Worm. Marine Drugs, 18(8). https://doi.org/10.3390/md18080407 Yacoub, T ., Rima, M., Karam, M., & Fajloun, J.-M. S. A. Z. (2020). Antimicrobials from Venomous Animals: An Overview. Molecules , 25(10). https://doi.org/10.3390/molecules25102402 Yoshino, T . P ., Brown, M., Wu, X.-J., Jackson, C. J., Ocadiz-Ruiz, R., Chalmers, I. W., Kolb, M., Hokke, C. H., & Hoffmann, K. F . (2014). Excreted/secreted Schistosoma mansoni venom allergen-like 9 (SmVAL9) modulates host extracellular matrix remodelling gene expression. .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 14, 2025. ; https://doi.org/10.1101/2025.01.09.632098doi: bioRxiv preprint 29 International Journal for Parasitology, 44(8), 551–563. Zancolli, G., Reijnders, M., Waterhouse, R. M., & Robinson-Rechavi, M. (2022). Convergent evolution of venom gland transcriptomes across Metazoa. Proceedings of the National Academy of Sciences of the United States of America, 119(1). https://doi.org/10.1073/pnas.2111392119 Zayas, R. M., Cebrià, F ., Guo, T ., Feng, J., & Newmark, P . A. (2010). The use of lectins as markers for differentiated secretory cells in planarians. Developmental Dynamics : An Official Publication of the American Association of Anatomists, 239(11), 2888–2897. Zhang, Q., Xu, J., Zhou, X., & Liu, Z. (2022). CAP superfamily proteins from venomous animals: Who we are and what to do? International Journal of Biological Macromolecules, 221, 691–702. Zhang, Y . (2015). Why do we study animal toxins? Dong Wu Xue Yan Jiu = Zoological Research, 36(4), 183–222. .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 14, 2025. ; https://doi.org/10.1101/2025.01.09.632098doi: bioRxiv preprint