Similarities in the Behavior of Transcription Factors in Young Oil Palm (Elaeis guineensis Jacq.) Plants under Salt or Water Stress Provide Strategies to Develop Stress-Tolerant Oil Palm Plants at Once

1. Introduction

Oil palm (Elaeis guineensis Jacq.) is the number one oilseed crop in the world, considering the amount of vegetal oil consumed yearly. The consumption of palm oil and palm kernel oil in the 2022/23 season was about 85 million metric tons, ahead of soybean oil (60 million metric tons), rapeseed oil (31), and sunflower-seed oil (18.5) [1].

Over 95% of the oil palm plantations in Brazil are in the Legal Amazon Area, specifically in the Northeast region of the Pará State. Although Brazil has millions of hectares with an aptitude for oil palm production, it currently has only about 200,000 hectares destined for this purpose [2]. Any initiative aiming to increase the area cultivated with oil palm in Brazil faces a dilemma, stay in the Amazon Biome, where this crop finds a natural environment of high rainfall throughout the year, or go outside this area. Both options have its pro and cons.

One of the main cons of going somewhere outside the Amazon Biome, more specifically to the Mid-West and Northeast regions in Brazil, is that there oil palm will face long periods of drought throughout the year when it will not receive the physiological water level required to maintain productivity and will demand the use of artificial irrigation. Accordingly to Corley et al. [3], oil palm requires ~2,000 mm/year of water and does not tolerate water deprivation for more than 90 days.

The climate in the Cerrado Biome of the Mid-West region of Brazil has two seasons, rainy summers and dry winters. Some areas in that region experience a period of about four months without rain. The Northeast region is a semi-arid region where precipitation is insufficient to meet the water demand of crops. Water scarcity in places like those leaches mineral salts from the root zone, promoting soil salinization, a global problem due to its negative impact on agricultural productivity and sustainability [4]. About 30% of the world's irrigated agricultural land suffers from secondary salinity. That is a problem experiencing constant growth and expecting to reach more than half of the land by 2050 [5,6,7].

Plant response to drought and salinity stresses is a complex and coordinated process that rapidly and effectively works to restore homeostasis. It consequently enables the plant to adapt to unfavorable conditions [8]. Such a process regulates the expression of several stress-responsive gene groups [9], especially at transcriptional and post-transcriptional levels, with attention given to transcription factors (TFs). There are plenty of studies available showing that transcription factors (TFs) play central roles in the plant's response to different abiotic stimuli, allowing them to respond and adapt promptly to certain conditions, such as heat, drought, salinity, cold, and others [10]. Approximately 7–10% of all plant genes code for TFs at distinct moments, and several TF gene families are known for precisely coordinating the spatial and temporal expression of downstream genes associated with abiotic stress [11,12]. So far, 60 distinct TF families got identified in plants, including heat shock factors (HSPs), which are a class of great importance that play crucial roles in protecting plants against the effects caused by various abiotic stresses through the regulation of heat shock proteins [9,13].

Studies show that overexpression of certain stress-inducible transcription factors, such as NF-Y, MYB, HD-ZIP, and DREB, can increase tolerance to salinity, drought, or extreme temperatures and that they play central roles in plant response to different abiotic stimuli [14,15]. The homeodomain leucine zipper (HD-Zip) family of transcription factors is one of the largest plant-specific superfamilies. The TFs from this family got linked to plant growth modulation and responses to environmental stresses [16,17,18,19], and studies show that their overexpression does promote tolerance to drought [18,20,21], and to drought and salt stresses [22].

In two previous studies, our research group reported comprehensive, large-scale miRNA and RNA-Seq analyses applied to characterize young oil palm responses to a high level of salinity [23] and water deprivation [24]. In those studies, we showed that TFs play a preponderant role in the early response of this oilseed crop to those two abiotic stresses. Here we describe an additional step toward analyzing those datasets, aiming at getting more insights into the TF families and their TF genes commonly affected by those two stresses. To do so, all TF genes differentially expressed got considered, not only those identified as putative miRNA-target genes. A group of 20 TF genes (and their promoter sequences) from four distinct TF families underwent structural and functional annotation, with subsequent qPCR analysis of their behavior in both stress conditions.

