Expression analysis of miRNAs and their putative target genes confirm a preponderant role of transcription factors in the early response of oil palm plants to salinity stress

Background Several mechanisms regulating gene expression contribute to restore and reestablish cellular homeostasis so that plants can adapt and survive in adverse situations. MicroRNAs (miRNAs) play roles important in the transcriptional and post-transcriptional regulation of gene expression, emerging as a regulatory molecule key in the responses to plant stress, such as cold, heat, drought, and salt. This work is a comprehensive and large-scale miRNA analysis performed to characterize the miRNA population present in oil palm (Elaeis guineensis Jacq.) exposed to a high level of salt stress, to identify miRNA-putative target genes in the oil palm genome, and to perform an in silico comparison of the expression profile of the miRNAs and their putative target genes. Results A group of 79 miRNAs was found in oil palm, been 52 known miRNAs and 27 new ones. The known miRNAs found belonged to 28 families. Those miRNAs led to 229 distinct miRNA-putative target genes identified in the genome of oil palm. miRNAs and putative target genes differentially expressed under salinity stress were then selected for functional annotation analysis. The regulation of transcription, DNA-templated, and the oxidation-reduction process were the biological processes with the highest number of hits to the putative target genes, while protein binding and DNA binding were the molecular functions with the highest number of hits. Finally, the nucleus was the cellular component with the highest number of hits. The functional annotation of the putative target genes differentially expressed under salinity stress showed several ones coding for transcription factors which have already proven able to result in tolerance to salinity stress by overexpression or knockout in other plant species. Conclusions Our findings provide new insights into the early response of young oil palm plants to salinity stress and confirm an expected preponderant role of transcription factors - such as NF-YA3, HOX32, and GRF1 - in this response. Besides, it points out potential salt-responsive miRNAs and miRNA-putative target genes that one can utilize to develop oil palm plants tolerant to salinity stress. Supplementary Information The online version contains supplementary material available at 10.1186/s12870-021-03296-9.

important in the transcriptional and post-transcriptional regulation of gene expression, emerging as a regulatory molecule key in the responses to plant stress, and the main components of miRNA show high conservation between species [2,3].
The synthesis of plant miRNAs happens in the nucleus, where the encoded plant miRNA genes are processed by polymerase II to form miRNAs [4]. Long sequences of miRNAs are folded into a hairpin structure, known as primary miRNAs or pri-miRNAs, and are then cleaved by DCL1 to form short and incomplete double-stranded structures called pre-miRNAs [5]. Pre-miRNAs are cleaved by DCL1 or DCL4 to form a double-stranded miRNA known as miRNA dimer. Subsequently, the methyltransferase HEN1 carries out methylation of the 3′ end, and the miRNA is then transported to the cytoplasm by the plant homologous protein exportin-5 (HASTY, HST). Then the miRNA single-chain and the AGO protein form the RISC complex, binding to the complementary target mRNA to cleave or inhibit translation, obtaining negative regulation of the target gene [6,7].
Several abiotic stressors, such as cold, heat, drought, and salt, affect the plant life cycle interfering with growth and productivity [8]. Several mechanisms regulating gene expression contribute to restore and reestablish cellular homeostasis so that plants can adapt and survive in adverse situations. miRNAs play a role important in regulating gene expression in response to stress conditions [1,9].
Abiotic stresses upregulate some genes and downregulate others, depending on the role played by the gene. miRNAs responsive to water stress are present in Oryza sativa [10,11], Arabidopsis thaliana [12], and Medicago truncatula [13,14]. Sunkar and Zhu [15] showed that miR319c is positively regulated in Arabidopsis when subjected to cold stress but did not change when subjected to dehydration, salt, or ABA. Several miRNAs, such as miR156, miR159, miR167, miR171, miR319, and miR396, showed differential expression levels during the response to salt stress in Arabidopsis sp. [12] and Zea mays [16]. Using state-of-the-art sequencing technology (NGS) to identify miRNAs responsive to salt, Dong et al. [17] identified 104 differentially expressed miRNAs in soybean nodules under salt stress.
Some studies have reported miRNAs in oil palm. Nasaruddin et al. [18] found five new potential miRNA encoding sequences in a collection of 7284 oil palm EST sequences by a combined homology and structural analysis approach, having roles in regulating the auxin response, floral development, and basal transcription. Low et al. [19] applied a homology approach to identify 14 miRNAs in contigs assembled from sequences generated from the hypomethylated or gene-rich regions of Elaeis guineensis and E. oleifera genomes. Silva et al. [20] identified 57 mature miRNA in E. guineensis and 52 in E oleifera, respectively, and miRNA-target prediction revealed that most of these miRNA-putative target genes are transcription factors involved in the plant development process, particularly the regulation of root development. Ho et al. [21] investigated microRNA expression in oil palm female inflorescences at two stages of floral development corresponding to the emergence of floral meristems and the formation of floral organs, identifying 15 oil palm-specific miRNA candidates. Zheng et al. [22] identified 52 miRNAs in a study aiming to gain insights into the regulatory mechanisms of lipid and fatty acid (FA) metabolism in oil palm.
Oil palm (Elaeis guineensis Jacq.) is a source of vegetable oil that has great importance in many economic sectors. In Brazil, over 95% of the oil palm plantations are in the Amazon region. Due to environmental restrictions imposed on the use of the Amazon rainforest, and the logistical difficulties to flow the production to the main industrial centers in the country, there is a crescent demand of growers for the cultivation of oil palm in other geographic regions in the country. One must use irrigation in oil palm plantations outside the Amazon region in Brazil, mainly due to long periods of drought observed in these alternative regions with potential for oil palm cultivation [23,24]. Between 25 and 30% of the irrigated land area in the World is affected by salt and is essentially commercially unproductive [25]. Because of that, Embrapa started working on a multi-omics approach to characterize morphophysiological and molecular responses of oil palm (E. guineensis) to drought and salinity stresses [26].
The current study is a follow-up to the characterization of the morphophysiological responses of oil palm plants to drought and salinity stress [24,26]. We carried out a comprehensive, large-scale miRNA analysis to characterize the miRNA population present in oil palm exposed to a high level of salt stress, to identify miRNA-putative target genes in the oil palm genome, and to perform an in silico comparison of the expression profile of the miR-NAs and their putative target genes.

