Expression Proﬁle of Selected Genes Involved in Na + Homeostasis and In Silico miRNA Identiﬁcation in Medicago sativa and Medicago arborea under Salinity Stress

: The accumulation of ions due to increased salinity in the soil is one of the major abiotic stressors of cultivated plants that negatively affect their productivity. The model plant, Medicago truncatula , is the only Medicago species that has been extensively studied, whereas research into increased salinity adaptation of two important forage legumes, M. sativa and M. arborea , has been limited. In the present study, the expression of six genes, namely SOS1 , SOS3 , NHX2 , AKT , AVP and HKT1 was monitored to investigate the manner in which sodium ions are blocked and transferred to the various plant parts. In addition, in silico miRNA analysis was performed to identify miRNAs that possibly control the expression of the genes studied. The following treatments were applied: (1) salt stress, with initial treatment of 50 mM NaCl and gradual acclimatization every 10 days, (2) salt shock, with continuous application of 100 mM NaCl concentration and (3) no application of NaCl. Results showed that M. arborea appeared to overexpress and activate all available mechanisms of resistance in conditions of increased salinity, while M. sativa acted in a more targeted way, overexpressing the HKT1 and AKT genes that contribute to the accumulation of sodium ions, particularly in the root. Regarding miRNA in silico analysis, ﬁve miRNAs with signiﬁcant complementarity to putative target genes, AKT1 , AVP and SOS3 were identiﬁed and served as a ﬁrst step in investigating miRNA regulatory networks. Further miRNA expression studies will validate these results. Our ﬁndings contribute to the understanding of the molecular mechanisms underlying salt-responsiveness in Medicago and could be used in the future for generating salt-tolerant genotypes in crop improvement programs.


