Integrated analysis of transcriptome, small RNA, and degradome sequencing provides insights into mango anthracnose resistance

Abstract

Mango anthracnose is a major biotic stress that limits mango production worldwide. An integrated transcriptome, small RNA (sRNA), and degradome sequencing analysis was conducted to determine the complex molecular mechanisms of mango anthracnose resistance. The transcriptome sequencing of 12 samples from a mango anthracnose resistant genotype under control and stress conditions on the 3rd and 5th day post inoculation identified 3,462 and 8,341 differentially expressed genes (DEGs), respectively, which were mainly involved in transcription factor activity, defense response, an obsolete oxidation − reduction process, and flavonoid biosynthetic processes. The sRNA sequencing of the samples identified 372 known and 104 novel miRNAs. A total of 81 differentially expressed miRNAs were identified, of which three were differentially expressed at both the 3rd and 5th day post-inoculation (dpi), including ath-MIR166e-p5_1ss15AC, mtr-miR156e, and csi-miR3954_L + 1_2ss17CG21CT. According to degradome sequencing, 2,274 targets were predicted for 341 miRNAs. The combined analysis of transcriptome and sRNA sequencing identified 257 miRNA-mRNA interaction pairs. In these pairs, csi-miR3954_L + 1_2ss17CG21CT was up-regulated at both the 3rd and 5th dpi under stress, which could cleave multiple sites of an NAC gene (LOC123212502) that was down-regulated under stress. Overall, these miRNAs and genes provide a molecular foundation for the miRNA-mediated response to mango anthracnose stress and can be regarded as promising candidates for mango improvement.

Introduction

MicroRNA (miRNA) is a short, non-coding single single-stranded RNA sequence that silences the complementary mRNAs through the RNAi mechanism, which plays an important role in RNA transcription and protein translation [1, 2]. In plants, miRNAs play a critical role in responding to the inhibition of pathogen infection by regulating the production of reactive oxygen species (ROS) [3], hormone signal transduction [4], and upregulating the expression of resistant genes [5]. The Sl-miR6022-SlRLP6/10 module regulates tomato resistance to Phytophthora infestans by inducing the accumulation of ROS and reducing the jasmonic acid (JA) content and evapotranspiration (ET) [6]. Recent studies have shown that plant-derived miRNA can be transferred to pathogenic fungi and then silence the virulence gene of fungi [7, 8]. Cotton can export miR166 and miR159 to Verticillium dahliae hyphae by targeting the Clp-1 and HiC-15 genes related to pathogenicity, respectively [9]. Compared with other plant fungal diseases, the RNA silencing in mango anthracnose has been less extensively studied.

Colletotrichum complexes that consist of several strains pose a significant threat to mango, including Colletotrichum alienum, Colletotrichum fructicola, Colletotrichum siamense, Colletotrichum tropicale, and Colletotrichum asianum [14]. Colletotrichum siamense causes mango anthracnose in China, Mexico, and Thailand [10,11,12]. It infects the mango fruit before harvest or during the storage period, causing massive yield losses and affecting fruit quality. Black spots emerge and then enlarge and fuse, followed by rotting. Synthetic pesticide applications are currently used to prevent mango anthracnose [13]. However, the drug resistance of Colletotrichum complexes has increased due to the long-term use of fungicides, highlighting the need for environmentally friendly and effective approaches to control mango anthracnose [14, 15]. Recent studies have suggested that the overexpression of miRNA precursors, which regulate the expression of disease-resistance genes in plants through the mechanism of RNA interference (RNAi), can effectively control plant diseases [3, 6]. Therefore, understanding the regulatory role of miRNAs and their target genes would enable the development of new approaches to control mango anthracnose.

In conclusion, numerous studies have demonstrated that miRNAs are involved in the plant disease resistance process; however, the role of miRNA in mango anthracnose is not well understood. Studying the role that miRNAs play in mango anthracnose will lead to a better understanding of the mechanism of mango's resistance to anthracnose. Therefore, this study conducted transcriptome, small RNA (sRNA), and degradome sequencing using a mango anthracnose-resistant genotype named Jinhuang. A comprehensive and integrated analysis of the resulting datasets identified the miRNA-mRNA regulatory network. The miRNA regulation of gene expression was described under mango anthracnose stress, resulting in a better understanding of mango anthracnose resistance mechanisms.

