Insights from a Multi-Omics Integration (MOI) Study in Oil Palm (Elaeis guineensis Jacq.) Response to Abiotic Stresses: Part One—Salinity

Oil palm (Elaeis guineensis Jacq.) is the number one source of consumed vegetable oil nowadays. It is cultivated in areas of tropical rainforest, where it meets its natural condition of high rainfall throughout the year. The palm oil industry faces criticism due to a series of practices that was considered not environmentally sustainable, and it finds itself under pressure to adopt new and innovative procedures to reverse this negative public perception. Cultivating this oilseed crop outside the rainforest zone is only possible using artificial irrigation. Close to 30% of the world’s irrigated agricultural lands also face problems due to salinity stress. Consequently, the research community must consider drought and salinity together when studying to empower breeding programs in order to develop superior genotypes adapted to those potential new areas for oil palm cultivation. Multi-Omics Integration (MOI) offers a new window of opportunity for the non-trivial challenge of unraveling the mechanisms behind multigenic traits, such as drought and salinity tolerance. The current study carried out a comprehensive, large-scale, single-omics analysis (SOA), and MOI study on the leaves of young oil palm plants submitted to very high salinity stress. Taken together, a total of 1239 proteins were positively regulated, and 1660 were negatively regulated in transcriptomics and proteomics analyses. Meanwhile, the metabolomics analysis revealed 37 metabolites that were upregulated and 92 that were downregulated. After performing SOA, 436 differentially expressed (DE) full-length transcripts, 74 DE proteins, and 19 DE metabolites underwent MOI analysis, revealing several pathways affected by this stress, with at least one DE molecule in all three omics platforms used. The Cysteine and methionine metabolism (map00270) and Glycolysis/Gluconeogenesis (map00010) pathways were the most affected ones, each one with 20 DE molecules.


Introduction
Oil palm (Elaeis guineensis Jacq.) has the highest productivity among the main oilseed crops, reaching 3-8 times more oil per area than any other crop [1]. In 2021/2022, approximately 82 million metric tons of palm oil and palm kernel oil was consumed worldwide, making oil palm the number one source of consumed vegetable oil in the world [2]. It is the raw source of a wide range of products used by many industries, such as the food and steel industries, the manufacture of cleaning products, the pharmaceutical and cosmetics industries, and the biofuels industry [3].
Several countries placed, in the equatorial belt, expanded oil palm plantations in tropical forests where this crop meets its natural condition of high rainfall throughout the year [4]. Despite the significant economic gains, this movement imposes environmental stresses, such as biodiversity loss, greenhouse gas emissions, land degradation, and air and water pollution [1]. In Brazil, over 95% of the oil palm plantations are in the Amazon rainforest, where only 2.14% of the total area destined for commercial plantations is currently in use [5]. This under-utilization status is due to many constraints, such as environmental legal restrictions imposed by the Brazilian Government and logistical difficulties, which together hinder the production flow to the main industrial areas in the country and also the occurrence of pests and diseases [6,7].
Outside the Amazon rainforest, there is an extensive area in Brazil with favorable conditions for cultivating oil palm [8]. There are many logistic advantages to producing oil palm outside the Amazon region, offering a window of opportunity to increase the area with oil palm plantations and, consequently, the total national palm oil yield. However, these areas experience long periods of drought throughout the year when oil palm does not meet the physiological water requirement to maintain productivity [8][9][10]. Consequently, the oil palm grower needs to irrigate the crop and must do so with proper management to avoid soil salinization.
Approximately 30% of the world's irrigated agricultural lands are damaged by salinity, negatively affecting the productivity of traditional crops [11]. Most crop plants have evolved under very low soil salinity, and, under high salt, their development is highly inhibited, even leading to death at a concentration ranging between 100 and 200 mM NaCl [12]. Salinity stress affects plants by decreasing the osmotic potential of the soil solution, making it harder for the root to absorb water from the soil and consequently experiencing drought stress, and by accumulating sodium and chloride ions in the cytoplasm, leading to the inhibition of many enzyme reactions due to ion toxicity [13]. Salt stress tolerance in plants involves many morphophysiological and biochemical changes, such as stomatal closing, osmolyte accumulation, and increased Na + /Cl − antiporter, governed by multigenic traits [14].
Considering those circumstances, it is clear that any initiative to promote oil palm cultivation outside the Amazon Forest in Brazil needs to take drought and salinity together when researching for knowledge and technology to empower breeding programs to develop superior genotypes for those regions. The first challenge is understanding the morphophysiological, biochemical, and molecular responses of oil palm to these two abiotic stresses. In doing so, our group has studied the response of young oil palm plants from different angles, applying different omics platforms, alone or in combination [6,15,16]. Vieira and colleagues showed that young oil palm submitted to a high concentration of NaCl reduces the rates of CO 2 assimilation, stomatal conductance to water vapor, and transpiration [6]. Then, Ref. [15] confirmed a preponderant role of transcription factors in the early response of oil palm plants to salinity stress, and [16] identified the pathways and the metabolites most affected by drought stress.
The current study is a new step on our research activities characterizing the biochemical and molecular responses of E. guineensis to salinity stress, where we carried out a comprehensive, large-scale, single-omics analysis (SOA), and Multi-Omics Integration (MOI) analysis of the metabolome, transcriptome, and proteome profiles on the leaves of young oil palm plants submitted by Vieira and colleagues [6] with repsect to very high salinity stress.
The set of 693 enzymes underwent gene ontology analyses, and only the ten most populated groups per GO term are shown in Figure 1. The biological process subgroups with the largest number of representatives were carbohydrate metabolic process, followed by protein phosphorylation and fatty acid biosynthesis process. For molecular function, the most populated subgroups were ATP binding, metal ion binding, and heme binding. Finally, for cellular component the integral component of membrane came in first, followed by cytoplasm and cytosol components.  Furthermore, the DE enzymes were also classified according to the Enzyme Commission (EC) number, a numerical classification scheme for enzymes based on the chemical reaction. At a first level of classification that involves a general type of enzyme-catalyzed reaction that ranges from one to six, enzymes were dominated by oxidoreductases (EC 1), transferases (EC 2) and hydrolases (EC 3) classes (Figure 2a). In the subclasses of oxidoreductases class (EC 1), DE enzymes were represented mainly by those acting on paired donors, with the incorporation of or reduction in molecular oxygen (EC 1.14), followed by enzymes acting on the CH-OH group of donors (EC 1.1) and acting on the aldehyde or oxo group of donors (EC 1.1). The most representative subclass of transferases class (EC 2) included those with a function of transferring phosphorus-containing groups (EC 2.7), acyltransferases (EC 2.3) and glycosyltransferases (EC 2.4). Finally, the hydrolases (EC 3) had subclasses with compounds involved and acting on ester bonds (EC 3.1), glycosylases (EC 3.2), and acting on peptide bonds (peptidases) (EC 3.4) subclasses that were standing out.