2. Results

2.1. Selection of transcription factors responsive to drought and salinity stress

The results from the GhostKOALA annotation tool, specifically the ones from the KO0300 – Transcription Factor group, revealed 19 and 36 groups of different families of transcription factors responsive to salinity and drought stresses, respectively. There were 18 families of differentially expressed TFs present in both sets of results (Supplementary Table 1). Based on information gathered from previous studies [23,24,25,26], 20 differentially expressed genes got selected, and they belong to five of those 18 TFs families. Among them, there were 13 genes from the MYB family, three genes from the homeobox-leucine zipper protein (HD-ZIP) family, two genes from the heat shock transcription factor (HSFF) family, one gene from the nuclear transcription family Y subunit beta (NFY-B), and one of the nuclear transcription factor Y gamma (NFY-C) family (Table 1).

The RNASeq analysis of those 20 genes selected showed that eleven expressed differentially under saline stress, with fold change values ranging from 1.8 to 10.4 times the expression level seen in the non-stressed plants. Among the nine genes negatively regulated under salinity, gene 17 (XP_010936941.1) had a 17.2 times reduction in the expression level found in the non-stressed plants. Under drought stress, nine out of the 20 TFs genes expressed positively; gene 1 (XP_010904805.1) led the rank with a fold change of 50.2. The remaining eleven genes negatively regulated showed reductions ranging from -1.9 to -21.2, compared to the control (Table 1).

2.2. Selection of transcription factors responsive to drought and salinity stress

The phylogenetic tree obtained using the stretch of 1,432 nucleotides upstream of the start codon of each gene separated the 20 TFs genes into four groups; however, when comparing the expression profiles of the genes under the two stress tested, there was no particular pattern at all (Figure 1A, Table 1). The MEME software revealed three conserved promoter motifs in all 20 genes, and, again, no particular pattern appeared within the four groups found in the phylogenetic tree (Figure 1B).

The use of the PlantCare database revealed seven cis-acting elements related to stress response present in the promoter sequences of those 20 genes, which were: ABRE (ABA-responsive element), STRE (activated by heat shock, osmotic stress, low pH and nutrient deprivation), CGTCA-motif (responsive to methyl jasmonate (MeJA)), TGACG-motif (involved in MeJA responsiveness), TATA-box, CAAT-box, and the salicylic acid and auxin-responsive element as-1 (Figure 1C).

The CAAT- and TATA- box elements were present in large amounts in all 20 promoter sequences, being directly related to the initiation of gene transcription and enhancing activity, improving the ability of genes to react to adversity to respond to stress abiotic process more efficiently and quickly [27]. Two genes - 3 and 5 - from the HSFF family presented the elements ABRE, as-1, CAAT-box, CGTCA – motif, TATA-box, TGACG – motif, while genes 2, 6, 7, 8, 9, and 10 showed all seven cis-acting elements (Figure 1C).

2.3. Structural and functional annotation of the coding region

Regarding the size of the genes and the numbers and disposition of introns and exons, there was considerable diversity in the coding sequence of the 20 TFs genes analyzed. The longest gene was gene 4, almost 7,500 nucleotides long. The number of introns varied from zero (gene 20) to nine (gene 12). When analyzed by family, one can observe that egMYB presented two introns in most genes, while egHSFF presented one, egNFY zero or one intron, and in egHD-ZIP ranging from two to nine introns (Figure 2).

A phylogenetic tree resulted from analyzing the protein sequences of the 20 genes to explore the evolutionary relationship between transcription factor families, besides characterizing the features of conserved motifs and domains (Figure 3). The 20 TFs did group into three clades, where one clade had only the protein from gene 9 - an MYB family individual, and the second one had eight proteins - from genes also belonging to the MYB family. At last, the third clade includes the proteins of the remaining four genes of the MYB family, three genes of the HD-ZIP family, two genes of the HSFF family, and the two genes of the NFY family; both being side by side with their respective groups.

As for the conserved domains, only the protein coded by gene six had no conserved domain. The remaining 19 proteins showed 23 distinct conserved domains. The PLN03091 and PLN03212 domains appeared in the proteins coded by genes 1, 2, 7-9, 11, 13, 14, 16-19. Those proteins belong to the MYB family of transcription factors; and had the REB1, Myb, and SANT domains randomly distributed in them (Figure 3b, Supplementary Table S2).

