Integration of Global Lipidomics and Gonad Histological Analysis via Multivariate Chemometrics and Machine Learning: Identification of Potential Lipid Markers of Ovarian Development in the Blue Mussel (Mytilus edulis)

1. Introduction

Mussels are important species for marine ecosystems as they are recognised as habitat engineers and provide valuable ecosystem services for coastal areas and their economies [1,2]. Mussels are valuable aquaculture resources characterised by an extremely low environmental footprint [3,4], high nutritional properties [5], and their demand as food products is projected to increase in the future net-zero society [6]. Unfortunately, despite their prospective importance, mussel production in Europe has stagnated since 1999 [7], and the expansion of the mussel aquaculture sector is hampered by various factors, including climate change [8]. Therefore, it is important to increase our knowledge of mussel physiology and biochemistry and to identify further useful tools to monitor and predict the ability of mussels to adapt and thrive in future environmental scenarios.

Gonadal development is the main physiological process that runs through the annual life cycle of temperate mussels [9]. Blue mussels (Mytilus edulis) have an annual reproductive cycle with an accumulation of energy resources in late summer and autumn, a gonadal maturation process that extends through winter, and one or more spawning events in spring and early summer [10,11,12,13,14]. During maturation, gonads proliferate on the mantle, which undergoes a series of ultrastructural and biochemical reorganisations. In mature mussels, the mantle can account for 50% of the total body dry weight [14]. In resting females, the mantle consists of a single-layered cuboidal epithelium that is ciliated on the inner surface and composed of adipogranular (ADG) and vesicular connective tissue (VCT) cells [15]. ADG and VCT cells support oocyte development by providing glycogen and lipids, respectively [12]. The oocytes increase in numbers and mature, accumulating lipid droplets that form the yolk reserves of the embryos [16]. Glycogen is consumed during the development process, while lipids are accumulated in the oocytes [17,18]. This process leads to the annual accumulation patterns of lipids in several bivalve species, including Crassostrea gigas [19,20,21], Ostrea edulis [22], M. edulis [17,23], Mytilus galloprovincialis [24], Modiolus barbatus [25], Chlamis nobilis [26], Tapes decussatus and Tapes philippinarum [27].

Detailed information about the maturation of the gonads is important for the industrial production of mussels in hatcheries, where the mussels are artificially conditioned (so-called broodstock conditioning) to maximise gamete production. Mussels do not store their oocytes for prolonged periods of time, as it is their physiology to partially or completely spawn within a few weeks of sexual maturity [28]. Missing the right time window for gonad development can lead to unnecessary resource consumption, uncontrolled spawning or the occurrence of pre-spawning atresia [28,29,30].

Various histological, biochemical and molecular methods have been used to assess gonadal maturation in bivalves. Historical approaches included the use of allometric ratios as a measure of mussel condition. The condition index, which measures the ratio of soft tissue to shell [31], and the gonadosomatic index, which measures the ratio of mantle to soft tissue [32], are two commonly cited examples in the literature. Traditional techniques also include gonadal smears and gonadal histological analysis (GHA) [33,34]. The fastest method to obtain information on the gonadal status of mussels is the gonadal smear, although the loss of tissue structure due to sample preparation makes it difficult to distinguish between different stages of maturity. The ultrastructure of the tissue is well preserved in the GHA, making it the most important technique for analysing the reproductive cycle in bivalves [10,12,30,34,35,36]. The limitations of GHA are the limited sample size (e.g. the number of samples that can be processed per time unit), the difficulty in distinguishing transitional stages of development, and a certain degree of subjectivity that requires experienced operators in classifying the degree of development of gonads [37,38]. GHA has also been coupled with stereological indices, leading to the development of several quantitative indices of gonadal maturation, such as gamete volume fraction (GVF) or oocyte cytoplasmic area, overcoming some of the limitations of classical GHA [11,39,40,41,42]. Hines and co-workers compared classical histological approaches with metabolic and molecular techniques and showed that the latter two methods are advantageous in distinguishing sex in dormant or spent mussels [43].

Lipids play important role in mussel gonad maturation, as several authors have reported direct changes in lipid composition as a result of the maturation process [18,21,44,45,46]. Triglycerides (TG) were the major lipid class that changed during gonadal development in the ovaries of both Pecten maximus and C. gigas, accumulating in mature individuals [21,44]. Martinez-Pita and co-workers found no changes in the composition of lipid classes between the different stages of ovarian development in M. galloprovincialis and Donax trunchulus [18,45]. When analysing fatty acids (FA), these authors observed an increase in FA18:4n-3 in mature females of both species and of FA18:0, FA18:3n-3 and docosahexaenoic acid (DHA, FA22:6n-3) specifically in mature female mussels. The complexity of observing a direct effect of gamete development on gonad lipid composition may be related to the seasonal behaviour of the gonad development in bivalves, which tends to span over several months [11,12,13,32]. Furthermore, gonads are also a non-uniform tissue, as several developmental stages (or regressions) can occur simultaneously in the same organism [18,47], which greatly increases the difficulty of such a task.