Introduction
Alfalfa (Medicago sativa L.) is the most extensively cultivated perennial forage legume with high ecological importance [1]. Tree medick (Medicago arborea L.) is a Mediterranean stress-tolerant leguminous fodder shrub [2] that could play an important role as a genetic resource for the improvement of Medicago sp. salt tolerance [3]. Salinity is a major environmental abiotic stress factor in crop productivity worldwide that impacts plant development and yield via two main mechanisms: hyperosmotic pressure, which is caused by low water availability, and ion toxicity (especially Na + ), arising from solute imbalances [4]. However, little is known regarding the molecular mechanisms involved in the responses of Medicago sp. to salinity except Medicago truncatula.
Salt tolerance in higher plants is a quantitative trait that is regulated by polygenes and numerous signaling pathways, and a complex interaction system exists between the genes Figure 1 displays the expression levels of six salt stress-related genes in M. sativa and M. arborea. Information and additional data regarding phenotypical and physiological differences among the two species have already been published in previous works [3,19].
From the SOS gene family, we studied two genes, namely SOS1 and SOS3, which are directly involved in the regulation of salt tolerance by eliminating Na + from plant roots. According to our results, M. arborea showed generally higher gene expression values than M. sativa for all treatments, for both shoots and roots, while M. sativa overexpressed the genes more targeted focusing on roots. However, the highest SOS1 and SOS3 gene expression values appeared in plants (both roots and shoots) that were subjected to 100 mM NaCl. (Figure 1(a1,b2)). According to Mäser et al. [20], overexpression of SOS genes in M. falcata and M. truncatula shows that Na + release from cells is an important factor in plant survival under saline conditions, and activation of this mechanism can improve plant resistance. expression values appeared in plants (both roots and shoots) that were subjected to 100 mM NaCl. (Figure 1(a1,b2)). According to Mäser et al. [20], overexpression of SOS genes in M. falcata and M. truncatula shows that Na + release from cells is an important factor in plant survival under saline conditions, and activation of this mechanism can improve plant resistance. Comparisons between genotypes differ significantly at p < 0.05. Means followed by a different letter differ significantly (p < 0.05) based on the Tukey analysis. Relative expression ratios were calculated using the 2 −∆∆CT method [21]. The first two bars are for control treatment, the next two bars are for the 50-100 mM NaCl treatment and the final two bars are for the 100 mM NaCl treatment.
They also argue that Na + /H + antiporter SOS genes had a role in counteracting root to shoot Na + transport. SOS genes are expressed in the root and leaf vasculature, predominately in the root tips, and have been considered to re-extract Na + from the root xylem under salt stress via a proposed reversal of the SOS transporter. Sandhu [10] highlighted the major role that AKT and HKT1 genes play in the transport of Na + from shoots to roots. In this study, the two genes were overexpressed in the roots of M. sativa while the expression values for both genes in leaves were very low for all treatments (Figure 1(c1,f1)). Specifically for the HKT1 gene, the most elevated expression level was observed in the roots of plants subjected to 100 mM NaCl. Similarly, ref. [20] highlighted that in Arabidopsis the HKT1 gene was strongly expressed in roots. Inversely, the overexpression of these two genes in M. arborea plants was very profound for both leaves and roots of the stressed plants, especially in plants which were subjected to salt shock ( Figure 1(c2,f2)). Additionally, it was [10] demonstrated by testing 12 genotypes of M. sativa that the two genotypes of highest salinity resistance overexpressed the AKT and HKT1 genes, mainly in roots compared to leaves, suggesting the contribution of these genes to salinity tolerance. Mäser et al. [21] also showed that in the Arabidopsis thaliana species, K + starvation induced HKT1 expression in root cortex cells and enhanced inward Na + currents that included HKT1-like properties in root cortex cells; therefore, disruption of HKT1 facilitates the net transfer of Na + to shoots. Furthermore, HKT1 could counteract shoot-ward Na + transport either by sequestering Na + in the roots and/or preventing it from being loaded into the xylem by re-extracting Na + that is being loaded into the xylem, or by removing Na + from leaves and transporting Na + to roots via the phloem.
The NHX2 and AVP genes are known to play an important role in isolating sodium in the vacuole to protect the cytoplasm from toxic effects. Overexpression of NHX2 gene in alfalfa transgenic plants [8,22] and the AVP gene in Arabidopsis transgenic plants [23,24] improved salinity tolerance by improving intracellular potassium homeostasis. Additionally, several studies in cultivated crops highlighted the importance of AVP overexpression for providing increased plant yield under salinity and drought stress conditions [25][26][27]. In this study, M. sativa demonstrated the highest gene expression value for NHX2 in both roots and leaves for plants subjected to gradually increased salt concentrations (Figure 1(d1)). On the contrary, regarding AVP expression profile in M. sativa, the highest up-regulation was observed in the roots of plants which were treated with 100 mM NaCl (Figure 1(e1)). In contrast to M. sativa, gene expression profiles for all M. arborea salt treatments were clearly elevated ( Figure 1(d2,e2)).
Considering all gene expression values under increased salinity conditions for both M. sativa and M. arborea, we can conclude that M. arborea seems to overexpress all examined genes in both shoots and roots, while M. sativa acts in a more targeted manner, overexpressing specific genes (HKT1 and AKT), which help the accumulation of Na + more specifically in roots. The same trend was observed in another study [3] in which gene expression related to ion transport was carried out under salinity stress. Both studies agree that M. sativa uses a more cost-efficient strategy to respond to salinity by activating specific genes.

miRNA In Silico Analysis
Given the importance of miRNAs in the regulation of gene expression during plant development and in response to abiotic stresses, we performed in silico analysis by use of the 'psRNATarget' tool, to uncover miRNAs that may target the genes under study (Table 1). The analysis used miRNA and mRNA sequences from M. truncatula since most Medicago miRNAs sequences in the miRNA databases have been derived from this model leguminous plant, as well as due to its high genetic relationship with related species of the Medicago genus, such as M. sativa and M. arborea. Significant complementarity was detected between the conserved miRNA mtr-miR156a and the mRNA of the gene MtAVP (Pyrophospate-energized vacuolar membrane proton pump) in 15 consecutive positions (1-15 bp) from the 5 end of the miRNA (Table 1). G-U wobbles are allowed in miRNA-target complementarities, and in plants, it seems that complementarity in at least 1-13 consecutive bases from the 5 end of the miRNA may be associated with actual cleavage of the mRNA of the target gene and further degradation, or translational inhibition. This implies putative binding of mtr-miR156a to the MtAVP mRNA target and subsequent cleavage or translational inhibition of the transcript [28,29], leading to possible downregulation of the MtAVP gene.