Results
As shown previously in Vieira et al. [24], the electrical conductivity (EC) of the saturation extract increased, and the water potential (Ψw) decreased in a NaCl dosedependent manner. At the 12th day after imposing the stress (DAT), the EC values ranged from ±2 dS m − 1 (control plants) to ±45 dS m − 1 (stressed plants in substrate treated with 2.0 g of NaCl per 100 g of the substrate), while the Ψw values varied from zero to − 1.42 MPa, respectively (data not shown). There were no differences in the average evapotranspiration between the groups on day zero; however, when subjected to stress, the plants started to show differences in the evapotranspiration rates, remaining until the end of the experiment (data not shown).
At 12 DAT, it is visible a reduction in the rates of CO 2 assimilation (A), stomatal conductance to water vapor (gs), and transpiration (E), which correlated with the amount of NaCl used (Fig. 1A, B, and D). On the other hand, the increase in intercellular CO 2 concentration (Ci) also correlated with the amount of salt used (Fig. 1C). Stressed plants at the highest NaCl dose were already showing senescence of the leaves at 12 DAT (Fig. 2). Based on the morphophysiological responses of young oil palm plants to salinity stress (Figs. 1 and 2, and Vieira et al. [24]), both the control and the 2.0 g of NaCl per 100 g of the substrate treatments -which will be from now on referred as control and stressed treatmentswere selected to the characterization of the microRNA and mRNA profiles.

