A dynamic degradome landscape on miRNAs and their predicted targets in sugarcane caused by Sporisorium scitamineum stress

Background Sugarcane smut is a fungal disease caused by Sporisorium scitamineum. Cultivation of smut-resistant sugarcane varieties is the most effective way to control this disease. The interaction between sugarcane and S. scitamineum is a complex network system. However, to date, there is no report on the identification of microRNA (miRNA) target genes of sugarcane in response to smut pathogen infection by degradome technology. Results TaqMan qRT-PCR detection and enzyme activity determination showed that S. scitamineum rapidly proliferated and incurred significant enzyme activity changes in the reactive oxygen species metabolic pathway and phenylpropanoid metabolic pathway at 2 d and 5 d after inoculation, which was the best time points to study target gene degradation during sugarcane and S. scitamineum interaction. A total of 122.33 Mb of raw data was obtained from degradome sequencing analysis of YC05–179 (smut-resistant) and ROC22 (smut-susceptible) after inoculation. The Q30 of each sample was > 93%, and the sequence used for degradation site analysis exactly matched the sugarcane reference sequence. A total of 309 target genes were predicted in sugarcane, corresponding to 97 known miRNAs and 112 novel miRNAs, and 337 degradation sites, suggesting that miRNAs can efficiently direct cleavage at multiple sites in the predicted target mRNAs. Gene Ontology (GO) annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis indicated that the predicted target genes were involved in various regulatory processes, such as signal transduction mechanisms, inorganic ion transport and metabolism, defense mechanisms, translation, posttranslational modifications, energy production and conversion, and glycerolipid metabolism. qRT-PCR analysis of the expression level of 13 predicted target genes and their corresponding miRNAs revealed that there was no obvious negative regulatory relationship between miRNAs and their target genes. In addition, a number of putative resistance-related target genes regulated by miRNA-mediated cleavage were accumulated in sugarcane during S. scitamineum infection, suggesting that feedback regulation of miRNAs may be involved in the response of sugarcane to S. scitamineum infection. Conclusions This study elucidates the underlying response of sugarcane to S. scitamineum infection, and also provides a resource for miRNAs and their predicted target genes for smut resistance improvement in sugarcane. Electronic supplementary material The online version of this article (10.1186/s12864-018-5400-8) contains supplementary material, which is available to authorized users.


Background
Sugarcane smut, caused by Sporisorium scitamineum, is a worldwide airborne fungal disease that affects sugarcane production [1]. The disease leads to severe losses in cane yield and reduces sucrose content and quality, thereby preventing further development of the sugarcane industry in China. A typical symptom of sugarcane smut involves cane tips of infected plants growing a black whip that points downwards and curls inwards around 120 d of planting. Simultaneously, mycelia invade the cane buds, and chlamydospores fall to the soil, thereby also infect sugarcane in the next growing season [2]. Compared to normal plants, the main stems of plants infected by S. scitamineum are small, the cane leaves are slender and light green in color, and tillers usually also grow black whips that lead to a sharp decline in sugarcane production [3]. Cultivating sugarcane varieties with excellent smut-resistance is the most economical and effective way to control the disease [4].
Due to the incomplete genome sequencing of Saccharum spp. hybrid, related genomics research is limited, thus hindering the progress of molecular improvement of sugarcane varieties. Current researches on the molecular mechanism of interaction between sugarcane and S. scitamineum mainly focus on genomics of smut pathogen [5][6][7], sugarcane molecular marker-assisted selection [8,9], transcriptome [10,11] and proteomic [12,13] analysis, and resistance-related gene mining [10,11]. The understanding of how miRNAs regulate the expression of their target genes in response to S. scitamineum infection is limited. The only earlier investigation relating to this matter was performed by our research group, which involved the identification of differentially expressed miR-NAs in sugarcane challenged with S. scitamineum by using high-throughput sequencing [14].
MicroRNAs (miRNAs) are single-stranded, non-coding RNA molecules of approximately 21-24 nt in length in vivo [15]. miRNAs were first reported in Caenorhabditis elegans by Lee et al. [16]. Plant miRNAs were first obtained from a small Arabidopsis thaliana library by Reinhart et al. [17,18]. miRNAs are encoded by endogenous miRNA genes and negatively regulate gene expression primarily through the degradation of target mRNAs during post-transcriptional gene silencing (PTGS) [19,20]. miRNAs are important regulators of organisms and are commonly involved in the response of plants to biotic stress [21,22]. Previous studies have shown that plant miRNAs enhance their resistance to pathogen infection by regulating the expression of key disease-resistance genes [22][23][24][25]. During the late stage of wheat stripe rust infection, wheat miR408 was downregulated, and the expression of target gene chemocyanin-like protein (CLP1) was induced, which in turn inhibited the growth of mycelia in leaf [22]. After inoculation of tomato stalks with cucumber mosaic virus (CMV) and tomato aspermy virus (TAV), miR156 was accumulated, resulting in hollow and fibrotic stalks [23]. In addition, tomato plants overexpressing miR156 showed phenotypic symptoms that are similar to that of pathogen infections, indicating that miR156 regulates the interaction of tomato with pathogens [24]. miRNAs were differentially expressed between maize varieties with high resistance and susceptibility to Rhizoctonia solani [25]. The expression patterns of zea-miR168a, miR-2, and miR-6 were generally upregulated in diseaseresistant plants, while miR-3 was only upregulated in pathogen-infected sites [25]. The expression level of miRNA (except for miR-4 and miR-5) in Zea mays resistant variety was significantly higher than that in susceptible one, which induced host defense mechanisms to resist pathogen infection [25].
In organisms, the most important regulatory mechanism of miRNAs involves the miRNA-directed target cleavage that regulate their life cycle [26,27]. Degradome sequencing, also called parallel analysis of RNA ends (PARE), is a high-throughput sequencing technique [28,29]. This method can determine pairing information with miRNA through high-throughput deep sequencing of miRNA-mediated target gene cleavage degradation fragments, and the target genes of miRNAs are screened to determine how they regulate plant life activities in specific environments [28,29]. Degradome sequencing has been extensively applied to the identification of miRNA target genes in crops such as wheat [28], maize [30], peanut [31], and other crops. However, due to different calculation methods and screening criteria, it is difficult to avoid false positive in degradome sequencing data. Therefore, further validation of the authenticity and reliability of the miRNAs and their target genes is needed. Previous studies have shown that qRT-PCR can effectively detect the abundance of target genes in samples and verify the reliability of sequencing results in the degradome [32]. In addition, the qRT-PCR method can also be used to detect the expression level of the corresponding miRNAs, and then verify the interaction mode between miRNAs and their target genes [33].
In China, the main sugarcane variety grown during the past 20 years is ROC22, which is a Saccharum spp. hybrid from ROC5 × 69-463 (high-sugar line) and encompasses approximately 60% of the total sugarcane cultivated area. ROC22 is susceptible to S. scitamineum and results in a poor ratoon performance. In this study, an intergeneric BC2 hybrid with smut resistant character named YC05-179 (YC01-134 × ROC20) and the smut-susceptible variety ROC22 were used as experimental materials. First, the TaqMan qRT-PCR technology and physiological enzyme activity determination were performed to analyze the proliferation of S. scitamineum and activity changes of key enzymes involved in reactive oxygen species metabolic and phenylpropanoid metabolic pathways at various time points in different sugarcane genotypes infected by S. scitamineum. And the critical time points during interaction between sugarcane and S. scitamineum were determined. Second, the RNA of YC05-179 and ROC22 at different stages of infection was used for degradome sequencing analysis to screen miRNAs and their predicted target genes. Third, the expression level of partial obtained miRNAs and their corresponding predicted target genes was verified by qRT-PCR. The purpose of this study is to understand the miRNA-mediated molecular mechanisms in sugarcane response to S. scitamineum stress.

