Identification of microRNAs and Their Expression in Leaf Tissues of Guava (Psidium guajava L.) under Salinity Stress

Abstract

Superfruit guava (Psidium guajava L.) is one of the healthiest fruits due to its high antioxidant dietary fiber and vitamin content. However, the growth and development of this plant are severely affected by salinity stress, mostly at the seedling stage. MicroRNAs (miRNAs) are small, noncoding, endogenous, highly conserved RNA molecules that play key regulatory roles in plant development, organ morphogenesis, and stress response signaling. In this study, applying computational approaches and following high stringent filtering criteria, a total of 40 potential microRNAs belonging to 19 families were characterized from guava. The identified miRNA precursors formed stable stem-loop structures and exhibited high sequence conservation among diverse and evolutionarily distant plant species. Differential expression pattern of seven selected guava miRNAs (pgu-miR156f-5p, pgu-miR160c-5p, pgu-miR162-3p, pgu-miR164b-5p, pgu-miR166t, pgu-miR167a-5p, and pgu-miR390b-5p) were recorded under salinity stress and pgu-miR162-3p, pgu-miR164b-5p as well as pgu-miR166t were found to be the most affected ones. Using the psRNATarget tool, a total of 49 potential target transcripts of the characterized guava miRNAs were identified in this study which are mostly involved in metabolic pathways, cellular development, and stress response signaling. A biological network has also been constructed to understand the miRNA mediated gene regulation using the minimum free energy (MFE) values of the miRNA-target interaction. To the best of our knowledge, this is the first report of guava miRNAs and their targets.

1. Introduction

MicroRNAs (miRNAs) are highly conserved, 20–24 nucleotides (nt) long non-coding RNA molecules that play critical roles in post-transcriptional gene regulation either by triggering transcriptional repression or via targeting mRNA degradation [1,2]. Plant microRNAs have been implicated in several biological functions, such as plant development [3], signaling [4], organ morphogenesis [5], secondary metabolite production [6], as well in the adaptation to abiotic and biotic stresses [7]. MicroRNA biogenesis in plants begins when miRNA genes are transcribed into long primary transcripts (pri-miRNAs) by the enzyme RNA polymerase II. Subsequently, the resulting pri-miRNAs are cleaved to generate stem-loop RNA precursors (pre-miRNAs) by ribonuclease III-like Dicer (DCL1) enzyme. Right away, DCL1 recognizes and cleaves further the hairpin loop of the pre-miRNAs and produces short double-stranded RNAs (dsRNAs) or duplexes, and finally, one strand of the mature miRNA duplexes associates with the RNA Induced Silencing Complex (RISC) to form a miRNA-ribonucleoprotein complex guided by Argonaute (AGO) protein to interact with the relevant mRNA targets [8,9].

The recent development of next-generation sequencing technology has led to the discovery of numerous miRNAs in several non-model plant species [8]; however, the overall procedure is expensive, time-consuming, and requires a high technical expertise. In the plant kingdom, multiple miRNAs are evolutionarily conserved, and this feature simplifies the process of characterization of novel miRNA orthologues in new plant species by identifying homologs [10,11,12]. Nonetheless, the only sequence-based in silico homology approaches to identify potential miRNAs in new plant species may yield false-positive results, and hence consideration must be given to the secondary structures as well as other parameters of the pre-miRNAs such as length, GC content, Minimum Folding Free Energy (MFE), and Minimum Folding Free Energy Index (MFEI) in order to increase the precision of the computational prediction by discriminating from other coding or non-coding RNAs [6,12,13]. However, experimental confirmation of the predicted miRNAs is strongly recommended [6,14,15].

Guava (Psidium guajava L.) is an important fruit crop belonging to the Myrtaceae family, originated from Mexico and distributed throughout Asia, Africa, Europe, and South America [16,17]. Guava is popularly known as the poor man’s apple due to its high nutritious value and low cost of cultivation. The compositional analysis showed that guava contains a variety of powerful health-promoting substances including flavonoids, phenols, tannins, saponins, triterpenes, lectins, carotenoids, essential oils, fatty acids, fibers, and vitamins [18,19]. However, salinity has been shown to be one of the major problems in guava cultivation, which impairs its growth and productivity [20]. Excess accumulation of salts in mature guava leaves has several detrimental consequences, such as quick leaf chlorosis, necrosis, and decreased photosynthetic activity [21].