Identification of known and novel miRNAs, and differential expression analysis of miRNAs
The small RNA raw sequence data (9 fastq files) used in this study have been uploaded in the Sequence Read Archive (SRA) database of the National Center for Biotechnology Information under Elaeis guineensis microRNA_Drought and Salinity Stresses -BioProject Fig. 1 Box plots of the changes in gas exchange parameters in oil palm plants grown under increasing concentrations of NaCl in the cultivation substrate. A Net CO 2 assimilation rate (A); (B) stomatal conductance rate to water vapor (gs); (C) intercellular CO 2 concentration (Ci); and (D) transpiration rate (E). Medians and interquartile ranges are shown. The values represent the average of four replicates. The significance of differences was calculated using the Kruskal-Wallis test with post-hoc Dunn (p < 0.05). nsP > 0.05, non-significant comparisons were not shown in the graph; *P ≤ 0.05; **P ≤ 0.01; ***P ≤ 0.001 number of PRJNA646488, BioSample SAMN12799239. All adapter-free small RNA sequences (from all nine fastq files) were concatenated into a single file and submitted for miRNA prediction using mireap version 0.2 and Shortstack version 3.4, generating 96 positive hits for potential miRNAs (data not shown). Concomitantly, all adapter-free small RNA sequences from control and salt-stressed samples (three replicates) were submitted to assemble and then mapped against the oil palm reference genome [27], generating 3384 positive hits (data not shown).
A search in the database with the 3384 hits to the oil palm genome, using the database with the potential miRNAs, led to a total of 79 miRNAs, being 52 known miRNAs and 27 new ones (Fig. 3). The length of the 27 new ones ranges from 21 (24 miRNAs) to 22 (3 miR-NAs) nucleotides (Supplementary Table 1). The genes of the 79 miRNAs identified in this study ranged from 68 to 285 bp in length and spread throughout all 16 chromosomes of the E. guineensis genome (Supplementary Table 2). Several miRNAs are present in more than one place in the genome, in different chromosomes, or at different positions in the same chromosome. Chromosomes 01, 04, and 08 had the highest miRNAs amount, 12, 11, and 11, respectively. Twenty-eight miRNAs got mapped to 28 unplaced scaffolds. The highest number of miRNAs in one unplaced scaffold was three, in scaffold NW_011551039.1 (Supplementary Table 2).
The new putative miRNA genes are between 68 and 267 bp in length (Supplementary Table 2). Regarding the location of these genes in the genome of E. guineensis, 20 of them were in intragenic and nine in intergenic regions (Table 1). Of the genes present inside genes, only the egu-miR16sds, egu-miR17sds, egu-miR22sds, and egu-miR26sds genes are in non-characterized ones ( Table 1).
Among the 79 known and novel miRNAs, 72 showed a significant (probability ≥0.95) different level of expression under saline stress; however, all were downregulated (Supplementary Table 2). These differentially expressed (DE) miRNAs had their expression level reduced in the range from 39.75 to 99.82%. In general, those DE miR-NAs with their genes located in different regions in the genome did not present very distinct Log 2 (FC) values; the only exception was ppe-miR397 (Supplementary Table 2).

Prediction and differential expression analysis of miRNA-putative target genes
The psRNA-Target online program, version 2, led to 425 positive hits as miRNA-putative target genes. When analyzing the mode of inhibition of these positive hits, the vast majority, 398, presented an mRNA cleavage mode and the remaining 27 a translation inhibition (Supplementary Table 3). It usually occurs because of some incompatibility around the center of the complementary region, as the central area is essential for cleavage [27].
Out of the 425 positive hits, there were 229 distinct putative target genes; based on the LOC Ids from the oil palm reference genome. Among these putative target genes, 150 were target to just one miRNA, and 79 were target to more than one miRNA -ranging from two to 36 miRNA per target gene (Supplementary Table 3).
The RNA-seq fastq files used in this study -from control and stressed plant samples, three replicates/ treatment -are part of a group of 18 fastq files that have been uploaded in the Sequence Read Archive (SRA) database of the National Center for Biotechnology Information under Elaeis guineensis Tran-scriptome_Drought and Salinity Stresses -BioProject PRJNA573093, BioSample SAMN12799239. The transcriptome analysis showed that over 90% of the raw read pairs survived the preprocessing stage requiring a minimum average quality of reads ≥30 and the minimum length of reads ≥75 nucleotides. Over 95% of the high-quality read pairs mapped to the reference genome available at NCBI [28]. The reference genome has 29,567 genomic features of type 'gene' retrieved from 2781 ref. sequences in GCF_000442705.1_EG5_ genomic.fna file; however, 4213 of these features had no aligned reads detected in any of the samples ( Table 2).
When comparing control against stressed plants, the pairwise differential expression analysis revealed that out of the 29,567 features from the E. guineensis genome [28], 5366 were DE at False Discovery Rate (FDR) < 0.05 (data not shown); being 2380 upregulated (Log 2 (FC) > 0) and 2986 downregulated (Log 2 (FC) < 0). By applying the same criteria for the differential expression analysis of the 229 distinct miRNA-putative target genes previously prospected (Supplementary Table 3), a group of 24 upregulated and 27 downregulated genes were identified (data not shown). These 51 DE putative target genes were integratively and functionally annotated (Table 3, Supplementary Table 4).