Plant materials and inoculation
The tested sugarcane varieties were YC05-179 (smut resistant) and ROC22 (smut susceptible). The pathogen strain was mixed spores of S. scitamineum collected from different sugarcane varieties and different planting areas in the field. They were all provided by the Key Laboratory of Sugarcane Biology and Genetic Breeding, Ministry of Agriculture, Fujian Agriculture and Forestry University (Fuzhou, Fujian, China). The plants were treated with S. scitamineum using the method of Su et al. [14]. The consistent and robust cane stems were selected, which were cut into single-bud stem, and then soaked in running water for 2 d. The cane stems were then placed on a tray in a 32°C incubator (65% relative humidity) with 12 h light/12 h dark conditions for germination. The water was replaced once in the morning and in the evening. When the cane buds had grown to about 2 cm in height, the treatment group were subjected to acupuncture inoculation with 5 × 10 6 spores/mL S. scitamineum spore suspension (0.01% Tween-20) in buds, and the control group underwent acupuncture inoculation with sterile water (0.01% Tween-20). The cane stems were then placed in a 28°C incubator (65% relative humidity) with 12 h light/12 h dark conditions. The sugarcane bud tissues were collected at 0, 1, 2, 3, 5 and 7 d after inoculation, fixed in liquid nitrogen, and then stored at − 80°C until analysis.
Sugarcane genomic DNA and total RNA extraction and quality testing Three sugarcane buds were taken from each sample. Genomic DNA of the samples was extracted according to the modified CTAB method of Yao et al. [34]. Total RNA was extracted using TRIzol (Invitrogen, Carlsbad, CA, USA). The integrity of total RNA was detected by 1.5% agarose gel electrophoresis. The concentration and purity of genomic DNA and total RNA were assessed using a Nano-Drop (Thermo Fisher, USA) and Agilent 2100 (Agilent Technologies, Palo Alto, CA, USA) systems.

Quantification of smut pathogen in sugarcane
The quantification of S. scitamineum in YC05-179 and ROC22 was detected using TaqMan qRT-PCR technology that was developed by Su et al. [35]. The detection primers were bEQ-F: 5'-TGAAAGTTCTCATGCAAGCC-3′ and bEQ-R: 5'-TGAGAGGTCGATTGAGGTTG-3′, and the TaqMan probe was 5'-FAM-TGCTCGACGCCAATTCG GAG-TAMRA-3′. A standard curve was established by a 10-fold gradient dilution with recombinant plasmid DNA containing the bE gene (b East mating type gene, GenBank Accession No. U61290.1). The concentration of DNA was 500 ng/μL. The amplification program was 50°C for 2 min; followed by 40 cycles of 95°C for 10 min, 95°C for 15 s, 60°C for 1 min; and a single-point fluorescence detection at 60°C. The blank control replaces the template DNA with an equal volume of ddH 2 O. The negative control is the DNA of pathogenfree FN41 (a newly released Saccharum spp. hybrid from Yuetang91-976 × ROC20 in China) 4-month-old plantlets and the positive control is the genomic DNA of S. scitamineum. Each run of TaqMan qRT-PCR contained three replicates. A standard curve was drawn using Microsoft Excel and GraphPad Prism softwares, and the quantification of S. scitamineum in each sample was calculated.
Determination the activity of key enzymes involved in reactive oxygen species metabolic and phenylpropanoid metabolic pathways At 0, 1, 2, 3, 5 and 7 d after inoculation, 0.1 mol/L borate buffer (pH 8.7) was used for the extraction of crude enzyme solution in cane bud tissues [36]. All samples (five buds per sample) were weighed and homogenized with ice-cold borate buffer at the ratio of 1 g/10 mL, as well as with a small amount of quartz sand and polyvinylpyrrolidone (PVP). The supernatant was centrifuged at 5000 rpm/min for 15 min at 4°C. The final supernatant was used as the crude enzyme solution for the determination of peroxidase (POD), superoxide dismutase (SOD), catalase (CAT), polyphenol oxidase (PPO), phenylalanine ammonia lyase (PAL), and tyrosine ammonia lyase (TAL) activities. The guaiacol method was used to determine POD enzyme activity [37]. SOD activity was measured using nitroblue tetrazolium (NBT) photoreduction [38]. The activities of CAT, PPO, PAL, and TAL were determined according to the methods of Beers and Sizer [39], Galeazzi et al. [40], Green et al. [41], and Kofalvi and Nassuth [42], respectively. Three biological replicates were prepared for each treatment. The data were analyzed using Microsoft Excel and SPSS softwares. GraphPad Prism software was used to produce figures. Duncan's method was used for statistical analysis.

