Small RNA profiling of Cavendish banana roots inoculated with Fusarium oxysporum f. sp. cubense race 1 and tropical race 4

Abstract

Fusarium wilt, caused by the soil-borne fungal pathogen, Fusarium oxysporum f. sp. cubense (Foc), is considered as one of the most threatening diseases of banana. The Cavendish variety, resistant to Foc race 1 (R1), is susceptible to tropical race 4 (TR4), an aggressive race of the pathogen which is of increasing concern worldwide. Previous studies have revealed that plant small RNAs (sRNAs) play crucial roles in the host response to pathogen infection. To investigate the roles of sRNAs involved in the interaction of the banana-Foc pathosystem, small RNA profiles of Cavendish banana roots inoculated with Foc TR4 and Foc R1 were obtained and analyzed in the present study using Next-Generation Sequencing (NGS) technology. A total of 112 discrete mature microRNAs (miRNAs) belonging to 22 known miRNA families were found across all constructed sRNA libraries. The expression of miR166, miR159 and miR156 was up-regulated in TR4-inoculated samples as compared to mock-inoculated samples, while the expression of these miRNAs was approximately the same in R1-inoculated and mock-inoculated samples. Consistent with the sequencing data, qRT-PCR results demonstrated up-regulation of these miRNAs and down-regulation of their target genes in TR4-inoculated samples, but not in R1-inoculated samples. Considering Cavendish banana is resistant to R1 and susceptible to TR4, it is possible that these sRNAs and their target genes are involved in particular plant defence pathways such as salicylic acid-based defence. The findings will pave way for future investigations of the defence mechanism and potential approaches of resistance improvement.

Background

Fusarium wilt is a devastating disease of banana, caused by the fungal soil-borne pathogen, Fusarium oxysporum f. sp. cubense (Foc). It infects the banana plants from the roots, then travels up and blocks the transport of water and nutrients, causing plant death (Pietro et al. 2003; Ploetz 2015b). The outbreak of Foc race 1 (R1) almost destroyed the banana cultivar Gros Mitchel, which was almost the sole banana cultivar traded worldwide in the 1900s (Ploetz 2015b). Even though the alternative option, Cavendish, is resistant to Foc R1, a relatively new Foc race, tropical race 4 (TR4) is pathogenic on all available commercially produced banana cultivars, including Cavendish. The longevity of F. oxysporum in soil is very high and is confounded by the lack of effective control methods (Hwang and Ko 2004; Ploetz 2006; Dale et al. 2017). In order to develop potential novel solutions for management of Foc, a greater understanding of the molecular mechanisms underlying the interaction between the pathogen and the host needs to be explored.

Small RNAs (sRNAs) are short non-coding RNA molecules of 20–26 nucleotides (nt) that regulate a range of eukaryotic biological processes and play crucial roles in responding to biotic and abiotic stresses in plants. Two major classes of sRNAs in plants are microRNAs (miRNAs) and small interfering RNAs (siRNAs). The class miRNAs are derived from endogenous genes and processed from single-stranded hairpin RNA (hpRNA), which can form a stem-loop secondary structure. With regards to the siRNAs, they generally originate from double-stranded RNA (dsRNA) either from exogenous viruses or from endogenous repetitive regions or other complementary transcripts (Axtell 2013; Kamthan et al. 2015). The intermediate dsRNA or hpRNA is then cleaved by Dicer-Like (DCL) endonucleases into 20–26 base-pair (bp) double-stranded duplexes, one strand of which is loaded into Argonaute (AGO) proteins and forms the RNA-Induced Silencing Complex (RISC). Then the RNAi mechanism is achieved by either the post-transcriptional silencing (PTGS) or transcriptional silencing (TGS) pathways. In PTGS, the RISC leads the miRNA or siRNA to its complementary sequences of the targeted mRNAs and represses expression by slicing the targeted mRNAs. In general, the action of siRNAs leads to cleavage of mRNAs with perfectly matching sequences; whereas miRNAs target transcripts that are not fully complementary. In TGS, the 24 nt siRNAs bind to genomic DNA and inhibit gene transcription by DNA methylation and chromatin remodeling (Vaucheret 2006; Castel and Martienssen 2013).