miRNA In Silico Analysis
Given the importance of miRNAs in the regulation of gene expression during plant development and in response to abiotic stresses, we performed in silico analysis by use of the 'psRNAΤarget' tool, to uncover miRNAs that may target the genes under study ( Table  1). The analysis used miRNA and mRNA sequences from M. truncatula since most Medicago miRNAs sequences in the miRNA databases have been derived from this model leguminous plant, as well as due to its high genetic relationship with related species of the Medicago genus, such as M. sativa and M. arborea. Significant complementarity was detected between the conserved miRNA mtr-miR156a and the mRNA of the gene MtAVP (Pyrophospate-energized vacuolar membrane proton pump) in 15 consecutive positions (1-15 bp) from the 5′ end of the miRNA (Table 1). G-U wobbles are allowed in miRNAtarget complementarities, and in plants, it seems that complementarity in at least 1-13 consecutive bases from the 5′ end of the miRNA may be associated with actual cleavage of the mRNA of the target gene and further degradation, or translational inhibition. This implies putative binding of mtr-miR156a to the MtAVP mRNA target and subsequent cleavage or translational inhibition of the transcript [28,29], leading to possible downregulation of the MtAVP gene.  Similarly, a significant match was observed between miR167b and MtAKT1 mRNA in 17 consecutives bases from the 5′ end of the miRNA. A gap has been introduced, indicating that there may be a bulge in that position in the miRNA-target dimer structure. MiR156a and miR167b both belong to well-conserved plant miRNA families and have been implicated in salt responsiveness in a variety of species [30][31][32]. Based on the above observations, we could speculate that if miR156a targets AVP mRNA in M. sativa saltstressed leaf tissue, the AVP target is downregulated, presumably due to high expression of miR156a. On the other hand, in root tissue miR156a, it is probably downregulated and, consequently, the AVP gene is upregulated ( Figure 1). Likewise, if miR167b targets AKT1, a negative correlation between miR167b and AKT1 messenger RNA levels might be expected in salt-stressed leaves and roots in a manner equivalent to the regulation of the miRNA156-AVP putative module. However, these hypotheses await further experimentation and validation. Thus far, miR167 and miR156 have been implicated in the salt response as well as in other abiotic stresses and developmental processes in a variety

miRNA In Silico Analysis
Given the importance of miRNAs in the regulation of gene expression during plant development and in response to abiotic stresses, we performed in silico analysis by use of the 'psRNAΤarget' tool, to uncover miRNAs that may target the genes under study ( Table  1). The analysis used miRNA and mRNA sequences from M. truncatula since most Medicago miRNAs sequences in the miRNA databases have been derived from this model leguminous plant, as well as due to its high genetic relationship with related species of the Medicago genus, such as M. sativa and M. arborea. Significant complementarity was detected between the conserved miRNA mtr-miR156a and the mRNA of the gene MtAVP (Pyrophospate-energized vacuolar membrane proton pump) in 15 consecutive positions (1-15 bp) from the 5′ end of the miRNA (Table 1). G-U wobbles are allowed in miRNAtarget complementarities, and in plants, it seems that complementarity in at least 1-13 consecutive bases from the 5′ end of the miRNA may be associated with actual cleavage of the mRNA of the target gene and further degradation, or translational inhibition. This implies putative binding of mtr-miR156a to the MtAVP mRNA target and subsequent cleavage or translational inhibition of the transcript [28,29], leading to possible downregulation of the MtAVP gene.  Similarly, a significant match was observed between miR167b and MtAKT1 mRNA in 17 consecutives bases from the 5′ end of the miRNA. A gap has been introduced, indicating that there may be a bulge in that position in the miRNA-target dimer structure. MiR156a and miR167b both belong to well-conserved plant miRNA families and have been implicated in salt responsiveness in a variety of species [30][31][32]. Based on the above observations, we could speculate that if miR156a targets AVP mRNA in M. sativa saltstressed leaf tissue, the AVP target is downregulated, presumably due to high expression of miR156a. On the other hand, in root tissue miR156a, it is probably downregulated and, consequently, the AVP gene is upregulated (Figure 1). Likewise, if miR167b targets AKT1, a negative correlation between miR167b and AKT1 messenger RNA levels might be expected in salt-stressed leaves and roots in a manner equivalent to the regulation of the miRNA156-AVP putative module. However, these hypotheses await further experimentation and validation. Thus far, miR167 and miR156 have been implicated in the salt response as well as in other abiotic stresses and developmental processes in a variety