Oil Palm Proteome under Salinity Stress
A global proteomics analysis led to the identification of 3234 and 2872 peptides with high confidence (FDR ≤ 0.01) in control and stressed samples, respectively, which infers up to 1809 protein entries from E. guineensis proteome (Uniprot) in both conditions-control and stressed (Table 2). Approximately 38% of proteins (688) were inferred from more than four peptides, and about 34% (622) had at least one proteotypic peptide observation. A list of all peptides and proteins confidently identified, as well as a simplified list of 1019 proteins according to the maximum parsimony criterion, is presented in Supplementary Tables S2-S4   This group of 171 DE protein sequences-including those found exclusively in control (62) and stressed (89) in at least two replicates and those 20 proteins that attended to the statistical criteria of PatternLab V software-was submitted to functional annotation and MOI analyses. The KEGG mapper reconstruction results revealed 131 proteins with 84 distinct K numbers, including 99 enzymes. Seventy-three enzymes belonged to known pathways and were used in the MOI analysis.
This set of 171 selected proteins was then submitted to gene ontology analyses, and again only the ten most populated groups per GO term are shown in Figure 1. The biological process subgroups with the largest number of proteins were translation, followed Twenty proteins showed statistically significant differences in their abundance between stressed and control samples ( Table 3). As shown in Figure 3b, 12 proteins (in blue) were significantly up-regulated while eight (in red) were significantly down-regulated between stressed and control samples. A group of 642 proteins did not meet the statistical criteria and was not considered for further analysis (Table 1). Our differential abundance analysis considered proteins identified at least in two replicates in each condition. This filtering process decreased the list of 316 and 380 proteins uniquely identified in stressed and control to 89 and 62, respectively. This group of 171 DE protein sequences-including those found exclusively in control (62) and stressed (89) in at least two replicates and those 20 proteins that attended to the statistical criteria of PatternLab V software-was submitted to functional annotation and MOI analyses. The KEGG mapper reconstruction results revealed 131 proteins with 84 distinct K numbers, including 99 enzymes. Seventy-three enzymes belonged to known pathways and were used in the MOI analysis.
This set of 171 selected proteins was then submitted to gene ontology analyses, and again only the ten most populated groups per GO term are shown in Figure 1. The biological process subgroups with the largest number of proteins were translation, followed by carbohydrate metabolic process, fatty acid biosynthetic process, proteolysis, pentosephosphate shunt, oxidative branch, and glucose metabolic process. For molecular function, the proteins were mainly distributed in the subgroups of ATP binding, structural constituent of ribosome, and metal ion binding. Finally, the cellular component of the cytoplasm came in first, followed by the nucleus and cytosol.
The prevalent chemical reactions by which proteins were classified according to EC were Oxireductases (EC 1), transferases (EC 2), and hydrolases (EC 3) classes ( Figure 2b). In the subclasses of oxidoreductases, the main groups were acting on the CH-CH group of donors (EC 1.3), acting on Ch or CH(2) groups of donors (EC 1.17), and acting on the aldehyde or oxo group of donors (EC 1.2). The most representative subclasses of transferases class included those with a transferring phosphorus-containing groups (EC 2.7) and transferring nitrogenous groups (EC 2.6). For hydrolases, those acting on carbonnitrogen bonds (EC 3.5) and acting on acid anhydrides (EC 3.6) came first.