Since miRNAs have been considered as crucial players in plants’ response towards salinity and other stress modulation, these molecules have been proposed as major genetic engineering targets to produce abiotic stress tolerant transgenic plants through loss-of-function or gain-of-function approaches [22,23]. Thus, profiling miRNAs in nonmodel plants is essential not only for understanding the regulation of various biological phenomena but also to explore their role in stress response signaling. To date, no scientific information about the guava miRNAs and their targets are available. Hence, using the recently published guava draft genome sequence (GenBank assembly accession GCA_002914565.1) several miRNAs and their corresponding targets have been characterized as well as their expression pattern under salinity stress were studied to gain a better understanding of the physiological role of miRNAs in guava.

2. Materials and Methods

2.1. Computational Prediction of Potential Guava miRNAs and Their Pre-miRNA Candidates

We performed an in silico analysis to predict and identify potential guava miRNAs using a reference set of plant miRNAs obtained from the miRbase database [24]. The reference set confined a total of 1580 known mature miRNAs from several plant species including Arabidopsis thaliana (428), Malus domestica (322), Glycine max (756), and Eugenia uniflora (74). The workflow is summarized in Figure 1. Briefly, the preceding set of known miRNAs were BLASTn against the guava genome; the sequences with exact matches were selected manually. The potential precursor sequences of approx 400 nt (200 nt upstream and 200 nt downstream to the BLAST hit region) were mined and protein-coding sequences were discarded. Stable secondary structures of the selected precursors were generated using the mFold web server [25] and the stability was evaluated by the previously demonstrated strict filtering criteria [26]: (i) the pre-miRNA must have a stem-loop structure containing the mature miRNA sequence within one arm; (ii) the potential mature miRNA should not be presented in the hairpin structures’ terminal loop, (iii) mature miRNA should have fewer than nine mismatches with the opposite miRNA* sequence [27], and (iv) the potential stem-loop candidate should have minimum negative folding free energy (MFE) or ΔG (−kcal/mol) and higher minimum folding free energy index (MFEI). The formula for calculating MFEI is as follows:

$MFEI=\frac{\left(MFE/length of RNA sequence\right) \times 100}{\%GC content}$

2.2. Phylogenetic and Conservation Analysis of Guava miRNA and Their Pre-miRNAs

To perform the conservation analysis of the predicted guava miRNAs and their precursor candidates, we retrieve the FASTA format of the miRNA precursor sequences of several plant species such as Amborella trichopoda (atr), Arabidopsis lyrata (aly), Arabidopsis thaliana (ath), Asparagus officinalis (aff), Brachypodium distachyon (bdi), Brassica napus (bna), Carica papaya (cpa), Citrus sinensis (csi), Cucumis melo (cme), Fragaria vesca (fve), Glycine max (gma), Linum usitatissimum (lus), Malus domestica (mdm), Manihot esculenta (mes), Medicago truncatula (mtr), Nicotiana tabacum (nta), Oryza sativa (osa), Populus trichocarpa (ptc), Prunus persica (ppe), Ricinus communis (rco), Theobroma cacao (tcc) and Vitis vinífera (vvi) available at miRbase and aligned accordingly. Multiple sequence alignment and phylogenetic tree construction (based on the Tamura-Nei model with 1000 boot-strapped replicates) were carried out using MEGA X software (version 10.0.5). The conservation analysis of the identified guava pre-miRNAs of miR160c, miR390b, miR396b was performed by the WebLogo tool [28] using their orthologs. Moreover, to elucidate the conserved nature of miRNAs across species and their cross-species transferability, a syntenic map was generated using potential guava pre-miRNAs against the well-annotated apple (Malus domestica) genome (GenBank assembly accession GCF_002114115.1), phylogenetically close species of guava [29].

2.3. Target Prediction of Guava miRNAs and Their Functional Annotations

For the prediction of the potential guava miRNA targets “Plant Small RNA Target Analysis Server” (psRNATarget) [30] was employed. The target transcript search was executed against the Malus domestica database due to the non-availability of the guava protein database on the psRNATarget list. The selection parameters were designated as follows: maximum expectation value of 3, translation inhibition ranges of 9 nt to 11 nt, number of top targets of 10, the penalty for G:U pair of 0.5, and number of mismatches allowed in the seed region of 1.5. Protein information of the matched sequences was obtained using UniProt BLAST. Consequently, gene ontology (GO) analysis of the potential guava target transcripts was executed, and the biological processes, cellular components, and molecular functions associated with each GO term were inferred using QuickGO [31]. Moreover, to know the coregulation of the potential targets, a biological network was generated using the MFE values of the miRNA-target interaction and visualized by Cytoscape 3.2 [32]. Finally, KEGG analysis [33] was performed to investigate the metabolic pathways and their networks regulated by the potential guava miRNAs with the Bi-directional Best Hit (BBH) method.