Numerous studies have demonstrated that plant sRNAs are actively involved in the defence response upon pathogen attack (Ruiz-Ferrer and Voinnet 2009; Sheeba et al. 2013; Gupta et al. 2014; Islam et al. 2017). Up-regulation of miR393 was observed following infection by Pseudomonas syringae in Arabidopsis thaliana, which consequently repressed the expression of auxin receptor genes so that growth of the pathogenic P. syringae was inhibited (Navarro et al. 2006). Research by Fahlgren et al. (2007) confirmed this phenomenon. They also demonstrated miR160 and miR167 were up-regulated after P. syringae infection, suppressing the expression of transcription factors of the ADP-ribosylation factor (ARF) family, which are key factors of auxin signaling. Overexpression of miR396a-5p was found to increase the susceptibility to Phytophthora nicotianae infection in tobacco (Chen et al. 2015), suggesting miR396a-5p plays a crucial role in biotic stress by targeting growth-regulating factors (GRFs). An early example of sRNAs involved in host-fungi interaction was in loblolly pine challenged by the rust fungus Cronartium quercuum f. sp. fusiforme (Lu et al. 2007). In the loblolly pine-rust system, ten miRNA families such as miR156, miR159 and miR160 were found to be significantly down-regulated during the infection. The expression of their target transcripts, ARF10, SQUAMOSA-PROMOTER BINDING PROTEIN (SBP), MYELOBLASTOSIS (MYB) and Aux/IAA protein genes, was correspondingly enhanced to defend against the infection. In tomato infected with Fusarium oxysporum, two miRNAs, miR482f and miR5300, targeting NB-LRR immune genes were repressed in resistant cultivars so that their target genes were activated and able to play their roles in the plant immune response (Ouyang et al. 2014). In addition to the endogenous sRNAs functioning in the host, miR166 and miR159 in cotton were found to be exported into the fungal pathogen Verticillium dahliae (Zhang et al. 2016), which inhibited the expression of two fungal virulence genes and increased disease resistance. In terms of the banana-Foc pathosystem, studies on banana sRNA identification have emerged in recent years (Ghag et al. 2015; Song et al. 2015; Song et al. 2016), however little is known concerning the roles of banana sRNAs in response to Foc invasion. In a recent report, members of the miR169 family showed higher expression level in resistant banana cultivars upon Foc inoculation than in susceptible cultivars, suggesting their involvement in Foc resistance in banana (Song et al. 2018).

In this study, we applied deep sequencing technology to obtain comparative sRNA profiles of Cavendish banana roots inoculated with Foc R1 and Foc TR4, aiming to identify the sRNAs that were differentially regulated by R1 and TR4 and thus gain more knowledge of the molecular mechanisms of the host’s defence response during pathogen invasion.

Results

An overview of the sequencing data

At 28 days post inoculation (dpi), banana root tissues were harvested for total RNA extraction at the two inoculation locations. Twelve small RNA libraries were constructed from four treatment groups, namely the TR4_C (TR4 mock-inoculated), TR4_I (TR4-inoculated) both from the Northern Territory (NT) site and the R1_C (R1 mock-inoculated), and R1_I (R1-inoculated) from the Brisbane site. In all cases three biological replicates were used for each treatment. From the twelve sRNA libraries, a total of 12.21 Gb of raw data was generated. Following adapter removal and high-quality read screening, an average of 2,869,477, 3,438,027, 2,112,927 and 2,262,496 processed reads from three replicates were obtained from TR4_C, TR4_I, R1_C and R1_I groups, respectively.

Identification of known miRNAs in response to Foc infection

To investigate the regulatory response of known conserved miRNAs to Foc infection, the MiRDeep2 package was used to align the processed reads from all four groups to the miRNA database, miRBase. A total of 112 miRNAs belonging to 22 known miRNA families were found with at least one raw read from the sRNA libraries mapping to the database. More than ten members from each of the miR156, miR166, miR167, miR395 and miR396 families were found in the twelve sRNA libraries, while fewer members from other miRNA families were present in these sRNA libraries (Fig. 1). With regard to the overall expression of different miRNA families (Fig. 2), miR166, miR396 and miR159 had the highest abundance regardless of their differential expression in each treatment group. Whereas miR398, miR529, miR444, miR2118, miR394, miR393 and miR5077 were only found in TR4 libraries but not in R1 libraries.

Fig. 1

Number of miRNA members belonging to the miRNA families found in this study. The miRNA families listed on the X-axis are the identified miRNA families in the twelve sRNA libraries. Among them, miR156, miR395, miR166 and miR167 possessed the most members present in the sRNA libraries