Oil Palm Metabolome under Salinity Stress
Statistical analysis on Metaboanalyst returned 2843 and 1855 peaks, respectively, in the polar-positive and polar-negative fractions (Table 1). Fifty-two peaks were differentially expressed, and eighteen were up-regulated and thirty-four were down-regulated in the polar-positive while seventy-seven were differentially expressed in the polar-negative, in which nineteen were up-regulated and fifty-eight were down-regulated.
All 129 peaks differentially expressed were then submitted to functional interpretation via analysis in the MS Peaks to Pathway module, and the combined mummichog and GSEA pathway meta-analysis resulted in a list of 19 differentially expressed metabolites (DEMs), which was then submitted to the pathway topology analysis module ( Table 4). The monobactam biosynthesis (map00261); arginine biosynthesis (map00220); beta-alanine metabolism (map00410); pentose phosphate pathway (map00030); carbon fixation in photosynthetic organisms (map00710); alanine, aspartate and glutamate metabolism (map00250); galactose metabolism (map00052); and glutathione metabolism (map00480) pathways came out as the one with a raw p ≤ 0.05 ( Figure 4).   The metabolome view resulted from the analysis in the Pathway Topology Analysis module using the Hypergeometric test, the relative betweenness centrality node importance measure, and the latest KEGG version of the Oryza sativa pathway library. Pathway impact takes into account both node centrality parameters-betweenness centrality and degree centrality-and represents the importance of annotated compounds in a specific pathway.

Integrating Oil Palm Transcriptome, Proteome and Metabolome
A total of 510 enzymes (436 from transcriptomics analysis and 74 from proteomics analysis) (Supplementary Table S8) and 19 metabolites from metabolomics analysis (Table  3), all selected as differentially expressed in the leaves of young oil palm plants (stressed/control), were submitted to MOI analysis.
By applying the Omics Fusion platform to perform the MOI analysis, results revealed a group of eleven pathways affected by salinity stress, and with at least one molecule dif- The metabolome view resulted from the analysis in the Pathway Topology Analysis module using the Hypergeometric test, the relative betweenness centrality node importance measure, and the latest KEGG version of the Oryza sativa pathway library. Pathway impact takes into account both node centrality parameters-betweenness centrality and degree centrality-and represents the importance of annotated compounds in a specific pathway.

Integrating Oil Palm Transcriptome, Proteome and Metabolome
A total of 510 enzymes (436 from transcriptomics analysis and 74 from proteomics analysis) (Supplementary Table S8) and 19 metabolites from metabolomics analysis (Table 3), all selected as differentially expressed in the leaves of young oil palm plants (stressed/control), were submitted to MOI analysis.
By applying the Omics Fusion platform to perform the MOI analysis, results revealed a group of eleven pathways affected by salinity stress, and with at least one molecule differentially expressed in each one of the three omics platforms used ( Table 5). The Cysteine and methionine metabolism (map00270) and the Glycolysis/Gluconeogenesis (map00010) pathways came tied first in this list, each one with 20 unique molecules from the transcriptome/proteome/metabolome integrative analysis (Supplementary Table S9).