2.4. Plant Materials, Stress Treatment, RNA Extraction, and miRNA Expression Analysis

Guava seeds after surface sterilization (with 70% ethanol, 2.5% Sodium hypochlorite, and water) were germinated in Petri dishes and subsequently transferred into a hydroponic system containing Hoagland nutrient solution (pH 6.5) [34] and allowed to grow under controlled environmental conditions (25 °C, 70% humidity, and artificial illumination 250 μmol m⁻² s⁻¹ with a 12 h photoperiod) for four weeks. Consequently, half of the guava seedlings were transferred to a nutrient solution containing 200 mM NaCl for salinity stress, while the rest were placed into a nutrient solution without NaCl as control. Leaves of stressed and control plants were collected at 24, 48, and 72 h. Small RNA (<200 nt) was isolated from leaf tissues using the mirVana^TM miRNA Isolation kit (Thermo Scientific, Wilmington, NC, USA) following the manufacturer’s instructions and pooled separately each for stressed and control samples. The quality and quantity of RNA samples were checked with Nanodrop One (Thermo Scientific, Wilmington, NC, USA), and 1 µg of RNA for individual samples was subsequently polyadenylated (using modified oligo dT primer) and reverse transcribed using mRQ Buffer and enzyme provided with Mir-X miRNA First-Stand Synthesis kit (Takara, Tokyo, Japan). Prior to the qRT-PCR experiment, a No-RT (-RT) control PCR was performed to monitor any genomic DNA contamination. The qRT-PCR experiment was performed by Step One Real-Time PCR System (Applied Biosystems, Carlsbad, CA, USA) and Mir-X miRNA TB Green qRT-PCR kit (Takara, Tokyo, Japan) using the entire predicted miRNA sequence as a forward primer and the adapter-specific mRQ3′ primer provided with the kit as the reverse primer. Each reaction was made in 12.5 µL volume containing 1× SYBR Advantage Premix, 1× ROX dye, 0.2 µM each of forward and reverse primers as indicated above, and 2 µL of the first-strand cDNA. Seven miRNAs (pgu-miR156f-5p, pgu-miR160c-5p, pgu-miR162-3p, pgu-miR164b-5p, pgu-miR166t, pgu-miR167a-5p, and pgu-miR390b-5p), previously reported to have key roles in both biotic and abiotic stress responses [35] were selected for the qRT-PCR experiment, and the qPCR conditions were as follows: initial denaturation at 95 °C for 10 s followed by 45 cycles of denaturation at 95 °C for 5 s and annealing at 63 °C for 20 s, and finally a dissociation curve 95 °C for 30 s, 55 °C for 20 s and 95 °C for 20 s. The relative fold change values were obtained using the comparative Ct method or Ct (2^−ΔΔCT). Recently, the biological averaging method, in which individual biological replicates were replaced with pooled biological replicates, has been employed in several real-time PCR, microarray, and RNA-Seq experiments [36,37], and hence in this study, all the qPCR experiments were carried out with two pooled biological replicates and three technical replicates.

3. Results

3.1. Characterization of Guava miRNAs and Their Candidate Precursors

In this study using a high stringent filtering method, a total of 40 potential guava miRNAs belonging to 19 families were identified (Table 1). The majority of the identified guava miRNAs were 21 nt long. The precursors of guava miRNAs showed large variability in their size ranging from 72 to 215 nt with an average of 112 ± 32 nt. Guava miRNA pgu-miR159d exhibited the longest precursor length of 215 nt, while pgu-miR395-3p showed the shortest one of 72 nt. In this study, MFE values of precursors ranged from −33.1 to −95.2 kcal/mol with an average of −49.64 ± 14.39, whereas the MFEI values fluctuated between 0.70 to 1.36 with an average of 0.94 ± 0.14. The predicted secondary structures of guava miRNA precursors with higher MFEI values (top 10) are shown in Figure 2.

