We selected six clownfish species spanning the full range of mutualistic behavior (see Litsios et al., 2012) from strict specialists (A. frenatus and A. biaculeatus), to specialists on a few anemone species (A. ocellaris and A. percula), to two generalists (A. clarkii and A. perideraion). We also sampled two damselfish species (Chromis viridis and Pomacentrus moluccensis) that are closely related to clownfishes but unable to interact within sea anemone. The number of individuals sampled ranged from 4 to 10 depending on the species. The adult fish were all imported from different areas of Indonesia (Supplementary Figure S1) and kept in sea water aquariums without sea anemones at Aqua-Home Aigle, Switzerland, for a quarantine of 10 days. The aquarium temperature was 25°C, salinity was 1.025, NH₃/NH₄ were kept at 0–0.2 ppm and nitrates were around 10 ppm. All fishes were under a 12 h day and 12 h night light environment.

We sampled the mucus on the skin surface of all individuals between September 2018 and May 2019. All fishes were anesthetized individually in a solution of sea water and Tricaine methanesulfonate (100 mg/l) for a few minutes (1–2 min). They were then removed from this solution and placed in a sterilized glass petri dish for mucus collection (3–4 scrapes per side) starting from the pectoral fin until the end of the tail with the help of a soft cell scraper (SARSTEDT, Nümbrecht Germany). Fishes were then put back into fresh sea water until they awakened and regained mobility. The mucus was washed off from the cell scraper using 500 μl of UHPLC water, transferred into an Eppendorf tube, snapped frozen on dry ice and stored at −80°C prior to metabolite extraction. All procedures were performed in accordance with protocols approved by the veterinary authority of the Canton de Vaud.

We lyophilized the mucus samples and extracted the metabolites using a biphasic extraction with the following solvent mixture: 540 μl of MTBE, 200 μl of MeOH and 360 μl of H₂O (Reverter et al., 2017). We vortexed the sample extracts before bath-sonication for 2 min to improve the metabolite extraction. We finally centrifuged the samples 15 min at 14,000 rpm at 4°C. The organic, upper phase (non-polar) was separated from the aqueous, lower (polar) phase and evaporated to dryness in a vacuum concentrator (LabConco, Missouri, US). The dry polar (aqueous phase) and lipid metabolite (organic phase) extracts were stored at −80°C prior to UHPLC-HRMS analysis. We normalized the sample amounts by measuring the total protein content in the protein pellet of each sample following metabolite extraction (see Supplementary Methods for Details).

We reconstituted the dry polar metabolite extracts in 100 μl MeOH:H₂O (4,1, v/v) and the lipid extracts in 100 μl Butanol:MeOH (1,1, v/v) spiked with the ISTD Splash Lipidomix mix (Avanti Polar Lipids, Birmingham, United States). Both reconstituted extracts were vortexed, sonicated for 30 s and centrifuged for 15 min at 14,000 rpm at 4°C to remove insoluble debris. We transferred the supernatants to LC vials for UHPLC-HRMS analysis.

We analyzed the aqueous mucus extracts by hydrophilic interaction liquid chromatography coupled with high-resolution mass spectrometry using electrospray ionization source operating in positive mode (HILIC (+) ESI HRMS) and in negative mode (HILIC (−) ESI HRMS) as described by Gallart-Ayala et al. (2018). The instrument was set to acquire over a m/z range between 30 and 1,200 Da, with the MS acquisition rate of 2 spectra/s.

The organic phase of clownfish mucus extracts was analyzed by reversed phase liquid chromatography coupled to high-resolution mass spectrometry using electrospray ionization source operating in positive mode (RP (+) ESI HRMS) and in negative mode (RP (−) ESI HRMS).

Analyses were performed on a 6,550 iFunnel Q-TOF mass spectrometer interfaced with the 1,290 UHPLC system (Agilent Technologies). In both modes, we analyzed samples using a pre-column ZORBAX Eclipse Plus UHPLC Guard 3PK, 1.8 μm, 2.1 mm × 5 mm coupled to a ZORBAX Eclipse plus rapid resolution HD, 1.8 μm, 2.1 mm × 100 mm (Agilent Technologies). The mobile phases were composed of A = 0.1% formic acid in water and B = 0.1% formic acid in IPA/ACN (90:10). For both modes, the linear gradient elution started at 15% B, rising to 30% B (2 min) then to 48–72% B (2.50–8.50 min) reaching 99% B (11.50–12 min) and back at 15% B (12.10–15 min).