Degradome library construction, sequencing and bioinformatics analysis
The samples for degradome sequencing were YC05-179 and ROC22 cane buds inoculated with sterile water at 0 d (control group) and with S. scitamineum at 2 and 5 d (treated group). These samples were named Y0, Y2, Y5, R0, R2, and R5, respectively. The quality of total RNA met the OD 260/280 of 1.8-2.2, with normal absorption peak at 260 nm, and 28S/18S ratio ≥ 1.5. Library construction, degradome sequencing, and data analysis were commissioned by Beijing BioMed Biotech Co., Ltd. (Beijing, China) using the Illumina HiSeq™ 2500 [43]. Raw reads were generated from degradome sequencing. Then the reads with adapters, the low-quality reads with mass values below 30 and bases in excess of 20%, the reads with unknown base N contents greater than or equal to 10%, and the reads less than 18 nt in length were filtered out to eventually obtain clean reads and cluster reads [31]. Non-coding RNAs (rRNA, scRNA, snoRNA, snRNA, and tRNA) were removed by aligning clean reads and cluster reads with the Rfam database [44]. Because whole genome sequencing of Saccharum spp. hybrid has not been completed to date, we mapped the remaining sequences using the bowtie software and targetfinder software to align with the sugarcane reference sequences (sugarcane transcriptome under S. scitamineum stress [11], GSS database, and EST database in NCBI) and the known miRNAs from miRBase 21.0 (http://www.mirbase.org) or miRNAs identified in our previous study [14] to predict miRNA target genes. When the score of miRNA that matches the mRNA was less than or equal to 7, the transcript sequence was considered as the miRNA target gene [45,46].

miRNA targets prediction and identification
Based on the depth statistics of mRNA target genes and the abundance of transcripts, the target genes were grouped into five categories, namely, 0, 1, 2, 3, and 4 [45,47]. Category 0 indicated that the position had a depth > 1, an abundance equal to the maximum of the transcript abundance, and the transcript had only one maximum value. Category 1 indicated that the position had a depth > 1, an abundance equal to the maximum value of the transcript abundance, and the transcript had two or more maxima. Category 2 represented the depth of the position was > 1 and the abundance was less than the maximum but higher than the mean of the transcript abundance. Category 3 represented that the depth was > 1 and the abundance was less than or equal to the mean of the transcript abundance. Category 4 represented the depth of the position equal to 1.
The degradation site of predicted target gene was analyzed by Cleaveland software at p-value < 0.05. The screened predicted target gene sequences with degradation sites were compared to the COG, GO, KEGG, NR, NT, and Swiss-Port databases to obtain the predicted target gene annotation information. The expression level of the predicted target gene was calculated using the fragments per kilobase of transcript per million mapped reads (FPKM) method which eliminated the influence of gene length and sequencing difference in high-throughput sequencing [48]. The treatment group and the control group were compared to analyze the differential expression of the predicted target genes, and the ratio of expression level was expressed as the fold-change [49]. A fold-change ≥2 and a false discovery rate (FDR) < 0.01 were used as screening criteria for differentially expressed predicted target genes. The FDR was obtained by correcting the difference in the significance of the p-value [50]. The p-value was corrected using the Benjamini-Hochberg calibration method, and the FDR was used as a screening index to ensure the quality of differentially expressed genes. DY2 and DY5 represent the differentially expressed predicted target genes at 2 and 5 d after YC05-179 was inoculated with S. scitamineum, whereas DR2 and DR5 represent the differentially expressed predicted target genes at 2 and 5 d after ROC22 was inoculated with S. scitamineum. DY2, DY5, DR2, and DR5 were respectively analyzed by COG, GO, and KEGG to investigate functional and related metabolic pathways of the differentially expressed predicted target genes. The continuously and non-continuously differentially expressed target genes predicted in YC05-179 and ROC22 at different time points were analyzed to reveal potential genes that were related to smut resistance.
qRT-PCR validation of the expression level of miRNAs and their predicted targets MiR168a-5p, miR5293, miR160a, nov-miR132, nov-mir-143, nov-mir-63, nov-mir-18, miR5368, nov-mir-10, miR858b, nov-mir-97, miR162a, miR529-3p, and their partial predicted target genes (Table 1) were verified by qRT-PCR. The qRT-PCR primers of the predicted target genes (Additional file 1: Table S1) were designed using Beacon Designer software. Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) was used as the internal reference gene [51]. Total RNA of YC05-179 and ROC22 inoculated with sterile water for 0 d and with S. scitamineum for 2 and 5 d were digested with RNase-Free DNase (Promega, USA) to remove DNA contamination, followed by reverse transcription to generate first-strand cDNA using a PrimeScript® RT reagent kit (Perfect Real Time) (TaKaRa, Dalian, China). The qRT-PCR reaction system was prepared using the SYBR Green dye method following the instructions of the FastStart Universal SYBR Green PCR Master (ROX) kit (Roche, Shanghai, China). qRT-PCR amplification was performed on an ABI 7500 instrument at 50°C for 2 min, followed by 40 cycles of 95°C for 15 s and 60°C for 1 min. Three replicates were used for each sample. The blank controls were used to replace the cDNA template with sterile water. Stemloop method was used to detect the expression level of candidate miRNA [52]. According to the method of Varkonyi-Gasic et al. [33], miRNA stem-loop primers (RT primer) and upstream primers (Additional file 1: Table S2) were designed. The downstream primer was a universal primer for anchoring the stem-loop region. GAPDH was used as the internal reference gene [53]. Total RNA was reverse-transcribed using an Applied Biosystems® TaqMan® MicroRNA Reverse Transcription kit (Invitrogen, Carlsbad, CA, USA) with a 15-μL reverse-transcription system. In the process of reverse transcription, when using cDNA template as internal control, the RT primers were replaced by random primers. The miRNA qRT-PCR reaction system and procedure were the same as that for the predicted target gene. qRT-PCR data were calculated using the 2 −ΔΔCT method [54].