In the coding regions of genes 3 and 5 (both of the HSFF family), the domains HSF, HSF1, HSF_DNA-bind, and HSF_DNA-bind superfamily were found. While in the protein coded by genes 4, 10, 12, 15, and 20, the domains HOX, homeodomain, homeodomain superfamily, Homeobox, HHT1, HALZ, H4, Histone superfamily, START, HAP5, HD-ZIP_N, START_ArGLABRA2_like, CBFD_NFYB_HMF were also found and SRPBCC (Figure 3b, Supplementary Table S2).

The MEME program predicted seven conserved motifs in 18 of the 20 proteins analyzed (Figure 3c). The numbers of conserved motifs vary from one to five per protein, with the proteins coded by genes 11 and 19 having the highest number of motifs, and the length of the motifs ranged from eight to 50 amino acids. Furthermore, in domains PLN03091 and PLN03212, one can see motifs 1, 2, and 4. Motif 5 was present in the proteins coded by genes 3 and 5, while motif 6 was present in the proteins coded by genes 10, 12, and 15. Motif 3 was present in nine proteins, while motif 7 was in three (genes 10, 11, and 19) (Figure 3c).

From the functional annotation of the TFs coded by the 20 genes responsive to drought and salinity stress, 87 activities for biological processes, nine for molecular functions, and four for cell components were identified (Figure 4, Supplementary Table S3). One should notice that the protein coded by gene 14 had the highest number of GO terms annotated, with 29 GO terms, followed by those coded by genes 9, 11, and 19, with 23 GO terms for each one (Supplementary Table S3).

2.4. Differential expression and correlation analyses

To further characterize the 20 genes selected, their expression levels got determined using qPCR (Figure 5, Supplementary Table S4). In water deprivation, genes 1, 5, 9, and 11 were significantly positively regulated, while the remaining genes were negatively regulated. Under salt stress, genes 1, 5, 9, 11-16, and 19 were positively regulated, while genes 2-4, 6-8, 10, 17, 18, and 20 were negatively regulated.

The differential expression of the 20 genes under salt and drought stress underwent a correlation analysis. Interestingly, almost all genes showed the same behavior when analyzed under drought and salinity stress conditions, with an R2 score of 0.7317. Only genes 12 and 13 showed opposite behavior when subjected to different stress conditions, being positively regulated under salt stress and negatively under water deprivation (Figure 6).

The expression levels of all genes, measured by RNA-Seq and qPCR, were compared under both abiotic stresses (Figure 7). Among the 20 genes tested, only gene 4 showed a different expression pattern, showing a positive regulation response to both conditions when using the RNA-Seq and a negative regulation response when using qPCR.

3. Discussion

Soil salinity and water deficit are two of the most severe abiotic stress factors, restricting the growth and development of many plant species and generating significant losses in crop productivity [28]. Plants respond to unfavorable conditions through two pathways, namely the abscisic acid (ABA)-dependent pathway and the ABA-independent pathway [29]. Transcription factors (TFs) play vital roles in mediating stress tolerance in both of them, stimulating osmotic homeostasis, and ionic balance, among other elements, in an attempt to reduce the damage caused by stress in plants [30,31].

There is a lack of knowledge about the functioning of transcription factors (TFs) in regulating gene expression in oil palm plants under drought and salt stresses. Recently, our group reported several transcription factors been affected by those two stresses in connection with miRNAs [23,24]. The goal of the present study was to characterize TFs' found as responsive to both stresses in young oil palm plants, independently of been or not a miRNA-target gene. To do so, a group of 20 TF genes got selected for further analysis. Those TFs belong to four distinct families - MYB, HD-ZIP, HSFF, and NF-Y - previously reported as directly involved in salt and drought resistance in some plant species.