We applied a 5-min post-run to ensure the re-equilibration of the columns at the end of every injection. The flow rate was 600 μl/min, column temperature 60°C and sample injection volume 2 μl. We set the source conditions for ESI as follows: dry gas temperature of 200°C and flow of 14 l/min, fragmentor voltage of 175 V, sheath gas temperature of 300°C and flow of 11 l/min, nozzle voltage of 1,000 V, and capillary voltage of 3,500 V. The instrument was set to acquire over a m/z range between 30 and 1,700 Da, with the MS acquisition rate of 3 spectra/s.

MS/MS data were acquired in data-dependent acquisition (DDA) mode on QC samples using iterative injections through PC-driven exclusion (Broadhurst et al., 2018). Parameters were set according to Agilent Technologies, Inc. Santa Clara, CA, United States (Koelmel et al., 2020). The additional MS/MS data were acquired in targeted mode on selected precursors with narrow isolation width (≈ 1.3 m/z), MS acquisition rate at 3 spectra/s and MS/MS acquisition at 3 spectra/s. The collision energy for polar metabolites was fixed at 25 eV, and for lipids at 20 and 40 eV.

We converted the raw UHPLC-HRMS metabolome data to mzXML files using ProteoWizard MS Convert. The mzXML files were uploaded to XCMS software for data processing including peak detection, retention time correction, profile alignment, putative identification and isotope annotation (Kuhl et al., 2012; Tautenhahn et al., 2012; Gowda et al., 2014). We processed the data as a multi-group experiment and the parameter settings were as follows: centWave algorithm for feature detection (Δm/z = 20 ppm, minimum peak width = 5 s and maximum peak width = 20 s, S/N threshold = 10, mzdiff = 0.01, integration method = 1, prefilter peaks = 3, prefilter intensity = 10,000, noise filter = 5,000); obiwarp settings for retention time correction (profStep = 1); and parameters for chromatogram alignment, including mzwid = 0.02, minfrac = 0.5 and = 5 (Libiseller et al., 2015).

We converted the raw UHPLC-HRMS lipidome data to ABF files using Reifycs Abf Converter. The ABF files were uploaded to MS-DIAL RIKEN software downloadable from the PRIMe website for data processing including peak detection, deconvolution, profile alignment, putative identification, and isotope annotation (Tsugawa et al., 2015). The following parameters were used: retention time begin: 1 min; retention time end: 13 min; mass range begin: 0 Da; mass range end: 1500 Da; MS1 (centroiding) tolerance: 0.01 Da; smoothing level: 6 scans; minimum peak height: 10,000 amplitude; mass slice width: 0.1 Da; accurate mass tolerance (MS1): 0.01 Da; accurate mass tolerance (MS2): 0.005 Da; retention time tolerance for alignment: 0.05 min; MS1 tolerance for alignment: 0.015 Da.

We exported the metabolite feature table (containing peak areas of all detected polar and lipid metabolites across all samples) to the R software (Chambers, 2008) and used the batchCorr package¹ for signal intensity drift correction and noise filtering (CV (QC features) > 30; Brunius et al., 2016).

We analyzed pooled quality control (QC) samples (representative of all analyzed samples per species) at regular intervals (every 5–6 samples) across the entire run to correct for the signal intensity drift and to filter chemical and/or informatic noise. Data was curated by applying the following filters: CV (QC) < 30%, peak intensity >10′000 ion counts, visual inspection of peak shape.

Metabolite levels were log transformed prior to statistical analysis. Positive and negative acquisition modes were merged to represent the global polar and lipid metabolome of fish mucus. We reduced the dimensionality of the data to better visualize the distribution of samples based on the molecular signature of the mucus found in the curated metabolome and lipidome using T-distributed stochastic neighbor embedding (t-SNE) as implemented in the R package “tsne”². We also used the UMAP algorithm on our data, but it led to the same pattern as with the t-SNE and we therefore do not show it. We also used PCA analyses to reduce the dimensionality of the data prior to hierarchical clustering based on Manhattan distances. The analyses were done using the R packages “FactoMineR”³ and “factoextra”⁴.

Finally, we selected the metabolite features whose levels varied between different groups (clownfishes vs. damselfishes and generalists vs. specialists) using pairwise t-tests. We did not correct for phylogenetic relationships because of the small number of species involved and the fact that the main comparison (clownfish vs. damselfish) is done between two monophyletic groups. We corrected the p-values obtained using an optimized False Discovery Rate (FDR) approach based on the Benjamini Hochberg procedure (Benjamini and Hochberg, 1995) and set the FDR threshold to 0.05. For the metabolites that were significant after FDR correction, we evaluated the difference in mucus profiles by calculating the log-fold change of each metabolite between the groups tested. We used the peak area of each metabolite feature as the measure to scale the data. Finally, we highlighted the metabolite features showing the highest log-fold change in a circular bar plot using a custom-made R script⁵.