Results
Smut pathogen proliferation and changes in the activity of key enzymes involved in reactive oxygen species metabolic and phenylpropanoid metabolic pathways during the early stage of infection The results of TaqMan qRT-PCR ( Fig. 1) showed that the Ct values of smut-resistant genotype (YC05-179) and -susceptible genotype (ROC22) inoculated with sterile water, negative control, and blank control were all higher than 37, indicating the absence of S. scitamineum, whereas the Ct value of samples inoculated with S. scitamineum was between 27 and 33. The quantification of S. scitamineum in YC05-179 and ROC22 increased with inoculation time, which were 323,995.15 ± 53,563.55 copies/μL-2,935,184.09 ± 36,789.33 copies/μL and 340,733.51 ± 29,137.42 copies/μL-7,525,544.93 ± 358,488.58 copies/μL, respectively. At 0 d after inoculation, the quantification of S. scitamineum was similar in YC05-179 and ROC22. The dynamic increase in DNA copy number of the S. scitamineum in ROC22 was more distinct at 1, 2 and 3 d after Cleavage site, nucleotide number from 5′ end of cDNA; Category, the "category" of this cleaveage site Fig. 1 The amount of smut pathogen in YC05-179 and ROC22 inoculated with Sporisorium scitamineum by TaqMan qRT-PCR analysis. The quantification of smut pathogen were calculated with the equation of the linear regression line. All data points were means ± standard error (n = 3). Different lowercase letters indicated a significant difference between resistant and susceptible genotypes, as determined by the Duncan's new multiple range test (p < 0.05). YC05-179, smut-resistant genotype; ROC22, smut-susceptible genotype inoculation, which was 3.27-, 4.31-, and 6.82-fold respectively of that at 0 d after inoculation. The quantification of S. scitamineum in YC05-179 was increased at 1 d but decreased at 2 d, which was 1.34-fold and 0.57-fold that at 0 d, respectively. At 5 and 7 d, the content of S. scitamineum in YC05-179 showed minimal change, which was 6.96-and 9.06-fold that at 0 d, whereas the content of S. scitamineum in ROC22 continued to increase and peaked at 7 d, which was 22.09-fold that at 0 d after inoculation. Figure 2 showed the activity changes of six key enzymes involved in reactive oxygen species metabolic and phenylpropanoid metabolic pathways in two sugarcane genotypes post inoculation with S. scitamineum. After inoculation, the activity of POD in YC05-179 and ROC22 was all increased. It peaked at 2, 3 and 5 d in YC05-179, whereas at 3 and 7 d in ROC22. The activity of SOD in YC05-179 was decreased at 1 d but gradually increased and reached a peak value at 7 d. In ROC22, the activity of SOD was stable at 1 d, and decreased at 2 and 7 d, whereas reached a peak value at 5 d. There was an opposite expression pattern in the activity of CAT between YC05-179 and ROC22. After inoculation with S. scitamineum, the activity of CAT in YC05-179 was significantly increased and peaked at 1 d, follow by 2, 5 and 7 d, whereas that in ROC22 was decreased at 1, 2 and 5 d. PAL and TAL, which are upstream of the phenylpropanoid metabolic pathway, inhibit lignin biosynthesis when their activities are inhibited [55]. Inoculation of YC05-179 with S. scitamineum resulted in a stronger activity of PAL at 3 d, followed by 2 and 5 d. The activity of PAL in ROC22 was increased at 3 and 7 d but decreased at 2 and 5 d. In YC05-179, the activity of TAL reached a peak value at 1 and 5 d, followed by 3 d, and remained stable at 2 and 7 d. The activity of TAL in ROC22 were all increased after inoculation and peaked at 1, 2, 5 and 7 d. PPO, which catalyzes the phenols to quinones, was assumed to be involved in plant defense against pathogens [56]. After inoculation, the activity of PPO in YC05-179 was significantly increased and peaked at 1, 2 and 7 d, followed by 3 and 5 d. The activity of PPO in ROC22 was significantly increased at 2 and 5 d but remained stable at the other time points. In summary, the activities of POD, SOD, CAT, PPO, PAL and TAL Fig. 2 Variations of the activity of key enzymes involved in reactive oxygen species metabolic and phenylpropanoid metabolic pathways during the early stage of Sporisorium scitamineum infection. All data points were means ± standard error (n = 3). Different lowercase letters indicated a significant difference, as determined by the Duncan's new multiple range test (p < 0.05). YC05-179, smut-resistant genotype; ROC22, smut-susceptible genotype; POD, peroxidase; SOD, superoxide dismutase; CAT, catalase; PPO, polyphenol oxidase; PAL, phenylalnine ammonialyase; TAL, tyrosine ammonia-lyase in YC05-179 were generally higher than those in the control (0 d) and ROC22 plants after inoculation with S. scitamineum. At 2 and 5 d, the metabolic levels of activated oxygen and phenylpropanoid in YC05-179 were stronger than those in ROC22. Furthermore, there was a significant difference in smut pathogen content between YC05-179 and ROC22 after inoculation ( Fig. 1). At the early stage of 2 d, the amount of S. scitamineum was decreased in YC05-179 but increased in ROC22 (Fig. 1). Therefore, we used 2 and 5 d after inoculation with S. scitamineum as the best sampling time points for degradome sequencing.