The first step in the characterization process focused on identifying features commonly found in the promoter sequences of those genes. Gene promoters contain essential regulatory cis-elements for regulating gene expression [32], and, in the present study, as expected, many important cis-elements were present in all 20 promoter regions. Numerous CAAT-box and TATA-box sequences appeared, presenting themselves possibly as the most relevant cis-elements of the ones identified. CAAT-box is a cis-acting element commonly found in promoter sequences with enhancing activity that may contribute to increased activity of promoters in response to salinity and drought stress [32]. Meanwhile, the genomic distribution of the TATA box tends to favor function-specific genes. Such as, TATA-box is rich in promoters of stress-related genes and is absent in promoters of essential genes, in addition to being associated with rapid and variable regulation under stress [27,33].

Besides CAAT – box and TATA – box, other cis-acting elements also reported as having links with response to abiotic stresses appeared, such as ABRE (response to ABA), STRE (activated by heat and osmotic stresses), as-1 (Oxidative Stress-Responsive Element activated by salicylic acid), and CGTCA – motif and TGACG – motif, both responsive to methyl jasmonate (MEJA). Although 18 out of the 20 genes presented similar qualitative profiles in the two scenarios - salt and drought stresses - evaluated, which got validated by RNA-Seq and qPCR analysis, not one specific pattern in terms of the number and position of those above cited cis-acting elements was evident.

The second step in the characterization process focused on identifying features in the sequence of the proteins coded by those genes. As expected, regarding the motifs and domains present in the proteins coded by the selected genes, it was possible to identify equivalent patterns among the representatives of each specific TFs family, corroborating the results from the GhostKOALA annotation tool, specifically the ones from the KO0300 – Transcription Factor group.

Thirteen of the selected genes belong to the MYB transcription factors family. Several studies have shown that MYB proteins perform multiple functions in response to stresses such as salt, drought, and cold [34,35,36]. MYB20 (gene 13) showed a positive regulation in the leaves of young oil palm plants under saline stress and a negative under drought. This gene is known for acting as a negative regulator of ABA-mediated stomatal closure under conditions of water stress, and its overexpression results in higher susceptibility to desiccation [37]. According to Wang et al. [36], it positively regulates plant salt tolerance by repressing key ABA repressor, PP2Cs. MYB30 (gene 16) had its expression increased in the leaves of young oil palm plants in response to both stresses. This FT is known for modulating plants' tolerance to salt through the regulation of mitochondrial alternative oxidase AOX1s; and the AOX1a overexpression rescued salt-sensitive phenotypes [36,38]. On the other hand, MYB30 positively regulates wax synthesis genes that have a direct link with the biosynthesis of the cuticle, one of the most important plant structures protecting plant leaves from desiccation [36,39].

Two of the selected genes belong to the NF-Y transcription factors family. Nuclear factor Y (NF-Y) proteins, also known as CCAAT Binding Factor (CBF) or Heme Activator Protein (HAP), modulate the expression of downstream genes via a heterotrimer – NF-YA-YB-YC – that binds through NF-YA to the CCAAT box present in the promoter region of those genes. Several members of the NF-Y family and their subunits are regulators of abiotic stress responses in plants [40,41,42,43,44]. NF-YA is a target of the miR169 family, and the overexpression of NF-YA in Arabidopsis increases its tolerance to salt stress [42].

Our previous studies have shown that members of the NF-Y family are regulated by microRNAs (miR169 and egu-miR02sds) in the leaves of young oil palm plants under salinity and drought stress [23,24]. Even though the miR169 gene got negatively regulated in oil palm plants under salinity and drought stress, the levels of NF-YA increased only under salt stress [23]. In the present study, NF-YB3 (gene 20) and NF-YC2 (gene 4) got differentially expressed in the leaves of oil palm plants under both abiotic stresses, with the former negatively regulated under salinity and drought stress and the latter positively regulated.

A hypothetical model presented by Leyva-González et al. [41] proposes that in plants exposed to abiotic stress, NF-YA levels increase due to their transcriptional activation and to the reduction in the miR169 levels, corroborating the results from Salgado et al. [23]. The higher amount of NF-YA represses early abiotic stress response genes probably by sequestering NF-YB-YC, creating a regulatory loop to arrest early responses that represent high energy and carbon costs. If that is the case, that might happen in the leaves of young oil palm plants only under salinity stress, considering that the reduction in the levels of miR169 due to drought stress does not show to affect the level of NF-YA, as seen in Salgado et al. [24].