Discussion
Soil salinization reduces plant growth and productivity of most terrestrial crops with economic importance, including oil palm [13]. In the case of oil palm, Refs. [6,18,19] reported the development of salinization protocols, which are necessary to study the response of E. guineensis to this abiotic stress in the search for some intraspecific trait variability. Such protocols are also needed to select tolerant oil palm genotypes developed via genetic engineering or genome editing strategies. These two studies generated not only morphophysiological, biochemical, and molecular insights into the response of this species to this abiotic stress but also reported on the ionic imbalance in the substrate, roots, and leaves of young oil palm plants under salinity stress.
Salinity tolerance is a multigenic trait that governs physiological, biochemical, and molecular mechanisms to facilitate water retention and/or acquisition, protect chloroplast functions, and maintain ion homeostasis [13]. Datasets in genomics, transcriptomics, proteomics, metabolomics, epigenomics, ionomics, and phenomics are accumulating everywhere, intending to gain insights into the mechanisms behind plant interaction with abiotic stresses; however, due to the molecular complexity of such interaction, single-omics analyses (SOA) will have limited power in delivering a more systemic and accurate picture of those responses. Multi-Omics Integration (MOI) strategies [20] are a new window of opportunity facilitating hypothesis generation, leading to the non-trivial challenge of unraveling the mechanisms behind this multigenic trait.
In the present study, SOA showed that carbohydrate metabolism and translation were the most affected biological process subcategories for differentially expressed genes and proteins, respectively. In addition to being a substrate for energy production, carbohydrates play a role in plant stress perception and signal transduction and can also mediate osmotic regulation and carbohydrate distribution [21].
From the translation point of view, the protein synthesis machinery is quite sensitive to salt since the production of new proteins is crucial for salinity tolerance [22]. Salinitytolerant species have a more efficient system for regulating transcription, synthesis, and protein processing when compared to sensitive species [23]. Genes encoding the plastid translation machinery in Arabidopsis thaliana are salt responsive, indicating a possible role in supporting chloroplast functionality [23]. Reaumuria soongarica (Pall.) Maxim., a salttolerant species, showed a complex pattern of protein expression, mainly those involved in translation, ribosomal structure, and biogenesis [24].
The dominant molecular function for DE enzymes identified in this study was ATPbinding proteins. These enzymes use the energy of ATP hydrolysis to catalyze a series of chemical reactions [25]. Most ATP-binding proteins are intracellular and extracellular transmembrane proteins, participating in the movement of various molecules and, under stress conditions, in intracellular osmotic balance maintenance [26]. ABC transporters that constitute one of the most populated families of proteins driven by ATP hydrolysis revealed a complex expression pattern in E. guineensis under drought stress, pointing to their role in controlling the influx and efflux of chemical molecules while in water scarcity [27].
Integral membrane components, the most affected subcategory of cellular components, include, in this category, proteins incorporated into cell membranes. Salt stress causes damage to the cell membrane, altering its permeability, lipid composition, and enzyme activity [28]. Several factors cause changes in the structure of cell membrane components during salt stress; among them, the excessive production of reactive oxygen species (ROS) is highlighted, which causes conformational changes in membrane proteins and lipid peroxidation, reducing the efficiency of transport systems and increasing membrane permeability [29].
In salinity-tolerant cultivars, an increase in the antioxidant defense system occurs, reducing lipid peroxidation and maintaining adequate levels of membrane permeability. On the other hand, there is an increase in the leakage of electrolytes from the membranes, which indicates a loss of membrane integrity in the sensitive plants [29]. The increase in electrolyte leakage has already been evident in oil palm leaves, indicating possible damage caused by salinity, with direct consequences in photosynthetic capacity reduction and biomass accumulation [18].
In the present study, an attempt to integrate three distinct omics platformstranscriptomics, proteomics, and metabolomics-was reported for the first time to gain further insights into the mechanisms behind the response of young oil palm plants to salinity stress. The MOI strategy used in the present study is a pathway-based approach for integrating omics datasets. Such integration was only possible due to the selection and characterization of salt-responsive genes coding for enzymes in the oil palm reference genome. Enzymes catalyzing reactions in a metabolic pathway are the bridges to connect transcriptomics, proteomics, and metabolomics datasets in such an integratomics study.
The present MOI study revealed eleven pathways affected by the salinity stress in the leaves of the young oil palm plant, with at least one molecule differentially expressed in all three platforms used. The Cysteine and methionine metabolism (map00270) and Glycolysis/Gluconeogenesis (map00010) pathways were the most affected ones. Even though this study identified other pathways, further discussion will concentrate only on these two.
Reactions that promote cysteine (Cys) biosynthesis are involved in the pathway of cysteine and methionine metabolism [30]. Cys acts as a sulfur donor for the biosynthesis of many essential bio-molecules, such as methionine, vitamins, co-factors, and Fe-S groups, and for the production of glutathione (GSH), considered the principal determinant of cellular redox homeostasis [30]. The enzymes serine O-acetyltransferase (EC 2.3.1.30) and cysteine synthase (EC 2.5.1.47) usually carry out the Cysteine biosynthesis in two steps. These enzymes are highly conserved in plants and are responsible for maintaining homeostasis between cysteine consumption and sulfate reduction [31].
Among the enzymes integrated into the cysteine and methionine metabolism pathway, seven and twelve were up-and down-regulated in the leaves of young oil palm plants under salinity stress, respectively. Serine O-acetyltransferase experienced an approximately 11-fold increase in expression, while L-lactate dehydrogenase experienced a decrease of about 90% in its original expression level. As cysteine is the first organic compound in the primary metabolism of sulfate, the elevated transcription of serine O-acetyltransferase may indicate that sulfate entry into the pathway plays a role in the saline stress response in oil palm. In tobacco, plants over-expressing bacterial serine O-acetyltransferase conferred resistance to high levels of oxidative stress with a four-fold higher cysteine expression [32]. Recently, Ref. [33] demonstrated that the exogenous application of nitric oxide (NO), a compound that regulates the response to different stresses in plants, increased the content of enzymes synthesizing Cys, helping maintain the cellular homeostasis in plants under the osmotic tension.
Amino acid methionine has nutritional value for plants, participating in the initiation of translation, in addition to being a precursor of S-Adenosyl methionine (SAM), the donor of the methyl group that regulates different essential cellular processes, such as cell division, synthesis cell wall, chlorophyll synthesis, and membrane synthesis [34]. SAM is synthesized from adenosine triphosphate (ATP) and methionine by the enzyme S-adenosylmethionine synthetase-SAMS (EC 2.5.1.6). The present study showed that this enzyme had a 2.6 fold increase in expression under salinity stress. Overexpression of the SsSAMS2 gene from the halophyte plant of Suaeda salsa L. in transgenic tobacco plants enhanced salt tolerance, as indicated by maintaining a higher photosynthetic rate and accumulation of more biomass [35].
GSH is a low molecular weight thiol crucial for maintaining the regulation of cellular redox homeostasis [30]. Two ATP-dependent enzymes, glutamate-cysteine ligase (EC:6.3.2.2) and GSH synthetase (EC:6.3.2.3), catalyze GSH synthesis from cysteine, glutamate, and glycine [36]. In the present study, the metabolite glutathione (C00051) upregulated 2.9 fold while the enzyme glutathione synthase downregulated to 70% of its original levels in the leaves of young oil palm plants under salinity stress. Meanwhile, all versions of glutamate-cysteine ligase found in the reference genome of oil palm [37,38] were non-DE. The exogenous application of GSH reversed the effects of salt stress on seedlings of tomatoes, as well as the expression and activities of enzymes related to the synthesis and metabolism of GSH, including gamma-glutamylcysteine synthetase (γ-ECS) and glutathione synthetase (GS), among others [39].
The glycolysis pathway directly supplies energy to plant cells from reactions that oxidize hexoses to produce ATP and pyruvate, the latter acting as a substrate for entry into the citric acid cycle (TCA). Conversely, the gluconeogenesis pathway synthesizes hexoses using low molecular weight compounds to meet energy needs under the conditions of reduced glucose supply [40]. The ATP-dependent 6-phosphofructokinase 2 (PKF) and pyruvate kinase (PK) enzymes from the glycolytic pathway did downregulate in the leaves of young oil palm plants under saline stress. Those two enzymes, together with hexokinase, are regulators of glycolysis, as they participate in irreversible reactions [40]. The energy production via glycolysis plays a role in the saline stress response in plants as it provides ATP to support the stress condition [41].
That was evident in the study by [41], where salt stress inhibited the growth of Cucumis sativus L. with a significant reduction in ATP production rates [41] and applied exogenous putrescine (Put), reversing the saline stress with positive modulation in the PFK and PK levels. In halophyte species Bruguiera sexangular, both PFK and PK enzymes increased expression in response to long-term salinity [42]. This suggests that increased PFK and PK activity increases the activity of the glycolytic pathway to maintain normal physiological metabolism under saline stress conditions in halophyte species. Salinity stress possibly promoted a reduction in ATP production as it negatively affected enzymes in the flow of glycolysis and TCA in oil palm.
Fructose-1,6-bisphosphatase (FBPase) and phosphate dikinase (PPDK), enzymes of the gluconeogenesis pathway, were negatively regulated in oil palm under salt stress. Under saline stress conditions, the active synthesis of sugars by this route contributes to mitigating the osmotic stress effect resulting from the submission of plants to a saline environment. In maize (Zea mays L.), the photosynthesis rate was similar between control plants and plants under neutral salt stress, suggesting that gluconeogenesis acted on the active synthesis of sugars and the maintenance of osmotic balance [43].
The overexpression of TaFBA-10 in A. thaliana (L.) Heynh increased FAB activity with positive effects on scavenging ROS under cold stress, whereas chlorophyll content was severely affected [44]. The PPDK enzyme, in turn, is an enzyme involved in the regulation of the C4 pathway in plants. PPDK enzyme activity increases in salinity tolerant accessions of Miscanthus sinensis Andersson [45]. In this manner, the activity of this enzyme compensates for the suppression of the Calvin Cycle by saline stress.
Another enzyme that participates in energy production is L-lactate dehydrogenase type B (LDH). This enzyme did downregulate in young oil palm plants under salinity stress. Lactate dehydrogenase (LDH) converts pyruvate to lactate that regenerates NAD + to maintain cellular respiration under anaerobic conditions. Under flood stress conditions, the initiation of fermentation responds to keeping energy supply in hypoxia conditions [46]. The downregulation of the LDH enzyme in palm oil at 12 DAT indicates a deficiency in response to salinity stress and may indicate a possible anaerobic condition caused by it.
Vieira and colleagues showed that young oil palm plants are sensitive to high concentrations of NaCl [6]. The present MOI study, which used datasets generated from leaf tissue collected by Vieira and colleagues, showed that enzymes competing for energy production in the glycolysis and gluconeogenesis pathways were negatively affected by salinity stress in the leaves of young oil palm plants. Concomitantly, gluconeogenesis, which involves the synthesis of glucose from non-carbohydrate substances, apparently does not represent an immediate response to reduced glucose supply in this oilseed crop under such stress.
The samples for transcriptome, metabolome, and proteome analyses were collected at once, exactly 12 days after the onset of the stress, using a split sample study design. In terms of data integration, accordingly to [47], the ideal situation is to have samples originating from the same biological source material and obtained at the same time-a piece of tissue may be cut into several sections and one used for a specific omics platform analysis, whilst the other is used to another one. In such design, the samples are more similar in that they all are assumed to produce data without batch effects between the different omics data sets [47].
It is clear that the small number of DE metabolites was the main limitation of the pathway-based MOI approach used in this study. The biggest number of DE metabolites in the eleven most affected pathways identified in the MOI analysis was three, while there were up to six proteins and 17 transcripts. The possible main reason to the fact that only 19 metabolites were differentially expressed in the leaves of young oil palm plants submitted to salinity stress was the mmetabolomics approach used. The untargeted metabolomics is an exciting technology for searching for novel metabolic perturbations in various biological systems, allowing the profile of many hundreds or thousands peaks with varying chemical properties at once; however, there are still various obstacles, such as the limited capability to identify novel compounds of interest and the need for advanced and more robust databases [48]. In the present study, we used the latest KEGG version of the O. sativa pathway library.