Integrating the expression profiles from DE miRNA-putative target genes and their respective DE miRNAs
By integrating the expression profiles of the upregulated miRNA-putative target genes and their respective miR-NAs, twenty-one of them showed just one DE miRNA, one showed two (LOC105047586), and one showed four (LOC105046708) linked to the miRNA-target gene, while one had a non-differentially expressed miRNA linked to it (LOC105059776) ( Table 3). On the other hand, by integrating the expression profiles of the downregulated miRNA-putative target genes and their respective miRNAs, twenty-five of them showed just one, and two showed two (LOC105054987, LOC105059511) DE miRNA linked to the miRNA-target gene, while one had a non-differentially expressed miRNA linked to it (LOC105058639) ( Table 3).
Among the 40 DE miRNAs with DE miRNA-putative target genes, 28 had only one target gene, nine had three, and three had three (Table 3). DE miRNAs vvi-miR171j, gu-miR18sds, and egu-miR03sds had three distinct DE putative target genes presenting different profiles when submitted to salinity stress. Vvi-miR171j downregulated to 15% of its expression level in the control plants, while its putative target genes downregulated to 38 (LOC105060969), 41 (LOC105059511), and 50% (LOC105054987). DE miRNA egu-miR18sds also downregulated to 15% of its expression level in the control plants, while one of its putative target genes upregulated to 189% (LOC105055689) and two downregulated to 39 (LOC105031985) and 67% (LOC105061136). On the other hand, while egu-miR03sds downregulated to 21% of its expression level in the control plants, one of the putative target genes downregulated to 74% (LOC105048718) and the other two upregulated to 137 (LOC105040914) and 143% (LOC105059001) of their initial expression level (Table 3). LOC109505530 was the miRNA-target gene that experienced the highest expression level increase in the leaf of oil palm plants due to saline stress. This gene is one of two found as targetted by aof-miR536, but the only one to upregulate due to this stress. The saline stress led aorf-miR536 to downregulate to less than 10% the level found in the control plants, while LOC109505530 upregulated to almost eight times its initial expression level (Fig. 4A). LOC105046708 is target of three distinct miRNAs (ata-miR166d-3p, sly-miR166c-3p, and osa-miR166i-3p) in the genome of oil palm. This gene experienced an increase of approximately 80% due to saline stress, while the miRNAs targetting it downregulated to between 13 and 18% the level found in the control plants (Fig. 4A).
All oil palm new miRNAs identified in this study downregulate due to saline stress. However, their differentially expressed putative target genes belonged to two groups according to their response to saline stress, seven upregulated and 13 downregulated (Fig. 4B). In the case of egu-miR13sds and egu-miR18sds, the fate of their respective putative target genes -each one had 3 -were completely distinct, with some upregulating and some downregulating due to saline stress (Fig. 4B).