3.2. Conservation Analysis of Guava miRNAs and Their Potential Precursors

The recently identified guava miRNAs displayed a high degree of sequence homology (≤1 mismatch) to their respective homologs (orthologs) from several other monocots and dicot plant species (Figure 3). Moreover, high sequence conservation among the pre-miRNA orthologs was also noticed (Figure 4).

The conserved nature of the pre-miRNAs as well as mature miRNAs offers the chance to explore their evolutionary relationships (Figure 5). Phylogenetic analysis of pre-miRNAs pgu-160c-5p suggested its closeness to a group of plant species including apple, papaya, and black cottonwood pre-miRNAs (mdm-miR160c, cpa-miR160c, and ptc-miR160c), while pre-miRNA pgu-miR390b-5p is closer to both muskmelon (cme-miR390b) and flax (lus-miR390b). On the other hand, pgu-miR396-5p is closer to a group of plant species that included apple (mdm-miR396b), sweet orange (csi-miR396b), papaya (cpa-miR396b), cassava (mes-miR396b), black cottonwood (ptc-miR396b), muskmelon (cme-miR396b), and cacao (tcc-miR396b) (Figure 5). Furthermore, the comparative synteny map indicated the widespread distribution of potential guava miRNA orthologs in the apple genome demonstrating their cross-species transferability during the course of evolution (Figure 6).

3.3. Predicted Targets for Guava miRNAs and Their Functional Annotations

In this study using the psRNAtarget tool, a total of 49 potential target transcripts of guava miRNAs were identified. Important targets include squamosa promoter-binding like proteins/SPBs/SPLs, auxin response factors (ARFs), NAC domain protein, nuclear transcription factor Y, WRKY, and myb-like, laccase, thioredoxin, cytochrome f, among others. However, GO analysis of the predicted targets was conducted to obtain a deeper understanding of the miRNA function, which helps to unravel the biological mechanism, molecular function, and cellular component regulatory network of the miRNA gene (Figure 7). GO enrichment analysis highlighted different targets with molecular functions such as binding activity (protein binding, DNA binding, metal ion binding, etc.), catalytic activity, kinase activity, oxidoreductase activity, and structural activity (Figure 7A) are involved in significant biological processes such as transcription regulation, metabolic processes, transport, and oxidation-reduction (Figure 7B) in guava. The cellular components involved were found to be the nucleus, plasma membrane, cytoplasm, Golgi apparatus, ribosome, apoplast, chloroplast, and proteasome (Figure 7C). By conducting KEGG analysis, a deeper insight into the related biosynthetic and metabolic pathways was obtained (Figure 8). The KEGG analysis revealed that the potential guava miRNA targets are involved in a total of 22 different metabolic pathways in plants and animals. “Plant hormone signal transduction” was the most significantly enriched, followed by “plant-pathogen interaction” and “MAPK signaling pathway“ (Figure 8). Moreover, the coregulation of several potential target genes was observed by gene network analysis (Figure 9).

3.4. Expression Analysis of Guava miRNAs under Salinity Stress

To address whether the salinity stress influences the expression of the selected guava miRNAs (pgu-miR156f-5p, pgu-miR160c-5p, pgu-miR162-3p, pgu-miR164b-5p, pgu-miR166t, pgu-miR167a-5p, and pgu-miR390b-5p) in leaves, a qRT-PCR experiment was carried out and the results showed the differential expression of all the 7 miRNA under salinity stress. The expression levels of pgu-miR156f-5p, pgu-miR160c-5p, pgu-miR162-3p, pgu-miR164b-5p, pgu-miR166t, pgu-miR167a-5p were downregulated, while pgu-miR390b-5p was upregulated (Figure 10). Pgu-miR162-3p and pgu-miR164b-5p expressions have been shown to be most influenced by salinity stress, whereas the expressions of others have been less affected (Figure 10).

4. Discussion

The fundamental role of miRNAs in plant development and adaptation to environmental stresses has positioned them as an interesting molecule of study. Although miRNAs have been widely studied in several important crops and model plant species such as Arabidopsis, soybean, wheat, rice, and Tobacco [38,39,40,41,42], none of the systematic studies has been performed on guava so far. The recently identified guava miRNAs, as well as their precursors, displayed sequence conservation with their orthologs from different monocot and dicot plant species. This result indicates that between monocotyledonous and dicotyledonous species, miRNAs are universally conserved and may perform the same physiological role [43,44].