Materials and methods

Plant materials, stress treatment, and RNA extraction

The Jinhuang mango (anthracnose resistant genotype [16]) was planted in a field at Baise, Guangxi Province, China. Healthy and uniformly-sized mango fruits were sterilized by soaking in 70% ethanol for 15 min twice, and the fruits were then washed with sterilized water. Colletotrichum siamense is a major pathogen responsible for mango anthracnose in China. The C. siamense strain HN10, isolated from anthracnose on mango fruit in Guangxi Province, China, was selected for subsequent inoculation studies. A 6 mm-diameter punch was used to make 2 − 3 mm-deep holes on the surface of mango fruits, then a 7-day C. siamense mycelium cake was placed into a hole. Mangoes inoculated with PDA cakes were used as the control. Inoculated mango fruits were sampled at 3 d post-inoculation (3rd dpi), referred to as T3 (treatment) and C3 (control), and 5 d post-inoculation (5 th dpi), referred to as T5 (treatment) and C5 (control). Two inoculated sites, from which the fungal mycelium cake was removed, were used as samples for a single biological replicate. There were three such biological replicate samples per treatment. All samples were immediately frozen in liquid nitrogen and then stored at −8℃ for RNA extraction. Twelve RNA libraries and twelve sRNA libraries of three biological replicates of T3, T5, C3, and C5 were constructed. The Trizol Reagent (Thermo Fisher Scientific, Waltham, MA, USA) was used for RNA extraction following the manufacturer’s instructions. The NanoDrop ND-1000 (NanoDrop, Wilmington, DE, USA) and Bioanalyzer 2100 (Agilent, Santa Clara, CA, USA) were used to analyze the integrity and quality of total RNA. The RNA integrity number (RIN) values of the RNA are listed in Table S1.

Transcriptome sequencing and analysis