Following data pre-processing (metabolite/lipid profile alignment) and statistical analysis, metabolite/lipid features were matched against Human Metabolome Database (HMDB) (Wishart et al., 2007), Metlin (Smith et al., 2005) and LipidMaps (Fahy et al., 2007) based on accurate mass to charge ratio (m/z) with Δppm = 20, using the online CEU mass Mediator tool⁶.

We validated the identifications of the metabolites by MS/MS data matching (Mylonas et al., 2009) against the in-house PCDL database as well as METLIN⁷ (Smith et al., 2005; Benton et al., 2008). The spectra from lipid profiles were matched against the LipidBlast database for lipid identification (Koelmel et al., 2020). The GNPS platform was used to perform feature-based molecular networking (Aron et al., 2020) to help us identify features that could not be recognized by other available databases.

Raw polar metabolome (HILIC-HRMS) and lipidome profiles (RPLC-HRMS) comprised 3,366 and 4,126 features in total, respectively. Following noise [using CV (QCs) < 30%] and redundancy filtering (by feature annotation and grouping), we reduced these complex datasets to 938 and 903 unique polar and lipid metabolite features, respectively. Within these comprehensive profiles, we were able to confirm the identity of 46 known polar metabolites and 232 lipid species.

The t-SNE visualization of polar metabolite and of lipid signatures (including unknown features) revealed two main clusters: clownfishes and damselfishes (Figures 1A,B). While both signatures allowed for the discrimination between the clownfishes and the damselfishes, the polar metabolome profile was also species-specific and allowed for narrower grouping of individuals of the same species. The exceptions were two samples of the species A. perideraion that clustered close to A. ocellaris and A. percula. Based on the lipid signature, two individuals of the strict specialist A. biaculeatus and three individuals of A. frenatus also remained isolated farther away from their main cluster in the center of the t-SNE (Figure 1A).

Moreover, the t-SNE visualization of the polar metabolome signature (Figure 1B) showed three additional, well-defined, clusters. The first one included all samples from the most generalist species A. clarkii and A. perideraion. The second grouped the two strict specialists A. frenatus and A. biaculeatus, which share the same sea anemone host species. The third cluster represented the specialists and sister species A. ocellaris and A. percula with two individuals (the same as in Figure 1A) of the generalist A. perideraion.

The hierarchical clustering analyses done on the PCA loadings for both the metabolome and lipidome features corroborated the t-SNE clusters (Supplementary Figure S2). It further showed that the clownfish mucus composition was more variable than the damselfishes for both datasets. While clownfishes were again grouped in clear clusters based on species identity (except for the two individuals of A. perideraion and a few other exceptions; Supplementary Figure S2), they were distributed across the deeper splits of the hierarchical clustering. In contrast, the two damselfish species remained grouped in narrow clusters regardless of the data type.

We identified the features that were specific to clownfishes versus damselfishes and, within the clownfishes, generalists versus specialists. In the lipidome, the abundance patterns of 92 lipid features were significantly different {FDR < 0.05 and abs [log₂(FC)] > 2} between clownfishes and damselfishes. Among them, 52 lipid species were more abundant in clownfish than in damselfish mucus, while 40 lipids showed the opposite pattern (Figure 2A). Most of the lipids present in higher levels in clownfish mucus belonged to the sphingolipid class of ceramides, whereas the more abundant species in damselfish mucus were mainly glycerophospholipids and glycerolipids (Figures 2B,C).

Distinct lipid signatures of clownfish generalists versus specialists mucus were characterized by 63 lipid features {FDR < 0.05 and abs [log₂(FC)] > 2}, out of which 16 were more abundant in generalists and 47 in specialists (Figure 3A). Among the lipid features identified, generalist mucus was richer in glycerolipids whereas specialist mucus was richer in phosphatidylcholines (Figures 3B,C).

