Genome-wide analysis of conserved and novel miRNAs in mango mesocarp reveals early regulatory networks involved in postharvest heat stress response
Mitzuko Dautt-Castro, Abraham Cruz-Mendívil, Libertad Ulloa-Álvarez, Jorge A. Enríquez-Domínguez, Miranda Ávila-Sánchez, Lourdes K. Ulloa-Llanes, Sergio Casas-Flores, Tiago D. Serafim, María A. Islas-Osuna
TL;DR
This study identifies miRNAs in mango fruit that respond to heat stress, revealing regulatory networks that could help improve fruit thermotolerance and quality.
Contribution
The first genome-wide analysis of miRNAs in mango mesocarp under postharvest heat stress, identifying novel regulatory modules and their roles in thermotolerance.
Findings
90 miRNAs were identified across 27 families, with most targeting transcription factors and hormone signaling regulators.
miR168, miR319, and miR482 were identified as early heat-responsive miRNAs involved in ROS homeostasis and thermotolerance.
A putative interaction between miR482 and a lncRNA suggests coordination of phasiRNA generation during heat stress.
Abstract
MicroRNAs (miRNAs) are key post-transcriptional regulators involved in plant development and stress responses. Although extensively studied in model species, their functional roles remain largely unexplored in non-model fruits, such as mango (Mangifera indica L.). In this study, we present, to our knowledge, the first genome-wide analysis of miRNAs in mango mesocarp under postharvest heat stress induced by hot water treatment (HWT), a widely used quarantine method. Using small RNA sequencing (sRNA-seq), we identified 90 miRNAs distributed across 27 families. Phylogenetic analysis revealed evolutionary trajectories shaped by whole-genome duplication. Most miRNAs were predicted to target transcription factors and regulators involved in growth and hormone signaling. Differential expression profiling across multiple time points post-HWT identified miR168, miR319, and miR482 as early…
Genes, proteins, chemicals, diseases, species, mutations and cell lines named across the full text — each resolved to its canonical identifier and authoritative record.
Click any figure to enlarge with its caption.
Figure 1
Figure 2
Figure 3
Figure 4
Figure 5
Figure 6
Figure 7
Figure 8- —SECIHTI
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
TopicsPlant Molecular Biology Research · Plant Gene Expression Analysis · Postharvest Quality and Shelf Life Management
Introduction
Small RNAs (sRNAs) are critical regulators of gene expression in eukaryotic organisms, typically ranging from 20 to 24 nucleotides (nt) in length. The main classes of sRNAs include microRNAs (miRNAs) and small interfering RNAs (siRNAs), which mediate transcriptional gene silencing (TGS) via RNA-directed DNA methylation (RdDM) or post-transcriptional gene silencing (PTGS), through mRNA degradation or translational inhibition^1–3^.
Among these, miRNAs form a highly conserved and functionally specific class of sRNAs that regulate gene expression with remarkable precision. They are transcribed by RNA polymerase II into primary transcripts (pri-miRNAs), which are processed in the nucleus by the DCL1-HYL1-SE protein complex into precursor miRNAs (pre-miRNAs), and subsequently into mature miRNA duplexes. One strand of the duplex is incorporated into the RNA-induced silencing complex (RISC), which guides either cleavage or translational repression of target mRNAs based on sequence complementarity^4^. Through this pathway, miRNAs modulate a wide range of biological processes, including organogenesis, hormone signaling, and response to biotic and abiotic stresses.
Several miRNAs have been functionally linked to heat-stress responses across various plant species. For instance, miR169, miR167b, and miR398 have been shown to modulate thermotolerance in Arabidopsis, tomato, and grapevine through regulation of heat shock transcription factors (HSFs) and NF-YA genes^5^, modulation of auxin signaling^6^, or regulation of reactive oxygen species (ROS) homeostasis^7^.
High-throughput small RNA sequencing (sRNA-seq) has enabled genome-wide identification and quantification of small RNAs across several plant species, revealing miRNA-mediated regulatory networks involved in heat-stress responses. In tomato, sRNA-seq and degradome analysis identified several heat-responsive miRNAs, including sly-miR482d-3p, sly-miR398b, sly-miR167h, sly-miR164b-3p, and sly-miR172d, which are predicted to regulate genes associated with auxin signaling, oxidative stress responses, and transcriptional regulation^8,9^.
In wheat, sRNA-seq and degradome profiling identified 84 miRNAs targeting over 200 transcripts, many of which are linked to redox homeostasis, mitochondrial function, and heat stress signaling, underscoring their crucial role in maintaining thermotolerance, particularly in reproductive tissues^10^.
Despite advances in model species and some fruit crops, non-model fruits, such as mango (Mangifera indica L.), remain largely underexplored for sRNA-mediated regulatory mechanisms. Mango is a climacteric fruit of significant economic importance, particularly in countries like Mexico and India. A common postharvest practice in mango export is hot water treatment (HWT), a phytosanitary procedure that involves immersing fruits in warm water to eliminate quarantine pests^11^. However, this treatment can induce heat stress, altering physiological, biochemical, and molecular processes in the fruit and thereby accelerating ripening, potentially reducing shelf life and overall consumer quality^12^.
Previous genome-wide studies conducted by our research group have identified components of the miRNA biogenesis pathway in mango and suggested their involvement in the heat stress response^13^. However, the specific miRNA repertoire expressed under HWT, along with their differential expression patterns and regulatory functions, remains largely uncharacterized. Addressing this knowledge gap is particularly relevant given the potential of miRNAs to regulate key processes associated with postharvest fruit quality. Therefore, this study aimed to comprehensively characterize the miRNA profile of mango subjected to HWT using sRNA-seq, to identify heat-responsive miRNAs, and to explore their potential regulatory roles in the fruit’s molecular response to heat stress.
Results
A total of ninety miRNAs were identified in the mango mesocarp
Twelve small-RNA libraries were sequenced, yielding an average of 12.2 million raw reads and 9 million clean reads per library, with 78.2% of the clean reads successfully mapped to the mango genome (Table 1). The samples used to construct the libraries comprised three biological replicates at each post-treatment time point (1, 3, 6, and 24 h post-treatment). However, principal component analysis (PCA) of the expression profiles revealed that one outlier replicate was detected in each treatment group and was removed for downstream studies (Table 1, marked with an asterisk), leaving two biological replicates per treatment for subsequent analyses (Fig. S1). Length distribution analysis showed that more than 90% of the 23,625 predicted small RNAs were from 24 nt (Fig. 1a). According to ShortStack parameters, ninety small RNAs met the criteria to be classified as miRNAs after mapping the trimmed reads to the mango reference genome (cv. Alphonso) and aligning the mature homologous miRNA sequences from miRBase. Among these, 77 were annotated as known miRNAs and 13 as novel candidates (Table 1). Of the total, 73 were 21 nt and 17 were 22 nt in length (Fig. 1b). Further analysis of precursor loci indicated that most miRNAs (54.44%) originated from intergenic regions, while 26.66% derived from exonic sequences, 6.66% from intronic regions, and 12.22% from mixed loci spanning multiple genomic contexts (Fig. 1c). Among miRNAs mapped to intronic or exonic regions, thirty-nine were associated with non-coding RNA loci. In contrast, three appeared to originate from protein-coding genes, including those encoding a G-type lectin S-receptor-like serine/threonine-protein kinase, a DNA replication complex GINS protein, and an uncharacterized RNA (Table S1). In terms of chromosomal distribution, miRNAs were mapped to 19 of the 20 mango chromosomes, except for chromosome 11 (Table S1). Furthermore, over 73% of mature miRNAs had a 5’ terminal Uracil (U), while 12.22%, 7.77%, and 6.66% began with Guanine (G), Cytosine (C), and Adenine (A), respectively (Fig. 1d). Interestingly, only 63 of the 90 identified miRNAs were unique mature sequences. The remaining 27 miRNAs shared 100% sequence identity at the mature miRNA level with one of these unique sequences but differed in their precursor and genomic origins (Table S1).
Table 1. Summary of small RNA-sequencing (sRNA-seq) libraries from Mangifera indica L. mesocarp under heat stress induced by hot water treatment (HWT).LibrariesRaw readsClean readsReads mapped to genomeKnown miRNAs identifiedNovel miRNAs identified1 h – R112,812,6419,135,8747,799,898 (85.37%)77131 h – R213,461,90411,797,3059,172,145 (77.74%)1 h – R313,320,03410,315,7328,193,234 (79.42%)3 h – R113,365,74210,124,5677,883,630 (77.86%)3 h – R211,247,1359,695,8127,676,611 (79.17%)3 h – R310,413,5716,463,0894,713,631 (72.93%)6 h – R113,385,24210,699,7148,514,972 (79.58%)6 h – R210,887,1467,538,8586,004,097 (79.64%)6 h – R311,038,8687,394,6855,958,123 (80.57%)24 h – R115,027,88310,192,6107,839,770 (76.91%)24 h – R210,589,5548,656,7126,727,525 (77.71%)24 h – R310,660,8016,820,4904,843,856 (71.01%)h: hours post-treatment; R: biological replicate.
Fig. 1. Characterization of miRNAs identified in Mangifera indica L. (a) Size distribution of small RNAs obtained from sRNA-seq, classified by nucleotide length; (b) Length distribution of known and novel miRNAs; (c) Genomic origin of miRNAs, categorized by the type of genomic region (intergenic, exonic, intronic, or mixed); (d) Nucleotide Composition at the 5’ end of mature miRNAs.
Mango mesocarp miRNAs are distributed across twenty-seven conserved families
According to miRNA annotations from sRNA-seq and supported by phylogenetic analyses, twenty-seven miRNA families were identified in mango mesocarp. The miR166 family was the most represented, with eleven members, followed by miR396 with nine members, and miR319 with seven (Fig. 2). To further investigate the evolutionary history of these miRNA families, a phylogenetic analysis was conducted using their pre-miRNA sequences. The results showed that most miRNA families are highly conserved across the analyzed species (Mangifera indica L., Arabidopsis thaliana,* Citrus sinensis*,* Fragaria vesca*, and Solanum lycopersicum), suggesting that they have functionally conserved roles within each family. The miR156 family appeared as the earliest diverging clade, possibly representing one of the oldest miRNA lineages in plant evolution. In contrast, the miR399 family exhibited the most recent divergence (Fig. 3). Notably, some members of the miR482 family displayed divergence across distant evolutionary timepoints. Interestingly, several miRNAs did not cluster with their canonical families, such as MimiR164 (cluster 3502), MimiR530 (cluster 3168), MimiR156 (cluster 17053), and MimiR166 (cluster 3389), suggesting that functional specialization have occurred in these miRNAs (Fig. 3). On the other hand, the only member of the miR828 family appeared to be poorly conserved across species, as no homologs were detected among the analyzed taxa. Furthermore, this phylogenetic analysis enabled the annotation of eleven out of the thirteen novel miRNAs identified in the sRNA-seq analysis, as they clustered within specific known families (Fig. 3). This included the classification of four as miR482, three as miR395, and one each member of the miR396, miR530, miR169, and miR390 families. For the two novel miRNA candidates that did not clearly cluster with known families (Fig. 3: Cluster_4783 and Cluster_12685), additional evidence supported their authenticity, including canonical stem–loop precursor structures, strand-specific and dominant Dicer-sized read accumulation across the precursor loci, and the absence of homology to annotated non-miRNA ncRNAs.
Fig. 2. The mango mesocarp contains twenty-seven conserved miRNA families. Bar graph showing the number of miRNA members identified per conserved family in mango mesocarp based on sRNA-seq data. Families such as miR166, miR396, and miR319 displayed the highest number of members, consistent with their broad conservation and functional diversification in other plant species.
Fig. 3. Phylogeny tree constructed using pre-miRNAs from M. indica and homologous sequences from other plant species. Each color represents a distinct miRNA family. Mango pre-miRNAs are shown in bold, and novel miRNAs identified through sRNA-seq are highlighted in blue. Asterisks (✻) indicate novel miRNAs that did not clearly cluster with any previously annotated families. The tree was generated using the maximum-likelihood method and the best-fit model inferred by IQ-TREE 2. Species abbreviations: Mi, Mangifera indica; At, Arabidopsis thaliana; Cs, Citrus sinensis; Fv, Fragaria vesca; Sl, Solanum lycopersicum.
The target genes of mango miRNAs are mostly transcription factors
To investigate the biological processes in which mango mesocarp miRNAs may be involved, we first predicted their target genes using the psRNAtarget server^14^. A total of 395 target transcripts were identified, of which 370 were predicted to be regulated via mRNA cleavage, and 24 via translation inhibition. No suitable targets were found for ten miRNAs under the defined prediction criteria (Table S1). Using all identified target sequences, we performed a Gene Ontology (GO) enrichment analysis. Most targets were involved in DNA-templated transcription regulation, followed by cell differentiation and the auxin-activated signaling pathway. The most enriched molecular function categories were DNA binding, DNA-binding transcription factor activity, and ATP binding, while the most prominent cellular component was the nucleus (Fig. 4). These findings suggest that the majority of mango miRNA targets are transcription factors or proteins associated with transcription, underscoring their potential regulatory role in gene expression. Interestingly, none of the enriched GO terms were directly linked to heat stress responses, such as “cellular response to heat” (GO:0034605) or “response to temperature stimulus” (GO:0009266). However, a comprehensive literature review revealed that many of the identified targets are functionally associated with heat stress in other plant species.
Fig. 4. Gene ontology (GO) classification of predicted target genes of mango miRNAs. (a) Biological processes (GO terms with ≥ 5 associated genes); (b) Molecular functions (GO terms with ≥ 3 associated genes); (c) Cellular components (GO terms with ≥ 2 associated genes). GO annotation was performed using the Blast2GO module in OmicsBox, based on functional enrichment of the predicted miRNA targets.
Some mango miRNAs respond to HWT and may regulate their predicted target genes to mitigate heat stress
To identify miRNAs responsive to heat stress at different time points (3, 6, and 24 h post-treatment, hpt) compared to the 1-hour time point, a differential expression analysis was performed using the read counts from ShortStack. At 3 hpt, 109 sRNAs were differentially expressed (71 up-regulated and 38 down-regulated) (Fig. S2a, b; Table S2); however, only one corresponded to a known miRNA, miR169. At 6 hpt, the highest number of differentially expressed sRNAs was observed, totaling 239 (183 up-regulated and 56 down-regulated) (Fig. S2c, d; Table S2), among which four miRNAs were identified: miR168, miR319, miR169, and miR319b. This transcriptional change at 6 hpt is consistent with the PCA results, in which samples from this time point showed the most significant separation from the other groups, reflecting a marked shift in the global expression profile (Fig. S1). At 24 hpt, 179 sRNAs were differentially expressed (128 up-regulated and 51 down-regulated), including three miRNAs: miR482, miR169, and miR319b (Fig. S2e, f; Table S2). The predicted precursors hairpins of all these differentially expressed miRNAs exhibited canonical stem-loop secondary structures, with minimum free energy values ranging from − 61.5 to -83.6 kcal/mol (Table 2; Fig. S3). From this set, three miRNAs were selected for expression validation by stem-loop RT-qPCR, based on the predicted duplex stability (ΔG values) of their duplex with target transcripts (Table 2). The selected miRNAs were: miR168 (up-regulated at 6 hpt) predicted to target MiAGO1 and MiPUM2; miR319 (up-regulated at 6 hpt) targeting MiGAMYB and MiTCP4; and miR482 (down-regulated at 24 hpt) targeting a long non-coding RNA (Milnc-1) (Table 2). miR169, although differentially expressed at all three time points, lacked a high-confidence predicted target. For miR319, which was down-regulated at 6 and 24 hpt, the predicted interaction with its target MiMYB123 was less stable than those of the validated miRNAs (Table 2).
Table 2. Heat-responsive miRNAs identified in mango mesocarp and their predicted target genes.miRNASequenceLocusChrGenomic locationDifferential expression (time/fold-change)Putative Target IDPutative targetsInhibitionΔG duplex formation (kcal/mol)ΔG precursor hairpin (kcal/molmiR168UCGCUUGGUGCAGGUCGGGAACNC_058144.1:10593810–10,594,0128Intergenic6 h (1.9)XM_044621157.1 MiAGO1 Cleavage-37.41-75.5XM_044649578.1 MiPUM2 Cleavage-31.67miR319UUGGACUGAAGGGAGCUCCCUNC_058139.1:5004416–5,004,8363Intergenic6 h (3.68)XM_044606613.1 MiGAMYB Cleavage-33.43-83.9NC_058140.1:1308917–1,309,3374MixedXM_044656421.1 MiMYB33 Cleavage-33.43NC_058156.1:830778–830,98220IntronXM_044609727.1 MiTCP4 Cleavage-37.7NC_058142.1:10476063–10,476,2756IntergenicXM_044605661.1 MiMNS3 Traslation-30.21NC_058139.1:1054542–1,054,7493IntergenicmiR482UCUUGCCAACUCCGCCCAUACCNC_058141.1:811685–811,8305Exon24 h (-29.15)c24640_g2_i10 Milnc-1 Cleavage-40.13-111c11593_g1_i2Uncharacterized proteinCleavage-43.92NC_058141.1:843322–843,4665IntergenicXM_044641238.1 MiRP Cleavage-37.35XM_044634119.1 MiRPP-13 Cleavage-61XM_044631012.1 MiNIP6-1 Cleavage-33.91miR169UGUAGGCUUAGAAAGUAUCGUGNC_058138.1:22168361–22,168,5552-3 h (2.35)- - ---82.86 h (3.87)24 h (4.43)miR319bUGGACUGAAGGGAGCUCCUUCNC_058146.1:5739952–5,740,16810Intron6 h (-8.13)XM_044616538.1 MiMYB123 Cleavage-26.79-83.612Intergenic24 h (-5.14)AGO1: protein argonaute 1; APUM2: pumilio homolog 2; GAMYB: similar to transcription factor GAM1; MYB33 and MYB123: similar to transcription factor MYB33 or MYB123; TCP4: similar to transcription factor TCP4; MNS3: mannosyl-oligosaccharide 1,2-alpha-mannosidase MNS3; lnc-1, long non- coding-1; RP, disease resistance protein; RPP13: putative disease resistance RPP13-like protein 1; NIP6-1: aquaporin. Bold letters: miRNAs and targets selected to validate.
miR168 was mapped to exon 1 of MiAGO1 and exon 2 of MiPUM2 (Fig. 5a, c), and both interactions showed highly thermodynamic stability (Fig. 5b, d). RT-qPCR revealed a significant up-regulation of miR168 at all evaluated time points (Fig. 5e), while MiAGO1 showed an inverse expression pattern, particularly evident at 1 hpt and maintained throughout the time course (Fig. 5e, f), supporting a regulatory interaction. In contrast, MiPUM2 was up-regulated at 3 and 6 hpt (Fig. 5g) and did not exhibit a consistent inverse trend with miR168.
Fig. 5miR168 potentially regulates MiAGO1 expression in mango mesocarp under hot water treatment (HWT). (a, c) Gene structures of the predicted targets MiAGO1 and MiPUM2, and their sequences alignment with miR168. (b, d) Predicted miRNA-target interactions using RNAup web server. (e) Relative expression of miR168 compared with the normalized sRNA-seq read counts. (f, g) Relative expression levels of MiAGO1 and MiPUM2, respectively. Expression levels of miRNA and mRNAs were quantified by Stem loop RT-qPCR and RT-qPCR, respectively, at 1, 3, 6, and 24 h post-treatment (hpt) and normalized against U6 (for miRNA) and Actin 7 (for mRNAs) using the 2^−Ct^ method. Untreated fruits were used as controls. Bars represent standard error. Different letters indicate statistically significant differences as determined by Tukey-Kramer test (p < 0.05).
miR319 was predicted to bind to exon 4 of MiTCP4 and exon 2 of MiGAMYB (Fig. 6a, c), forming thermodynamically stable duplexes (Fig. 6b, d). Among the nine pre-miR319 sequences identified, seven shared an identical mature miR319 sequence. As a result, the oligonucleotide used for stem-loop RT-qPCR could hybridize to all seven variants (Fig. 6e, red box), indicating that the results reflect cumulative expression. Stem loop RT-qPCR showed consistently elevated expression of miR319 in treated samples compared to controls (Fig. 6f). Interestingly, the expression profiles of its predicted targets MiTCP4 and MiGAMYB followed a similar pattern -downregulated at 1 and 6 hpt, upregulated at 3 hpt, and returning to baseline by 24 hpt (Fig. 6g, h), suggesting a potential regulatory activity at early time points.
Fig. 6miR319 potentially regulates the expression of MiTCP4 and MiGAMYB in mango mesocarp under hot water treatment (HWT). (a, c) Gene structures of the predicted targets MiTCP4 and MiGAMYB, and their sequence alignments with miR319. (b, d) miRNA-target interaction analysis performed using the RNAup web server. (e) Alignment of the seven pre-miR319 identified in mango by sRNA-seq, visualized in Jalview, v2.11.4.1. The red box indicates the mature miR319 variant evaluated in this study; green boxes denote other mature miR319 sequences derived from the same precursors. (f) Relative expression of miR319 compared with the normalized sRNA-seq read counts. (g, h) Relative expression levels of MiTCP4 and MiGAMYB, respectively. Expression levels of miRNA and mRNAs were quantified by Stem loop RT-qPCR and RT-qPCR, respectively, at 1, 3, 6, and 24 h post-treatment (hpt) and normalized against U6 (for miRNA) and Actin 7 (for mRNAs) using the 2^−Ct^ method. Untreated fruits were used as controls. Bars represent standard error. Different letters indicate statistically significant differences as determined by Tukey-Kramer test (p < 0.05).
Furthermore, within the same pre-miR319 group, three additional isomiRs with sequence variants were identified based on precursor sequence alignments (Fig. 6e, green boxes). However, no corresponding reads were detected for these isomiRs in the sRNA-seq dataset, suggesting tissue and condition-specific expression.
For miR482, the predicted target was Milnc-1 (Fig. 7a, b). The non-coding nature of Milnc-1 was confirmed by a coding potential score below zero (-0.79) and the absence of a signal peptide. Two pre-miR482 sequences shared the same mature miR482 form (Fig. 7c), thus RT-qPCR represents their combined expression. miR482 was up-regulated in response to heat stress (Fig. 7d), while Milnc-1 was down-regulated at 1 hpt, showing an inverse expression trend across all time points (Fig. 7e), supporting a potential regulatory interaction.
Fig. 7miR482 potentially regulates the expression of a long non-coding RNA (Milnc-1) in mango mesocarp under hot water treatment (HWT). (a) Gene structure of the predicted target Milnc-1 and its alignment with miR482. (b) Predicted miRNA-Milnc-1 interaction analyzed using the RNAup web server. (c) Alignment of two pre-miR482 sequences identified in mango by sRNA-seq, visualized in Jalview, v2.11.4.1. The red box indicates the mature miR482 variant evaluated in this study. (d) Relative expression of miR482 compared with normalized sRNA-seq read counts. (e) Relative expression levels of Milnc-1. Expression levels of miRNA and mRNAs were quantified by Stem loop RT-qPCR and RT-qPCR, respectively, at 1, 3, 6, and 24 h post-treatment (hpt) and normalized against U6 (for miRNA) and Actin 7 (for mRNAs) using the 2−Ct method. Untreated fruits were used as controls. Bars represent standard error. Different letters indicate statistically significant differences as determined by Tukey-Kramer test (p < 0.05).
Fig. 8miR319 specifically targets and represses MiTCP4 from Mangifera indica L. in a transient expression system. Relative expression levels of MiTCP4 in Nicotiana benthamiana leaves agroinfiltrated with the following constructs: empty vector + MiTCP4 (EV), miR319 + wild-type MiTCP4 (*MiTCP4-*WT), and miR319 + resistant MiTCP4 version (MiTCP4-mut). Leaf tissue was collected 48 h post-infiltration for RNA extraction and RT-qPCR analysis. Expression values were normalized to N. benthamiana PP2A using the 2^−ΔΔCt^ method. Bars represent mean ± SE from three biological replicates. Representative data from two independent experiments are shown. Different letters indicate statistically significant differences (Tukey-Kramer test, p < 0.05).
Finally, RT-qPCR expression data were compared with normalized read counts from sRNA-seq, showing high concordance for all three miRNAs (Figs. 5e, 6f and 7d), confirming the reliability of the sequence results. Notably, for miR319 and miR482, expression values reflect the average expression across multiple pre-miRNA variants.
miR319 effectively represses Mangifera indica L. MiTCP4 expression in Nicotiana benthamiana
To functionally validate one of the predicted miRNA-target modules, we selected miR319 and its putative target MiTCP4 for heterologous expression assays in N. benthamiana. Three treatments were tested: (i) an empty vector plus MiTCP4 (EV + MiTCP4, control); (ii) miR319 co-infiltrated with wild-type MiTCP4 (miR319 + MiTCP4 WT); and (iii) miR319 co-infiltrated with mutant version of MiTCP4 lacking the miR319-binding site (miR319 + MiTCP4 mut). This regulatory module was chosen due to its conserved role in plant development and fruit ripening across species.
Quantitative RT-PCR analysis at 48 h post-infiltration revealed a significant reduction in MiTCP4 transcript levels in leaves co-expressing miR319 and MiTCP4 WT compared to both the EV control and the miR319 + MiTCP4 mut treatment (Fig. 8). Notably, the mutated MiTCP4 construct exhibited a partial recovery in transcript abundance, confirming resistance to miR319-mediated cleavage (Fig. 8). Together, these results demonstrate that miR319 directly downregulates mango MiTCP4 expression in N. benthamiana, supporting the regulatory relationships inferred by our sRNA-seq analysis.
Discussion
MicroRNAs (miRNAs) are essential post-transcriptional regulators involved in a wide range of biological processes in plants, including reproduction, growth, development, ripening, and stress responses. In model species such as Arabidopsis and tomato, miRNA characterization has substantially advanced our understanding of plant responses to adverse environmental conditions, including heat stress. However, in non-model species such as mango (Mangifera indica L.), our knowledge of miRNA-mediated regulation remains limited, despite its high economic, commercial, and nutritional value. Notably, mangoes exported from Mexico are subjected to quarantine hot-water treatment (HWT), a phytosanitary practice that induces thermal stress.
To better understand the regulatory landscape in which mango miRNAs function under heat stress, small RNAs from mesocarp tissue were sequenced using sRNA-seq. The predominance of 24-nt siRNAs observed in mango is consistent with findings in other species, including Arabidopsis^15^and sweet orange (3). Similarly, the low relative abundance of miRNAs (0.4%), compared to other sRNA types, mirrors findings in wheat^16^. Consistent with previous observations in other plants, most mango miRNA genes are in intergenic regions^17,18^. Additionally, the 5’ nucleotide composition of mango miRNAs suggests preferential loading into AGO1, which aligns with previous evidence indicating that canonical miRNAs are primarily incorporated into AGO1 in plants^3^.
In mango, miRNAs are distributed across twenty-seven families, a number comparable to those reported in Fragaria vesca^19^, Cucumis melo^20^, and Pisum sativum L^21^. Families such as miR166, miR396, miR319, miR171, and miR156 were among the most represented, consistent with findings in other plant species and underscoring their evolutionary conservation and functional significance^17^. Phylogenetic analysis of pre-miRNAs from Mangifera indica L. revealed numerous highly conserved families, with clustering patterns reflecting shared evolutionary relationships with other eudicots, such as A. thaliana, C. sinensis, F. vesca, and S. lycopersicum. These results corroborate previous reports on the ancestral origin of core plant miRNAs^22,23^. In several families, including miR156, miR166, miR396, miR171, and miR164, multiple clusters were identified in mango. While some grouped within the same phylogenetic branch, others appeared in distinct clades, suggesting a complex evolutionary trajectory, potentially shaped by whole-genome duplication (WGD) events reported in mango approximately 65 million years ago (mya)^24^or 33 mya^25^. These events were likely followed by differential retention, subfunctionalization, or structural divergence of miRNA precursors^22,26^. A particularly notable case was miR828, for which no conserved homologs were detected in other species analyzed. This may indicate lineage-specific evolution in mango or different losses in other lineages, supported by its distinct phylogenetic placement. The presence of identical mature miRNA sequences arising from different precursors located in distinct genomic regions suggests selective pressure to conserve mature miRNA function, while allowing precursor diversification that could influence transcriptional regulation, processing efficiency, or spatiotemporal expression^27,28^. Additionally, miRNAs that do not cluster into any known conserved families may represent species-specific regulators uniquely adapted to mango’s developmental programs or stress responses. Consistent with this, the predicted mango pre-miRNA sequences exhibited canonical stem-loop secondary structure and minimum free energy values comparable to those reported for conserved plant miRNA, including those characterized in Arabidopsis^29^.
On the other hand, the lack of GO terms directly associated with heat stress likely reflects the indirect mode of miRNA-mediated regulation, whereby miRNAs fine-tune upstream regulatory and signaling components rather than directly targeting classical heat-responsive genes. This observation highlights a potentially intriguing regulatory feature that merits further investigation using functional genomic approaches.
The relative abundance of miRNAs under heat stress varies across plant species and developmental stages. In this study, we analyzed the expression of miR168, miR319, and miR482, which were differentially expressed in mango mesocarp following HWT, as identified in the sRNA-seq data. The results suggest that these miRNAs may function as early regulators of the heat stress response, potentially contributing to the modulation of processes such as miRNA homeostasis (miR168/AGO1), thermotolerance (miR319/GAMYB), ROS balance (miR319/TCP4), and phasiRNA biogenesis (miR482/lnc-1), as depicted in the integrative model (Fig. S4).
The role of miR168 and AGO1 in heat stress regulation has been previously reported in Solanum pimpinellifolium and Malus hupehensis, where miR168 levels decrease, and AGO1 expression increases in response to heat^30,31^. In mango, MiAGO1 was upregulated at 3 hpt, possibly reflecting increased miRNA biogenesis and a conserved activation of gene silencing pathways under heat stress, as described in tomato, papaya, guava, and grapevine. In this mechanism, miR168 restricts AGO1 accumulation, thereby preventing excessive silencing activity^32,33^, as evidenced by the marked reduction at 1 hpt in our study. Another predicted target of miR168 is MiPUM2, a member of the PUF-family of RNA-binding proteins involved in post-transcriptional regulation through 3′-UTR interactions. Although no inverse correlation between miR168 and MiPUM2 expression was observed, MiPUM2 may be regulated through miRNA-independent pathways. Given its reported role in maintaining meristematic cell activity under stress conditions^34^, MiPUM2 might act independently to protect sensitive tissues in mango during heat exposure.
In mango, miR319 is predicted to target MiGAMYB, MiMYB33, MiTCP4, and MiMNS3. The expression patterns of miR319 and MiTCP4 suggest a functional regulatory interaction, consistent with findings in rose^35^and tomato^36^. Our transient expression assay in N. benthamiana further validates this interaction, demonstrating that miR319 significantly downregulates mango MiTCP4 expression. These results support the regulatory relationship inferred from expression profiles and align with findings in other species, such as Gossypium hirsutum^37^, reinforcing the conservation of this miRNA–target interaction in heat stress response and fruit physiology.
On the other hand, TCP4 has been implicated in maintaining ovule identity under heat stress^38^, and its downregulation may help limit inappropriate differentiation and modulate ROS responses. Notably, STRING-based protein interaction analysis (data not shown) predicted an interaction between TCP4 and Squamosa promoter-binding protein-like 9 (SPL9). The miR156/SPL9 regulatory module has been shown to influence ROS accumulation via the salicylic acid pathway^39^, suggesting that in mango, miR319 may fine-tune ROS homeostasis through a TCP4-SPL9 axis, thereby preventing oxidative damage^40^. Supporting this role, overexpression of ZmTCP14 in maize enhances drought tolerance by promoting ROS accumulation^41^, reinforcing the importance of TCP-like proteins in abiotic stress adaptation. Furthermore, MiTCP4 and MiGAMYB displayed temporally coordinated expression patterns, suggesting that both may be co-regulated by miR319. In other species, such as Solanum habrochaites, overexpression of sha-miR319d has been linked to improved heat stress tolerance. However, silencing of GAMYB-like 1 did not result in a significant phenotypic change under heat stress conditions^36^, indicating possible functional redundancy or context-specific activity. In contrast, miR159, a highly conserved miRNA family related to miR319^42^, has been shown to negatively regulate heat tolerance. In rice, miR159 overexpression increased plant sensitivity to heat, highlighting the role of GAMYB repression in thermotolerance^43^. Although our findings suggest that miR319 may participate in stress response regulation in mango, further experimental validation is needed to clarify its specific role, particularly concerning MiGAMYB regulation and thermotolerance.
We also identified a member of the miR482 family in mango, with predicted targets including a long non-coding RNA and several resistance genes, such as PR and RPP-13, members of the NBS-LRR family. In other plant species, miR482 is known to regulate NBS-LRR genes and their associated lncRNAs. For instance, in tomato, lncRNA23468 acts as an endogenous target mimic (eTM) of miR482b; its overexpression induced downregulation of miR482b and subsequent upregulation of NBS-LRR defense genes^44^. In litchi (Litchi chinensis), the miR482/2118 family induces phasiRNAs production from lncRNA loci, through mechanisms involving alternative splicing and polyadenylation, thereby broadening the post-transcriptional regulatory landscape^45^.
Although such mechanisms have been well documented in the context of plant–pathogen interactions, their roles under abiotic stress, particularly heat stress, remain poorly understood. Our data could represent initial support for the involvement of the miR482/lncRNA regulatory axis in the heat stress response. Beyond its role in defense against biotic stressors, miR482 may contribute to broader stress adaptation networks under abiotic stress. The repression of the lncRNA target observed in mango suggests a potential prioritization mechanism, whereby the plant could redirect resources to activate alternative pathways, such as phasiRNA production, that may serve as secondary post-transcriptional regulators. This hypothesis is supported by López-Virgen et al. (2024), who reported a significant upregulation of MiDCL3 and MiAGO4 at 24 hpt, two key components of the RdDM pathway. This suggests a temporal shift in regulatory dynamics, in which miRNAs initiate the early stress response, followed by the accumulation of siRNAs, including phasiRNAs, at later stages. These findings expand the functional understanding of miR482 in mango, highlighting its potential role in coordinating coding and non-coding gene RNA regulatory networks to enhance fruit resilience under heat stress.
Overall, the expression patterns observed for the analyzed miRNA-target pairs indicate that heat stress in mango triggers finely tuned, transient post-transcriptional regulation rather than strong inverse transcript-level correlations. This mode of regulation aligns with mechanisms described for miR168/AGO1 and miR482 in other plants [32, 46-47].^32,^^46^^,47^. Thus, the subtle changes detected likely represent homeostatic control mechanisms that allow mango to rapidly adjust to heat stress while avoiding the detrimental effects of over-suppressing key regulatory genes.
Within this framework, our study integrates sequencing data, bioinformatic prediction, and expression profiling to infer miRNA-target interactions, with functional validation achieved specifically for the miR319-MiTCP4 module. While this approach provides robust support for dynamic miRNA-mediated regulation, it does not directly resolve genome-wide miRNA-guided cleavage events. Thus, incorporating degradome sequencing represents an important future direction for further refining and validating the miRNA-target relationship in mango.
Conclusion
This study presents a comprehensive characterization of miRNAs expressed in mango mesocarp under heat stress, uncovering a complex regulatory landscape involving both conserved and novel miRNAs. Our findings suggest that specific miRNAs, such as miR168, miR319, and miR482, function as early regulators of the heat stress response by modulating the expression of genes associated with development, stress signaling, and post-transcriptional regulation. Through target prediction and expression analysis, we identified key regulatory interactions, including the miR168/AGO1 feedback loop, the miR319/TCP4 axis linked to ROS homeostasis, and, notably, the miR482-mediated regulation of a long non-coding RNA. This last interaction may represent a novel regulatory module under thermal stress, suggesting an expanded role for miR482 beyond its well-established function in biotic stress responses.
Although RT-qPCR analysis did not reveal strong inverse expression patterns for all predicted targets, the functional validation of the miR19/TCP4 module in a heterologous system supports the robustness of the regulatory relationships inferred in this study.
Altogether, our results highlight the evolutionary conservation and functional diversification of miRNA-mediated regulation in mango, supporting the idea that miRNAs play a central role in the fruit’s molecular adaptation to postharvest heat treatment. This work provides a valuable foundation for future functional studies and opens avenues for the development of biomarkers or molecular tools to enhance mango stress resilience and postharvest quality.
Methods
Plant material and HWT application
Mango (Mangifera indica L.) fruit of the cultivar Ataulfo were harvested from the orchard “La Aviación,” located in Escuinapa, Sinaloa, Mexico (22°48′48″ N, 105°31′ W). Selected fruits were uniform in shape and size, free from visible damage, and at physiological maturity (approximately 120 days after flowering). The fruits were transported to the laboratory at CIAD, A.C. (Culiacán, Sinaloa), disinfected with 200 ppm chlorinated water, and stored at 20 °C for 12 h. The HWT was applied according to USDA-APHIS guidelines^11^: fruits were immersed in water at 115 °F (46.1 °C) for 75 min, followed by hydrocooling at 77 °F (25 °C) for 30 min. Mesocarp tissue was sampled at 1, 3, 6, and 24 hpt, immediately frozen in liquid nitrogen, ground into a fine powder, and stored at -20 °C until further use. Untreated fruits were used as controls. Three biological replicates were collected for both treatment and control groups, with each replicate consisting of three individual mango fruits.
RNA extraction, small RNA library construction and sequencing
Total RNA was extracted from the samples using the method reported by López-Gómez and Gómez-Lim^48^, followed by DNase I treatment (Roche) to remove residual genomic DNA. RNA quantity was measured using a NanoDrop ND-1000 UV-Vis spectrophotometer (NanoDrop Technologies Inc., Wilmington, DE, USA), and RNA integrity was verified by electrophoresis on a 1% agarose gel under denaturing conditions. Samples collected at 1, 3, 6, and 24 hpt were selected for library preparation. 2 µg of total RNA was sent for library preparation and sequencing (Novogene Co.). Briefly, 3’ and 5’ adaptors were ligated to the 3’ and 5’ ends of small RNA, respectively. Then, the first strand cDNA was synthesized after hybridization with the reverse transcription primer. The double-stranded cDNA library was generated through PCR enrichment. After purification and size selection, libraries with insertions between 18 ~ 40 bp were ready for sequencing with NovaSeq SE50 with 10 M reads per sample. A total of twelve small-RNA libraries (three biological replicates per time point) were constructed. The resulting raw sequencing data have been deposited under NCBI Sequence Bioproject PRJNA1310765 and SRA accessions SAMN50782523, SAMN50782524, SAMN50782525, SAMN50782526.
Identification of known and novel miRNAs
Raw reads were trimmed to remove small RNA Illumina adapters using TrimGalore v0.6.10 (https://github.com/FelixKrueger/TrimGalore), followed by Trimmomatic v0.39^49^ to retain reads with a minimum quality of 20, a minimum length of 18 bp, and a maximum length of 27 bp. Quality control reports for both raw and trimmed reads were generated using FastQC v0.12.1 (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) and summarized with MultiQC v1.14^50^. Trimmed reads were aligned to sequences of non-coding RNA (rRNA, tRNA, snRNA, and snoRNA) from Rfam v15.0 (https://ftp.ebi.ac.uk/pub/databases/Rfam/CURRENT/fasta_files/) using Bowtie v1.3.1 with the option --un to keep the reads that were unaligned to Rfam-selected sequences. Known and novel miRNAs were identified with ShortStack v4.0.2^51^, using the following parameters: --genomefile to align the filtered reads to the mango cv. Alphonso reference genome v2.1 (GCF_011075055.1), --known_miRNAs to align mature miRNA sequences from miRBase to the mango genome, and --dn_mirna to activate a de novo search for miRNA loci.
The fasta file for known miRNAs was compiled from species phylogenetically close to M. indica^25^and available at miRbase v22.1 (https://www.mirbase.org/)^52^, including Arabidopsis thaliana,* Carica papaya*,* Citrus clementina*,* Citrus reticulata*,* Citrus sinensis*,* Citrus trifoliata*,* Fragaria vesca*,* Prunus persica*,* Solanum lycopersicum*,* Theobroma cacao*, and Vitis vinifera.
Precursor sequences and their corresponding secondary stem-loop structures were obtained as part of the ShortStack miRNA annotation output and subsequently used for their analysis, visualization, and minimum free energy analysis.
Differential expression analysis of miRNAs
Raw counts generated by ShortStack were imported into the R environment and analyzed using the DESeq2 package (v1.39.8) to identify differentially expressed small RNAs. Raw counts were normalized using the variance-stabilizing transformation (VST) method to perform principal component analysis (PCA) and identify global expression patterns and atypical replicates. One atypical replicate from each condition was removed for further analysis. The remaining eight libraries were normalized by the median of ratios method and used for pair-wise comparisons. sRNAs were considered differentially expressed genes (DEGs) if they met the criteria of an adjusted p-value ≤ 0.05 and |log2 fold change| > 1with respect to the 1 hpt. DEGs were visualized using volcano plots and heatmaps generated with the R/EnhancedVolcano v1.18 and R/pheatmap v1.0.12 packages in R.
Phylogenetic analysis of miRNA families of mango mesocarp
To gain deeper insights into the evolutionary history of mango miRNAs, a comprehensive phylogenetic analysis was performed using all pre-miRNA sequences identified in this study, along with homologous pre-miRNAs from Arabidopsis thaliana, Citrus sinensis, Fragaria vesca, and Solanum lycopersicum, obtained from the miRBase database (v22)^52^. Sequences were first aligned using Clustal Omega method in Geneious Prime 2025.0.3 (https://www.geneious.com). The resulting FASTA alignment file was used to construct a phylogenetic tree using the maximum likelihood method with the IQ-TREE program (version 3.0), applying the best-fit model predicted for our sequences and 100,000 bootstrap replicates to assess branch support^53^. The final tree was visualized and edited using iTOL (Interactive Tree of Life, v7)^54^.
Prediction, annotation, and functional analysis of target genes
Target prediction was conducted using the psRNAtarget tool^14^, utilizing the mature sequences of mango miRNAs identified in this study, along with transcript sequences from the mango transcriptome^12^ and genome^25^. Only genes that met plant-specific target prediction criteria established by Schwab et al. (2005) were considered^55^. The maximum expectation value was set to 3 to ensure specificity.
Target gene annotations were validated by cross-referencing with existing annotations in both the mango transcriptome and genome, as further confirmed via BLAST analysis against the NCBI non-redundant protein database^56^. Gene ontology (GO) analysis was performed to classify the predicted target genes based on biological process, molecular function, and cellular component. This functional annotation was carried out using OmicsBox software with Blast2GO mapping^57^ using the complete set of predicted target genes associated with all miRNAs identified by sRNA-seq.
Selection and validation of miRNAs and target genes by RT-qPCR
To select miRNAs for experimental validation, duplex formation energies between differentially expressed miRNAs and their predicted targets were calculated using the Vienna RNA Websuite^58^ miRNAs-target pairs with the most stable interactions (lowest free energy values) were selected: miR168 with MiAGO1 and MiPUM2; miR319 with MiTCP4 and MiGAMYB; and miR482 with a long non-coding RNA (Milnc-1). Alignments of miR319 and miR482 precursors were performed in Jalview version 2.11.4.1^59^.
For cDNA synthesis, 1 µg of total RNA from 1-, 3-, 6-, and 24-h post-treatment samples, as well as from control samples, was reverse-transcribed with SuperScript III kit reagents (Invitrogen), according to the manufacturer’s instructions. For miRNA quantification, cDNA was synthesized using a universal stem-loop RT primer based on the design reported by Chen et al. (2005)^60^with an additional N-tail at the 3’ end (Table S4). For target gene quantification, cDNA was synthesized using oligo (dT) primers. Specific forward primers were designed for each miRNA^61^, along with a universal reverse primer derived from the stem-loop sequence. For target genes, both forward and reverse primers were designed in non-conserved regions to avoid off-target amplification (Table S4).
RT-qPCR was performed on a StepOne Real-Time PCR system (Applied Biosystems, USA) under the following cycling conditions: 1 cycle at 95 °C, 40 cycles at 95 °C for 15 s and 60 °C for 1 min, followed by melt curve analysis. Relative expression levels were calculated using the 2^− ΔCt^ method^62^. U6 small nuclear RNA was used as the internal reference for miRNA quantification, while Actin 7 (Manin06g006480.1) served as the reference gene for target mRNA normalization (Table S4). Data visualization and statistical analyses were performed using GraphPad Prism (http://www.graphpad.com). A one-way ANOVA followed by the Tukey-Kramer test was applied to assess statistically significant differences (p < 0.05).
To compare qPCR results with sequencing data, the normalized mean read counts from the sRNA-seq dataset were used as a reference, as only a single delta value was applied for relative expression calculation.
Functional validation of miR319-MiTCP4 regulation
To assess the regulatory interaction between miR319 and its target MiTCP4, transient expression assays were performed in Nicotiana benthamiana following the protocol of Li (2011)^63^. Plants were grown under controlled conditions (23 ± 2 °C, with a 16/8 h light/dark cycle and 65% relative humidity) until they reached five weeks of age.
Constructs pEG100::miR319, pEG101::MiTCP4-WT, and pEG101::MiTCP4-mut, all driven by the CaMV35S promoter, were synthesized and cloned by GenScript (NovoPro bioscience Inc.). The mutated version (MiTCP4-mut) carries seven synonymous point mutations within the predicted miRNA binding site to prevent miRNA-mediated cleavage without altering the amino acid sequence. These plasmids, along with the empty vector pEG100 (EV), were introduced into Agrobacterium tumefaciens strain GV3101 via electroporation.
For co-infiltration, bacterial cultures were adjusted to OD_600_ = 0.4 and mixed at a 1:1 ratio. Infiltration was carried out on the abaxial leaf surface using needleless syringes. Each treatment included three biological replicates (independent plants), with two leaves infiltrated per plant. Leaf tissue was harvested 48 h post-infiltration for RNA extraction using the cetyltrimethylammonium bromide (CTAB) method^64^. 1 µg of DNA-free RNA was used for cDNA synthesis with the SuperScript III reverse transcriptase kit (Invitrogen), following the manufacturer’s protocol. Quantitative RT-PCR was performed using MiTCP4- specific primers, with PPA2 as the reference gene for N. benthamiana (Table S4). Relative transcript levels were calculated using the 2^− ΔΔCt^ method. Statistical significance among treatments was assessed by one-way ANOVA followed by the Tukey-Kramer pos hoc test (p < 0.05).
Supplementary Information
Below is the link to the electronic supplementary material.
Supplementary Material 1
Supplementary Material 2
Supplementary Material 3
Supplementary Material 4
Supplementary Material 5
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1USDA-APHIS. Treatment manual. 1–128. at (2016). https://www.aphis.usda.gov/sites/default/files/treatment.pdf.