miRNA In Silico Analysis
Given the importance of miRNAs in the regulation of gene expression during plant development and in response to abiotic stresses, we performed in silico analysis by use of the 'psRNAΤarget' tool, to uncover miRNAs that may target the genes under study ( Table  1). The analysis used miRNA and mRNA sequences from M. truncatula since most Medicago miRNAs sequences in the miRNA databases have been derived from this model leguminous plant, as well as due to its high genetic relationship with related species of the Medicago genus, such as M. sativa and M. arborea. Significant complementarity was detected between the conserved miRNA mtr-miR156a and the mRNA of the gene MtAVP (Pyrophospate-energized vacuolar membrane proton pump) in 15 consecutive positions (1-15 bp) from the 5′ end of the miRNA (Table 1). G-U wobbles are allowed in miRNAtarget complementarities, and in plants, it seems that complementarity in at least 1-13 consecutive bases from the 5′ end of the miRNA may be associated with actual cleavage of the mRNA of the target gene and further degradation, or translational inhibition. This implies putative binding of mtr-miR156a to the MtAVP mRNA target and subsequent cleavage or translational inhibition of the transcript [28,29], leading to possible downregulation of the MtAVP gene.  Similarly, a significant match was observed between miR167b and MtAKT1 mRNA in 17 consecutives bases from the 5′ end of the miRNA. A gap has been introduced, indicating that there may be a bulge in that position in the miRNA-target dimer structure. MiR156a and miR167b both belong to well-conserved plant miRNA families and have been implicated in salt responsiveness in a variety of species [30][31][32]. Based on the above observations, we could speculate that if miR156a targets AVP mRNA in M. sativa saltstressed leaf tissue, the AVP target is downregulated, presumably due to high expression of miR156a. On the other hand, in root tissue miR156a, it is probably downregulated and, consequently, the AVP gene is upregulated (Figure 1). Likewise, if miR167b targets AKT1, a negative correlation between miR167b and AKT1 messenger RNA levels might be expected in salt-stressed leaves and roots in a manner equivalent to the regulation of the miRNA156-AVP putative module. However, these hypotheses await further experimentation and validation. Thus far, miR167 and miR156 have been implicated in the salt response as well as in other abiotic stresses and developmental processes in a variety

miRNA In Silico Analysis
Given the importance of miRNAs in the regulation of gene expression during plant development and in response to abiotic stresses, we performed in silico analysis by use of the 'psRNAΤarget' tool, to uncover miRNAs that may target the genes under study ( Table  1). The analysis used miRNA and mRNA sequences from M. truncatula since most Medicago miRNAs sequences in the miRNA databases have been derived from this model leguminous plant, as well as due to its high genetic relationship with related species of the Medicago genus, such as M. sativa and M. arborea. Significant complementarity was detected between the conserved miRNA mtr-miR156a and the mRNA of the gene MtAVP (Pyrophospate-energized vacuolar membrane proton pump) in 15 consecutive positions (1-15 bp) from the 5′ end of the miRNA (Table 1). G-U wobbles are allowed in miRNAtarget complementarities, and in plants, it seems that complementarity in at least 1-13 consecutive bases from the 5′ end of the miRNA may be associated with actual cleavage of the mRNA of the target gene and further degradation, or translational inhibition. This implies putative binding of mtr-miR156a to the MtAVP mRNA target and subsequent cleavage or translational inhibition of the transcript [28,29], leading to possible downregulation of the MtAVP gene.  Similarly, a significant match was observed between miR167b and MtAKT1 mRNA in 17 consecutives bases from the 5′ end of the miRNA. A gap has been introduced, indicating that there may be a bulge in that position in the miRNA-target dimer structure. MiR156a and miR167b both belong to well-conserved plant miRNA families and have been implicated in salt responsiveness in a variety of species [30][31][32]. Based on the above observations, we could speculate that if miR156a targets AVP mRNA in M. sativa saltstressed leaf tissue, the AVP target is downregulated, presumably due to high expression of miR156a. On the other hand, in root tissue miR156a, it is probably downregulated and, consequently, the AVP gene is upregulated (Figure 1). Likewise, if miR167b targets AKT1, a negative correlation between miR167b and AKT1 messenger RNA levels might be expected in salt-stressed leaves and roots in a manner equivalent to the regulation of the miRNA156-AVP putative module. However, these hypotheses await further experimentation and validation. Thus far, miR167 and miR156 have been implicated in the salt response as well as in other abiotic stresses and developmental processes in a variety