Based on the polar metabolome signature, family specific clustering of clownfishes versus damselfishes was characterized by differential patterns of 224 metabolite features {FDR < 0.05 and abs [log₂(FC)] > 2}: 175 were more abundant in clownfish mucus compared to damselfishes whereas 49 showed the reverse pattern (Figure 4A). The main class of polar features identified as more abundant in clownfish mucus were carnitines (polar lipids) followed by a few amino acids and derivatives, and peptides (Figure 4B). Different amino acids and peptides were, however, also present in higher levels in damselfish mucus compared to clownfish (Figure 4C). Among the unidentified polar features higher in damselfish mucus (i.e., “unknowns” that did not return any matches against databases), many eluted at 2 min, which suggested that these may be a distinct class of polar molecules specific to damselfish mucus (Figure 4C). Within clownfishes, 100 polar metabolite features significantly differed between clownfish generalists and specialists {FDR < 0.05 and abs [log₂(FC)] > 2} with the majority present in higher levels in specialists (Figure 5). However, very few polar features differing amongst clownfish generalists and specialists could be identified (See Supplementary Table S1; Supplementary Table S2 for lipid and polar metabolite feature IDs in fish mucus respectively).

It should be noted that because of the nature of our samples, which come from non-model species, and incomplete metabolite/lipid databases available, there are still many lipids and polar metabolites significantly differing in our fish groups that could not be identified.

We detected several different classes of metabolite features in clownfish and damselfish skin mucus such as amino acids, peptides, carnitines, sphingolipids and glycerophospholipids (Supplementary Tables 1, 2; Supplementary Figure S1).

While the chemical characterization of polar and lipid metabolites in samples from non-model species remains challenging, we successfully characterized the presence of 46 polar metabolites and 232 lipid species in the global molecular profile of clownfish and damselfish mucus. Using pattern recognition approaches, we showed that the acquired metabolic signatures (including annotated and unknown features) are highly discriminative of the clownfish and damselfish skin mucus. These results suggest that the clownfishes have a specific chemical signature that could be associated with their adaptation to thrive within sea anemones. Although the current data do not allow us to fully understand the mechanisms involved and the metabolite function, we found that the signatures of the two damselfish species are very similar and far less variable compared to the signatures of closely related clownfish species. This was unexpected because the divergence time between the two damselfish species is much larger than between the clownfishes (50-million-year time span within the damselfish family whereas clownfishes diversified within the last 15 million years) (McCord et al., 2021). Instead, social structure and living behavior of damselfishes are more similar compared to those within the clownfish clade. It suggests that the host-specific association of clownfishes with sea anemones could lead to a larger variation in their mucus composition because of the interaction and adaptation to this new environment.

Among a multitude of differences in metabolic signatures, we revealed one major lipid class which was significantly more abundant in clownfish mucus. Clownfish mucus was characterized by the presence of higher (chemical) diversity and relative abundance of ceramides. Interestingly, most of the detected ceramides had an even number of carbons in their chains. Ceramides are precursors of most sphingolipids and one of the enzymes involved in their biosynthesis, serine palmitoyltransferase (SPT) (Merrill Jr, 2002), is highly conserved in both eukaryotes and bacteria (Yard et al., 2007). Eukaryotes, however, mainly synthesize even-chained sphingoid backbones and side chains, while bacterial sphingolipids consist mainly in odd-chain sphingoid backbone and fatty acyl chains (Zitomer et al., 2009). The roles of sphingolipids are many, ranging from plasma membrane constituents and inflammation regulators to important messengers in apoptosis and cell growth (Bikman and Summers, 2011; Olsen and Færgeman, 2017). The ceramides identified in clownfish mucus range from 32 to 50 carbon atoms and the majority has an even number of carbons in their acyl chains suggesting that they are synthesized by the clownfish itself. Ceramides could be responsible for the mediation of the interaction with sea anemones, and this particularly concerns species with a carbohydrate backbone containing glucose, galactose, N-acetylglucosamine, N-acetylgalactosamine, fucose, glucuronic acid or sialic acid (Miller-Podraza, 2000). This specific part of ceramides can be recognized by receptors of human and mouse immune cells thereby activating the production of several cytokines (Morales-Serna et al., 2007). The presence of sensory and chemoreceptors binding N-acetylated sugars triggers the nematocyst discharge (Ozacmak et al., 2021) although it has been shown that the co-presence of proline can counteract this effect (Watson and Roberts, 1995). The large number of ceramides found in the skin mucus of clownfishes compared to damselfish may act both as a physical and immune protective layer but also as the mediator of chemical communication between the fish and the host sea anemone. It further suggests that investigating the N-terminal sugars of glycolipids in clownfishes and damselfishes could shed light on the molecules triggering or inhibiting the discharge on these fish.