Functional annotation of the differencially expressed putative target genes
Among the 51 DE miRNA-putative target genes selected for functional annotation analysis, twenty had positive hits for biological process, 33 for molecular function, and eight for cellular component (Supplementary Table 4). The regulation of transcription, DNA-templated (GO:0006355), was the biological process with the highest number of hits, six, followed by the oxidation-reduction process (GO:0055114) with five. Protein binding (GO:0005515) and DNA binding (GO:0003677) were the molecular functions with the highest number of hits, six, five, respectively. The cellular component with the highest number of hits was the nucleus (GO:0005634), four, followed by the membrane (GO:0016021). Four genes had hits for domains from the GRAS family, a player important in gibberellin signaling [29].
The functional annotation analysis led to 18 known proteins, besides three lncRNA genes and three uncharacterized/unknown proteins. Among the 18 proteins there are six different kinds of transcription factors (LOC105054175, LOC105056468, LOC105046708, LOC105039459, LOC105048659, and LOC105043768). The remaining protein are: SPX-MFS proteins (LOC105043377,

Discussion
Among the few studies reporting miRNA in oil palm [18][19][20][21][22], none of them has studied the role of this type of macromolecule in the response of this species to salinity stress. So, to the best of our knowledge, this is the first report on the expression profile of oil palm miRNAs and their putative target genes when subjected to saline stress. The 52 orthologous miRNAs identified in this study belong to 28 different families previously reported in oil palm as expressed in floral meristems -miR156, miR160, miR166, miR167, miR168, miR172, miR396, miR528, and miR535, by Ho et al. [21]; in the development of the mesocarp -miR156, miR395, and miR528, by Fang et al. [30]; and in shoot apical meristem, immature and mature flowers -miR156, miR159, and miR160, by Nasaruddin et al. [18] using ESTs from a study by Ho et al. [31].
The expression profile of the 79 miRNAs found in this study revealed that 72 of them downregulate in the saltstressed plants, and the remaining had no significant differential expression (Supplementary Table 2). The genotypic specificity of the miRNAs behavior is usually evident since different genotypes express the same miRNA but at different levels [8,32]. According to Sunkar et al. [8], a miRNA that presents negative regulation during stress probably targets positive regulators of stress tolerance, generating an accumulation of gene products.
Some studies show that a single miRNA can selectively regulate its targets in a non-linear dose-dependent manner, so preferred mRNA targets may vary over the developmental stages, depending on the level of expression of regulatory miRNAs [21,37]. In Ho et al. [21], validation using RNA degradome data supports that a single miRNA may regulate multiple targets, and an mRNA may be regulated by more than one miRNA, suggesting a complex and fine-tuned interaction network between miRNAs and their targets at the post-transcriptional level. Ho and colleagues' results corroborate with ours in oil palm subjected to salt stress, as miR171, egu-miR03sds and egu-miR18sds showed differential expression, and each one of them had three putative target genes also differentially expressed; the same is true to a putative target gene regulated by multiple miRNAs.
Many miRNAs play a role in the response of plants to abiotic stresses, such as salinity, through posttranscriptional regulation, which has been the focus of several studies [1,8,22,38]. Different reports demonstrate that several miRNA-target genes are transcription factor mRNAs (TFs), indicating miRNAdependent post-transitional regulation during the development and response to the environment [39,40]. In plants, approximately 7-10% of the genes code for TFs at distinct moments, and dozen plant TF gene families precisely coordinate the spatial and temporal expression of downstream genes associated with abiotic stress [41,42]. The present study showed that miRNAs miR166, miR169, miR319, miR396, miR529, and egu-miR24sds showed altered expression profiles in young oil palm plants subjected to salt stress; and, through functional analysis, that they regulate TFs transcript levels, which in turn affect the protein levels of the TFs. The squamosa promoter-binding-like protein 17 homolog gene in oil palm (LOC105039459) is the putative target gene to miR529. While this miRNA underwent an 80% decrease in expression, its target gene showed a 70% increase in the apical leaf of young oil palm plants under very high salt stress. Squamosa promoter-binding (SBP) and SBP-Like (SPL) proteins are putative transcription factors having a plant-specific SBP domain consisting of 76 amino acids in length that regulates several biological processes, including salinity stress. Hou et al. [43] overexpressed the VpSBP16 gene from grape (Vitis vinifera) in A. thaliana and observed an enhancement in the tolerance to salt and drought stress during seed germination, as well in seedlings and mature plants, by regulating SOS and ROS signaling cascades. On the other hand, the OsSPL10 and CaSBP12 genes negatively control salt tolerance in rice and pepper, respectively [29,44].
The transcription factor bHLH143 homolog gene in oil palm (LOC105048659) is the putative target gene of egu-miR24sds. While this miRNA underwent a 65% decrease in expression, its putative target gene showed a 63% increase in the apical leaf of young oil palm plants under very high salt stress. The basic helix-loop-helix (bHLH) TFs are a large gene family in the plant genome, and some of these TFs regulate plant responses to abiotic stresses, including salt stress. The overexpression of the SlbHLH22 gene in tomato plants increased the tolerance to drought and salinity stress by improving the ROS scavenging system, increasing osmotic potential, and enhanced the accumulation of secondary metabolites [45]. Qiu et al. [46] overexpressed the MfbHLH38 gene -from the resurrection plant Myrothamnus flabellifoliain Arabidopsis and observed enhanced tolerance to both drought and salinity stresses through increasing water retention ability, regulating osmotic balance, decreasing stress-induced oxidation damage, and possibly participated in ABA-dependent stress-responding pathway.
The GATA transcription factor 27 homolog gene in oil palm (LOC105043768) is the putative target gene to miR319. While this miRNA underwent a 74% decrease in expression, its target gene showed a 44% increase under very high salt stress. GATA TFs belong to one of the most conserved families of zinc-finger TFs [47]. Transcript abundance analysis using salt-sensitive and salttolerant rice genotypes indicated differential expression of GATA TF genes in response to various abiotic stresses such as salinity, drought, and exogenous ABA, suggesting inherent roles of diverse GATA factors in abiotic stress signaling [48]. Nutan et al. [49] have shown that the overexpression of the OsGATA8 gene results in salinity tolerance in rice seedlings, as it maintains ion homeostasis and restricts membrane damage.
The homeobox-leucine zipper protein HOX32 homolog gene in oil palm (LOC105046708) is the putative target gene to four distinct miRNAs from the miR166 family. While these miRNAs underwent a 70-87% decrease in expression, its target gene showed a 78% increase in the apical leaf of young oil palm plants under very high salt stress. Arabidopsis thaliana has four distinct classes of homeodomain leucine zipper (HD-ZIP) transcription factors -HD-ZIPI to HD-ZIPIV -organized in multi-genes families [50]. Bhattacharjee et al. [51] carried out a functional analysis of two candidates stressresponsive HD-ZIP I class homeobox genes from rice, OsHOX22, and OsHOX24, and showed that these genes were highly upregulated under various abiotic stress conditions, including salinity stress, at different stages of development, including seedling, mature and reproductive stages. Besides that, Bhattacharjee and colleagues also overexpressed the OsHOX24 gene in Arabidopsis plants showing that its overexpression does not result in a detectable difference in the phenotype and various growth parameters compared to the wild type under normal growth conditions; however, it does result in higher sensitivity to salinity stress.
The growth-regulating factor 10 homolog gene in oil palm (LOC105054175) is the putative target gene to miR396. While this miRNA underwent an 83% decrease in expression, its target gene showed a 131% increase in the apical leaf of young oil palm plants under very high salt stress. microRNA miR396 controls the expression of several growth-regulating factors (GRFs), and the GRF-miRNA396 regulatory module appears to be central to several developmental processes, including flower and seed formation, root development, and the coordination of growth processes under adverse environmental conditions, including salt stress [52][53][54]. Genetically modified creeping bentgrass (Agrostis stolonifera) overexpressing Osa-miR396c, a rice miRNA396 gene, showed enhanced salt tolerance associated with improved water retention, increased chlorophyll content, cell membrane integrity, and Na + exclusion during high salinity exposure; however, they exhibited altered development [53]. RNA-sequencing analysis revealed that GRF1 and GRF3 regulate the expression of many clock core genes and genes with stress-and defense-related functions [55]. AtGRF7 -a repressor of stress-responsive genes under non-stress conditions -suppresses DREB2A expression to preserve plant growth rate [56]. DREB2A is a TF whose transcriptional and post-translational activation increases osmotic stress tolerance in Arabidopsis [55]. atgrf7 lost function mutants are more tolerant to drought and salinity stresses [52].
Nuclear factor Y (NF-Y) proteins are widespread in plants, animals, and other eukaryotes and are also known as CCAAT Binding Factor (CBF) or Heme Activator Protein (HAP); and they modulate the expression of downstream putative target genes via two main mechanisms [57]. The heterotrimer -NF-YA-YB-YC -binds to the CCAAT box present in the promoter region of the downstream putative target genes through NF-YA and regulates the expression of the putative target genes. The idea of NF-YA competing with TFs, and suppressing the formation of the NF-YB-YC-TF complex, was postulated [58]; however, according to Zhao et al. [57], there is still no direct molecular evidence to support it.
Different members of the NFY gene family, including NF-YA, are targets of the miR169 family, and studies have shown that overexpression of NF-YA in Arabidopsis increased the plant's tolerance to salt stress, increasing the expression of abscisic acid [59]. A hypothetical model presented by Leyva-González et al. [58] proposes that in plants growing under non-stress conditions, NF-YA expression is low due to high levels of miR169 but sufficient to activate the transcription of genes which promoters contain the CCAAT box. In plants exposed to abiotic stress, NF-YA levels increase due to their transcriptional activation and to the reduction in the miR169 levels. Increased NF-YA levels repress 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, and participating in the activation of a late one.
Our results showed that the NF-YA3 homolog gene (LOC105056468) expression level in the apical leaf of salt-stressed young oil palm increased 158%, while miR169 had its expression level decreased to 17%. The gene that expresses the miR169 targetting the NF-YA3 homolog gene in oil palm is at two places in the oil palm genome, chromosomes 08 and 13, but only the one in the former chromosome differentially expressed under salinity stress (Supplementary Table 2). As the salinity stress reduces the amount of miR169 in the leaves of young oil palm plants, we postulate that more NF-YA3 would be available to compete with any NF-YB-YC-TF complex, resulting in more of the NF-YA-YB-YC complex, which could restore some of the main biological functions of this complex, such as drought tolerance.
A high concentration of soluble salts in the soil can directly affect plant growth in two distinct phasesosmotic and ionic -, whose duration and intensity vary according to the plant species and salt levels [60]. There is a rapid reduction in the osmotic potential in the osmotic stress phase that restricts water absorption and, therefore, reduces transpiration rates [60,61]. Salinity in its first phase of salt stress is much similar to that of drought stress, and many common responses between salinity and drought stresses are also expected [62]. In the present study, the stressed plants showed a rate of evapotranspiration about half of the one in the control ones, which shows that the young oil palm plants were experiencing the osmotic stress at 12 DAT [24]. The ionic phase, on the contrary, occurs more slowly and depends not only on the saline concentration but also on the exposure time and on the plant's capacity to accumulate or expel toxic ions [60,63].
The young oil palm plants used in this study had been for 12 days under salinity stress when collecting leaves for the transcriptome characterization [24], which can be considered a short period when dealing with a perennial crop. Those plants had already shown -in the highest level of NaCl used -premature senescence, chlorosis, and necrosis of adult leaves, and consequently a reduction in the photosynthetic area available to support continued growth (Fig. 2). Such symptoms result from a high Na + level in the plant that disrupts protein synthesis and interferes with enzyme activity [25]. In the case of these young oil palm plants, one can see an increase of almost 4X in Na + and 2X in Cl − in the absorption roots in the highest level of NaCl used, but not in the apical and basal leaves; showing that these plants were already starting to experience ionic stress [24].

Conclusion
This comprehensive and large-scale miRNA analysis characterized the miRNA population present in the leaves of young oil palm plants exposed to a high level of salt stress, to identify miRNA-putative target genes in the oil palm genome, and to perform an in silico comparison of the expression profile of the miRNAs and their putative target genes, resulting in: a) The identification of 79 miRNAs, 52 known miRNAs, and 27 new ones; 72 of them differentially expressed under salinity stress. The new ones received the names egu-miR(01to27)sds, where egu is the abbreviation of Elaeis guineensis and sds stands for salinity and drought stress; b) The prediction of 229 distinct genes as the targets to these 79 miRNAs in the oil palm genome; 150 of them were target to just one miRNA and the remaining 79 to two or more. Fifty-one miRNA-putative target genes differentially expressed under salinity stress; c) The functional annotation of 24 putative target genes upregulated under salinity stress. Among these genes, there were six that code for transcription factors and three for lncRNA; and.
d) The identification of potential targets genes -based upon evidence of a target gene-miRNA interaction under salinity stress -that can be tested as candidate genes to develop salinity stress tolerant oil palm plants. The development of salt-tolerant oil palm genotypes can come from overexpression or knock out of some of these miRNA or their respective putative target genes, either by a CRISPR/Cas genome editing strategy or by employing classic Agrobacterium-or biolistic-mediated genetic modification.