Lipidomics and the use of liquid chromatography coupled with mass spectrometry (LC-MS) could offer several advantages over conventional lipid analysis techniques when analysing a complex tissue such as bivalve gonads. First, the study of lipids at molecules level combines information from FA and lipid class composition analysis, revealing patterns that are otherwise lost when relying on only one of these two traditional lipid profiling techniques [48]. Second, modern MS platforms are characterised by higher sensitivity than traditional detectors used in lipid class analysis, providing an opportunity to evidence and extend existing knowledge on the role covered by minor, uncommon and poorly understood lipid classes in ovarian maturation of mussels. To date, rare are the studies investigating gamete development in marine invertebrates using LC-MS. One of the first examples of this type of study comes from Chansela and coworkers who investigated the lipid composition of Paeneus marguirensis ovaries using a mass spectrometry imaging approach [49]. The authors showed how follicles from the same histological preparation, but with different maturation stages, differed in their composition of major phosphocholine (PC) and TG molecular species. The results of this study demonstrate the presence of low unsaturated PC species (PC 32:1) and highly unsaturated TG (58:9) especially in the early developing gonadal areas. However, PC and TG associated with essential polyunsaturated FA such as eicosapentaenoic acid (EPA, FA20:5n-3), arachidonic acid (AA, FA20:4n-6) and DHA were only abundant in mature follicles.

More recently, the lipidome of the gonads of M. edulis has been characterised in detail, focusing on sex-specific differences between testis and ovaries [50]. In the present study, we aim to build on this previous work and continue here to analyse the blue mussel ovarian (BMO) lipidome by integrating GHA with LC-MS global untargeted lipidomics through a comprehensive statistical approach, incorporating unsupervised and supervised multivariate chemometrics and machine learning (RandomForest – RF [51]), applied here to highlight potential lipid markers that distinguish between histologically ripe (stage III) and non-ripe BMOs (stages I-II).

2. Materials and Methods

2.1. Mussels Collection and Sample Preparation

Commercial size blue mussels (Mytilus edulis, average shell size 5.6 ± 0.6 cm) were obtained from local growers on the West Coast of Scotland (Inverlussa Marine Services - www.inverlussa.com, Loch Spelve, Isle of Mull, UK) in February 2017. Mussels rearing and processing (sizing, feeding, aquarium rearing conditions) used in this study has been described in a previous publication [50] and details are provided in Methods Supplement 1. Each adult was opened “book-like” by cutting the hinge and abductor muscle with a scalpel and carefully removing the gills with bow scissors. The left valve mantle lobe was dissected and placed in a 1.5 mL tube (Eppendorf), snap-frozen in liquid nitrogen and stored at -80 °C for lipid analysis. The right mantle lobe was excised and stored in 7 mL Bijoux tubes (VWR) in 4% buffered seawater formalin for GHA.

2.2. Histological Analyses

The sex of each adult animal was determined by GHA after haematoxylin and eosin (H&E) staining, the standard method for sex determination and assessment of gonadal maturation [30,35,43,52]. H&E-stained slides were photographed with an Axioshop 2 inverted microscope (CarlZeiss) equipped with an AxioCam MRc5 camera (CarlZeiss) and acquired with the corresponding AxioVision 4.7 software. The GHA scores were assigned according to the ultrastructure of the gonads using a 4-staging classification system [10,12,53]. The reference samples used in the gonad classification system are shown in Methods Supplement 2a. Three different sections (anterior, medial and posterior) of the mantle of each mussel were scored to distinguish between ripe (stage III, N=25) and non-ripe (stages I-II, N=26) ovaries. The distinction between vitellogenic and atretic/degenerating gonads was made according to the criteria reported by [28] and [16] as described in Methods Supplement 2b. When different stages occurred simultaneously in a single section, the decision on staging was based on the condition of the largest part of the preparation [18]. A stereological variable to quantify the percentage of mantle area covered by gametes (gamete volume fraction – GVF) was calculated by image analysis of the GHA preparations using ImageJ software (www.imagej.nih.gov) via the colour threshold method (see Methods Supplement 2c for details). GVF was calculated as the percentage between the total mantle surface area and the mantle surface area covered with gametes. The relationship between GHA and GVF is shown in Appendix A (Figure A1). All GHA and GVF scores are provided in Supplementary Data 1.

2.3. Lipid Extraction