In this study, precursors of guava miRNAs displayed large variability in the size ranging from 72 to 215 nt as well as formed stable stem-loop secondary structures corroborating with the data reported in several other species including maize, cotton, soybean, flax, and passion fruit [12,45,46,47]. Moreover, all the precursors showed high MFEI values (0.70–1.36) with an average of 0.94 which is much higher than that of mRNAs (0.62–0.66), tRNAs (0.64), or rRNAs (0.59) [48] thus ruling out the possibilities of being other noncoding RNAs, furthermore, some experimental studies also confirmed that plant miRNAs were found to be correlated with high valued MFEI precursors [8,49]. Similarly, in our analysis, the prevalence of uracil (70%) at the first position of predicted miRNAs improved the authenticity of the findings agreeing with the reports that proved that miRNA-mediated regulation is highly dependent on the uracil present at the first position of the mature miRNA [26].

Consistent with previous studies, it was also observed that the majority of the predicted targets of guava miRNAs were transcription factors that are mostly involved in plant growth, developmental patterning, or cell differentiation. For example, transcription factors SBP/SPL have a crucial role in plant growth, vegetative phase transition, and root development and those are the main targets of miRNA family 156 [50]; while the miR160 family targets ARFs and plays a vital role in root development and auxin signaling pathways [51]. Likewise, the miR828 family has been found to target MYB-like transcripts which significantly participate in stress response signaling and plant development [52]. Moreover, transcription factors NAM and NAC which were involved in fruit ripening and shoot development are mostly targeted by the miR164 family, and the miR171 family targets the GRAS transcription factors which participate in nodule morphogenesis and floral development [53,54,55]. Therefore, the current study re-established the fact that the majority of miRNAs act on transcription factors controlling plant development and organ formations.

Numerous genomic and proteomic studies have elucidated that plants’ response to saline and other stresses comprises a broad spectrum of processes, such as protein biosynthesis, membrane trafficking, and signal transduction [56,57], and it has been well established that miRNAs and their targets influence directly on plant stress tolerance [58,59]. In this context, in Arabidopsis, maize, and cowpea an upregulation of miR156, miR159, miR160, miR162, miR168, miR169, and downregulation of their corresponding targets such as SBPs/SPLs, TCP family transcription factor, ARFs, RNaseIII CAF protein, AGO1, and CBF during salinity stress has been documented [1,60,61]. In addition, during salinity stress, three members of the miR169 family-miR169g, miR169n, and miR169o, as well as miR393, have been upregulated in rice, which precisely cleaves the NF-YA transcription factor gene transcript [62,63]. Similarly, a microarray experiment on two cotton cultivars (salt resistance SN-011 and salt-sensitive LM-6) revealed that miR156, miR169, miR535, and miR827 were substantially upregulated in LM-6, while miR167, miR397, and miR399 were downregulated [64]. Nevertheless, the current research on the crucial role of miRNAs in salt stress responses is largely based on the expression profiling in plant species with varying salt sensitivities under variable salt levels [65]. For example, recently a report demonstrated the downregulation of miR164 under salinity stress (which is consistent with our results) in safflower significantly increased NAC expression [66]. Moreover, it has been shown that miR164 is a negative regulator of lateral root development (auxin-mediated) by controlling NAC1 levels in Arabidopsis thaliana [67]. Thus, it is possible that the downregulation of miR164 in guava under salinity stress, may contribute to salt stress adaption and with root/shoot formation as previously reported in sweet potato [68,69]. In the same way, miR166 is downregulated in guava in response to salinity stress which agrees with the previous results reported in maize and chickpea [61,70], while an upregulation was evidenced in both sweet potato leaves and roots [69]. Interestingly, miR166 family members are proven to regulate auxiliary meristem initiation and leaf morphology by controlling the expression of the HD-ZIP III protein [70,71,72,73]. It has been demonstrated that the downregulation of miR166 leads to salt tolerance in safflower, suggesting the same response in guava [66]. Similarly, it has been widely reported that miR156 has a significant role during salt stress conditions in many plants [70,74,75]. In this study, one of the potential targets of guava miR156 was the squamosa promoter binding protein-likes (SPLs) which is a key regulator of plant abiotic stress tolerance [69]. Similar expression pattern with safflower and Japanese white birch, downregulation of this miRNA in guava may enhance the inhibition of translation or cleavage SPL mRNAs which participate in controlling trichome patterning on the inflorescence stems and floral organs since SPL transcription factors suppress trichome formation by activating TCL1 and TRY gene expression [66,76,77]. Furthermore, the downregulation of both miR160 and miR167 during salt stress in guava is consistent with results published for maize, B. napus, beet, and tomato but contradict with the outcome shown in green cotton, foxtail, and Chinese tamarisk [61,64,78,79,80,81,82]. Interestingly, these miRNAs have been demonstrated to target ARFs that play significant roles in the lateral and adventitious root formation but at high salt concentrations, a diminution of the root growth rate has been evidenced [17]. It has also been reported that miR162 is involved in miRNA processing by targeting the DCL1 gene [83] and corroborating with our results downregulation of miR162 has been reported in maize, cotton, and radish under salinity stress [61,84,85,86]. Moreover, in many plant species, salt responsive differential expression of miR390 was noticed that might target different gene families [86]. For example, it targets the noncoding TAS3 precursor RNA to trigger the biogenesis by cleaving ARFs transcripts leading to the regulation of lateral root growth [87]. In addition, miR390 targets different protein kinases such as AtBAM3 kinase-like receptor that regulates both floral and shoot meristem formation [88]. In this study, the downregulation of guava miR390 under salinity stress showed consistency with the previous reports in chickpea and rice while contradicted the report from Caragana intermedia under salinity stress [70,89,90]. Nonetheless, all the reports indicated that a systematic study of miRNA response to salt stress in closely related genotypes with conflicting stress sensitivities would provide deeper insights into miRNA-guided gene regulation.