Previous studies have reported that AtNF-YB1 improves plant performance in drought conditions, and its ortholog in maize, ZmNF-YB2, conferred resistance to water stress in the field [40]. That suggests that this would be a starting point for identifying possible markers for genetic improvement. Gene 4, an NF-YC2, was overexpressed in the leaves of oil palm plants under drought and salinity stress; this behavior was reported in conifers (Picea wilsonii) as well, where overexpression of PwHAP5 (Arabidopsis homolog NF-YC2) partially rescued the increased sensitivity of NF-YC2 to salt, drought, and ABA treatments, through its positive regulation [42].

Three of the selected genes belong to the HD-Zip transcription factors family, one positively (gene 15) and one negatively (gene 10) regulated in both scenarios, and one positively regulated under salinity stress and negatively under drought stress (gene 12). The homeodomain leucine zipper (HD-Zip) transcription factor family is one of the largest plant-specific superfamilies and includes genes with roles in modulating plant growth and response to environmental stresses, and categorized into four subfamilies - HD-Zip I, HD-Zip II, HD-Zip III, and HD-Zip IV [16,18,19].

The overexpression in Arabidopsis of HaHB4, a drought/ABA-inducible HD-Zip family I gene from sunflower, upregulated transcript levels of several genes involved in JA biosynthesis and defense-related processes, resulting in plants with shorter stem and internode distances, more round leaves and denser hairs; altogether, such phenotypes made the plants tolerant to drought [18,20]. Overexpression of ATHB-6, a drought-induced HD-Zip I gene from maize, has activated the expression of critical genes in the ABA-dependent and ROS signals pathways and improved tolerance to drought in maize [21]. Overexpression of Zmhdz10, also an HD-Zip I gene from maize, in rice and arabidopsis led to enhanced tolerance to drought and salt stresses [22]. At last, previous research using class I HD-ZIP homeobox genes showed that they get highly regulated under salt stress and that their overexpression in transgenic rice and cotton plants improved salt stress sensitivity [17,45]. The three HD-Zip genes characterized in this present study belong to the HD-Zip I (gene 10), HD-Zip II (gene 15), and HD-Zip IV (gene 12) subfamilies.

The remaining two of the selected genes belong to the HSF family. Heat shock transcription factors (HSF) are a TF class of great importance that play crucial roles in protecting plants against the effects caused by various abiotic stresses [46]. HSFs operate by inducing the accumulation of heat-shock proteins (Hsps) and by mediating the activation of other heat-responsive genes involved in cell protective mechanisms and the homeostasis of reactive oxygen species (ROS) [47]; and are reported as responsive to heat, drought, cold, and salt stresses [48]. There are three classes of HSFs, denominated A, B, and C. The HSFs from class-A have multiple acidic motifs (AHAs) at the C-terminus and function as transcriptional activators, and those from class-B contain a tetrapeptide -LFGV- in the C-terminus, which is assumed to function as a repressor motif by interaction with an unknown co-repressor, while the function of those from class-C is still unknown [47,49].

In Arabidopsis, Pérez-Salamó et al. [50] showed that estradiol-dependent induction of HSFA4A confers enhanced tolerance to salt and oxidative agents, whereas inactivation of that gene results in hypersensitivity to salt stress. The overexpression of LlHsfA3A or LlHsfA3B from lily (Lilium longiflorum) in Arabidopsis showed hypersensitivity to salt stress, and a lack of sucrose exacerbated this salt sensitivity [49]. In the present study, gene 3 - coding for an HSFA3 protein - is negatively regulated in both scenarios, while gene 5 - coding for an HSFC2b protein - was positively regulated.

Regarding the expression profile of the selected genes measured by RNA-Seq and qPCR techniques, 18 out of 20 showed similar (qualitatively) expression profiles in the two scenarios; only genes 12 and 13 differed. Such similarity in behavior in the initial responses of plants to salt stress and drought got reported previously [48,51,52,53], corroborating our results.