Frozen dissected BMOs were freeze-dried (ẞ18 LO plus, Christ) and ground to a fine powder in liquid nitrogen; ≈ 0.05 g of powdered dry BMOs were individually homogenised in HPLC water (Fisherbrand) for 1 min in glass potters (1:8 weight:volume) on ice. Aliquots of the homogenate (20 µL) from each gonad were pooled and formed the biological quality control (BQC) samples. Lipid extraction was performed according to [54] as described in [50]. Total lipids were extracted from 50 µL of BMO homogenate and BQC with 6 mL of 2:1 chloroform:methanol (both Fisherbrand) in pre-combusted (450°C 5 hr) 20 mL screw-cap vials (Pyrex). The mixture was shaken vigorously, then each vial was flushed with nitrogen and extracted on ice for 1 hr under dark conditions. The organic extract was separated from the pellet by centrifugation (350 g, 5 min, 4 °C) and collected in a clean tube. The pellet was re-extracted with a further 3 mL of chloroform:methanol. A 0.9 % KCl solution (VWR) was added to the organic extract to obtain a final ratio of 2:1:0.8 chloroform:methanol:KCl, and the polar and organic phases were separated after a second centrifugation step (350 g, 5 min, 4 °C). The organic layer was recovered and evaporated under nitrogen flow (NVap, Organomation). The dried lipid extracts were weighed (0.0001 g accuracy) and resuspended in 0.5 mL chloroform, yielding the BMO total lipid extract (TLE).

2.4. LC-MS Untargeted Lipidomics

Global lipidomics analyses of BMO TLE were performed by liquid chromatography coupled with high-resolution mass spectrometry (LC-MS) via a binary ultra-high pressure liquid chromatography (UHPLC) system (1250 Accela, ThermoFisher) connected to an electron spray ionization (ESI) and an Orbitrap mass analyser (Exactive, ThermoFisher). Separation was performed on a C18 Hypersil Gold 100 × 2.1 mm 1.9 nm particle size (ThermoFisher) maintained at 50 °C. The binary solvent system comprised a constant flow rate of 400 μL min^-1 with a gradient as described in Supplementary Table S1. Water and acetonitrile were HPLC grade and purchased from Fisherbrand, isopropanol was LC-MS grade (Hypergrade LiChrosolv, Merck), while ammonium formate and formic acid were both LC-MS grade and purchased from Sigma-Aldrich. TLEs were resuspended in 3:1 methanol:chloroform at a concentration of 1 mg mL^{- 1} for 3 μL injection volume.

All samples were analysed in positive (POS) and negative (NEG) ion mode over the m/z range of 250–2000 Da at a resolution of 100,000 FWHM and a scan speed of 1 Hz. Mass error was kept below 5 ppm by routine calibrations for both polarities using a calibration solution (Pierce™ LTQ ESI calibration solutions, ThermoFisher). Chromatograms and mass spectra were reviewed and integrated using Xcalibur software (ThermoFisher).

2.5. Data Processing