It is widely documented that plants face dynamic environmental challenges during their lives that have a substantial effect on their resilience, enlargement, maturity, and yield [91,92,93]. Since miRNAs are involved in plant responses to various environmental stimuli, recent research has identified miRNAs as powerful targets for improving plant stress tolerance [94,95], and hence, some important plant miRNAs have already been engineered to achieve better abiotic stress tolerance. For example, transgenic creeping bentgrass overexpressing any of the following miRNAs osa-miR319a, osa-miR528, or osa-miR393a exhibited enhanced salt tolerance [96,97,98]. More specifically, miR393 was reported to influence abiotic stress responses through directly repressing the expression of its targets AsTIR1 and AsAFB2 (auxin receptors), which leads to the enhanced tolerance to salinity, heat, and drought stress [98]. On the contrary, a work with transgenic rice and Arabidopsis thaliana that overexpressed osa-miR393 demonstrated that these transgenic lines are more sensitive to salt as compared to wild-type plants [63]. More recently, it has been established that in soybean the overexpression and knockdown of miR172c activity resulted in substantially increased and reduced root sensitivity to salt stress, respectively [99]. Nevertheless, all these case studies confirm miRNAs as crucial players and fine tuners in plant responses to stresses. Therefore, it is important to explore the possible role of miRNAs in nonmodel plant species since we assume that a better understanding of these molecules will provide an excellent platform to comprehend numerous physiological and developmental processes in plants that can help to develop authentic stress-tolerant transgenic lines.

5. Conclusions

The post-transcriptional function of miRNAs has a great influence on the overall gene regulatory network in plants. Due to the availability of advanced software tools and sequence resources in public databases, there has been a growing interest in computer-based miRNA identification in the last few years. Nevertheless, this is the first report of guava microRNAs and their targets. In this study using homology-based analysis and strict filtering criteria, we have identified 40 potential guava microRNAs belonging to 19 families as well as 49 corresponding targets and investigated the influence of salinity stress on selected miRNAs which may contribute to further understanding the miRNAs function and regulatory mechanism in guava. Moreover, recently, artificial microRNA mediated gene silencing technology has been utilized successfully for crop improvement. Hence, we believe this study will not only provide considerable aid in miRNA research on economically important fruit crops but also help to initiate an artificial miRNA-based salinity stress tolerance study on guava.

Author Contributions

Data analysis, writing, editing A.S. (Ashutosh Sharma); Experimental procedures, data collection and analysis L.M.R.-M., F.I.S.-C., P.R.R.-P., C.K.T.A., Y.E.B.A., A.K.H.A.; Bioinformatics, and statistical analysis A.S. (Aashish Srivastava); Conceptualization, experimental design, critical writing, reviewing S.P. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Table