Degradome library construction and data summary
Six cane-bud samples, including YC05-179 and ROC22 inoculated with sterilized water for 0 d and those inoculated with S. scitamineum for 2 and 5 d, were sequenced. Each sample was generated with no less than 10 M reads, and a total of 122.33 M raw reads in six samples was gained. After data evaluation, clean reads and cluster reads were obtained with a length of 47 nt.  (Table 2). After comparing the clean reads and cluster tags with Rfam database to exclude rRNAs, tRNAs, snoRNAs, and snRNA (except for miR-NAs), the remaining sequences were aligned with the reference sequence of sugarcane to obtain the fully mapped data, which consisted of 2,748,695 (Y0), 3,333,316 (Y2), 3,446,802 (Y5), 2,490,391 (R0), 3,151,567 (R2), and 3,001,075 (R5) reads ( Table 2) that were then used in the analysis of degradation sites.

Degradation site specificity and diversity
According to sequence homology, each miRNA can simultaneously target two or more target genes belonging to the same type or having similar conserved domains, and a target gene can also be cleaved by multiple miRNAs [57,58]. In this study, the degradome sequencing results showed that the predicted target gene could be cleaved by different miRNAs at a specific cleavage site. For example, Sugarca-ne_Unigene_BMK.74449 could be cleaved simultaneously by miR165a, miR166a, and miR166g-3p at position 4321 (Figs. 3a, b, c). The predicted target gene could also be cleaved by different miRNAs at different cleavage sites, e.g., Sugarcane_Unigene_BMK.61043 could be cleaved by nov-mir-84 and nov-mir-41 at positions 324 and 785, respectively (Figs. 3d, e). In addition, multiple target genes could be cleaved by the same miRNA. For example, Sug-arcane_Unigene_BMK.74449 and Sugarcane_Unigen-e_BMK.72615 could be cleaved by miR165a (Figs. 3a, f).

Target gene identification, annotation, and classification
Based on the miRNA database [14] and Unigene database [11] of sugarcane after smut pathogen infection, 2922 miRNA-mRNA pairs were screened from six libraries using targetfinder software. A total of 337 degradation sites were detected by Cleaveland software, corresponding to 219 miRNAs (97 known miRNAs and 112 new miRNAs) and 309 predicted target mRNAs (Additional file 2: Table S3). The predicted target mRNAs were all mainly classified as Category 0 in all six libraries, without Category 4 (Fig. 4). We found that the predicted target genes cleaved by miRNAs, such as squamosa promoter-binding-like protein (SPL), no apical meristem (NAM), v-myb avian myeloblastosis viral oncogene homolog (MYB), auxin response factor, extensin-like protein, somatic embryogenesis receptor kinase (SERK), CAT, ethylene-insensitive 3-like 3 protein (EIL3), and miRNA precursor encoding genes, may be involved in various life-controlling processes of sugarcane.
qRT-PCR validation of the expression level of miRNAs and their corresponding predicted target genes qRT-PCR analysis of miRNAs and their corresponding target genes will not only verify the accuracy of our degradome sequencing results, but also determine the miRNA-mediated regulatory role of miRNAs and their predicted target genes in sugarcane responses to S. scitamineum stress. The results showed that the expression level of 13 predicted target genes in qRT-PCR analysis and degradome sequencing was similar, but not completely consistent (Fig. 6). There was also a certain deviation in the differences of gene expression fold. The predicted target genes and their corresponding miRNAs were expressed in opposite patterns in at least one sugarcane variety. However, the expression patterns of different predicted target genes and their corresponding miRNAs in two sugarcane varieties (YC05-179 and ROC22) varied as follows: (i) Ubiquitin carboxyl-terminal hydrolase isozyme L5-like (UCH-L5), protein argonaute 1B (AGO 1B), and ARF8 followed a negative miRNA-mediated regulatory mode in two sugarcane varieties. The expression patterns of predicted target genes UCH-L5, AGO 1B, and ARF8 were opposite to that of their corresponding miRNAs, i.e., UCH-L5, AGO 1B, and ARF8 (slightly decreased but not significant) were downregulated, whereas the corresponding miR529-3p, miR168a-5p, and miR160a were upregulated in YC05-179, and the expression pattern in ROC22 was the opposite (Fig. 6A). (ii) Auxin-induced protein (AIP), cinnamoyl-CoA reductase (CCR), and S-adenosylmethionine decarboxylase (SAMDC) fit the negative miRNA-mediated regulatory mode only in ROC22. While the expression trends of predicted target genes and their corresponding miRNA in YC05-179 were consistent, i.e., AIP, CCR, and SAMDC were upregulated, the corresponding miR5293, nov-mir-132, and miR162a were downregulated in ROC22, whereas these predicted target genes and their corresponding miR-NAs were upregulated in YC05-179 (Fig. 6B). (iii)EIL3 and HIR1 were upregulated and followed the negative miRNA-mediated regulatory mode only in ROC22. Whereas the expression level of EIL3 and HIR1 in YC05-179 showed little change, the expression amount of their corresponding nov-mir-143 and miR5368 varied greatly and all was upregulated (Fig. 6C). (iv) Growth-regulating factor 8 (GRF8), glycerol kinase (GK), PP2C, mildew resistance locus o (MLO), and Myb-related protein Hv33 (MYB2) fit the negative miRNA-mediated regulatory mode only in YC05-179 (Fig. 6D). GRF8 was downregulated in both YC05-179 and ROC22, and the corresponding nov-mir-18, nov-mir-66, and miR396e-5p were all upregulated in YC05-179, but upregulated in ROC22 only at 5 d. GK was upregulated at 2 d and downregulated at 5 d after inoculation with S. scitamineum in YC05-179, whereas the corresponding nov-mir-63 showed the opposite pattern. The expression level of GK in ROC22 decreased with prolongation of inoculation time.
The expression pattern of nov-mir-63 was the same as that in YC05-179. PP2C was downregulated and the corresponding nov-mir-97 was upregulated in YC05-179 at 2 d, but both of them were stable in ROC22. MLO was downregulated in YC05-179 and ROC22 at 2 d, and both of them recovered at The above results indicated the complexity of miRNA regulation of S. scitamineum that infects sugarcane. There was no obvious negative linear regulation between the miRNAs and their predicted target genes, and there were also significant differences in the regulation of expression patterns between different genotypes of sugarcane. The expression level of AGO 1B, UCH-L5, AIP, CCR, EIL3, HIR1, and SAMDC significantly changed at 2 d after inoculation, indicating that they responded earlier to S. scitamineum infection. ARF and GRF are associated with plant growth and development [59][60][61]. The significantly decreased expression level of GRF8 and the slightly decreased expression level of ARF8 after inoculation suggests that the infection of S. scitamineum might inhibit the growth of sugarcane to a certain extent. MLO is a calcium-binding protein whose expression is negatively correlated with plant disease resistance [62,63]. The MLO gene in sugarcane was downregulated at 2 d after inoculation with S. scitamineum compared to ROC22, with a significant decrease in expression level in YC05-179, suggesting that MLO responded earlier to smut pathogen infection and that the intensity of the response was greater in the resistant variety than the susceptible one. MYB plays a more important role in the defense response of plants to stress [64]. MYB2 was upregulated in YC05-179 at 5 d after inoculation with S. scitamineum, indicating a delay in response. GK is a rate-limiting enzyme in the metabolic pathway of glycerol, and its expression level is closely related to the innate immune response in plants [65,66]. The expression level of GK in the resistant variety YC05-179 was increased and higher than that in the susceptible variety ROC22 after infection with S. scitamineum, which may corroborate the positive correlation between the expression level of GK and smut resistance in sugarcane varieties.