miRNA In Silico Analysis
Given the importance of miRNAs in the regulation of gene expression during plant development and in response to abiotic stresses, we performed in silico analysis by use of the 'psRNAΤarget' tool, to uncover miRNAs that may target the genes under study ( Table  1). The analysis used miRNA and mRNA sequences from M. truncatula since most Medicago miRNAs sequences in the miRNA databases have been derived from this model leguminous plant, as well as due to its high genetic relationship with related species of the Medicago genus, such as M. sativa and M. arborea. Significant complementarity was detected between the conserved miRNA mtr-miR156a and the mRNA of the gene MtAVP (Pyrophospate-energized vacuolar membrane proton pump) in 15 consecutive positions (1-15 bp) from the 5′ end of the miRNA (Table 1). G-U wobbles are allowed in miRNAtarget complementarities, and in plants, it seems that complementarity in at least 1-13 consecutive bases from the 5′ end of the miRNA may be associated with actual cleavage of the mRNA of the target gene and further degradation, or translational inhibition. This implies putative binding of mtr-miR156a to the MtAVP mRNA target and subsequent cleavage or translational inhibition of the transcript [28,29], leading to possible downregulation of the MtAVP gene.  Similarly, a significant match was observed between miR167b and MtAKT1 mRNA in 17 consecutives bases from the 5′ end of the miRNA. A gap has been introduced, indicating that there may be a bulge in that position in the miRNA-target dimer structure. MiR156a and miR167b both belong to well-conserved plant miRNA families and have been implicated in salt responsiveness in a variety of species [30][31][32]. Based on the above observations, we could speculate that if miR156a targets AVP mRNA in M. sativa saltstressed leaf tissue, the AVP target is downregulated, presumably due to high expression of miR156a. On the other hand, in root tissue miR156a, it is probably downregulated and, consequently, the AVP gene is upregulated (Figure 1). Likewise, if miR167b targets AKT1, a negative correlation between miR167b and AKT1 messenger RNA levels might be expected in salt-stressed leaves and roots in a manner equivalent to the regulation of the miRNA156-AVP putative module. However, these hypotheses await further experimentation and validation. Thus far, miR167 and miR156 have been implicated in the salt response as well as in other abiotic stresses and developmental processes in a variety Similarly, a significant match was observed between miR167b and MtAKT1 mRNA in 17 consecutives bases from the 5 end of the miRNA. A gap has been introduced, indicating that there may be a bulge in that position in the miRNA-target dimer structure. MiR156a and miR167b both belong to well-conserved plant miRNA families and have been implicated in salt responsiveness in a variety of species [30][31][32]. Based on the above observations, we could speculate that if miR156a targets AVP mRNA in M. sativa salt-stressed leaf tissue, the AVP target is downregulated, presumably due to high expression of miR156a. On the other hand, in root tissue miR156a, it is probably downregulated and, consequently, the AVP gene is upregulated (Figure 1). Likewise, if miR167b targets AKT1, a negative correlation between miR167b and AKT1 messenger RNA levels might be expected in saltstressed leaves and roots in a manner equivalent to the regulation of the miRNA156-AVP putative module. However, these hypotheses await further experimentation and validation.
Thus far, miR167 and miR156 have been implicated in the salt response as well as in other abiotic stresses and developmental processes in a variety of species through the targeting and regulation of transcription factor genes, mainly ARF (AUXIN RESPONSE FACTOR) and SQUAMOSA LIKE PROTEIN (SLP), respectively [33][34][35][36][37][38]. Research on miRNA regulation of gene expression has focused on transcription factor target genes, and modules such as miR156-SPL, miR159-MYB, miR164-NAC, miR167-ARF and miR172-AP have been widely investigated. However, to date, little information exists regarding the engagement of other target genes in these regulatory processes through miRNA-mRNA modules, including the genes under study. An ever-increasing number of miRNAs and putative target genes identified in recent years [33,34] will facilitate the study of other miRNAtarget gene pairs and their function in salt-tolerance. Our analysis of AKT1, AVP and SOS3 will provide additional information to these studies and contribute towards this goal. Furthermore, candidate miRNAs targeting SOS3 seem to be the miR2605, mir2670f and miR2670g. Nevertheless, there is no information thus far regarding these sequences, and they may represent species-specific miRNAs that need to be studied further.
This work highlights the importance of Na + compartmentation and the action of the SOS1, SOS3, NHX2, AKT, AVP and HKT1 genes towards salinity tolerance. MiRNAs are also known to contribute to salinity resilience. Expression studies via qPCR experiments for the identified miRNA to validate their significance will possibly enable their future use into breeding programs.