Plant Material, Experimental Design and Saline Stress
The oil palm plants used in this study were clones regenerated out of embryogenic calluses obtained from the leaves of an adult plant-genotype AM33, a Deli × Ghana from ASD Costa Rica, as previously reported by [6]. Before starting the experiments, plants were standardized accordingly to the developmental stage, size, and number of leaves. They were in the growth stage known as bifid saplings, and the experiment was performed in March 2018 in a greenhouse at Embrapa Agroenergy in Brasília, DF, Brazil (S-15.732 • , W-47.900 • ). The main environmental variables (temperature, humidity, and radiation) fluctuated according to the weather conditions and underwent monitoring throughout the entire experimental period using the data collected at a nearby meteorological station (S-15.789 • , W-47.925 • ).
The experiment consisted of five salinity levels (0.0, 0.5, 1.0, 1.5, and 2.0 g of NaCl per 100 g of substrate (a mixture of vermiculite, soil, and the Bioplant commercial substrate (Bioplant Agrícola Ltda, Nova Ponte, MG, Brazil), in a 1:1:1 ratio, on a dry basis), with four replicates in a completely randomized design (for additional details, see [6]). The substrate mixture was fertilized using 2.5 g L −1 of the N-P2O5-K2O formula (20-20-20). For the omics (transcriptomics, metabolomics, and proteomics) analysis described in the present study, we collected the apical leaves from control and stressed plants (0.0 and 2.0 g of NaCl per 100 g of substrate) 12 days after imposition of the treatments (DAT).

Transcriptomics Data Analysis
Leaves harvested from control and stressed plants were immediately immersed in liquid nitrogen and stored at −80 • C until RNA extraction; three plants for treatments. Details regarding total RNA extraction and quality analysis, library preparation, and sequencing are in [15,19]. RNA-Seq raw sequence data are in the Sequence Read Archive (SRA) database of the National Center for Biotechnology Information (NCBI) under BioProject number PRJNA573093.
All the transcriptomics analysis was performed with OmicsBox platform-version 2.0.36 [49], as previously described by [15,17]. The oil palm genome [19,20]-downloaded from NCBI (BioProject PRJNA268357; BioSample SAMN02981535) in September 2021-was the reference genome for RNA-Seq data alignment. The pairwise differential expression analysis between experimental conditions (Stressed Plants-12 DAT X Control-12 DAT) was performed through edgeR software version 3.28.0 [50], applying a simple design and an exact statistical test without the use of a filter for low counts genes.

Proteomics Data Analysis
Leaves samples for proteomics analysis were harvested, immediately immersed in liquid nitrogen, and then stored at −80 • C until protein extraction; three plants for control and three from stressed plants. Approximately 5.0 g of ground tissue-with 0.02 g/g of PVP (polyvinylpolypyrrolidone) added to it-was weighed and mixed with 3.0 mL of buffer (50 mM Tris HCl + 14 mM β-mercaptoethanol, pH 7.5) and 30 µL of protease inhibitor. After gently stirring the suspension on ice for 10 min, it was centrifuged at 10.000 G at 4.0 • C for 15 min. Then, 1.0 mL of the supernatant was transferred to 2.0 mL microtubes, mixed with 1.0 mL of 10% TCA (trichloroacetic acid) solution in acetone, kept at −20 • C for 2 h for protein precipitation, and then centrifuged at 10,000 G at 4.0 • C for 15 min. The protein pellet underwent wash with ice-cold 80% acetone, followed by centrifugation under the same conditions as above. After washing twice, we stored it at −80 • C until protein quantification [51] and visualization in an SDS-PAGE Gel.
After protein quantification, all samples went to the GenOne company (Rio de Janeiro, RJ, Brazil) fort protein preparation and LC-MS/MS analysis. After undergoing treatment with 10 mM DTT at 56 • C for 30 min, followed by 40 mM iodoacetamide (IDA) at room temperature in the dark and also for 30 min. Then, samples were incubated for 20 h at 37 • C with trypsin (1:50) in a thermomixer at 800 rpm. At last, after adding 50 µL of 95% acetonitrile and 5% TFA, samples were stirred three times at 1000 rpm for 15 min for tryptic peptides extraction, vacuum dried, and dissolved in 20 µL of 0.1% formic acid in water.
For a global proteomics analysis, we adopted a label-free quantitation approach using spectral counting by LC-MS/MS passing the samples through a nano-high performance liquid chromatography (EASY 1000; Thermo Fisher, Waltham, MA, USA) coupled to Orbitrap Q Exactive Plus (Thermo Scientific, Waltham, MA, USA) mass spectrometer. An MS scan spectra ranging from 375 to 2000 m/z were acquired using a resolution of 70,000 in the Orbitrap. We used the Xcalibur software (version 2.0.7) (Thermo Scientific, Waltham, MA, USA) to obtain the data in biological triplicates.
The PatternLab for Proteomics V software [23] was used to process the raw files. We used the Comet algorithm [52], the E. guineensis Uniprot reference database (30,667 entries), and 123 common contaminant proteins (Proteome ID: UP000504607) to perform peptide sequence matching (PSM) and employed a target-reverse strategy to increase confidence in protein identifications [53]. The search considered semi-specific candidates and allowed a maximum of two missed cleavages. Fixed cysteine carbamidomethylation and variable methionine oxidation were applied; the Comet search engine used a precursor mass tolerance of 40 ppm and a fragment compartment tolerance of 0.02.
We employed the SEPro-Search Engine Processor-module of PatternLab [54] to validate the peptide spectrum matches and, subsequently, grouped identifications by enzymatic specificity (semi-specific), resulting in two distinct subgroups. Then, we applied XCorr, DeltaCN, Spectral Counting Score, and Peaks Matched values to generate a Bayesian discriminator. SEPro automatically establishes a cutoff score to accept a 1% false discovery rate (FDR) based on the number of decoys performed independently on each subset of data, resulting in a false positive rate independent of the triptych status. We chose a minimum sequence length of six amino acid residues and discarded similar proteins that represent an identical sequence and consist of a fragment of another one. At last, a final list of mapped proteins was composed only of PSMs with less than five ppm.

Metabolomics Data Analysis
Leaves harvested from control and stressed plants were immediately immersed in liquid nitrogen and stored at −80 • C until metabolite extraction: four plants for treatments. Before solvent extraction, all samples underwent grounding in liquid nitrogen. The solvents used were methanol grade UHPLC, acetonitrile grade LC-MS, formic acid grade LC-MS, sodium hydroxide ACS grade LC-MS, all from Sigma-Aldrich, and water treated in a Milli-Q system from Millipore. We employed a well-established protocol [16,55,56] to extract the metabolites in three phases (polar, non-polar, and protein pellet). Aliquots of 50 mg of ground sample were transferred to 2 mL microtubes, and then 1 mL of a mixture of 1:3 (v:v) methanol/methyl tert-butyl ether (MTBE) at −20 • C was added. Homogenization on an orbital shaker at 4.0 • C and ultrasound treatment in an ice bath were each performed for 10 min. As the next step, we added 500 µL of a mixture of 1:3 (v:v) methanol/water to each microtube. After centrifugation (15,300× g at 4.0 • C for 5 min), an upper non-polar (green) and a lower polar (brown) phase and a protein pellet remained in each microtube. After transferring both fractions separately to 1.5 mL microtubes, they were submitted to a Speed vac system (Centrivap, Labconco) to be vacuum dried. Finally, the dry-fraction, resuspended in 500 µL of 1:3 (v:v) methanol and water mixture and transferred to vials, were now ready for analysis.
High-resolution mass spectrometry was used for detection (MaXis 4G Q-TOF MS, Bruker Daltonics) equipped with an electrospray source in positive (ESI-(+)-MS) and negative (ESI-(−)-MS) modes. The settings of the mass spectrometer were as follows: capillary voltage, 3800 V; dry gas flow, 9 L min −1 ; dry temperature, 200 • C; nebulizer pressure, 4 bar; final plate offset, 500 V. For the external calibration of the equipment, we used a sodium formate solution (10 mM HCOONa solution in 50:50 v:v isopropanol and water containing 0.2% formic acid) injected through a six-way valve at the beginning of each chromatographic run. Ampicillin ([M+H] + m/z 350.1186729 and [M-H]-m/z 348.1028826) was the internal standard for later peak normalization on data analysis. DataAnalysis 4.2 software (Bruker Daltonics, Bremen, Germany) was the first used to analyze the raw data from UHPLC-MS, as mzMXL files. Pre-processing of data was performed using XCMS Online [57,58], including peak detection, retention time correction, and alignment of the metabolites. CentWave was used for peak detection (∆m/z = 10 ppm; minimum peak width, 5 s; maximum peak width, 20 s). For the alignment of retention times, the parameters were as follows: mzwid = 0.015; minfrac = 0.5; bw = 5. The unpaired parametric t-test (Welch t-test) was used for the statistical analysis at the pre-processing stage. Then, a data set was created from control (0.0 g) and stressed plants subjected to NaCl/100 g of the fresh substrate at 12 DAT. All with four biological repeats.
The pre-processed data (csv file) underwent analysis in the Statistical Analysis module of the MetaboAnalyst 5.0 [59,60]. The scaling used was the Pareto method [61]. Afterward, the differentially expressed peaks (DEPs) selected were those passing the criteria of false rate discovery (FDR) ≤ 0.05 and Log 2 (fold change (FC)) = 1. When using the MS Peaks to Pathway module to analyze the selected DEPs, we employed the following parameters: molecular weight tolerance of 5 ppm; mixed ion mode; joint analysis using both the mummichog [62] and Gene Set Enrichment Analysis-GSEA [63] algorithms; the latest KEGG version of the O. sativa pathway library. The p-value cutoff from the mummichog algorithm was at 1.0 × 10 −5 .
When two or more matched forms were observed as DEP (in the case of isotopes), the mass error was the criteria for the feature selection for the comparison with metabolite databases, keeping the smallest [56]. The mass error was also the criteria in the case of a single matched compound relative to two or more DEPs. The mass spectra of all DEPs underwent analysis for more information about the adduct forms obtained from the database comparison. Subsequently, we performed the putative annotation of the metabolites of interest by applying the filtered exact mass data to the molecular formula from KEGG.
Finally, the KEGG IDs of the matched compounds were submitted to the pathway analysis module for visualization through integrating enrichment and pathway topology analysis [64]. Parameter sets were as follows: the hypergeometric test and the latest KEGG version of the O. sativa pathway library.

Functional Annotation and Itegratomics Analysis
The results obtained using OmicsBox and PatternLab V underwent a functional classification. Distinct multiFASTA files generated were submitted to the functional classification in the BlastKOALA platform [64].
The approach used to integrate the three omics was pathway mapping, and the analysis was performed using the Omics Fusion platform [65]. Previously to the integration of multi-omics data, the NCBI accession of transcripts related to enzymes was converted to UniProt ID. Thus, the input data used were the IDs of each omics, which include UniProt Accession for transcriptomics and proteomics, and KEGG ID for metabolomics. Firstly, the data were enriched through several databases (EMBL, KEGG, NCBI, and UniProt), and then the module "KEGG feature distribution" was used to map these omics data in known pathways.

Conclusions
Previously, in addition to showing that young oil palm submitted to a high concentration of NaCl reduces the rates of CO 2 assimilation, stomatal conductance to water vapor, and transpiration, our group also confirmed a preponderant role of transcription factors in the early response of oil palm plants to salinity stress [6,15]. Data from ionomics, phenomics, and transcriptomics (mRNA and miRNA) were employed to show that. Currently, two new omics platforms joined this list-metabolomics and proteomics-and a first MOI study was performed. For phenomics-morphophysiological characterization, data came from two salinity stress experiments carried out in November 2017 and March 2018, and all the transcriptome, metabolome, and proteome data came from one of the experiments at once, 12 days after the onset of the stress, using a split-sample study design. Extensive leaf necrosis was already visible when the samples from the stressed treatment (electrical conductivity of~40 dS m −1 ) were collected, and one must consider that when analyzing these omics data sets.
The SOA and MOI studies here reported generated new insights on the response the early response of oil palm plants to salinity stress, pointing out genes, proteins, metabolites, and pathways directly affected by this stress. The eleven pathways identified by MOI analysis definitely appear at the top of the list as priorities for further studies. However, it is clear that two factors limited the accomplishments of the MOI study-the small number of differentially expressed metabolites identified via an untargeted metabolomics approach and the lack of data regarding the Log 2 (FC) from the proteins found exclusively in the control and stressed treatments when using the global proteomics analysis. No Log 2 (FC) from most of the DE proteins was identified, and only 19 DE metabolites limited the use of correlation studies.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/plants11131755/s1, Table S1: Proteins from the oil palm genome differentially expressed (FDR ≤ 0.05) in the leaves of young oil palm plants under salinity stress. FC-Fold Change; CPM-Count Per Million; FDR-False Discovery Rate; Kegg Orthology number; E.C. number; Table S2: a-List of all peptides confidently identified in oil palm conditions; b-list of all proteins confidently identified in oil palm conditions; c-all proteins confidently identified in oil palm conditions with maximum parsimony criterion; Table S3: a-list of the proteins detected only in control biological condition; b-list of the proteins detected only in stressed biological condition; c-list of all differentially expressed proteins detected in both biological conditions (Stressed and Control) with statistical significance (FDR ≤ 0.05); Table S4: a-list of 510 differentially expressed enzymes (436 from transcriptomics and 74 from proteomics) prospected in the transcriptome and proteome of young oil palm plants under salinity stress and submitted to integrative analysis; b-The top two pathways with the largest number of components and its description, obtained from transcriptomics, proteomics, and metabolomics integrative analysis using the Omics Fusion platform. Table S5: List of the proteins detected only in control biological condition; Table S6: List of the proteins detected only in stressed biological condition; Table S7: List of all proteins detected in both biological conditions (Stressed and Control) with statistical significance (FDR ≤ 0.05); Table S8: List of 510 differentially expressed enzymes (436 from transcriptomics and 74 from proteomics) prospected in the transcriptome and proteome of young oil palm plants under salinity stress and submitted to integrative analysis; Table S9: The top two pathways with the largest number of components and its description, obtained from transcriptomics, proteomics, and metabolomics integrative analysis using the Omics Fusion platform.