Discussion
Plant miRNAs mainly regulate target gene expression by mediating target mRNAs cleavage or repressing gene translation during plant development [26,27]. Degradome sequencing is a high-throughput method to identify miRNAs and their predicted target genes at a certain developmental stage or under specific stress in plants, which in turn may reveal the target genes of miRNAs that are related to plant development or response to stress [28,29]. This technique has now been successfully applied to various plants such as Populus tomentosa [57], cotton [58], rice [67], and peanut [31]. In the present study, the expression of the predicted target genes regulated by miRNA-mediated cleavage in smut-resistant and -susceptible varieties of sugarcane under the infection of S. scitamineum was analyzed using degradome sequencing. The results showed that an initial data amount of 122.33 M was obtained on six sugarcane samples. The data amount of each sample was not less than 10 M and the Q30 was > 93%, suggesting that the sequencing quality was relatively high. In addition, the sequence of degraded fragments obtained by degradome sequencing was the same as that of transcriptome sequencing [11], which demonstrates the accuracy of degradome sequencing results.
The target genes were annotated with GO, KEGG, NR, NT, Swiss-Prot and COG databases to obtain their basic information, functional classification, and involved metabolic pathways [71,72]. In this study, predicted target genes involved in various regulatory processes of life activity such as signal transduction, ion transport, translation and posttranslational modification, energy production and transduction, and metabolism of glycerides. In addition, a miRNA precursor, namely, Sugarcane_Unigene_BMK.40037 (miR171e-3 precursor miRNA), cleaved by nov-mir-219, was found in the YC05-179 at 2 d after inoculation with S. scitamineum. The interaction process between sugarcane and S. scitamineum is regulated by a multi-gene network system. Correspondingly, the process of sugarcane in response to S. scitamineum infection is involved in the regulation of multiple metabolic pathways [73,74]. The differentially expressed predicted target genes in YC05-179 and ROC22 were basically the same in the classification of targets, mainly playing a catalytic and binding role in response to stimulation, signaling pathway and immune process. KEGG analysis showed that the differentially expressed predicted target genes involved in plant hormone signal transduction, plant-pathogen interactions, oxidative phosphorylation, and other disease-related metabolic pathways, in which plant-pathogen interaction pathways appeared only in the differentially expressed predicted target genes of YC05-179, namely Sugarcane_Unigene_BMK.75694. Sugarcane_Unigene_ BMK.75694 is cleaved by nov-mir-63, encodes glycerol kinase, which is involved in energy production and transduction. After inoculating with S. scitamineum, the expression of Sugarcane_Unigene_BMK.75694 in YC05-179 was increased with the prolongation of inoculation time, which was opposite to that in ROC22. Sugarcane_Unigene_BMK.62557 (catalase), which participated in the peroxisomal pathway and was regulated by miR858b and nov-mir-88, was only differentially expressed in YC05-179 and ROC22 that were inoculated for 5 d. The expression level of Sugarcane_Unigene_BMK.62557 in YC05-179 was increased with the elongation of inoculation, but was opposite to that in ROC22. This does not agree with our finding that catalase activity was higher at 2 d compared to that at 5 d after inoculation (Fig. 2), which may be because the expression level of other catalase family members had greater changes than that of Sugarcane_Unigene_BMK.62557, thus fitting the trend of changes in the activity of catalase. In addition, YC05-179 and ROC22 also responded to S. scitamineum infection through their own differential expression of genes. Therefore, link the changes in transcript levels of predicted target genes to observe biology of infection will provide a better insight on sugarcane in response to smut pathogen attack.
Pathogen invasion triggers various plant immune responses [75,76]. When infected by pathogens, the responses of plant genes may differ such as early or late stress responses or high or low expression level [11]. miRNAs mainly regulate target gene expression by cleavage, which in turn influences the life processes of plants [26,27]. At the same time, the response of plants to pathogen infection involves multiple metabolic pathways and several genes [10,13]. This study showed that inoculation with S. scitamineum induces changes in multiple resistance-related metabolic pathways in YC05-179 and ROC22 as following:

Lignin biosynthesis pathway
Lignin biosynthesis is a branch of the phenylpropanoid metabolic pathway that plays an important role in plant disease resistance [77]. CCR is a key enzyme that catalyzes lignin biosynthesis and promotes the formation of lignin [78,79]. CCR first catalyzes the formation of corresponding aldehydes by coumaryl-CoA, feruloyl-CoA and sinapoyl-CoA, which are then catalyzed by cinnamyl alcohol dehydrogenase (CAD) to form lignin monomers [78,79]. Previous studies have found that the reduced expression of CCR can cause a significant decrease in lignin content [80]. In addition, CCR gene expression in plants is upregulated and the accumulation of lignin monomer is increased during fungal infections [81,82]. PPO has often been suggested to participate in plant defense against pathogens and pests by promoting the formation of lignin and quinones [56,83]. Li and Steffens reported that overexpression of PPO in tomato plants results in enhanced resistance to Pseudomonas syringae pv. tomato [56]. In this study, degradome sequencing analysis showed that nov-mir-132 targets the CCR gene and regulates its expression by cleavage. After inoculation with S. scitamineum, the expression level of CCR in YC05-179 and ROC22 was upregulated (Fig. 6B), which coincided with the increase in PPO activity at 2 and 5 d after inoculation (Fig. 2), suggesting that high expression of CCR transcript and PPO activity may promote the synthesis and accumulation of lignin to actively respond to S. scitamineum infection.

Interaction pathways between plants and pathogens
Interaction between plants and pathogens will stimulate various defense mechanisms in plants such as the formation and accumulation of phytoalexins and disease-related proteins, second messenger production, hypersensitive reactions (HRs), and programmed cell death (PCD). GK catalyzes the phosphorylation of glycerol, which is a key rate-limiting enzyme in the glycerol metabolic pathway [90]. Nonhost resistance 1 (NHO1), a member of GK, is involved in the JA and SA signal transduction pathways and plays a role in plant disease resistance [66,91]. The NHO gene is necessary for the R-gene related pathway and its expression can be activated by flagellin [65]. After inoculation with Xanthomonas oryzae, the expression level of OsNHO1 in rice rapidly increased within 3 h, peaked at 9 h, then gradually decreased and dropped to that of the control at 1 d, indicating that NHO1 is highly responsive to X. oryzae in rice [92]. In this study, nov-mir-63 targets GK. KEGG pathway analysis showed that this GK has NHO1 activity and participates in the plant-pathogen interaction pathway and glycerolipid metabolism. Moreover, after inoculation with S. scitamineum, the GK expression level was upregulated in the resistant variety YC05-179 but downregulated in the susceptible variety ROC22 at 2 d (Fig. 6D), suggesting that the GK gene may play a positive role in sugarcane resistance to S. scitamineum.
Ca 2+ is a second messenger involved in plant-pathogen interactions. MLO is a recessive susceptible gene. Ca 2+ -mediated MLO binds to calmodulin and is negatively correlation to plant disease resistance [62]. Plants often show susceptibility when the MLO gene normally expresses the MLO protein, whereas exhibit broad-spectrum disease resistance when the MLO gene does not express or express non-functional proteins [63]. Previous studies have shown that mol mutations can improve plant resistance to bacteria and fungi [93,94]. Therefore, in plants, the MLO gene is equivalent to a susceptible gene. The MLO gene has been associated with resistance to powdery mildew of wheat [95], as well as resistance to leaf blight of barley [96], necrosis of infected parts, and suppression of pathogen expansion [96]. In this study, we found that nov-mir-10 could target MLO. After inoculation with S. scitamineum, MLO was downregulated in YC05-179 at 2 d, and its expression level was lower than that in ROC22 (Fig. 6D). It has been suggested that nov-mir-10 in YC05-179 could more efficiently cleave MLO, thereby reducing the inhibition of defense response by the MLO protein and improving sugarcane resistance to smut pathogen.