The raw LC-MS data were processed using Progenesis QI (Nonlinear Dynamics, Waters). Lipids were identified by the exact mass of the precursor ion (MS', 5 ppm Δm/z), and with the exception of CAEP (ceramide aminoethyl phosphonates [55]), lipid IDs are reported according to the LIPID MAPs nomenclature standard system [56,57]. Raw data from LC-MS lipid adducts for POS and NEG were experimentally evaluated using the lipid standard mixture (Supplementary Table S2) purchased from Avanti Polar Lipids (www.avantilipids.com) and Matreya (www.matreya.com). Lipid IDs were obtained by searching for exact masses in the lipid mass spectra libraries LIPID MAPS (www.lipidmaps.org), Human Metabolome Data Base – HMDB (www.hmdb.ca) and Metlin (www.metlin.scripps.edu), as well as from an in-house lipid database derived from lipidomics studies on bivalves and other marine invertebrates, which served as a tool for the identification of unusual marine lipids [55,58,59,60,61,62,63,64,65,66]. Isomeric lipid molecules resolved by reverse phase separation (same exact mass 5 Δppm, but different retention time) are distinguished with a lowercase letter after the lipid ID (i.e. PC 38:5a and PC 38:5b). The full set of lipids identified in the BMO and used in this study is included in Supplementary Data 1. All LC-MS raw data files and experimental metadata have been deposited in the Metabolights [67] data repository (www.ebi.ac.uk/metabolights/MTBLS8492) and are available under the study identifier MTBLS8492.

2.6. Data Analysis

The statistical analysis was performed with the statistical platform R (version 4.2.0). Data are presented as mean ± standard deviation unless otherwise stated. Statistical differences were considered significant at p<0.05.

The peak intensity tables (PIT), created from LC-MS data obtained by POS and NEG ionisation modalities, were manually sorted to remove contaminant masses. Further data filtering included removing features with >30% missing values, with the remaining missing values replaced by a small value (half of the minimum intensity value) and characteristics with low repeatability removed by filtering out using BQC and the interquartile range. The filtering of the raw data was performed using the filtering functions of the R-based package "MetaboAnalystR" [68]. At this stage, POS and NEG data were merged into a single dataset to avoid redundancy in the analysis. Duplicate lipids (i.e., observed in the POS and NEG data) were manually checked for their retention time between the different samples; if they matched, the ionization polarity that gave the best results in terms of signal for that specific lipid class was retained, assessed using as reference the lipid standard mixture.

The PIT was autoscaled (mean centered and divided by each lipid standard deviation), and the preliminary approach used to highlight lipids that differ between ripe and non-ripe BMOs involved principal component analysis (PCA) obtained from the R package "factoextra" [69] as unsupervised statistical model [70]. This approach was used to highlight the presence of outliers; a brief discussion of the results can be found in Appendix A. Other important markers were identified by supervised statistical methods on Pareto-scaled data (mean centered and divided by the square root of the standard deviation of each lipid). Pareto-scaling of the data was used to include information regarding molecule abundances that would otherwise be lost by autoscaling [71]. The first approach used in biomarker discovery involved the identification of significant variables via volcano plots. The volcano plot was calculated via the R package "MetaboAnalystR" [68] and plotted via the package "EnhancedVolcano" [72]. For the volcano plot, the lipids were filtered by a fold change (FC) >1.3 and following Wilcoxon rank sum test with a false discovery rate (FDR, [73]), adjusted p-value <0.05. Orthogonal partial least squares discriminant analysis (OPLS-DA, [74]), from the "MetaboAnalystR" package was the second supervised method used to identify lipids associated with the two gonad maturation conditions. The performance of OPLS-DA was assessed by 10-fold leave one out cross-validation (LOOCV) and model fit was tested by a 2000-fold permutation test (see Appendix A [75]). Important lipids were identified in the Sigma plot, and lipids characterised by a p(corr)<-0.6 and >0.6 were retained as considered important in distinguishing ripe and non-ripe BMOs. RandomForest – RF [51] was the last method used here. RF is a machine learning algorithm that was trained as shown in [76] and was obtained from the package "RandomForest" [77]. In this study, RF was used as a classification method (RFclass) using the grouping variables (ripe/non-ripe) and as a regression method (RFreg) considering the GVF. In both cases, the best 10% of RF classification (lowest out-of-bag error rate – OOB-ER) and RF regression (highest percentage of variability explained – PVE) were considered as important markers. A backward purging approach for the top 5% of lipids for both RFclass and RFreg was used to identify the most important molecules detected by the two approaches (see Methods Supplement 2d for details on model training). The results of all four statistical approaches were compared using a Venn diagram and ranked according to their ROC (Receiving Operating Characteristics) curves and AUC (Area Under the Curve) values. ROC and AUC values were calculated using the biomarker discovery function of the MetaboAnalyst package. Details of the full ROC and AUC analyses of all biomarkers identified in this chapter can be found in Table C1 (Appendix C).

The Lipid Ontology (LiOn) enrichment analysis terms were calculated and analysed via the LiOn web platform [78]. PCA analysis of LiOn terms was visualised via GO-PCA [79] and presented as a heatmap plot using Euclid as the distance matrix and Complete as the clustering algorithm. Enriched LiOn terms were ordered by p-value (calculated with a t-test), and terms were considered significantly different if the FDR-adj p-value was <0.05 [73].

3. Results and Discussion

In this study, a total of 51 BMOs were analysed using LC-MS global untargeted lipidomics. Ovaries were categorised as ripe (N=25) or not ripe (N=26) based on their GHA scores and GVF (Supplementary Data 1). Stereological analyses such as GVF are less prone to subjectivity bias than GHA, especially when assessing histological preparation, which indicates the presence of different gonadal stages in the same tissue [80]. In contrast, gonadal atresia and gamete degeneration are hardly recognised when using GVF alone [28]. Therefore, the two approaches offer the best performance in combination. GVF increased significantly in ripe BMOs, providing higher confidence in the criteria used for GHA classification (Figure A1).

During ovarian maturation, the area covered by follicles increases due to oocyte proliferation, with the volume covered by energy support cells decreasing [10,15,16]. As the oocytes proliferate, lipids are accumulated in the gonads at the expense of glycogen content [17,18,41]. This usually leads to the reported strong seasonal fluctuations in the major lipid classes of bivalves [19,23,27,44,81]. BMOs lipidomics profiles here are integrated with GHA and GVF through various chemometric and machine learning approaches to identify specific lipid molecules associated with the maturation process of BMOs.

3.1. Multivariate Chemometrics

Our investigation started with the application of a multivariate unsupervised chemometric approach such as PCA (Figure A2). As unsupervised approach (e.g. samples are plotted without “a priori” consideration of grouping variables), PCA does not suffer from problems of overfitting, however, PCA can be difficult to interpret to distinguish highly similar and heterogeneous groups with small multivariate differences, a large number of features and a small amount of variability explained by the main axes [82]. Selecting the correct orthogonal axes that best explain the differences between ripe and immature BMO was problematic through PCA, because a high percentage of variability is explained by the first four pairs of principal components of the PCA model (>10%, Figure A2A). Indeed, in all the PCA plots showing the combination of main 5 principal components, the two BMO groups did not evidence marked differences on their lipidomic profiles, observation expected given the heterogeneous origin of the different gonad samples considered in this study (see Method Supplement 1).

To highlight specific lipids linked to the two BMO groups, we continued the biomarker discovery approach by using more powerful supervised statistical models such as volcano plot, OPLS-DA and RF. Through volcano analysis important lipids distinguishing ripe from undeveloped BMOs are highlighted based on their fold change (FC, here 1.3 is considered as cut-off value) and their statistical significance (FDR p<0.05, calculated according to a non-parametric Wilcoxon rank test). The resulting plot is shown in Figure 1A and led to the selection of 25 significant features (Table A1). Volcano plots provide a quick and reliable option for identifying markers that define a particular condition, although this statistical approach is limited to pairwise comparisons and often lacks statistical power, especially when analysing typical -omics datasets where the number of variables is far greater than the number of observations. In these cases, correction factors for the false discovery rate are usually applied, which reduce the statistical power of the test [83,84], and, indeed, volcano analysis resulted in the most restrictive (lowest number of important feature highlighted) of the methods discussed in this work.

OPLS-DA is an extension of partial least squares regression and an approach characterised by a fast computation time and a decomposition of the variance of the samples into uncorrelated (orthogonal to the groups, orthogonal T-score) and correlated (predictive, T-score) to the grouping factor [74]. The results of OPLS-DA can be found in Figure 1B and show that 10.2% of the variance within our observation is correlated with the two BMO maturation classes (T-score). OPLS-DA as a supervised method may suffer from overfitting problems, especially when the number of observations is small. Therefore, an assessment of model performance (cross-validation) and permutation tests (time-consuming) are required to obtain a measure of model fit [75]. Model performances and fitting were evaluated respectively by a 10-fold LOOCV (predicted: R²x:0.102, R²y: 0.5, Q²:0.363, orthogonal: R²x: 0.167, R²y: 0.17, Q²: 0.07, Figure A3A), and 2000-fold permutation (Q²: 0.436 p<0.001 and R²y: 0.669 p<0.05, Figure A3B). From the resulting S-Plot, a set of 40 important lipids (Table A1) correlated with the grouping variable (ripe and non-ripe BMOs) were selected, setting an arbitrary threshold for the p(corr) value <-0.6 and >0.6 (Figure 1C).

The final approach used in this study is RF, a machine learning method that creates a series of classification (or regression) trees (a “forest”) by continuously bootstrapping the sample set. For each bootstrapping, the data set is split into a test sample (66% of the samples) and an out-of-bag sample (training data with 33% of the observations). As a bootstrapping method, RF is also robust to several attributes that may constrain other statistical approaches (e.g., zero values, normality and homoscedasticity assumptions), while in classification mode it assumes a balanced number of observations [51,76], and compared to the previous statistical approaches, RF is far more computationally intensive [76]. From the training data, the algorithm recursively selects a random set of variables whose ability to predict a categorical response (ripe or non-ripe BMO) or its correlation with a continuous response variable (as GVF) is evaluated using training data and tested on OOB samples.

Partitioning is based on the identification of the predictor that minimises either the within-group variance (categorical response variable) or the error in the regression against the response, in the case of a continuous response variable [51,76]. The optimal predictor becomes the best node of the forest, and the process is repeated for all the number of trees specified in the function. In RF classification (RFclass), the variables with the lowest OOB error rate (OOB-ER) are the best predictors for the categorical variables, with the lower OOB-ER indicating the higher predictive significance of a given molecule. In RF regression (RFreg), the importance of each predictor is evaluated as the proportion of variation explained (PVE), a measure of the reduction in the predictive ability of the tree when that variable is removed. A backward purging (BP) approach is used to select the most important features in the RF. BP consists of performing multiple RFs on the most informative subset of variables (here the top 5% of features are selected from both RFclass and RFreg, Figure 2) and recording the performance of each variable (in terms of OOB-ER or PVE). The results of the BP process of the RFclass data showed that the saturated PC 28:0 (OOB-ER 0.12) was the best predictor of ripe/non-ripe condition (Figure 2A). As for RFreg, the BP process yielded unsaturated phosphatidic acid (PA) 44:4 (PVE 0.70, Figure 2B), followed by saturated ether-linked PC such as PC(O-30:1/P-30:0), PC(O-30:0) and PC(30:0), all with a PVE of 0.69 as most important molecules. All the saturated plasmanyl/plasmenyl and diacyl PC molecules and PA(44:4) were included in the panel of 35 lipids reported in Table 1, and will be discussed later in the text. In the next sections, we will compare the different roles of important lipids (through LiOn enrichment analysis) and the results of the four different statistical approaches here described. Aiming to maintain the maximum amount of meaningful information [76], it was convened to select the top 10% of molecules (threshold over which model performances substantially improved, see Methods Supplement 2d) as important for RF (≈83 molecules for RFclass and RFreg, Table C1).

3.2. Comparison of the Outputs Obtained by the Different Statisatical Approaches

The four different statistical models applied in this study allowed the selection of a total of 123 relevant lipid molecules (respectively 25, 40, 72 and 92 via Volcano plot, OPLS-DA, RFclass and RFreg, Table A1) differentiating ripe from not ripe BMO from 838 molecules present in the PIT (Supplementary Data 2). It should be noticed that none of the PC, glycerophosphoethanolamines (PE), glycerophosphoinositols (PI) and triglycerides (TG) detected by the present statistical approach was a major component of the blue mussel gonad lipidome [50] supporting also the previous observation made on the PCA analysis plots.

Key features highlighted by the different approaches are compared using a Venn diagram (Figure 3). RFreg and RFclass yielded the largest number of unique molecules (respectively 41 and 12 molecules, i.e. 33.3% and 9.8% of the total), while only 4 variables were exclusive to OPLS-DA (3.3% of the total) and 2 variables were detected only by volcano analysis (1.6%). Ten variables were common for all four approaches (8.1%). The comparison between classification (Volcano, OPLS-DA and RFclass) and regression approaches (RFreg) shows that 51 (41.5% of the total) molecules (out the 92 considered important and highlighted by RFreg) was common with one of more of the models here evaluated. This, again, indicates that there was a large consensus between GVF and the classification criteria used to define the two BMO groups.

The functional role of the lipids that differentiated undeveloped from ripe BMOs was investigated via LiOn enrichment analysis (Figure 4). Figure 4A shows the PCA analysis for LiOn terms, which reveals a clear distinction of the function of those molecules since 84% of the samples (43/51) could be correctly clustered according to the functional role of the highlighted lipids. Terms related to storage lipids, lipid bodies and TG were enriched in ripe BMOs, which could be attributed to vitellogenesis and accumulation of lipid bodies in the yolk granules [16]. Interestingly, TG species, identified by the various statistical methods, showed inverse pattern according to their unsaturation levels: TG molecules bearing highly unsaturated acyl residuals (DB >10) resulted inversely correlated with GVF, and also more abundant in undeveloped BMOs, opposite pattern was observed for TG with <10 DB (Table A1). This observation could be the result of rearrangment and/or resorbition of PUFA during ovarian maturation/regression [85]. Although increased content of unsaturated TG was reported in ripe ovarian follicles of P. marguensis [49], others reported reported dominance of saturated and monounsaturated acyl moyeties in TG of mud crab ovaries, Scylla paramamosain [86] and sea sponges [87].

In the not-ripe BMO group was instead observed a significant up-regulation of ‘membrane components’, an enrichment terms associated with membrane lipids (Figure 4B). The previous term is mainly linked to the larger abundance in undeveloped ovaries of diacyl and ether-linked (associated with the term ‘plasmalogen’) PC and PE molecules. This observation could be related to the greater surface area of the gonad covered by support cells, such as ADG and VCT, that are characteristic of early maturation stages of bivalve ovaries [15], which is also supported by the reduction of GVF found in non-ripe BMOs compared to mature ones (Figure A1).

At the same time many of such membrane lipids included saturated and short chained molecules, features connected with the up-regulation of terms including ‘below average transition temperature’ and ‘below average bilayer thickens’. On the contrary, longer and more unsaturated PC and PE were more abundant in ripe BMOs, leading to the up-regulation of ‘PUFA’ and ‘fatty acids with more than 18C’ terms for this group.

In addition, terms related to ‘lysosomes’, the ‘Golgi apparatus’ and ‘mitochondrion’ were also upregulated in the immature BMO group. The former two of these terms are connected with sphingolipids (Ceramides - Cer, CAEP and CerPE), since these classes of lipids are produced in the Golgi apparatus and act in autolysis and atretic processes [88]. The larger content of such lipids in undeveloped BMOs (Table A1) could be related to gonadal degeneration and atresia, as the presence of lysosomes and vacuoles together with macrophage infiltration into follicles characterises atresia in mussel ovaries [28]. Interestingly, opposite patterns could be evidenced between Cer, higher on ripe BMOs and correlated with GVF, CAEP and CerPE. Such pattern could be related with the role of Cer as stabilizers in adipose tissue, and hence could be again related to accumulation of yolk reserves in ripe ovaries [89]. The mitochondrion term is associated to the detection by RFreg of a cardiolipin (CL) 88:21 which resulted also weakly inversely correlated with GVF (Spearman rank r=-0.05, Table A1). The reason behind this observation is unclear, since the same molecule did not evidence a significantly different content between the two BMO groups (FDR p>0.05). CL are fundamental lipids found in the inner mitochondrial membrane, participating in cellular respiration [90]. Highly unsaturated CL were found predominant in mussel gonads [50] and generally distinguish mussels from other groups of bivalves such as clams [91]. CL(88:21) must contain four 22C acyl residuals, since 24C FA detected in blue mussel gonads were saturated or monounsaturated species [50]. Therefore, it is likely that this molecule included three DHA and one FA22:3 non-methylene interrupted dienoic (NMID). Nevertheless, the connection of the aforementioned molecule with mussel ovarian development is unclear, if we exclude the structural role of NMID FA in mussel gonads [92]. An increase of CL(88:24) was observed in quiescent scallops (Placopecten magellanus), whereas larger extents of CL containing 18C FA (FA18:0 and FA18:1) were found in ripe (pre-spawning) individuals [93]. Given the different biology of the two species (dioecious for mussels and hermaphrodite for scallops) it becomes more difficult to directly compare our result and the previous study.

3.4. Ranking the Panel of Potential Biomarkers

ROC curves were generated for all important variables selected by each statistical method to distinguish ripe from not ripe BMOs and these were ranked using the AUC (Table A1). ROC curves are the graphical representation between the sensitivity (true positive) and specificity (true negative) of a biomarker [94]. AUC can range between 0.5 (diagonal line, random effect, 50% specific and 50% sensitive) and 1 (perfect biomarker, 100% specific and sensitive). The optimal cut-off value (maximum specificity and sensitivity) of a specific biomarker is calculated using the maximum distance between the ROC curve and the diagonal line (Youden index [95]). Molecules with AUC values between 0.8 and 0.9 have been considered as potential biomarker candidates for various human diseases in the past [96,97]. Out of all the important variables highlighted by the statistical approach used in this study, 35 had an AUC value >0.8 and are listed in Table 1. From the AUC>0.8 variables (35/123), 10 were detected by all four statistical approaches, 17 by a combination of three of these approaches, 6 by two models and just 2 by a single model. This observation emphasises the importance of combining multiple orthogonal statistical approaches, as any biomarker discovery method will inevitably have some limitations if unaccompanied.

In the present study, no lipid was found to be a perfect marker of the two different gonadal stages (e.g. AUC =1). One molecule showed an AUC >0.9: CerPE(40:2), which had and AUC of 0.905 (95%CI: 0.799-0.981, Figure 5A) and was more abundant in undeveloped ovaries (Figure 5B). CerPE(40:2) was characterised by a sensitivity of 80% (true positive rate) and a specificity of 100% (true negative rate). Three other CerPE ranked among AUC>0.8 molecules (Table 1): CerPE(35:3/2), CerPE(34:3)a and CerPE(40:1), with AUC values of 0.874, 0.845 and 0.806, respectively. Other ceramide lipids important for distinguishing between ripe and immature BMO included CAEP 42:2 and 42:3 with AUC values of 0.808 and 0.804, respectively. All the previous CerPE and CAEP molecules were minor components of these classes detected in Mytilus gonads [50] and were all found mainly in undeveloped BMO; the role of these sphingolipid classes in marine invertebrates is not yet fully understood. CerPE and CAEP are mainly found in invertebrate membranes [55,98,99,100,101] and are synthesised from ceramide precursors by the action of a specific CDP ethanolamine:ceramide phosphotransferase localised in the Golgi apparatus [102]. In the scientific literature there are scarce reports on the role of CerPE and CAEP in bivalve ovarian maturation. Pauletto and coworkers [103] reported a differential expression for genes involved in complex ceramide metabolism between ovarian oocytes and spawned oocytes in the clam Venerupsis decussatus. Other authors hypothesised a similar function for CerPE and CAEP of the vertebrate analogue sphingomyelin (SM), which, can form specific domains of the cell membrane defined as membrane rafts [101]. The decrease in CerPE and CAEP observed in ripe females could also be a consequence of the gamete maturation process, as a decrease in SM was observed in Xenopus laevis during oocyte maturation under the influence of progesterone [104].

Glycerophosphoserines (PS) represented the lipid class that, overhall, changed mostly between the two BMO groups (Table 1). PS(38:5) was the second most important biomarker of BMO maturity (AUC 0.894) and several of the other nine PS molecules with AUC>0.8 (36:2, 38:5, 38:6a, 39:6, 40:2, 40:5, 40:6 and O-38:5/P-38:4 and O-38:6/P-38:5) are between the main diacyl and plasmanyl/plasmenyl PS species found in Mytilus gonads [50]. All AUC>0.8 PS were more abundant in undeveloped BMOs. PS are secondary components of cell membranes, found mainly on the cytosolic side of the lipid bilayer [105]. Their relevance as markers of BMO development could be linked to different reasons. PS may be converted into different glycerophospholipids, by the action of phosphotransferases, during ovarian maturation, as observed in largemouth gudgeon (Coreius guichenoti) ovaries [106]. On the contrary, the larger content of PS in undeveloped BMOs could be also linked to apoptotic processes in degenerating gonads, since the exteriorisation of PS is considered to be one of the “triggers” for phagocytosis of apoptotic cells [107,108], as observed in UV exposed O. edulis haemocytes [109] and in apoptotic ovarian cells of Drosophila melanogaster [110,111].

The third molecule for relevance as biomarker evidenced in this study is the ether-linked PE O-36:2/P-36:1c (AUC: 0.875; Table 1). All PE molecules with AUC>0.8 were plasmanyl/plasmenyl forms (O-34:3/P-34:2, O-36:1/P-36:0b, O-36:2/P-36:1b, O-36:3/P-36:2, O-37:2/P-37:1, O-37:3/P-37:2, O-37:4/P-37:3 and O-39:3/P-39:4) but did not include any of the principal PE forms described in blue mussel gonads [50], increasing the complexity of explaining they physiological motivation behind our observation. All the AUC>0.8 PE resulted mainly abundant in undeveloped BMOs, and this could possibly be related with synthesis of PC, during ovarian development, thorugh sequential methylation of PE [112]. A total of seven PC molecules had an AUC >0.8; six of which were saturated molecules (28:0, 30:0 and the plasmanyl/plasmenyl O-30.0, O-31:0 and O-32:0; Table 1) and resulted distincly abundant in not ripe BMOs. The long chained and highly unsaturated ether-linked PC(46:10/P-46:9)a characterised, instead, ripe BMOs. This trend is in agreement with what reported by other authors on the PC composition of ovaries in the shrimp P. marguiensis: saturated PC molecules were found in undeveloped regions of the ovary, whereas highly unsaturated PC (containing 20-22C PUFAs) had an higher coverage ripe ovarian parts [49]. PC and PE are the two major membrane lipids in bivalves [21,50,66,113,114]. Given the low unsaturation of PC and PE with AUC>0.8 and abundant in not ripe BMOs, it is likely that these molecules carried saturated or monounsaturated FA rather than long-chain PUFA. Acyl residues of membrane lipids represent energy sources during embryonic development in invertebrates [115,116], so their enrichment with PUFA at later developmental stages may be important for the proper development of the offspring. To such extent, PUFA enrichment patterns have been observed for PC and PE fractions in maturating scallop ovaries [44].

Lastly, four PA molecules (38:6, 40:6, 44:4 and the plasmalogen O-38:6/P-38:5), one glycerophosphoglycerol (PG) 34:1 and the lyso-glycerophosphocholine (LPC) 22:2 completed the panel of 35 molecules with an AUC >0.8 in distinguishing ripe from not ripe BMOs (Table 1). All these lipid species were more abuntand in undeveloped ovaries, with the exception of PA(44:4) and LPC(22:2). There are not many reports on the role of such minor glycerophospholipid classes in bivalves or invertebrates ovarian maturation, due also to their low abundance which is hardly observable through traditional lipid analysis methods [48]. Therefore, we can only speculate about their role in the development of the BMO. PA and PG are intermediates of glycerophospholipid biosynthesis, each forming the precursors of the major membrane lipid classes and of cardiolipins [105]. PA concentration increases in Xenopus oocytes as response to insulin treatment, as a result of phospholipase D, which produces lyso-glycerophosphate (LPA) which participate in lipid signalling of oocyte maturation [112]. A reduction in most of PG was instead observed in sea sponges in reproductive individuals [87], observation that could agree with data discussed in the present study. LPC are metabolites of PC formed by the Lands cycle by the action of a phospholipase that cleaves one of the FA residues [117]. LPC(22:2) is a molecule that is likely to bear a 22:2 NMID FA [50]. NMID FA are endogenous lipids of mussels which are preferentially incorporated in the mantle, competing with C20 and C22 PUFA under reproductive demands [92], which can suggest the a role for LPC(22:2) in restructuring of membrane lipids in ripe BMOs. Nevertheless, especially for these minor lipids, further evidences are needed to obtain a definitive overview of the system.

4. Conclusions

Rapid and informative tools for measuring the maturity stage of female mussels are useful both for industrial applications, as they can be used to monitor gamete maturation in hatcheries to reduce problems with uncontrolled spawning and pre-spawning atresia, and for monitoring the reproductive status of native mussel populations. The high resolution and ability to work at the lipid molecule level make LC-MS lipidomics an important platform for studying lipid biomarkers that can assess maturation status in BMOs. To fulfil this task, we applied different supervised statistical approaches that were more powerful than non-supervised approaches (PCA) in highlighting the small differences between the BMO groups studied. By combining the four supervised statistical approaches used here (Volcano, OPLS-DA, RFclass and RFreg), it was possible to highlight a panel of 35 lipid species that differ to a high degree (AUC>0.8) between ripe and undeveloped BMOs. These potential markers include CerPE, CAEP, LPC, PA, PC, PE, PG and PS. In many cases, however, the exact role and reason for the importance of AUC>0.8 lipids for classification could only be speculated upon similarities reported by others in different model organisms. Future studies should aim to better understand the role played by lipids during gamete development in marine invertebrates such as bivalves, focussing on less studied lipid classes.

The set of potential lipid biomarkers should be validated and possibly quantified using further independent gonad samples, and if confirmed, the protocol should be improved and refined to develop a non-lethal assay. GHA is considered a traditional approach for monitoring gonadal development in marine mussels [118] and is an important tool for investigating the reproductive physiology, toxicology, life cycle and population structure of bivalves [10,11,16,28,32,34,37,119]. In the future, GHA will continue to be a necessary diagnostic tool, but technological advances and the application of state-of-the-art high-throughput omics techniques should support and complement traditional techniques to increase our knowledge of the physiology of marine organisms and their ability to adapt to a changing environment.