Fig. 2

MiR166, miR396 and miR159 possessed the highest expression in the treatment groups. The miRNA families listed on the X-axis are the identified miRNA families in the twelve sRNA libraries. The normalized read count on the Y-axis was calculated as “reads per million reads” (RPMR) for each sRNA library. The logarithmic scale was used to display the normalized read count. Error Bars represent Standard Error of the Mean (SEM)

To investigate the role of miRNAs in root responses to Foc inoculation, differential expression analysis was performed for miRNAs which had at least one normalized read (RPMR, reads per million reads) in each sRNA library (Fig. 3). Fold change of TR4_I versus TR4_C and R1_I versus R1_C is shown in Additional file 1: Figure S1. Direct comparison of miRNA expression of TR4 groups with R1 groups was not conducted due to trials being at different locations. However, comparison patterns of the inoculated group versus control group in both the TR4 and R1 groups were obtained. Different expression patterns of the TR4 groups and R1 groups were observed. Overall, less difference was observed between R1_I and R1_C than between TR4_I and TR4_C. The expression of miR156a, miR159f and miR166c in TR4_I was 5, 3.9 and 2.9 times higher than in TR4_C, while in the R1 samples, the fold change for these same miRNAs between R1_I and R1_C was lower (Additional file 1: Figure S1). In addition to the above miRNAs with different expression patterns in response to TR4 and R1 inoculation, miR396 expression was reduced in TR4 inoculated samples, while miR408 expression was reduced in both TR4 and R1 inoculated samples.

Fig. 3

Expression of miRNAs in each treatment group. Members of miR156, miR159 and miR166 families showed increased expression in TR4_I as compared to TR4_C. Not as much increase was induced by R1 inoculation. The normalized read count on the Y-axis was calculated as “reads per million reads” for each sRNA library. The logarithmic scale was used to display the normalized read count. Error Bars represent Standard Error of the Mean (SEM)

Identification of novel miRNAs in four sRNA libraries