Phytohormone signal transduction pathway
Phytohormones play a role in plant growth and development and responses to environmental stresses [97]. Phytohormones mainly include auxin, gibberellin acid (GA), cytokinins (CKs), ethylene (ETH), SA, JA, polyamines, and brassinosteroids (BRs). Among these phytohormones, SA, JA, ETH, and BR are disease resistance signaling molecules, whereas ETH and auxin biosynthesis promote each other and participate in plant responses to stress [97]. Ethylene insensitive 3 (EIN3) and EIN3-like 1 (EIL1) are transcription factors in the ethylene signaling pathway that promote the expression of ethylene response factor 1 (ERF1) and thereby regulate defense genes in response to pathogen infection [98]. In our previous proteomics research on sugarcane at 2 d after inoculation with S. scitamineum by iTRAQ analysis, four proteins involved in the ETH pathway were observed, including two 1aminocyclopropane-1-carboxylate oxidases (ACOs) that were responsible for ETH biosynthesis, as well as one EIN3 and one ERF1 that was responsible for ETH signaling [12]. One ACO (gi35014290) was upregulated in both sugarcane genotypes (YC05-179 and ROC22) and the other one (gi41615358) was downregulated in ROC22 only, but remained unchanged in YC05-179. One EIN3 (Sugarcane_Unigene_BMK.65773) and one ERF1 (gi35045219) were both upregulated in YC05-179, whereas it remained stable in ROC22 [12]. In this study, another EIL3 gene (Sugarcane_Unigene_BMK.64656) was targeted by nov-mir-143. After inoculation with S. scitamineum, the expression of EIL3 (Sugarcane_Uni-gene_BMK.64656) in YC05-179 was stable, but in ROC22 it was rapidly increased at 2 d and 5 d (Fig. 6C). It is speculated that smut pathogen attack enhances sugarcane ethylene metabolism, which is involved in EIL3 in ROC22 and in turn favors resistance to S. scitamineum infection [99]. However, YC05-179 responded to smut pathogen infection by other family genes in the ETH pathway. In addition, we also identified an AIP gene. KEGG analysis showed that AIP has small auxin up RNAs (SAUR) activity. SAUR plays a negative regulatory role in auxin biosynthesis and transport [100]. After inoculation with S. scitamineum, the expression level of AIP was upregulated in YC05-179 at 2 and 5 d and increased at 2 d in ROC22 by qRT-PCR (Fig. 6B). However, since this result differs from the AIP expression patterns obtained from degradome sequencing, further experiments are needed to verify how AIP responds to the infection by S. scitamineum. SAMDC and S-adenosylmethionine synthetase (SAMS) are important genes for polyamine biosynthesis [101,102]. Previous studies have shown that overexpression of the SAMDC gene can increase wilt resistance that is induced by Verticillium dahliae and Fusarium oxysporum in transgenic tobacco plants [103]. The expression of a SAMS gene in sugarcane increases after infection with S. scitamineum [104]. We found that miR162a targets SAMDC after inoculation with S. scitamineum and promotes the expression of SAMDC in both YC05-179 and ROC22 at 2 d (Fig. 6B), indicating that infection by S. scitamineum may enhance polyamine metabolism pathway to improve sugarcane smut resistance at the early stage.

Resistance-related transcription factors
Infection by pathogenic bacteria can stimulate plant transcription factors to participate in defensive responses. MYB is a typical transcription factor in plants that regulates phenylpropanoid metabolism [105] and is involved in hormonal signal transduction pathways such as auxin [106], JA [107], and ETH [107]. MYB is also involved in the systematic acquired resistance (SAR) and HR reactions [108]. Previous studies have found that MYB regulates PAL synthesis and is a positive regulator of phenylpropanoid anabolism [105]. In addition, MYB can be induced by exogenous stresses such as JA and SA, or TMV infection, and it can also activate disease-resistant defense responses involving the PR gene [109]. In this study, we found that MYB2 could be targeted by miR858b. After 5 d of inoculation with S. scitamineum, MYB2 gene was upregulated in YC05-179 and remained stable in ROC22 (Fig. 6D), which coincided with the expression pattern of PPO activity (Fig. 2), suggesting that miR858b may have different cleavage effects on MYB2 then regulate PAL synthesis in different resistant varieties of sugarcane, ultimately leading to different levels of resistance to smut in YC05-179 and ROC22. Similarly, Yang et al. demonstrated that a novel anther-specific myb gene (NtMYBAS1) from tobacco was a functional anther-specific transcription factor, which was likely to be a positive regulator of PAL synthesis in sporophytic [105].

miRNA feedback regulation
Previous studies have found that miRNAs can regulate target gene expression and involved in plant growth, development, and stress response, but also feedback regulate their metabolic synthesis [27]. Xie et al. found that Arabidopsis dicer-like1 (DCL1), which plays an important role in the formation of mature miRNAs, is subject to negative feedback regulation through the activity of miR162 [110]. AGO protein is an important part of RNA-induced silencing complex (RISC) with the function of cleaving miRNA target gene or inhibiting translation [111]. miR168 can target AGO protein and thus regulate the miRNA-regulated target genes in plants through changes in AGO protein expression [111]. In this study, AGO 1B was targeted by miR168a-5p. After inoculation with S. scitamineum, the upregulated expression of miR168a-5p in YC05-179 causes a decrease in the expression level of target AGO 1B and may promote miRNA-mediated accumulation of disease-related target genes, which in turn resists further infection of S. scitamineum. Meanwhile, after inoculation with S. scitamineum, the expression patterns of AGO 1B and miR168a-5p in ROC22 were opposite to that in YC05-179, i.e., miR168a-5p was downregulated and the expression level of AGO 1B was upregulated, suggesting that the miRNA self-feedback pathway may involve in sugarcane responses to S. scitamineum.

Conclusions
In the present study, the S. scitamineum was rapidly proliferated and the enzyme activities of POD, SOD, CAT, PPO, PAL, and TAL in the reactive oxygen species metabolic pathway and phenylpropanoid metabolic pathway were significantly changed at 2 and 5 d in the compatible and incompatible interactions between sugarcane and S. scitamineum. Furthermore, 97 known miRNAs and 112 novel miRNAs with 309 predicted target genes were identified by degradome sequencing. GO and KEGG pathway analyses showed that many predicted target genes enriched in regulation and metabolism. qRT-PCR validation demonstrated that there was no obvious negative regulatory relationship between miRNAs and their target genes. This study elucidates the underlying response of sugarcane to S. scitamineum infection and provides useful information on the interplay between miRNAs and their predicted targets. In the future, genetic transformations can be done to further our understanding on whether these genes enhance smut resistance in sugarcane.

Additional files
Additional file 1: Table S1. The qRT-PCR primers of the predicted target genes. Table S2. The qRT-PCR primers of miRNAs. Figure S1.