Plant material and growth conditions
The oil palm plants used in this study were clones regenerated in our lab out of embryogenic calluses obtained from leaves of an adult plant belonging to the E. guineensis genotype AM33, a Deli x Ghana from ASD Costa Rica (http:// www. asd-cr. com). The protocols and procedures implemented to regenerate the plants are described in Corrêa et al. [64]. Plants were kept in black plastic pots (5 L), containing 1700 g of a mix of vermiculite, soil, and a commercial substrate (Bioplant ® ), in a 1:1:1 ratio on a dry basis, and fertilized using 2.5 g/L of the formula 20-20-20. Before starting the experiments, plants were standardized accordingly to the developmental stage, size, and number of leaves. The experiment was performed in a greenhouse at Embrapa Agroenergy (www. embra pa. br/ en/ agroe nergia) in Brasília, DF, Brazil (S-15.732°, W-47.900°). The main environmental variables (temperature, humidity, and radiation) measured at a nearby meteorological station (S-15.789°, W-47.925°) fluctuated according to the weather conditions. The oil palm plants used in this study were in the growth stage known as "bifid saplings" when subjected to salt stress.

Experimental design and saline stress
The experiment was carried out in March 2018 and consisted of five treatments (0.0, 0.5, 1.0, 1.5, and 2.0 g of NaCl per 100 g of substrate), with four replicates in a completely randomized design. For details regarding moisture content, field capacity, and electric conductivity in the substrate, determined preliminarily, see Vieira et al. [24].
To salinize the substrate, the amount of NaCl corresponding to the level to be applied to each treatment was dissolved in an amount of tap water standardized and calculated by the difference between the amount of water previously present in the fresh substrate and the amount of water retained for the substrate to reach field capacity. Applying the right amount of water to get the substrate field capacity was a means of ensuring that there was no extravasation of the solution and loss of Na + or Cl − . Thus, the amount of salt added would remain in the substrate.
Plants were under stress for 12 days, with daily water maintenance by replacing the lost volume with tap water. The difference between total weight (TW) (container, soil, water added to reach field capacity, and plant weights, altogether) and the daily weight (DW) is equal to the amount of water necessary to replace daily water losses due to evapotranspiration. Such a procedure was essential to allow the same level of electric conductivity and water potential accordingly to the dose of salt added to the substrate.

Gas exchange measurements
Gas exchange was measured on the middle third of the apical leaf, in a previously marked area, between 9:00 and 11:00 a.m. [24]. The parameters of leaf gas exchange [net CO 2 assimilation rate (A), transpiration rate (E), stomatal conductance to water vapor (gs), and intercellular CO 2 concentration (Ci)] were measured by a portable infrared gas analyzer LI-COR Mod. 6400XT (LI-COR, Lincoln, NE, USA) equipped with a measuring chamber (2 × 3 cm) with artificial light system LI-COR Mod. 6400-02B. The extracted data was provided by the OPEN software version 6.3. The block temperature was 25 °C, PAR was 1500 μmol/m 2 /s, the relative humidity of the air inside the measuring chamber was between 50 and 60%, the airflow index was 400 μmol/s, and the CO 2 concentration was 400 ppm in the reference cell, using the model 6400-01 CO 2 mixer with cylinder CO 2 (7.5 g). After submitting the gas exchange data to the Kruskal-Wallis test, we applied the Dunn's test (p < 0.05) to those data with significant differences between treatments.

Transcriptomics
Apical leaves from three control and stressed plants (0.0 and 2.0 g of NaCl per 100 g of substrate), collected 12 days after imposition of the treatments (DAT), were immediately immersed in liquid nitrogen and then stored at − 80 °C until RNA extraction, library preparation, and sequencing.

Total RNA extraction and quality analysis, library preparation and sequencing
Total RNA was isolated from oil palm leaves using the Qiagen RNeasy ® Plant Mini kit (QIAGEN, CA, USA) following the manufacturer's protocol. RNA quantity and quality were measured using a Nanodrop Qubit 2.0 Fluorometer (Life Technologies, CA, USA) and an Agilent Bioanalyzer Model 2100 (Agilent Technologies, Palo Alto, CA). The GenOne Company (Rio de Janeiro, RJ, Brazil) performed the RNA-Seq using an Illumina HiSeq platform and the paired-end strategy. The Functional Genomics Center / ESALQ-USP (Piracicaba, SP, Brazil) performed the small RNAs sequencing using an Illumina HiSeq platform.

RNA-Seq data analysis
The OmicsBox version 1.3 [65] was employed to perform all RNA-Seq analyses. We used FastQC [66] and Trimmomatic [67] for quality control, filter reads, and remove low-quality bases. The oil palm reference genome [27] files downloaded from NCBI (BioProject PRJNA192219; BioSample SAMN02981535) on October 2020 -was used to align the RNA-Seq data using default parameters from OmicsBox version 1.3 through software STAR [68]. The default parameters from OmicsBox version 1.3 through HTSeq version 0.9.0 were employed to quantify expression at the gene or transcript level [69]. The pairwise differential expression analysis between experimental conditions (Control vs. Stressed) was performed through edgeR version 3.28.0 [70], applying a simple design and an exact statistical test without a filter for low counts genes.

miRNAs data analysis
The small RNA raw data was submitted to the cutadapt software version 2.7 [71], generating adapter-free small RNA reads 20-24 nucleotides long. The Rfam version 12.0 database was used to remove contaminants, followed by mapping to the oil palm reference genome [27] using Bowtie2 [72].
All adapter-free small RNA sequences (stressed and control) were concatenated into a single file for miRNA prediction. The prediction was then made using mireap version 0.2 (https:// sourc eforge. net/ proje cts/ mireap) and Shortstack version 3.4 (https:// github. com/ MikeA xtell/ Short Stack), independently or in an association. Both programs generate clusters of sequences lined up in genomic regions. Ideally, these clusters indicate the genomic location and the miRNA precursor, mature miRNA, and miRNA* sequences. Shortstack also analyzes precursor and hairpin metrics formed according to parameters established by Axtell and Meyers and classifies them in Y (confirmed miRNA) or N1-N15, where N15 means that the candidate has all the correct metrics, but the miRNA* is absent [73]. The clusters formed by the mireap were analyzed by Shortack to obtain the classifications of each miRNA. StrucVis (https:// github. com/ MikeA xtell/ struc Vis) was used in sequences classified as Y or N15 by ShortStack and/or ShortStack-mireap for structural evaluation of miRNA. Finally, manual curation was made of all miRNAs classified as Y and N15. The length of the strings, the predicted structure of the hairpin, and the annotation by homology were evaluated (miRBase -http:// www. mirba se. org/ search. shtml).
The prediction of miRNA-putative target genes was performed using the psRNA-Target online program, version 2 (https:// bio. tools/ psrna target), with the following parameters: 5 of top targets, 5 expectation, 1 Penalty for other mismatches. For the analysis of differential expression of miRNAs, we used the NOISeq R package [74]. For this, the individual counts of each sample were used as input. The genes that showed p values ≥0.95 were designated as differentially expressed.
To functionally annotate the differentially expressed miRNA-putative target genes we used the LOC id to get to the protein sequence at NCBI, and then submitted it to the InterProScan search at InterPro (http:// www. ebi. ac. uk/ inter pro/) [75].