Figure 1. Diagrammatic representation of guava miRNA search procedure (workflow).

Figure 2. Secondary stem-loop structures of the potential guava miRNA precursors/pre-miRNAs. Respective miRNAs are represented with red font.

Figure 3. Conserved and nonconserved potential guava miRNA families (dark green boxes) and their homologs in other plant organisms. Identical color shades reflect closely related species, while filled and empty boxes indicate the presence or absence of miRNA families, respectively.

Figure 4. Web logo displaying conserved nucleotide sequences of (A) pre-miRNA160c (B) pre-miRNA390b and (C) pre-miRNA396b sequences.

Figure 5. Phylogenetic analysis of the identified guava miRNAs (marked with a red box) pgu-miR160c-5p, pgu-miR390b-5p, and pgu-miR396b-5p was performed using their potential precursor sequences.

Figure 6. Comparative synteny map of potential guava miRNAs against the well-annotated genome of phylogenetically close species apple (Malus domestica).

Figure 7. Gene ontology analysis of potential targets in Guava. (A) Molecular function, (B) Biological process. (C) Cellular component.

Figure 8. KEGG pathways mapped using the KAAS tool among the predicted targets of the identified miRNAs in guava.

Figure 9. Minimum free energy (MFE) based network interaction of potential guava miRNAs and their corresponding targets. Different guava miRNAs and their targets are marked with orange and blue circles, respectively. Targets marked with the green circle are shared by two miRNAs, while targets marked with the pink circle are shared by 3 or more miRNAs.

Figure 10. Evaluation of relative fold change of the selected guava miRNAs under salinity stress. The delta-delta CT method was used to determine the fold change and the U6 RNA was used as a normalization control. The values were further normalized with respect to the control condition that was set to 1.

Table 1

Potential miRNAs in guava.