Plant Material and Seed Pretreatment
The M. sativa population used in this study was a hybrid of two M. sativa parents [3]. Seeds from both M. sativa and M. arborea were obtained from E. Bingham, Agronomy Department, University of Wisconsin-Madison, Madison, WI, USA. The germplasm of M. arborea is a collection from Greece [3]. The seeds were scarified with absolute sulfuric acid for 10 min and then rinsed thoroughly with sterile distilled water. Bleach solution (3%) was added for 1.5-2 min and then rinsed thoroughly. The scarified seeds were placed on Petri dishes containing 0.5% agar at 4 • C overnight, and then transferred to 20 • C (dark) for 3-4 days.

Growth Conditions and Salt Stress Treatments
From each of the Medicago species, 60 seeds were planted, after their germination, in pots of 8.5 cm height and 10 cm in diameter. Each pot contained peat, Kronos N 50-300 mg/L, P 2 O 5 80-300 mg/L, K2O 80-300 mg/L, pH 5-6.5, salinity <1.75 g/L and one seed. The pots were then placed in a growth chamber at a constant temperature of 25 • C, and a photoperiod of 16 h of light and 8 h of darkness, where they were regularly watered and fertilized every week. For each treatment, 20 plants in seedling stage were used for each species in a completely randomized design. The following treatments were applied: (1) salt stress, with initial treatment of 50 mM NaCl and gradual step acclimatization 10 days after, with treatment of 100 mM NaCl, (2) salt shock, with continuous application of 100 mM NaCl concentration and (3) no application of NaCl. Plants were watered every five days after salt treatment initiation. Once per week, Hoagland solution [39] was added to the salt solution. The growth chamber experiment lasted 34 days. After the completion of the experiment, whole shoots, excluding cotyledons and the roots, free of pathogens, were cut and immediately transferred into liquid nitrogen in the middle of the light period. When harvested, all plants were in the vegetative stage, and roots did not show nodules.

RNA Isolation and cDNA Synthesis
For plant RNA isolation, the TRIzol reagent kit (Sigma-Aldrich, St. Louis, MO, USA) was used and RNA concentration was measured by Nanodrop Spectrophotometer (Thermo Scientific NanoDrop™ 1000 Spectrophotometer). For cDNA synthesis, Superscript II enzyme (Invitrogen, Carlsbad, CA, USA) was used. The quality of the cDNA was evaluated

Real-Time PCR Experiments
The StepOnePlus™ Real-Time PCR Systems (Applied Biosystems, Foster City, CA, USA) machine was used to detect relative transcript levels using SYBR Select Master Mix (Applied Biosystems, Foster City, CA, USA). The threshold cycle CT method was used, in which the relative amount of a target gene in the samples can be determined. In our experiment, six target genes were studied (SOS1, SOS3, NHX2, AKT, AVP and HKT1) and the primers that were used were the same as Sandhu et al., 2017. Primer selection was based on Sandhu et al., 2017 [10]. All reactions were performed in triplicates. Melting curves of the reaction products were generated and fluorescence data were collected at a temperature above the melting temperature of non-specific products. Relative expression levels of the studied genes were calculated according to the 2 −∆∆Ct method [20] (for calculating ∆∆Ct value, the control sample of each entry was used). Actin-2 was used as an internal control for normalization.

Statistical Analysis
The effect of variables on gene expression was tested using the one-way analysis of variance (ANOVA) followed by mean comparisons through Tukey's HSD post hoc test with p-values of p < 0.05. Values are expressed as mean ± SEM (standard error of mean). The statistical analysis was performed using the IBM SPSS Statistics 23 statistical package for Windows (SPSS Inc., Chicago, IL, USA).

In Silico miRNA Analysis
For the prediction of potential miRNAs that may target and regulate the expression of the genes under study, the "psRNATarget: A Plant Small RNA Target Analysis Server" tool (https://www.zhaolab.org/psRNATarget/-version, 20 December 2022) was employed. The default parameters were used and were as follows: maximum expectation value of 5, penalty for G:U pair of 0.5, penalty for other mismatches of 1, seed region of 2-13 nt, number of mismatches allowed in the seed region of 2.0 and translation inhibition ranges of 10 nt to 11 nt. Under these parameters, and following the analysis of the results, the most significant score which designates the best possibility for miRNA-mRNA pair formation was selected.