Last, and based on the study by McCoy et al. [54], we applied GENIE3 to mine 306 public oil palm transcriptome datasets and to construct a gene regulatory network (GRN) featuring epigenetic regulators and transcription factors from the oil palm genome. Such a study [55] used 1,333 unique regulators and 27,642 target genes from the oil palm reference genome [56]. A closer look at this oil palm GRN, and focusing on the 20 TFs genes characterized in the present study, allowed the identification of the genes regulating the TFs and those under regulation by them. Gene 10 (LOC105048648), a candidate gene concerning its potential capability to promote resistance to both salinity and drought stresses by heterologous overexpression, is a target for 16 epigenetic regulators and transcription factors while targeting 278 genes (Supplementary Table S5). A gene expressing a protein belonging to the Transcription activator GLK1/2-like (IPR044825) family is the one with the highest score among those potentially regulating gene 10, followed by another expressing a protein from the ADP-ribosylation factor GTPase-activating protein AGD11/12/13 (IPR044518) family.

4. Materials & Methods

4.1. RNA-Seq data generation

The RNA-Seq raw sequence data used in this study are available in the Sequence Read Archive (SRA) database of the National Center for Biotechnology Information under Elaeis guineensis Transcriptome_Drought and Salinity Stresses - BioProject PRJNA573093, BioSample SAMN12799239.

All information regarding the plant material, growth conditions, and abiotic (salt and drought) stress conditions used in the studies that generated RNA-Seq data are in Salgado et al. [23,24]. In summary, all plants used in this study were clones regenerated from embryogenic callus obtained from the leaves of an adult plant belonging to the AM33 genotype, and all experiments were carried out in a greenhouse at Embrapa Agroenergia (www.embrapa.br/en/agroenergia), in Brasília, DF, Brazil (S-15.732 °, W-47.900 °), using young oil palm plants in the growth phase known as bifid saplings.

In the salinity stress study, on the 12th day after the imposition of the stress, we collected the apical leaf from three control (0.0 g of NaCl per 100 g of the substrate) and three salt-stressed plants (2.0 g of NaCl per 100 g of the substrate) [23]. In the water deprivation study, apical leaves from control and drought-stressed plants were collected 14 days after the imposition of treatments, with three replicates per treatment [24]. In both studies, the plant tissue was immediately immersed in liquid nitrogen and then stored at -80 °C until RNA extraction, library preparation, and sequencing. Total RNA isolation, and RNA quality and quantity analysis, were performed as described by Salgado et al. [23,24]. The GenOne Company (Rio de Janeiro, RJ, Brazil) performed the RNA-Seq using an Illumina HiSeq platform and the paired-end strategy.

4.2. RNA-Seq data analysis

Details on the RNA-Seq analyses performed using OmicsBox version 1.3 [56] can be found in Salgado et al. [23,24]. The high-quality RNA-Seq data got aligned to the oil palm reference genome [57] - files downloaded from NCBI (BioProject PRJNA192219; BioSample: SAMN02981535) in March and October 2020 - using default parameters from OmicsBox version 1.3 through software STAR [58].

To quantify expression at the gene or transcript level, we used the default parameters of OmicsBox version 1.3 through HTSeq version 0.9.0 [59]. The pairwise differential expression analysis between the experimental conditions used edgeR version 3.28.0 [60], applying a simple design and an exact statistical test; without using a filter for low-count genes [23,24].

4.3. Genes selection and annotation

The differentially expressed genes identified in Salgado et al. [23,24] underwent annotation using the GhostKOALA annotation tool. Based on the results of the KO0300 – Transcription Factor group, we carried out a manual curation using the following conditions: 1°) Selection of the most relevant families according to the literature; 2°) Selection of differentially expressed genes within these previously selected groups; 3°) Genes that are present in both stresses; and 4°) The genes that showed a higher level of differential expression, both Up and Down.