Compared to clownfish mucus, the lipidome part of the damselfish mucus is richer in phospholipid content. We identified several abundant phosphatidylcholine (PC) and phosphatidylethanolamine (PE) species containing polyunsaturated long chain fatty acids, such as C20:4 arachidonic acid and C22:6 docosahexaenoic acid. The glycerophospholipids detected in the damselfish mucus could either come from the fish itself or from the damselfish-specific skin microbiome (LIPID MAPS® Lipidomics Gateway, n.d.).

Among polar lipid-like metabolites, we identified several acylcarnitines that were more abundant in clownfish mucus compared to damselfishes. In mammals, and probably in fish too, carnitines play an important role in fatty acid oxidation and constitute the markers of mitochondrial function (Vissing et al., 2019). Unfortunately, there are, to our knowledge, no studies discussing the role of carnitines in fish skin mucus. It therefore remains challenging to explain this difference observed between the two groups. Certain polar features such as ectoine and ergothioneine found at higher abundance in clownfish mucus are usually produced by bacteria, which could therefore mean that certain features detected may be coming from the microbiome present on fish skin surface. Ectoine has a role in protecting against osmotic stress and in skin health (Kocherovets, 2019). Ergothioneine is an antioxidant produced by bacteria and fungi (Kitsanayanyong et al., 2021). This compound may be produced at higher levels after contact with sea anemones by the bacterial flora on clownfish skin, but this will have to be tested and cannot be shown with our data. Trigonelline was also found to be significantly higher in clownfish mucus. This compound is a toxin that has already been detected in sea urchins and cnidarians (Arumugam et al., 2013) and was defined as a fear inducing cue in mud crabs (Poulin, 2017). The role of this compound in fish mucus is not known, but it may keep certain predators away from sea anemone locations if the chemical cue is released by clownfishes.

In damselfishes, the amino acids GABA and L-glutamate were found to be present at higher levels compared to clownfish mucus, but no specific role of these amino acids at the skin surface has been described. Amongst the unknown polar features higher in damselfish mucus, there are numerous molecules eluting at 2 min with respect to the applied gradient. This group accounted for most of the difference between clownfish and damselfish mucus and may therefore be an important portion of the moderately polar signature of damselfish mucus to be further explored.

The acquired metabolic signatures of mucus also allowed for the discrimination between generalists and specialists within the clownfish clade. The main lipid classes contributing to the differentiation between generalists and specialists mucus were sphingomyelins, phosphatidylcholines and some glycerolipids. Sphingomyelins are sphingolipids, considered as inactive precursors of ceramides, which suggest again their implication in the mediation of the mutualistic interaction with sea anemones.

The acquired metabolic signatures are highly species-specific. This was the case even for some species (e.g., the specialist A. biaculeatus and the generalist A. perideraion) whose individuals were sampled from different geographical locations in Indonesia. Metabolomic patterns are often strongly influenced by external environmental cues (Cachet et al., 2015), but the strong species-specific patterns that we detected could be due to either genetically induced responses (Ivanišević et al., 2011) or to the specific interaction with the species of host sea anemones from which they were reared in their natural habitat. Unfortunately, the latter hypothesis cannot be confirmed because we lack this information; further experiments should be done to test host influence on mucus molecular patterns.

The two strict specialists A. frenatus and A. biaculeatus which associate only with E. quadricolor in their natural environment and we observed that they clustered in proximity on the t-SNEs, suggesting that host preference may impact the metabolites present in their mucus (Lindon et al., 2007). It has already been shown that, in plant - pathogen and plant - fungi interactions, a metabolic shift occurred in the plant after penetration of the pathogen or after association with fungal endophytes (Jones et al., 2011; Nieva et al., 2021). This topic is scarcely investigated in marine non-model organisms, but suggests that for clownfishes, the external environmental factor weighing the most in mucus chemical modification could be the physical interaction with the host anemone. Since the study of clownfish and anemone microbiomes showed that after several days of symbiotic mutualism the clownfish microbiome shifts towards a similar composition as the anemone (Roux et al., 2019), it is plausible that this shift also occurs at a metabolite level. The metabolites present in anemone mucus should also be characterized to test this hypothesis since little is known about the anemone mucus chemical composition.

Considering the phylogenetic distance between the clownfish species in our study (Litsios et al., 2012), the mucus of A. frenatus, a strict specialist, clusters near one of the other specialists of this study. In addition, although the two generalist species are not very close phylogenetically, they also had a very similar chemical composition of their mucus. These results suggest that the chemical signature of fish mucus is species-specific, as a product of species-specific gene activity; there may also be an imprinted signature reflecting the choice of habitat such as the species of host sea anemone in natural environments.