The mRNA was specifically captured using oligo(dT) magnetic beads (Dynabeads Oligo (dT), cat.25–61,005, Thermo Fisher Scientific), which were then fragmented using a NEBNextR Magnesium RNA Fragmentation Module (Thermo Fisher Scientific). The fragmented RNA was reversed into complementary DNA (cDNA) using Invitrogen SuperScriptTM II Reverse Transcriptase (Thermo Fisher Scientific), following second-strand DNA synthesis using DNA polymerase I and RNase H (Thermo Fisher Scientific). After the heat-labile UDG enzyme treatment of the second-stranded DNAs, the ligated products were amplified with a polymerase chain reaction (PCR). The average insert size for the final cDNA libraries was 300 ± 50 bp. The 2 × 150 bp paired-end sequencing was performed on a Novaseq 6000 sequencing system (Illumina). The raw reads were filtered by Cutadapt v1.9 (https://cutadapt.readthedocs.io/en/stable/) to obtain high-quality clean reads. The clean reads were aligned to the mango genome (https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_011075055.1/), which was indexed beforehand, using HISAT2 v2.2.1 (https://daehwankimlab.github.io/hisat2/). The expression levels of all transcripts were estimated using StringTie v2.1.6 (http://ccb.jhu.edu/software/stringtie/) and ballgown (http://www.bioconductor.org/packages/release/bioc/html/ballgown.html). A differential expression analysis was performed using DESeq2 with a false detection rate (FDR) < 0.05 and |log2 FC|≥ 1. The gene ontology (GO) enrichment analysis of differentially expressed genes (DEGs) was conducted using R package of GOseq.

Small RNA sequencing and analysis

For sRNA sequencing, the sRNA libraries were constructed using TruSeq Small RNA Sample Preparation (Illumina) following the manufacturer’s instructions. The 50 bp single-end sequencing was performed on an Illumina Hiseq2500 sequencing system. The inhouse program ACGT101-miR v4.2 of LC Bio Technology CO. Ltd, (Hangzhou, China) was used to filter raw reads by removing adapter dimers, junk (i.e., fragments that did not map to the mango genome), low complexity reads (i.e., fragments with repetitive patterns), common RNA families (ribosomal RNA (rRNA), transfer RNA (tRNA), small nuclear RNA (snRNA), small nucleolar RNA (snoRNAs)), and repeats.

Identification of known and novel miRNAs

Unique sequences with lengths of 18 − 25 nt were mapped to known species precursors in miRbase 22.1 by a BLAST search to identify known miRNAs. The unmapped sequences were subjected to a BLAST analysis against the mango genome (https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_011075055.1/) and the resulting hairpin RNA structures containing sequences were predicted as novel miRNAs using RNAfold (http://rna.tbi.univie.ac. at/cgi-bin/RNAfold.cgi). This was based on the following criteria: (1) number of nucleotides in one bulge in the stem (≤ 12), (2) number of base pairs in the stem region of the predicted hairpin (≥ 16), (3) cutoff of free energy (kCal/mol ≤ −15), (4) length of hairpin (up and down stems + terminal loop ≥ 50), (5) length of hairpin loop (≤ 200), (6) number of nucleotides in one bulge in the mature region (≤ 4), (7) number of biased errors in one bulge in the mature region (≤ 2), (8) number of biased bulges in the mature region (≤ 2), (9) number of errors in the mature region (≤ 4), (10) number of base pairs in the mature region of the predicted hairpin (≥ 12), and (11) percentage of mature regions in the stem (≥ 80). The GSTAr v1.0 software was used to predict target genes by identifying binding sites with default parameters. The differential expression of miRNAs based on normalized counts was determined using a Student’s t-test with p < 0.05 and |log2 FC|≥ 1.

Degradome sequencing and analysis

Equivalent amounts of RNA from the controls (C3 and C5) and treatments (T3 and T5) were combined to generate two degradome libraries. The purified RNA with polyA was mixed with biotinylated random primers and ligated with 5’ adaptors. The cDNA was synthesized with 3’ adaptors followed by size selection with AMPureXP beads and then used as a template for the PCR. The Illumina Hiseq2500 was used to sequence the cDNA library. Raw data were filtered using Trimmomatic v0.35, then the reads mapped to rRNA, tRNA, snoRNAs, and repeats were removed. The clean reads were mapped to the mango transcriptome, and a degradome density file was generated. GSTAr v1.0 was used to predict the target genes of the miRNA. By combining the predicted target genes of the miRNA and genes in the degradome density file, the t-plots were generated.

Reverse-Transcription Quantitative PCR (RT-qPCR)

The RNA samples used for the RNA-seq were also used in an RT-qPCR analysis. The qPCR primers were designed using the primer-BLAST tool (https://www.ncbi.nlm.nih.gov/tools/primer-blast/). Reverse transcription and qPCR were conducted using ToloScript All-in-one RT EasyMix for qPCR (22,107, Tolobio, Shanghai, China) and 2 × Universal Blue SYBR Green qPCR Master Mix (G3326-01, Servicebio, Wuhan, China) according to each manufacturer’s instructions. The housekeeping gene Actin was selected as an internal standard in the qPCR analysis. Each sample was replicated three times. Expression data were analyzed using the 2^−ΔΔCt method. All primers are listed in Table S2.

Results

Transcriptome sequencing

In the study, 12 samples representing Jinhuang (anthracnose resistant genotype) under control (non-inoculated) and stress (C. siamense inoculated) conditions at two time points (3rd and 5 th dpi) (Fig. S1) were sequenced. A total of 481.60 million reads were generated from the pair-end sequencing of these 12 samples. After applying quality filters, 473.26 million clean reads (98.26%) were obtained for further analysis. These reads were mapped to the mango genome sequences, with an average of 402.77 million clean reads (85.10%) mapped to the mango genome. The data statistics indicated high-quality transcriptome sequencing (Table S3).

Differentially expressed genes under C. siamense infection

A differential expression analysis was performed to identify DEGs. A total of 19,557 genes were found to be expressed in at least one of the samples. At the 3rd and 5 th dpi, 3,462 and 8,341 genes were found to show a significant differential expression between the control and stress samples (3rd dpi: 1,943 up-regulated and 1,519 down-regulated; 5 th dpi: 4,325 up-regulated and 4,016 down-regulated), respectively (Fig. 1A). According to a Venn analysis, 2,606 common DEGs were identified in T3_vs_C3 and T5_vs_C5 (Fig. 1B). The time trend of the gene expression analysis revealed four clusters, with cluster 4 containing 1,285 genes that showed a low expression level in the control and a high expression level under stress (Fig. 1C). Overall, compared with the 3rd dpi, mango showed a greater response to C. siamense infection at the 5 th dpi. Three genes were selected to validate the RNA-seq data. The RNA-seq and RT-qPCR data were significantly correlated with R values of 0.943, 0.905, and 0.848 (p < 0.05) for each of the three genes (Fig. 2), indicating that the RNA-seq data were reliable.

Fig. 1

An overview of the differentially expressed genes (DEGs) in the mango response to Colletotrichum siamense infection. A The number of DEGs at the 3rd and 5 th day post inoculation (dpi) under control and stress conditions. B Venn analysis of DEGs at the 3rd and 5 th dpi under control and stress conditions. C Time trend of the gene expression analysis of DEGs at the 3rd and 5 th dpi under control and stress conditions

Fig. 2

Validation of the gene expression by RT-qPCR. Three randomly selected genes, XM_044654894.1, XM_044611652.1, and XM_044631664.1, were selected for an RT-qPCR analysis. Each treatment had three biological and three technical replicates. The RNA-seq (right-axis) and RT-qPCR (left-axis) data representing the five genes were significantly correlated (p < 0.05)

There were 113 and 69 enriched GO terms (Q value < 0.001) for the DEGs for T3_vs_C3 and T5_vs_C5, respectively (Tables S4 and S5). In T3_vs_C3, chloroplast thylakoid membrane, thylakoid, chloroplast thylakoid, chloroplast, obsolete oxidation–reduction process, chloroplast envelope, flavonoid biosynthetic process, chloroplast stroma, photosynthesis, and extracellular region were the top10 enriched GO terms (Fig. 3A). In T5_vs_C5, rRNA transcription, transcription factor binding, regulation of cellular respiration, obsolete oxidation–reduction process, flavonoid biosynthetic process, defense response, DNA-binding transcription factor activity, cell population proliferation, quercetin 3-O-glucosyltransferase activity, and extracellular region were the top10 enriched GO terms (Fig. 3B). Obsolete oxidation–reduction process and flavonoid biosynthetic process were the common GO terms in both T3_vs_C3 and T5_vs_C5.

Fig. 3

Enriched gene ontology (GO) terms for the differentially expressed genes (DEGs) at the 3rd (A) and 5 th (B) day post inoculation (dpi) under control and stress conditions. The size of the circles indicates the number of hit genes. BP, CC, and MF represent biological process, cellular component, and molecular function, respectively. A total of 113 and 69 GO terms were significantly enriched (FDR-value < 0.001) at the 3rd and 5 th dpi, respectively. The details of the enriched GO terms are shown in Tables S4 and S5

Small RNA sequencing and miRNA identification

A total of 12 sRNA libraries were constructed, and sequences for identifying C. siamense infection-related miRNAs in mango were obtained. A total of 178.55 million reads with an average of 14.87 million reads per sample were obtained. After filtering reads, a set of 127.96 million clean reads was obtained for the following analysis. A total of 28.98 million validated reads were obtained by removing repeat reads and clean reads mapping to rRNA, tRNA, snRNA, and snoRNA (Table S6). The most abundant length distribution of the unique sRNA reads was 21 nt (52.39%) (Fig. S2), followed by 22 nt (14.89%), 20 nt (7.71%), and 19 nt (7.45%).

To identify the known miRNAs in mango, the filtered reads were mapped to the miRNAs of plant species from miRBase (https://www.mirbase.org/), and 372 unique known miRNAs were identified from all samples. The reads that could not be mapped to miRBase, were mapped to the mango genome. Additionally, the precursor sequences were obtained by extending the full-match mapped reads, which could fold into a potential stem-loop, and a total of 104 unique novel miRNAs were obtained.

Differentially expressed miRNAs under C. siamense infection

To identify miRNA that responded to C. siamense infection, the expression patterns of all known and novel miRNAs were studied. All differentially expressed miRNAs are shown in Fig. 4 and Table S7. At the 3rd dpi, 15 miRNAs showed a significant differential expression between stress and control samples, of which two exhibited absolute log2 FC values greater than 2. At the 5 th dpi, 69 differentially expressed miRNAs were found under C. siamense infection, of which 13 had absolute log2 FC values greater than 2 (Table S7). A total of three miRNAs were identified in T3_vs_C3 and T5_vs_C5 according to a Venn analysis, including ath-MIR166e-p5_1 ss15 AC, mtr-miR156e, and csi-miR3954_L + 1_2 ss17 CG21 CT. Specifically, ath-MIR166e-p5_1 ss15 AC was down-regulated, while mtr-miR156e and csi-miR3954_L + 1_2 ss17 CG21 CT were up-regulated under stress. Based on a GO enrichment analysis of the computationally predicted targets, eight and 26 GO terms were significantly enriched (Q value < 0.001) in T3_vs_C3 and T5_vs_C5, respectively, including defense response and DNA-binding transcription factor activity (Tables S8 and S9). These miRNAs were involved in the stress response under C. siamense infection in mango.

Fig. 4

Differential expression of stress-response miRNAs at the 3rd and 5 th day post inoculation (dpi) under control and stress conditions. The expression levels are standardized by row. The miRNAs with asterisks were differentially expressed at both the 3rd and 5 th dpi

Target prediction according to the degradome

A degradome analysis was conducted to identify the targets of plant miRNAs. A total of 26.21 million reads were generated from control and treated samples, and 11.76 million unique reads were obtained, which were used for the identification of cleavage sites (Table S10). A total of 257 and 269 non-redundant targets with p ≤ 0.05 and category ≤ 4 were identified in control and treated samples, respectively, including a total of 127 (16 novel and 111 known) miRNAs (Table S11), which were represented in target plots (T-plots).

According to the degradome sequencing of treated samples, a total of 257 significant targets were identified, mainly belonging to transcription factors (TFs), transporter genes, disease resistance genes (e.g., putative disease resistance RPP13-like protein 1 and NBS-LRR resistance protein), and cellular enzymes (e.g., casein kinase and 3-ketoacyl-CoA thiolase). The TF encoding genes represented a considerable proportion of the targets, including the NAC (52), HD-ZIP (24), TCP2 (18), ARF (16), MYB (3), ZF (2), and GRAS (2) families, as shown in Fig. 5. A GO enrichment analysis was further conducted to reveal the potential role of miRNA targets in response to C. siamense stress in mango. The most significant GO terms were transcription factor binding, cytosol, double-stranded DNA binding, rRNA transcription, and regulation of cellular respiration (Fig. S3).

Fig. 5

A network representing the relationships between miRNAs and their target genes associated with the Colletotrichum siamense response. The red eclipses represent the miRNAs and the purple rectangles represent the target genes. HP, hypothetical protein; GRAS, scarecrow-like protein; RAP, ethylene-responsive transcription factor; PR, pathogenesis-related protein; HTP, heptahelical transmembrane protein; KI, kinase; 3-KAT, 3-ketoacyl-CoA thiolase; TANGO2, transport and Golgi organization 2; HyD, Hyoscyamine 6-dioxygenase; RPP, disease resistance protein; PPR, pentatricopeptide repeat-containing protein; GRF, growth-regulating factor; ARF, auxin response factor

Correlation analysis of miRNAs and their target’s expression profiles

According to the degradome results, the expression profiles of both the C. siamense response miRNA and their target genes were integrated and analyzed to reveal the mediatory role of miRNA under C. siamense infection in mango. Through a Pearson correlation analysis, a total of 34 and 257 miRNA interaction pairs were identified under stress conditions at the 3rd and 5 th dpi, respectively. These correlations can be either positive or negative. Of the 34 pairs, 19 were positive correlations and 15 were negative correlations at the 3rd dpi. At the 5 th dpi, 141 were positively correlated and 116 pairs were negatively correlated (Fig. 6, Tables S12 and S13). These results indicated the miRNA-mediated C. siamense response at the 5 th dpi was stronger than that at the 3rd dpi.

Fig. 6

Integrated correlation networks of the identified miRNAs and their target genes at the 3rd (A) and 5 th (B) day post inoculation (dpi). Ellipses indicate the target genes. Triangles indicate miRNAs. Red lines indicate a positive correlation and green lines indicate a negative correlation

To determine the mechanism by which miRNA degraded tmRNA, the negative pairs were further analyzed. At the 3rd dpi, 15 negative pairs consisting of nine miRNAs and 10 genes were identified. Similarly, at the 5 th dpi, 116 negative pairs consisting of 46 miRNAs and 74 genes were identified (Fig. 6, Tables S12 and S13). The number of genes was higher than that of miRNAs, which indicated that a single miRNA could cleave multiple targets. We further focused on three miRNAs (ath-MIR166e-p5_1 ss15 AC, mtr-miR156e, and csi-miR3954_L + 1_2 ss17 CG21 CT), which were significantly differentially expressed at both the 3rd and 5 th dpi. Both ath-MIR166e-p5_1 ss15 AC and mtr-miR156e had no significantly negatively expressed target genes, while csi-miR3954_L + 1_2 ss17 CG21 CT was significantly up-regulated. Additionally, its target gene LOC123212502, coding a NAC domain-containing protein, was significantly down-regulated at the 5 th dpi. The csi-miR3954_L + 1_2 ss17 CG21 CT could cleave several transcripts of LOC123212502 with multiple sites (Fig. 7).

Fig. 7

Expression profiles of the csi-miR3954_L + 1_2 ss17 CG21 CT-LOC123212502 interaction pairs and their validation. A Expression profile of csi-miR3954_L + 1_2 ss17 CG21 CT at the 3rd and 5 th day post inoculation (dpi) under control and stress conditions. B Expression profile of its target LOC123212502, an NAC gene, at the 3rd and 5 th day post inoculation (dpi) under control and stress conditions. C-F csi-miR3954_L + 1_2 ss17 CG21 CT cleaves different transcripts of LOC123212502 at multiple sites. The red dots indicate the cleavage nucleotide positions on the target gene

Discussion

High-throughput sequencing technologies enable the molecular basis of plant responses to pathogens to be better understood. In this study, we used three high-throughput approaches, namely, transcriptome sequencing, sRNA sequencing, and degradome sequencing, to determine the genetic and molecular responses behind mango anthracnose resistance. A comprehensive transcriptome analysis of the mango anthracnose-resistant genotype was performed. The obsolete oxidation–reduction process and flavonoid biosynthetic process were enriched in both T3_vs_C3 and T5_vs_C5. It has been reported that ROS play a crucial role in the biotic stress response [17], and the mango showed improved resistance to Colletotrichum gloeosporioides by activating ROS [18]. Flavonoid metabolites could scavenge ROS, thus contributing to biotic stress response [19]. In this study, the mango genotype Jinhuang may have scavenged ROS by synthesizing flavonoids to improve resistance to mango anthracnose. Furthermore, the number of DEGs at the 5 th dpi was higher than that at the 3rd dpi (Fig. 1A), suggesting that the response at the 5 th dpi was more evident. We focused on the GO enrichment terms at the 5 th dpi, and four TF-related GO terms were obtained from the top 10 enriched GO terms, implying that TFs played a crucial role in mango resistance. Transfer factors can regulate gene networks in the plant biotic stress response [20]. It has been reported that NAC coding genes negatively regulate defense responses against pathogens by suppressing pathogenesis-related (PR) coding genes [21, 22]. In this study, many NAC coding genes were down-regulated at the 5 th dpi, suggesting their potential roles in mango anthracnose resistance.

Small non-coding RNAs, especially miRNAs, are important modulators of gene expression at the post-transcription level and are potential tools for crop improvement [23]. Several miRNAs have been identified in the stress response of various plants using an sRNA sequencing approach [24, 25]. In mango, temperature-responsive miRNAs have been reported [26], but the role of miRNA in C. siamense resistance has not yet been characterized. In this study, we employed sRNA sequencing of the C. siamense resistant mango at the 3rd and 5 th dpi, and identified 476 miRNAs, of which 104 were novel miRNAs. Novel miRNAs might be involved in species-specific gene regulation. Similar to the transcriptome results, the number of differentially expressed miRNAs at the 5 th dpi was more than that at the 3rd dpi, indicating a stronger response at the 5 th dpi. The GO enrichment analysis highlighted that these targets were mainly involved in TF-related GO terms, which was consistent with the transcriptome results. These results indicated that the miRNAs of mango Jinhuang responded to C. siamense stress mainly by regulating TF coding genes. Similarly, the miRNAs of chickpea regulate the resistance to Ascochyta blight by targeting the WRKY, MYB, NAC, and ARF TF families according to high-throughput transcriptome and sRNA sequencing [24].

Degradome sequencing provides reliable insights into miRNA-mRNA regulation networks. Based on the current findings, a hypothetical scheme was postulated based on the transcriptome, sRNA, and degradome sequencing results (Fig. 8). Notably, the LRR coding genes involved in pathogen recognition were up-regulated and targeted by bol-MIR9410-p5_2 ss5 TG18 TA (Fig. 8). Similarly, miR482 targets NBS-LRR genes under fungal pathogen stress in cotton [27]. Of the three differentially expressed miRNAs at both the 3rd and 5 th dpi under stress (Fig. 4), csi-miR3954_L + 1_2 ss17 CG21 CT could cleave several transcripts of an NAC gene (LOC123212502) with multiple sites. The csi-miR3954_L + 1_2 ss17 CG21 CT was up-regulated under stress, while LOC123212502 was down-regulated (Fig. 8). Similarly, miR164 regulates the expression of NAC21/22 genes by controlling JA signaling under apple leaf spot disease stress [4]. In mulberry, Mno-miR164a and its target TF gene MnNAC100 manipulate the defense response to Botrytis cinerea stress [28]. Additionally, bol-MIR9410-p5_2 ss5 TG18 TA regulates the flavonol 4'-sulfotransferase gene (Fig. 8), which is involved in the synthesis of the secondary metabolite, flavonoid. The results were consistent with the GO enrichment analysis of the transcriptome, implying that flavonoid metabolites played an important role in C. siamense resistance. Similarly, the miR172a-SNB module controls multiple diseases via the accumulation of the secondary metabolite lignin in rice [29]. The functional characterization of the miRNA-mRNA regulatory role under C. siamense stress will be determined in a future study.

Fig. 8

A postulated model of various physiological and biochemical events incurred during the interaction of Colletotrichum siamense with the anthracnose-resistant mango genotype Jinhuang. The red arrows indicate up-regulation and green arrows indicate down-regulation under stress at the 5 th day post inoculation

In summary, this study represents the first attempt to integrate a transcriptome, sRNA, and degradome sequencing analysis to identify key miRNA-mRNA regulation networks in the mango response to C. siamense. A total of 3,462 and 8,341 DEGs were identified under control and stress conditions in mango at the 3rd and 5 th dpi, respectively, and they were mainly involved in transcription factor activity, defense response, an obsolete oxidation–reduction process, and flavonoid biosynthetic processes. A total of 81 differentially expressed miRNAs were identified, of which three were differentially expressed at both the 3rd and 5 th dpi, including ath-MIR166e-p5_1 ss15 AC, mtr-miR156e, and csi-miR3954_L + 1_2 ss17 CG21 CT. According to the degradome sequencing, 2,274 targets were predicted for 341 miRNAs. The combined transcriptome and sRNA sequencing identified 257 miRNA-mRNA interaction pairs. Notably, csi-miR3954_L + 1_2 ss17 CG21 CT was up-regulated at both the 3rd and 5 th dpi under stress. This miRNA could cleave multiple sites of an NAC gene (LOC123212502) that was down-regulated under stress. Overall, these miRNAs and genes provide a molecular basis for the miRNA-mediated C. siamense stress response in mango, and could serve as candidates for developing C. siamense resistance in mango using a genetic engineering approach.