To predict novel miRNAs present in the twelve sRNA libraries, processed reads were mapped to the banana genome sequence (http://banana-genome-hub.southgreen.fr/download) using the MiRDeep2 package. Reads that mapped to the miRBase and Rfam databases were removed (Griffiths-Jones et al. 2008). The folding free energy of the predicted miRNA precursors was determined and those with randfold P-value higher than 0.05 were removed (Bonnet et al. 2004). As a result, 29 putative novel miRNAs were obtained (Additional file 2: Table S1).

MiRNA target gene prediction and GO analysis

Target genes of the putative known and novel miRNAs were predicted using the online tool psRNATarget. With a maximum expectation score of 3.5, 279 and 365 target genes were predicted for the putative known miRNAs and novel miRNAs found above, respectively (Additional file 3: Table S2). To further understand their function, gene ontology (GO) analysis of the miRNA target genes was performed using Blast2Go. Overall, seven putative biological processes, three molecular functions and six cellular components were predicted for the putative target genes of the identified known miRNAs (Fig. 4). Binding and cell part comprised the most genes in molecular function category and cellular component category respectively. Regulation of transcription showed the highest percentage in biological processes (Fig. 5a), suggesting a high proportion of these target genes were involved in regulating gene expression during the pathogen invasion. GO analysis was also performed for predicted target genes of the putative novel miRNAs. Seven biological processes, three molecular functions and seven cellular components were predicted for the novel miRNA target genes. In the gene distribution of biological processes, “response to stimulus” accounted for 30% (Fig. 5b), which indicated that these miRNAs and their targeted genes could be actively involved in the host response to F. oxysporum.

Fig. 4

Gene ontology analysis of putative target genes of the identified known miRNAs. A total of six terms in the cellular component category, three terms in the molecular function category and seven terms in the biological process category were obtained. “Binding” and “Cellular process” possessed the most genes among all the GO terms found in this study

Fig. 5

Distribution of the predicted target genes in gene ontology terms of the biological process category. a Predicted target genes of the identified known miRNAs. b Predicted target genes of the putative novel miRNAs. “Regulation of transcription” and “Response to stimulus” comprised the most abundant target genes of the identified known miRNAs and putative novel miRNAs, respectively

qRT-PCR validation of miRNAs and target genes

Quantitative real-time PCR (qRT-PCR) was conducted to validate the expression level of the miRNAs and their target genes. Three miRNAs of the most interest, miR159, miR156 and miR166, were chosen for validation. For miRNA qRT-PCR, miR399 was used as the reference miRNA for normalization due to its constant expression level in different banana tissues, as well as under stress conditions (Muthusamy et al. 2015). L2 was used as the reference gene for target gene expression normalization (Podevin et al. 2012).

As shown in Fig. 6, miR166, miR159 and miR156 showed higher expression in TR4_I than in TR4_C. In contrast, expression of their predicted target genes, HOX9, GAMYB and SPL12 was significantly reduced in TR4_I as compared to in TR4_C. This negative correlation of expression between miRNAs and target genes is consistent with the sequencing data (Fig. 3). Thus, the up-regulation of these miRNAs upon TR4 inoculation caused downstream down-regulation of their target genes, suggesting the involvement of these genes during the interaction.

Fig. 6

qRT-PCR validation of selected miRNAs and their target genes. Negative-correlation of the expression between miR166, miR159, miR156 and their target genes was observed in TR4_C and TR4_I. a miR166 and its target gene, HOX9 (Homeobox-leucine zipper protein 9) (GSMUA_Achr5T18530_001). b miR159 and its target gene, GAMYB (GSMUA_Achr7T21560_001). c miR156 and its target gene, SPL12 (Squamosa Promoter Binding Protein-like 12) (GSMUA_Achr10T02970_001). One-way analysis of variation (ANOVA) was used to determine the significant difference of the target gene expression between groups. One asterisk represents significant difference with P-value< 0.05; two asterisks represent significant difference with P-value< 0.01

Fig. 7

Increased abundance of the reads mapping to some banana genes was observed in inoculated samples as compared to control. a 21 nt mapped reads. b 24 nt mapped reads. Banana genes with annotated protein name were listed in the left column. The four columns on the right represent the abundance of mapped reads from respective treatment group. The darkest blue color represents the highest mapped read count, and the darkest pink color represents the lowest mapped read count. Exact read count is listed in Additional file 4: Table S3

Exploration of banana siRNAs in response to Foc infection

To investigate whether any banana endogenous siRNAs played roles in responding to Foc invasion, the processed reads of the sRNA libraries were mapped to the banana genome sequence (with intron and UTRs). The abundance of 21 nt, 22 nt and 24 nt perfectly mapped reads was then compared between inoculated and mock-inoculated samples using the SCRAM program (Fletcher et al. 2018). Figure 7 is a list of representative banana genes with various 21 nt and 24 nt mapped read abundance in the four treatment groups (Additional file 4: Table S3). Ma02_g06960 had the highest abundance of 21 nt mapped reads and both TR4 and R1 inoculated samples had approximately doubled mapped reads as compared to control samples (Additional file 1: Figure S2). Ma08_g19220 also showed a similar increase pattern in TR4-inoculated and R1-inoculated samples as compared to TR4-control and R1-control samples respectively. However, both genes encode uncharacterized proteins. For the 24 nt sRNAs, the abundance of reads mapped to Ma11_g15080, encoding alkaline/neutral invertase B, was increased upon both TR4 and R1 inoculation, with R1 inoculation causing a more significant increase.

Some genes showed different comparative patterns between R1 groups and TR4 groups, which might provide an insight into how Cavendish responded differently to the two Foc races. Among these genes, Ma06_g22630 (Histone-lysine N-methyltransferase EZ1) and Ma07_g22960 (Putative wall-associated receptor kinase 5) showed 1.6-fold and 1.7-fold increase, respectively, of 21 nt mapped reads in TR4-inoculated samples than control samples, while no difference was observed between R1-inoculated and control samples. Similarly, 24 nt sRNA abundance of Ma00_g03320 (Uncharacterized protein), Ma08_g16280 (Protein kinase domain containing protein, expressed) and Ma04_g33910 (Putative S-norcoclaurine synthase 2) showed an increase of 2.6 times, 3.7 times and 3.5 times in TR4-inoculated samples over the control samples respectively, whereas the sRNA abundance in R1-inoculated samples and R1-control samples were parallel. Ma10_g20970 (Uncharacterized protein) showed 12-fold increase in TR4_I samples than TR4_C samples, while less than one mapped RPMR was found in R1_I and R1_C. The mapping profiles of these genes showed only one or two hotspots of the antisense-mapped sRNAs. Instead of being randomly distributed along the whole gene, which could be caused by degradation, it was likely that these sRNAs were playing roles in response to the pathogen inoculation. Apart from Ma10_g20970 encoding an uncharacterized protein, the sRNA mapping hotspots were mostly located in the intron regions of the genes of interest and mapped to the antisense strand. It is possible that the abundant 21 nt sRNAs were intron-derived siRNAs that can bind to their host DNA sequence and cause DNA methylation in cis manner (Chen et al. 2011), and the 24 nt sRNAs were hc-siRNAs (heterochromatin-siRNAs). Further work is required to confirm this speculation.

Discussion

Small RNAs are important regulators of plant development and the plant’s immunity system. Since the genome sequence of banana has been available, studies on banana miRNAs have been accumulating (D'Hont et al. 2012; Davey et al. 2013; Chai et al. 2015; Ghag et al. 2015). A number of banana miRNA families have been identified. MiRNAs responding to drought stress and salt stress were investigated by Muthusamy et al. (2015) and Lee et al. (2015). However, research on banana miRNAs in response to biotic stress remains limited. This study investigated the small RNA profile of Cavendish banana inoculated with the Foc R1 and TR4 with an attempt to understand the roles of the small RNA regulators in banana upon Foc invasion.

By means of high throughput sequencing technology, this study has identified 112 previously-identified miRNAs belonging to 22 miRNA families from the sRNA libraries of both R1-inoculated and TR4-inoculated banana roots. Consistent with a previous finding (Ghag et al. 2015), miR166 showed the highest abundance among all the known miRNAs in banana. By challenging the plants with Foc, we found the miRNA profiles were altered by the inoculation of Foc TR4, whereas not as much change occurred during Foc R1 inoculation. In particular, miR156, miR159 and miR166 were greatly up-regulated by TR4 inoculation. Because Cavendish banana is resistant to R1 and susceptible to TR4, it is possible that these differentially expressed miRNAs has facilitated the infection of TR4. On the other hand, this could also be a defence mechanism of the host to limit or/and delay the damage caused by the virulent Foc race.

Among the very few studies on small RNAs in pathogen inoculated banana, a study on Banana streak Mysore virus (BSMyV)-infected banana also reported up-regulation of miR166 during viral infection (Sheeba et al. 2013). In a range of studies miR156 has been shown to be a critical regulator of phase change of plant growth by targeting the SQUAMOSA PROMOTER BINDING PROTEIN-LIKE (SPL) genes (Wu and Poethig 2006; Chuck et al. 2007; Wu et al. 2009; Wang et al. 2011). In Arabidopsis, miR156 was reported to be up-regulated by virus infection (Kasschau et al. 2003) as well as in tobacco (Bazzini et al. 2007), though down-regulation of miR156 was also revealed in response to the infection of Cronartium quercuum in pines (Lu et al. 2007) and powdery mildew in wheat (Xin et al. 2010). Cui et al. (2014) disclosed that the induction of miR156 improved adaptation and tolerance to environmental stresses including salt and drought stress by delaying the phase change of the plant. This could also be a possible explanation for the up-regulation of miR156 in this study. Both miR159 and miR166 target the transcription factors from the MYB and HD-Zip families respectively. Both of these gene families are involved in abscisic acid (ABA) signaling pathways which has linkage with plant defence pathways. It was previously speculated that a series of miRNAs including miR159 and miR166 are negative regulators of the abscisic acid (ABA) signaling pathways and subsequently positive regulator of the salicylic acid (SA)-mediated defence pathways (de Torres-Zabala et al. 2007; Zhang et al. 2011). In this scenario, it appears the up-regulation of miR159 and miR166 may be a mechanism to activate the SA-mediated defence pathway to facilitate host immunity upon TR4 infection. Considering that SA primarily mediates resistance to biotrophic and hemibiotrophic pathogens rather than necrothophs (Mengiste 2012), and that Fusarium oxysporum is considered as a hemibiotroph or a necrothoph (Thatcher et al. 2009; Ploetz 2015a; Dita et al. 2018), we think this speculation is still premature at this stage. More recently, miR159 and miR166 of cotton were found to be exported into V. dahliae and reduce the virulence during the host-pathogen interaction (Zhang et al. 2016). Considering F. oxysporum and V. dahliae are both fungal pathogens infecting plant vascular tissues, it is possible that the increased miR159 and miR166 in this study could also act as effectors and travel to F. oxysporum through an unknown defence pathway against aggressive Fusarium strains, although this requires further investigation.

In addition to miRNAs, we also found several 21 nt and 24 nt sRNAs that mapped to particular antisense loci of some banana genes and with different mapped read abundance between TR4-inoculated samples and control samples. It is possible that they are banana endogenous siRNAs involved in the plant defence mechanisms in response to the Foc infection, although further investigation is needed to confirm their derivation and function.

Conclusions

This study carried out the global small RNA profiling of Cavendish banana roots inoculated with two Foc races with different virulence upon banana plants. Our findings facilitate the understanding of the small RNA regulatory mechanisms during the banana-Foc interaction, which are useful resources for developing RNAi-based approaches to increase banana tolerance against aggressive Foc races.

Methods

Plant materials and fungal culture

Tissue culture plantlets of Cavendish cv. ‘Williams’ were initially transplanted into University of California (UC) soil mix (Baker 1957) and grown in a misted shade house for approximately 4 weeks before being transferred to a glasshouse for an additional 8 weeks until the plantlets had a pseudostem height of approximately 15 cm.

Monoconidial cultures of accessions of Foc, BRIP 62880 (VCG 01213/16) and BRIP42161 (VCG 0124) were regenerated from lyophilized mycelia stored at -80 °C onto ½ strength potato dextrose agar (PDA, Difco Laboratories) from the Northern Territory Department of Primary Industries and Brisbane Herbarium collection (Department of Agriculture and Forestry, Dutton Park, Australia), respectively. The cultures were incubated in dark conditions for 10 days at 24 °C.

Inoculation and sample collection

Plants were inoculated using the Millet carrier method (Smith et al. 2008), where every 50 g sterile Japanese millet was inoculated with one 5 × 5 mm² mycelial plug of Foc. Due to quarantine restrictions on the handling of TR4 in Australia, inoculation of plants with TR4 BRIP 62880 was conducted at the Coastal Plains Research Station, Northern Territory (NT), Australia, whilst R1 BRIP 42161 was conducted at the University of Queensland, Brisbane, Queensland, Australia. Two mock-inoculated groups were operated in NT and Brisbane, respectively, to eliminate any influence resulting from the location difference. Plants were transferred into 40 mm ANOVA pots containing a bottom layer of 150 g of soil mix topped with 50 g of millet and 150 g of soil mix on top of that. Plants were transferred into these pots and filled with the soil mix. Plants were grown in the green house for the remainder of the trial. Control plants were inoculated using the same method described above but using uninoculated millet. Plant material was harvested at 28 days post-inoculation. Roots of the plants were harvested by removing excess soil, flash frozen and stored at -80 °C for RNA extractions. Three biological replicates consisting of five plants per replicate were used for each treatment group. For each biological replicate, root tissue was pooled prior to RNA extractions.

RNA extraction and deep sequencing

Pooled root tissue was snap frozen in liquid nitrogen and homogenized with a mortar and pestle. Total RNA was extracted using the EpiBio MasterPure Plant RNA Purification kit (Epibio Inc.) using the manufacturer’s protocol. Small RNA library construction was conducted by the Australian Genome Research Facility (AGRF) using the Truseq Small RNA Library Preparation Kit (Illumina, California, USA). Sequencing (50 bp, single-end) was then conducted using the Illumina HiSeq2500 platform by AGRF.

Bioinformatic analysis

The FASTX-Toolkit package (http://hannonlab.cshl.edu/fastx_toolkit/) (Gordon and Hannon 2010) was used to remove the adapter sequences from 50 bp single-end reads. Only reads with lengths between 18 and 32 nt were retained for further analysis. The SCRAM software package (Fletcher et al. 2018) was used to align processed reads to reference sequences and compare alignments of different treatments, with aligned reads normalized as reads aligned per million reads. The genome sequence of Musa acuminata (Musa) DH-Pahang (accessed at http://banana-genome-hub.southgreen.fr/) (D'Hont et al. 2012; Martin et al. 2016) was used as the banana reference. The genome assemblies Foc1_1.0 and FO_II5_V1 (accessed at http://fungi.ensembl.org/index.html) (Guo et al. 2014) were used as Foc R1 and Foc TR4 references respectively.

miRNA identification and target gene prediction

The MiRDeep2 software package (Friedländer et al. 2012) was used to identify known and novel miRNAs. Known miRNAs were determined by aligning sequencing reads to published mature miRNA and miRNA precursor sequences available in miRBase (Release 21, http://www.mirbase.org/) (Griffiths-Jones et al. 2008), allowing zero mismatches in the first 18 nucleotides and up to two mismatches thereafter. Reads that were not mapped to miRNAs from miRBase were aligned to the banana genome (D'Hont et al. 2012; Martin et al. 2016). Those with significant randfold P-value (< 0.05) in the MiRDeep2 pipeline were considered as putative novel miRNAs.

Target genes of identified miRNAs were predicted using the psRNATarget program (https://plantgrn.noble.org/psRNATarget/) (Dai and Zhao 2011; Dai et al. 2018). Version 1 of all gene coding sequences of Musa acuminata (banana) released from the Banana Genome Hub was used as the target transcript library. The default settings of Scoring Scheme V2 was put as predicting criteria. The predicted target genes of identified known miRNAs and putative novel miRNAs were then submitted to the Blast2Go program (http://www.blast2go.com) (Conesa et al. 2005) for gene ontology (GO) analysis.

qRT-PCR validation of differentially expressed miRNAs and their predicted targets

Total RNA extracted for deep sequencing was also used for qPCR analysis. RNA samples were treated with Invitrogen’s Turbo RNase-free DNase I (Invitrogen, California, USA) to eliminate any genomic DNA remaining in the samples. A 500 ng quantity of total RNA was used to synthesize cDNA of miRNAs using the miScript Plant RT Kit (QIAGEN, Limberg, Netherlands). The following conditions were used for reverse transcription: 37 °C for 2 h, 95 °C for 5 min with a final hold step at 4 °C. Quantitative real time PCR was performed with Rotor-Gene Q platform (QIAGEN, Limberg, Netherlands) in a 20 μL reaction mixture, which consisted of 2 μL diluted cDNA, 2 μL of 10× miScript universal forward primer, 2 μL of 8 μM reverse primer, 10 μL of 2× QuantiTect SYBR Green PCR Master Mix (QIAGEN, Limberg, Netherlands) and 4 μL RNase-free water. The following conditions were used for PCR amplification: 95 °C for 15 min, and then 40 cycles of 94 °C for 15 s, 55 °C for 30 s and 70 °C for 30 s. To validate the expression of predicted target genes, 1 μg of total RNA was used to synthesize cDNA using SensiFAST cDNA Synthesis Kit (Bioline, London, UK). Conditions for reverse transcription were as follows: 25 °C for 10 min, 42 °C for 15 min, 85 °C for 5 min and hold at 4 °C. The 20 μL reaction mixture for qRT-PCR consisted of 2 μL of diluted cDNA, 1 μL of 10 μM forward primer, 1 μL of 10 μM reverse primer, 10 μL of 2× SensiFAST SYBR® No-ROX Mix (Bioline, London, UK) and 6 μL RNase-free water. The reaction was performed using the Rotor-Gene Q platform (QIAGEN, Limberg, Netherlands) following the conditions below: 95 °C for 2 min, then 40 cycles of 95 °C for 5 s, 60 °C for 15 s and 72 °C for 10 s.

A melt-curve was generated in the Rotor-Gene Q Series Software to analyze the specificity of the miRNA and target gene PCR products. Data obtained from this software was imported to the LinRegPCR Software (University of Amsterdam, The Netherlands) to determine the primer efficiency. Ct is the cycle threshold of the PCR reaction. Primer efficiency was averaged from all reactions using the same primer. Ct value was averaged from two technical replicates. miR399 was used as an endogenous reference for miRNA expression analysis. L2 was used as an endogenous reference gene for target gene expression analysis. Relative expression of miRNAs and target genes was calculated using the following formula:

$$ \mathrm{Relative}\ \mathrm{Expression}\ \left(\mathrm{RE}\right)={{\mathrm{PE}}_{\mathrm{target}}}^{\left(-\mathrm{Ct}\ \mathrm{target}\right)}/{{\mathrm{PE}}_{\mathrm{reference}}}^{\left(-\mathrm{Ct}\ \mathrm{reference}\right)} $$

The RE of each treatment was then compared. One-way analysis of variation (ANOVA) was used to determine the significance of differences. Primers were designed using Primer 3.0 plus online system listed in Additional file 4: Table S4.