Identified miRNAs	LM * (nt)	Query miRNAs	miRNA Sequences	Accession	Strand	Location	LP * (nt)	MFEs (ΔG)	MFEI
pgu-miR156ac	21	mdm-miR156ac	UUGACAGAAGAUAGAGAGCAC	MI0022999	+/−	5′	86	−47	1.18
pgu-miR156f-5p	20	ath-miR156f-5p	UGACAGAAGAGAGUGAGCAC	MI0000183	+/+	5′	86	−50.8	1.1
pgu-miR159-5p	21	eun-miR159-5p	AGCUGCUGGUCUAUGGAUCCC	MI0033013	+/−	5′	134	−48.5	0.77
pgu-miR159d	21	mdm-miR159d	UUUGGAUUGAAGGGAGCUCUA	MI0035594	+/−	3′	215	−95.2	0.93
pgu-miR160a-3p	21	ath-miR160a-3p	GCGUAUGAGGAGCCAUGCAUA	MI0000190	+/+	3′	91	−49.7	0.97
pgu-miR160c-5p	21	ath-miR160c-5p	UGCCUGGCUCCCUGUAUGCCA	MI0000192	+/+	5′	86	−51.4	1.01
pgu-miR162-3p	21	eun-miR162-3p	UCGAUAAACCUCUGCAUCCAG	MI0033015	+/−	3′	111	−38.4	0.7
pgu-miR164b-5p	21	ath-miR164b-5p	UGGAGAAGCAGGGCACGUGCA	MI0000198	+/+	5′	91	−43.4	0.9
pgu-miR166-3p	21	eun-miR166-3p	UCGGACCAGGCUUCAUUCCCC	MI0033016	+/−	3′	123	−51.3	0.99
pgu-miR166t	20	gma-miR166t	UCGGACCAGGCUUCAUUCCC	MI0021696	+/+	5′	111	−47.7	0.75
pgu-miR167a-5p	21	ath-miR167a-5p	UGAAGCUGCCAGCAUGAUCUA	MI0000208	+/+	5′	111	−50.4	1.01
pgu-miR167c-5p	22	eun-miR167c-5p	UGAAGCUGCCAGCGUGAUCUCA	MI0033019	+/−	5′	86	−33.1	0.74
pgu-miR169b-5p	21	ath-miR169b-5p	CAGCCAAGGAUGACUUGCCGG	MI0000976	+/+	5′	116	−37.6	0.7
pgu-miR169f-5p	21	ath-miR169f-5p	UGAGCCAAGGAUGACUUGCCG	MI0000980	+/+	5′	96	−41.5	0.86
pgu-miR169k	21	mdm-miR169k	UAGCCAAGGAUGACUUGCCUG	MI0035667	+/+	5′	103	−45.5	1.08
pgu-miR169w	21	gma-miR169w	CAAGGAUGACUUGCCGGCAUU	MI0033468	+/+	5′	106	−38.2	0.8
pgu-miR171b	21	mdm-miR171b	UUGAGCCGCGUCAAUAUCUCC	MI0023043	+/+	3′	116	−46.1	0.85
pgu-miR171g	21	mdm-miR171g	UGAUUGAGCCGUGCCAAUAUC	MI0023048	+/+	3′	88	−36.7	1.05
pgu-miR171j	21	mdm-miR171j	UUGAGCCGCGCCAAUAUCACU	MI0023051	+/+	3′	106	−41.8	0.91
pgu-miR172a	21	ath-miR172a	AGAAUCUUGAUGAUGCUGCAU	MI0000215	+/+	3′	111	−48.4	0.99
pgu-miR172b-5p	21	eun-miR172b-5p	GCAGCAUCAUCAAGAUUCACA	MI0033022	+/−	5′	111	−48.4	0.99
pgu-miR172l	21	mdm-miR172l	GGAAUCUUGAUGAUGCUGCAG	MI0023067	+/−	3′	176	−65.1	0.97
pgu-miR319a	21	ath-miR319a	UUGGACUGAAGGGAGCUCCCU	MI0000544	+/+	3′	201	−92.4	0.9
pgu-miR319e	20	mdm-miR319e	GAGCUUUCUUCAGUCCACUC	MI0035621	+/+	5′	176	−85.7	0.88
pgu-miR390b-5p	21	ath-miR390b-5p	AAGCUCAGGAGGGAUAGCGCC	MI0001001	+/+	5′	84	−51.3	1.09
pgu-miR393a-5p	22	ath-miR393a-5p	UCCAAAGGGAUCGCAUUGAUCC	MI0001003	+/+	5′	96	−40.7	0.93
pgu-miR393c	22	mdm-miR393c	UCCAAAGGGAUCGCAUUGAUCU	MI0023081	+/+	5′	111	−44.4	0.91
pgu-miR393i	22	gma-miR393i	UUCCAAAGGGAUCGCAUUGAUC	MI0021710	+/+	5′	126	−56.4	0.97
pgu-miR395-3p	21	eun-miR395-3p	AUGAAGUGUUUGGGGGAACUC	MI0033024	+/+	3′	72	−47.5	1.36
pgu-miR395a	21	ath-miR395a	CUGAAGUGUUUGGGGGAACUC	MI0001007	+/+	3′	106	−37.2	0.89
pgu-miR396b-3p	21	ath-miR396b-3p	GCUCAAGAAAGCUGUGGGAAA	MI0001014	+/+	3′	126	−36.1	0.82
pgu-miR396b-5p	21	eun-miR396b-5p	UUCCACAGCUUUCUUGAACUG	MI0033026	+/+	5′	118	−52.4	1.01
pgu-miR396d	21	mdm-miR396d	UUCCACAGCUUUCUUGAACUU	MI0023096	+/+	5′	96	−50.6	0.9
pgu-miR399c	21	mdm-miR399c	UGCCAAAGGAGAAUUGCCCUG	MI0023107	+/−	3′	96	−50.4	1.12
pgu-miR399f	21	ath-miR399f	UGCCAAAGGAGAUUUGCCCGG	MI0001025	+/+	3′	111	−50	0.98
pgu-miR399g	21	mdm-miR399g	UGCCAAAGGAGAUUUGCUCGG	MI0023111	+/+	3′	81	−36.5	1.11
pgu-miR530-5p	21	eun-miR530-5p	UCUGCAUUUGCACCUGCACCU	MI0033032	+/+	5′	161	−73.6	0.89
pgu-miR535b-5p	21	eun-miR535b-5p	UGACAACGAGAGAGAGCACGC	MI0033034	+/+	5′	81	−41.6	1.01
pgu-miR535d	21	mdm-miR535d	UGACGACGAGAGAGAGCACGC	MI0023131	+/+	5′	101	−49	1.02
pgu-miR828b	22	mdm-miR828b	UCUUGCUCAAAUGAGUAUUCCA	MI0023134	+/+	5′	99	−33.8	0.89

* LM length of mature miRNAs, * LP length of precursors.