All information regarding intron and exon regions for each gene came from the NCBI (https://www.ncbi.nlm.nih.gov/). To identify the cis-elements in the promoter region, we got a stretch of 1,432 base pairs before the ATG start codon of each gene from the oil palm reference genome [56], and the known elements characterized using the database PlantCARE [61]; the statistical sieve was Score ≤ 5. To analyze the structure of proteins and promoters related to transcription factors, we applied the MEME software (https://meme-suite.org/meme/tools /meme), with a maximum of 10 found motifs. Motifs that did not pass the statistical sieve (E-value ≤ 0.05) got excluded from the analysis.

After obtaining the protein sequences of the selected genes from the oil palm reference genome [56], the proteins domains got characterized through the NCBI Conserved Domain Search (www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi). Those domains that did not pass the statistical sieve (E-value ≤ 1) got excluded from the analysis. For the functional analysis of proteins, we used the EggNOG platform [62]. For the multiple sequence alignment, we applied the ClustalW program [63]; for both the protein and promoter region data. Subsequently, the phylogenetic trees got built using the maximum likelihood method, using a bootstrap of 1000 replicates, using the MEGA11 software [64].

The figures referring to the results of the cis-element analysis; motifs, protein domains, and intron and exon regions got generated through the Python library “DNA Features Viewer” (https://edinburgh-genome-foundry.github.io/DnaFeaturesViewer/).

4.4. RNA extraction, reverse transcriptase-PCR and quantitative real-time PCR analyses

Total RNA isolation, and RNA quality and quantity analysis, were performed as described by Salgado et al. [23,24]. After quantified, the extracted total RNA served as a template for reverse transcription using the commercial kit SuperScriptR (Invitrogen) following the manufacturer's protocol. The positive control gene was the EgEfMPOB00119 40S mRNA of ribosomal protein S23, a constitutive gene in E. guineensis.

For primer design, fasta files from the selected genes got downloaded from NCBI and exported to PerlPrimer software to design primers (Table 2). qPCRs were performed in 96-well optical plates on a 7500 real-time PCR system (Applied Biosystems, Foster City, CA, USA) with a PowerUpTM SYBRTM Green Master Mix, REF A25741 (Thermo Fisher Scientific), following the manufacturer's instructions. The qPCR reaction conditions were as follows: 50°C for 2 min, 95°C for 2 min, followed by 40 cycles of 95°C for 5 seconds and 60°C for 30 seconds. A melting curve got created to analyze the specificity of the reactions. Then, three biological replicates got performed for each sample. The expression level of the untreated sample got set to 1 for normalization.

4.5. Pair-by-pair comparison and correlation analysis

We performed a pair-by-pair comparison between the different stresses and methods used to analyze and validate the genes; the input data were Log₂ (Fold Change). The correlation analysis got carried out using the data set referring to the differential expression of genes (DEGs) common to the two scenarios studied (salinity and drought); obtained by qPCR.

5. Conclusions

To the best of our knowledge, no report shows that oil palm (E. guineensis) has salinity or drought stress tolerance. All transcriptome data collected and analyzed in the present study and in Salgado et al. [23,24] derived from the apical leaf of young oil palm plants submitted to very high salinity levels in the soil or water deprivation. When we collected those leaves, the plants showed severe symptoms (browning, scorched, dying) as a consequence of being for almost two weeks under stress. Consequently, we assume that those plants were not showing any signs of tolerance; on the contrary, they were affected by high salinity in the soil or water deprivation and would die if kept longer under those conditions.

Starting from this premise and remembering that the initial goal of this study was to get more insights on the role of transcription factors in the response of young oil palm plants to those two stresses, it is necessary to state that the development of plants tolerant/resistant to those abiotic stresses was also in the horizon during this study. That is why a search for molecular symptoms that might allow us to set up strategies for genetic engineering/editing got performed.

In conclusion, the results presented here shed light on the role of 20 TFs from four different TF families during the response of young oil palm plants to severe levels of salinity in the soil or water deprivation. The behavior of such a group of TFs made it very clear that there is a high level of similarity in the molecular response to both stresses. That similarity allows the search for genes that might grant strategies aiming to develop oil palm plants with a tolerance to both salinity and drought at once. As an example, and based on the results obtained by our group - as well as on the scientific literature available regarding some of the genes here characterized and the success of some of them in achieving tolerance to those abiotic stresses through overexpression or silencing - we postulate that the overexpression of gene 10, an HD-Zip I gene from oil palm, is a valuable strategy to obtain plants with increased tolerance to drought and salinity. However, one must state that once this gene seems to be regulated by 16 and regulating 278 genes, such behavior might trigger several pleiotropic effects once overexpressed in the oil palm genome or heterologously.