microRNAs and Their Targets in Apple (Malus domestica cv. “Fuji”) Involved in Response to Infection of Pathogen Valsa mali

miRNAs are important regulators involving in plant-pathogen interactions. However, their roles in apple tree response to Valsa canker pathogen (Valsa mali, Vm) infection were poorly understood. In this study, we constructed two miRNA libraries using the twig bark tissues of apple tree (Malus domestica Borkh. cv. “Fuji”) inoculated with Vm (IVm) and PDA medium (control, BMd). Among all detected miRNAs, 23 miRNAs were specifically isolated from BMd and 39 miRNAs were specifically isolated from IVm. Meanwhile, the expression of 294 miRNAs decreased; and another 172 miRNAs showed an increased expression trend in IVm compared with that in BMd. Furthermore, two degradome sequencing libraries were also constructed to identify the target genes of these miRNAs. In total, 353 differentially expressed miRNAs between IVm and BMd were detected to be able to target 1,077 unigenes with 2,251 cleavage sites. Based on GO and KEGG analysis, these genes were found to be mainly related to transcription regulation and signal transduction. In addition, we selected 17 miRNAs and 22 corresponding target genes to screen the expression profiles when apple twigs were infected by Vm. The expression trends of most miRNAs/target genes were consist with the results of deep sequencing. Many of them may involve in the apple twig-Vm interaction by inducing/reducing their expression. What's more, miRNAs and their target genes regulate the apple twig-Vm interaction by forming many complicated regulation networks rather than one to one model. It is worth that a conserved miRNAs mdm-miR482b, which was down regulated in IVm compared with BMd, has 14 potential target genes, most of which are disease resistance related genes. This indicates that mdm-miR482b may play important roles in apple twig response to Vm. More important, the feedback regulation of sRNA pathway in apple twig is also very complex, and play critical role in the interaction between apple twig and Vm based on the results of expression analysis. In all, the results will provide insights into the crucial functions of miRNAs in the woody plant, apple tree-Vm interaction.


INTRODUCTION
As a class of conserved small RNAs in eukaryotes, microRNAs (miRNAs), with ∼21-24 nucleotides (nt) in length, are generated from single-stranded RNA transcript with a self-complementary hairpin structure, and act as specific repressors of target gene expression through DNA methylation, histone modification, cleavage of the target transcript or inhibition of translation Bartel, 2004;Ramachandran and Chen, 2009;Voinnet, 2009;Seo et al., 2013;Clancy et al., 2014). Since the first plant miRNA was discovered in 2002, miRNAs have been confirmed to act as powerful endogenous regulators in various biological processes, including developmental patterning, resistance to abiotic stresses and biotic stresses (Llave et al., 2002;Palatnik et al., 2003;Fujii et al., 2005;Navarro et al., 2006). Increasing evidences demonstrate that plant miRNAs are key players in plant-microbe interaction mechanisms (Jin, 2008;Padmanabhan et al., 2009;Weiberg et al., 2014;Thiebaut et al., 2015). However, most studies have focused on herbaceous plants, in particular, Arabidopsis, rice, wheat, etc. Regulation by sRNAs, particularly in response to pathogen infection, in economically important woody plants is poorly understood.
Apple (Malus domestica Borkh.) is one of the most widespread and economically important fruit trees worldwide. Specifically, in China, the apple planting area and production are nearly 2.26 million hectares and 40 million tons annually. The apple industry has been a pillar for rural economic development in northern China. However, the apple industry in China is suffering from huge economic losses due to the prevalence of various diseases. Among them, apple tree Valsa canker, caused by Valsa mali (Vm), is the major one . The pathogen causes decomposition of bark tissues of main trunk, lateral branch, and even leads to the death of twigs, branches, trunks, and eventually the entire tree. However, no immune or highly resistant variety to Vm has been reported. To date, the disease still cannot be effectively controlled, mainly because the information about the mechanism of apple tree in response to Vm was limited. Thus, accelerating studies on the mechanism of apple-Vm interaction will provide opportunities to increase the resistance to apple tree Valsa canker.
In the past few years, many studies about the pathogenic mechanism of Vm have been conducted at the genome level, transcriptome level, and functional genomics level (Ke et al., 2013;Hu et al., 2014;Yin et al., 2015). However, resistance mechanism of host in response to pathogen is not very clear. In 2016, Yin et al. (2016) reported that 2,713 apple genes were significantly up-regulated upon Vm infection, and they mainly focused on chitin, hormone and cell death biological processes. Because miRNAs are essential gene expression regulators, we hypothesize that they could be involved in the resistance response of apple tree to Vm. In fact, many miRNAs of apple have been isolated using both computational and high-throughput sequencing methods from different varieties and tissues, and these miRNAs have been confirmed to be related to tissue development (Gleave et al., 2008;Varkonyi-Gasic et al., 2010;Xia et al., 2012;Szczesniak and Makalowska, 2014;Qu et al., 2016;Xing et al., 2016). Moreover, many miRNAs were found to be involved in the biotic stress responses. In 2014, Yu et al. (2014) found that five apple miRNAs isolated from seedlings stems were responsive to apple ring rot infection. Ma et al. (2014) reported that the novel apple miRNA Md-miRLn11, isolated from a mixed sample of leaves, phloem, flower, and fruits, could affect the resistance to apple leaf spot disease by regulating a Md-NBS gene expression. In the next year, four apple miRNAs were isolated from the shoot tip tissues of four distinct rootstocks bench grafted with "Crimson Gala" scions, and presumed to be potentially important for resistance of apple trees to fire blight (Kaja et al., 2015). However, the specific miRNAs in the apple tree twig bark tissues, especially those involved in the interaction between the apple tree and Vm, are still not known.
In this study, high-throughput and degradome sequencing analyses were used to comprehensively identify the miRNAs and corresponding target genes from the apple twig bark tissues in response to Vm; these results will lay the foundation for exploration of the apple-Vm interaction mechanism at posttranscription level.

Plant Materials and Samples Collection
"Fuji" apple plants (M. domestica Borkh.) were grown in a greenhouse with 25 • C and natural daylight conditions. The twigs were inoculated with Vm as described by Li et al. (2016). The samples of lesion borders were collected at 0, 12, and 48 h post-inoculation (hpi). Three individual biological repeats were prepared.

Library Construction and RNA Deep Sequencing
Three repeats of control samples and inoculation samples were mixed for small RNA library construction, respectively. Three micrograms of total RNA per sample was used as input material for the small RNA library. Small RNA library construction and RNA deep sequencing were proceeded following the detailed protocol provided by the genome sequencing company (Novogene, China).

Bioinformatic Analysis of Sequence Data
Raw data (raw reads) were first processed through custom Perl and Python scripts to obtain clean data. The clean data were mapped to the reference sequence in miRBase21.0 by Bowtie (Langmead et al., 2009), without mismatch to look for known miRNAs. Then, the other reads were integrated to predict novel miRNAs using the available software miREvo (Wen et al., 2010) and mirdeep2 (Friedlander et al., 2011). The miRNA counts as well as base bias were identified using custom scripts. The miFam.dat (http://www.mirbase.org/ftp.shtml) was used to look for families of known miRNAs. The novel miRNA precursor was submitted to Rfam (http://rfam.sanger.ac.uk/search/) to look for Rfam families.

Venn Diagrams of Known miRNAs and Novel miRNAs
Normalization formula (Normalized expression = Mapped read count/Total reads * 1,000,000) was used to estimated miRNA expression levels (Zhou et al., 2010). DEGseq (2010) R package was used to analyze the differential expression of two samples with the criterion of Q < 0.01 and |log2(foldchange)| > 1 (Storey, 2003).

Construction of Degradome Libraries
Target genes of candidate miRNAs were verified by degradome sequencing using total RNA same to the RNA used for small RNA sequencing library construction following the published parallel analysis of RNA Ends (PARE) protocol (German et al., 2009). The data analysis was processed following the procedure instructions (LC Sciences, Hangzhou, China).

Target Gene Prediction and Annotation for Known and Novel miRNAs
The psRobot_tarin psRobot was performed to predict target genes of miRNA (Wu et al., 2012). To further explore the detailed molecular mechanism of miRNAs in apple twig response to Vm, the target transcripts of differentially expressed miRNAs were analyzed by Gene Ontology Functional Annotation Suite (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG). Subsequently, Revigo tool (http://revigo.irb.hr) was implemented for enrichment analysis of the target genes. The KOBAS software was used to test the statistical enrichment of the target gene candidates in the KEGG pathways (Mao et al., 2005).

Expression Analysis of miRNAs and Their Corresponding Target Genes
The expression of miRNAs in the apple twig bark tissues inoculated with Vm was determined using stem-loop RT-PCR technology (Feng et al., 2012). The expression analysis of the corresponding target genes was conducted as described by Feng et al. (2012). The apple translation elongation factor 1 alpha-subunit gene (EF-1a, Gene ID/Accession no. DQ341381) was selected to normalize the expression data of each miRNA and target gene. The relative expression of each miRNA and target gene was analyzed using the software of HemI (Heatmap Illustrator, version 1.0; Deng et al., 2014). All primers used in this work were listed in Table S12.

Overview of the sRNA Sequencing Results
Two sRNA libraries of apple twig (BMd and IVm) were constructed using Illumina high-throughput sequencing technology. In BMd, 13,065,359 raw reads (including Junk reads, Rfam, Repeats, valid reads, etc.) were generated, and 5,232,937 unique reads were isolated. Meanwhile, in IVm, 7,753,801 unique reads were isolated from 27,234,866 raw reads (Table S1). The proportion of different sRNAs for Rfam in BMd and IVm are shown in Figure S1. When the low-quality and junk sequences were removed, the left 3,491,357 valid reads in BMd and 4,949,663 valid reads in IVm were further analyzed, respectively. Among all sRNAs in both two libraries, 24 nt-length sRNAs were the most abundant, followed by 21, 22, and 23 nt-length. Despite this, an obvious difference existed between mock and treatment sample. For instance, the proportion of 24 nt-length sRNAs in BMd is 33.8%, whereas 24.9% in IVm. This difference was also found in other sRNAs with different length (Figure 1).

Identification of miRNAs in Apple Twig Bark
To explore all miRNAs in the bark of apple twig, we first compared all the valid reads in both BMd and IVm with known plant miRNAs in the miRBase to identify the conserved miRNA homologs. In total, 187 miRNAs of 46 families, generated from 138 miRNA precursors, were identified by mapping to the known miRNAs of apple. Meanwhile, a total of 37 miRNAs from 32 families, which could be mapped to miRNAs/pre-miRNAs of the selected species and the extended genome sequences from the mapped genome loci could form hairpin structures, were also preliminarily determined as candidate apple miRNAs. Moreover, we found 153 distinct reads could be mapped to the miRNAs/pre-miRNAs of the selected species in miRbase and the corresponding genome, but the extended sequences at the genome loci could not form hairpins, and another 51 reads could also be mapped to miRNAs/pre-miRNAs of the selected species in the miRbase, but the reads could not be mapped to the corresponding genome. In addition, 948 candidate miRNAs from 944 precursors, which could not be mapped to the pre-miRNAs of the selected species in miRbase were further determined as novel candidate miRNAs in apple twig bark, because they could be mapped to the corresponding genome, and the extended genome sequences from the genome may form hairpins (Table 1). Thus, we isolated 187 known miRNAs and 1189 candidate novel miRNAs from the apple twig bark tissues (Table S2). Most known miRNAs were conserved compared with those in other plants (Table S3; Figure S2). Moreover, the first nucleotide of these candidate miRNAs generally tended to be uracil (U) ( Figure S3). All the raw data and clean data of these two libraries have been submitted to the GEO repository of NCBI (GEO accession: GSE104752), with the sample number GSM2807179 for BMd and GSM2807180 for IVm.

Bio-informatic Analysis of Differentially Expressed miRNAs in Apple Twig Response to Vm
Specifically expressed miRNAs in apple twig inoculated with Vm were identified by comparing the miRNAs in IVm with those in BMd. Most miRNAs were co-expressed, with only 23 and  gp1a, Reads map to specific miRNAs/pre-miRNAs in miRbase and the pre-miRNAs further map to the genome & EST. gp1b, Reads map to selected (except for specific) miRNAs/pre-miRNAs in miRbase and the pre-miRNAs further map to the genome and EST. gp2a, Reads map to selected miRNAs/pre-miRNAs in miRbase. The mapped pre-miRNAs do not map to the genome, but the reads (and of course the miRNAs of the pre-miRNAs) map to genome. The extended genome sequences from the genome loci may form hairpins. gp2b, Reads were mapped to miRNAs/pre-miRNAs of selected species in miRbase and the mapped pre-miRNAs were not further mapped to genome, but the reads (and of course the miRNAs of the pre-miRNAs) were mapped to genome. The extended genome sequences from the genome loci may not form hairpins. gp3, Reads map to selected miRNAs/pre-miRNAs in miRbase. The mapped pre-miRNAs do not map to the genome, and the reads do not map to the genome. gp4, Reads do not map to selected pre-miRNAs in miRbase. But the reads map to genome and the extended genome sequences from genome may form hairpins.
39 miRNAs expressed exclusively in BMd or IVm, respectively (Figure 2). The relative expression of miRNAs in Vm-infected apple twig compared with the mock-inoculated control was also analyzed with the reads' abundance. The samples showed various expression profiles when apple twig was challenged with Vm ( Figure 3). Among them, 294 miRNAs were down-regulated when apple twigs were challenged with Vm, including 48 known apple miRNAs and 246 novel candidate miRNAs. Conversely, 172 miRNAs were up-regulated, including 36 known apple miRNAs and 136 novel candidate miRNAs (Table S4). In fact, only 113 and 40 candidate miRNAs showed down-and up-regulation in Vm-infected barks compared with the mock-inoculated control when the p ≤ 0.05 and |logfoldchange| ≥ 1 ( Figure S4).

Overview of the Degradome Sequencing Results
To better understand the functions of miRNAs in apple twigs responding to Vm infection, the target genes of miRNAs in this study were determined using degradome sequencing technology. Two high-throughput degradome sequencing libraries (TBMd and TBVm) were constructed to determine the target genes of miRNA. A total of 2,251 cleavage sites from 1,077 unigenes were identified for 353 candidate miRNAs, including 158 known apple miRNAs. What's more, all of these genes have potential roles in pathogen responses (Table S5). For example, the target of a known apple miRNA, mdm-miR397a_L+1R+2_1ss22AG, was identified to be a nitrite reductase. Moreover, another important gene for disease resistance, Ca 2+ -transporting ATPase, was identified to be the target of another known apple miRNA, mdm-miR3627a_R-1. The target genes for some novel candidate miRNAs were also determined. The gene ppe-MIR172c-p5 was found to regulate the expression of monodehydroascorbate reductase, and PC-5p-154495_9 could affect the ABC transporters pathway (KEGG) by targeting ATP-binding cassette ( Table 2). However, there were still large numbers of miRNAs whose targets remained undetected. All the results of these two libraries have also been submitted to the  GEO repository of NCBI (GEO accession: GSE104752), with the sample number GSM2807181 for TBMd and GSM2807180 for TBVm.

Target Annotation of Differentially Expressed miRNAs in Apple Response to Vm
The GO annotation showed that these target genes could participate in various cellular processes, such as apoptosis, transcription regulation, auxin-mediated signaling pathway, ATP binding, DNA binding and protein binding activity (Table S6). To better explore the biological meaning, we further processed an enrichment analysis of target genes with p ≤ 0.01 using the Revigo tool. 106 genes were classified into three categories, and most them were linked into a complex net. Based on the results, the two most basic "Cellular Component" categories are "nucleus" and "SCF ubiquitin ligase complex" (Figure S5). In "Biological Process, " the top two are "auxin mediated signaling pathway" and "regulation of transcription, DNA-dependent" (Figure S6). And in "Molecular Function, " the two most categories are "protein dimerization activity" and "translation elongation factor activity" (Figure S7). Meanwhile, 124 KEGG pathways were classified, and the most common pathways were "Cysteine and methionine metabolism, " "Phenylpropanoid biosynthesis, " "Ribosome, " "Regulation of autophagy, " "basal transcription factors, " "Aminoacyl-tRNA biosynthesis, " "Nitrogen metabolism, " and "Glutathione metabolism" (Figure S8). We also found some pathways related particularly to disease response, such as plant-pathogen interaction, the MAPK signaling pathway, and ABC transporters ( Table S7).

Vm Infection Triggers the Expression Change of miRNAs and Target Genes
Based on the results of deep sequencing, lots of miRNAs were speculated to play roles in interaction between apple tree and Vm. And 17 miRNAs, which showed an obvious induced/reduced expression in IVm compared with BMd (p ≤ 0.05), together with more than 10 normalized reads number in IVm or BMd, were selected for qRT-PCR analysis to explore their expression in apple twig tissues response to Vm infection. Our results showed that mdm-miR156a, mdm-miR858b, mdm-miR391, mdm-miR477a, and gma-miR6300 were up-regulated during Vm infection process, and mdm-miR391 showed noticeable induction in the early stage of infection (hpi). Compared with the control, the expression levels of mdm-miR319a, mdm-miR403a, ppe-miR530, and gma-miR160-p3 were decreased. Moreover, mdm-miR482b, mdm-miR535d, bol-miR9410, peu-miR2916-p3, PC-5p-409096, PC-3p-102462, and PC-3p-166024 expressed stably during the entire infection progress of Vm (Figure 4; Table S8).
Meanwhile, 22 identified targets, which were predicted to be important to stress response of plants and showed an obvious induced/reduced expression in TBVm compared with TBMd (p ≤ 0.05), were also selected for transcript analysis. Eleven genes were up-regulated when apple twigs were inoculated with Vm, and most of them were related to plant disease resistance, such as the homologs of ATL2, SNC1, WRKY75, RPPL1, and RGA3. Compared to the control, the expression of the homologs of RPM1, At1g62630, At1g12290, and NAC021 decreased, and the  other seven genes showed stable expression during the entire infection progress of Vm ( Figure 5; Table S9).

miRNAs and Corresponding Target Genes Function by Forming Complex Networks
The regulatory mechanisms of miRNAs in the apple twig bark response to Vm were very complex. Based on the detected targets, we found that the one-to-one regulation model was not prevalent. Majority of miRNAs could target more than one transcript, and most miRNAs and their targets constructed a complex regulatory network. For example, ppe-MIR169i-p5_2ss1CT19GT was found to cleave 80 target transcripts, and the transcript XM_008348120.1 could be regulated by 29 candidate miRNAs (Tables S10, S11). Notably, the expression of atlastin GTPase 2 (ATL2) could be regulated by nine novel candidate miRNAs, and the conserved miRNA mdm-miR482b in apple could target 14 genes; most target genes were disease resistance related genes, such as the Rho-type GTPase-activating protein and the disease resistance protein family members (Figure 6). Moreover, these miRNAs/targets on the cross point may be essential in the disease response signal network.

Feedback Regulation of the sRNA Pathway in Apple-Vm Interaction
Dicer and AGO genes are the core components of the sRNA pathway; however, the role of miRNAs in their regulation remains unknown. Based on the results of degradome sequencing, we found four transcripts of the "switch" gene controlling miRNA generation (Dicer1), which could be cleaved by three different but highly homologous miRNAs from the same family of miR162. Moreover, two transcripts coding AGO2 were identified to have a cleavage site of three similar miRNAs from the miR403 family. Moreover, the important component, AGO1, which mainly binds miRNAs to cleave the corresponding targets, could be regulated by four miRNAs. Among them, two FIGURE 6 | Disease resistance-related regulation network of miRNAs and corresponding targets in the interaction of apple twig bark and V. mali. Several miRNAs and their corresponding targets were selected to show the complex regulation net, and most of them were associated with plant disease resistance. (A) One miRNA was found to be able to regulate several targets. (B) One target could also be regulated by various miRNAs. ATL2, atlastin GTPase 2; ASA2, anthranilate synthase; At1g12290 and At1g59620, probable disease resistance protein; At3g14460 and RPM1, putative disease resistance protein; At1g62630 and At4g27190, disease resistance protein; ATJ8, J proteins; RGA2 and RGA3, Rho-type GTPase-activating protein; RPP13L2, putative disease resistance RPP13-like protein 2; RPP8L3, disease resistance RPP8-like protein 3; RPPL1, putative disease resistance RPP13-like protein 1; SNC1, TIR-NBS-LRR class disease resistance protein.
are from the miR168 family, and another two are from miR159 and miR530 families (Figure 7). To further verify the feedback regulation of the sRNA pathway in apple twig-Vm interaction, we detected the expression profiles for all these miRNAs and target genes (Figure 8). The gene mdm-miR162 was down-regulated at 12 hpi, while the expression of the corresponding target gene DCL1 began to increase and peaked at four-fold of the control. mdm-miR403a and AGO2 showed a similar expression trend during the Vm infection. However, mdm-miR159c and mdm-miR168b were increased, and the corresponding target gens, AGO1 and AGO1B, were down-regulated at 12 hpi, particularly the expression of mdm-miR168b and AGO1B, with ∼40-and 10-fold changes compared with control, respectively. These results indicate that the feedback regulation of the sRNA pathway in apple twig was much more complex and critical. The detailed mechanism will be another interesting scientific question remaining to be explored.

DISCUSSION
As an abundant class of endogenous, non-coding small RNAs, miRNAs play important regulatory roles in organ development, phase change and defense responses, etc. (Bartel, 2004). Some apple miRNAs were identified using computational and/or sequencing approaches, and were found to be responsive to tissue development and biotic stresses response (Varkonyi-Gasic et al., 2010;Xia et al., 2012;Yu et al., 2014). However, there is no report of miRNAs involved in the disease response to Vm. Thus, the identification of apple miRNAs and the corresponding target genes in response to the Vm infection will help us to understand the regulatory mechanism of apple twig bark and Vm interaction.
In this study, we performed Illumina deep sequencing to construct two sRNAs libraries (BMd and IVm). The amount of 24 nt-length sRNAs in apple twig tissues was biggest in both BMd and IVm, which is consistent with the length distribution of sRNAs found in other apple tissues (Xia et al., 2012;Qu et al., 2016). Actually, 21 and 24 nt-length sRNAs were the main two sRNA classes among unique sequences based on previous reports on many other plants, such as wheat (Feng et al., 2017), rice (Morin et al., 2008), grapevine (Wang et al., 2011), tomato (Cao et al., 2014), and so on. According to the previous studies, 21 nt-length sRNAs are mainly related to mRNA cleavage (Hamilton and Baulcombe, 1999), while the 24 nt-length sRNAs are primarily involved in RNA-directed DNA methylation (Hamilton et al., 2002;Qi et al., 2006). For instance, it has been demonstrated that 24 nt-length sRNAs are associated with transcriptional gene silencing by targeting DNA methylation of three endogenous transposable elements in Arabidopsis thaliana (Lewsey et al., 2015). In wheat, 24 nt-length sRNAs were found to mainly match class I and class II transposable elements, which were also predicted to be more likely to be methylated (Cantu et al., 2010). But these results don't mean that 21 nt-length sRNAs were unrelated with DNA methylation. In subsequent studies, 21 nt-length sRNAs were confirmed to be involved in the establishment of DNA methylation, whereas the 21 ntlength species contributed to the amplification and maintenance of DNA methylation (Kim and Zilberman, 2014;Bond and Baulcombe, 2015). However, whether the roles of 24 nt-length miRNAs in plant-disease interactions are conducted through DNA methylation is still unknown. What's more, the proportion of same length sRNAs in BMd and IVm is also different in this study. In fungi and animals, there are more than two pathways to generate various sRNAs, with the main ways of Dicer-dependent and Dicer-independent (Aravin et al., 2006;Lee et al., 2010). In FIGURE 7 | Feedback regulation of main components of RNAi pathway in apple twig bark-V. mali interaction. Three main components of RNAi pathway were found to be feedback regulated by several miRNAs. (A) Four transcripts of Dicer1 were detected to be cleaved at one common site by three miRNAs from the same family. (B) Three transcripts of AGO1 were detected to be cleaved at distinct sites by four miRNAs from three families. (C) Two transcripts of AGO2 were detected to be cleaved at one common site by three miRNAs from same family.
plant, there is no report to illustrate whether there is another way to generate sRNA like that. However, it has confirmed that many plant miRNAs expressed specially in different tissues or developmental stages, even response to the stresses response, thereby contributing to the specific target gene regulation (Inui et al., 2010). Actually, in our results, when the apple twig was challenged with Vm, the expression of key proteins in RNAi pathway changed a lot. Whether the difference of the proportion of same length sRNAs is associated with the expression change of these genes will still be explored in further study.
Plant miRNAs mainly regulate gene expression by cleaving mRNAs directly, and the cleavage normally occurs at the tenth nucleotide of the complementary region Varkonyi-Gasic et al., 2010). Thus, the determination of corresponding target genes of miRNAs in apple twig bark is critical for functional exploration of the miRNAs. Degradome sequencing, as a high-throughput technology for target detection, has been used in several studies in the past few years (Addo-Quaye et al., 2008;Li et al., 2010;Yang et al., 2014;Yao et al., 2015). In this study, the target genes were also detected using this technology. To get more accurate results, gene mapping and annotation were referred to the transcriptomic information of apple twig bark challenged with Vm. As reported, many plant miRNAs play roles by targeting transcription factors or other regulatory/stress-response proteins (Ehrenreich and Purugganan, 2008;Guo et al., 2008;Varkonyi-Gasic et al., 2010;Zhang et al., 2013;Niu et al., 2016). In this study, the target genes of 84.5% of known apple miRNAs were determined, and many of them were associated with disease responses, such as Ca 2+ -transporting ATPase, aarF domain-containing kinase, E3 ubiquitin-protein ligase, and the mlo protein. What's more, most miRNAs and the corresponding target genes function by forming a complex interaction network. That means one miRNA could target several genes in different families, and one target gene could also be regulated by distinct miRNAs (German et al., 2009;Li et al., 2010;Lin et al., 2010). For example, mdm-miR482b could target 14 genes, and most these genes are CC-NBS-LRR class disease resistance proteins (At1g12290, At1g59620, At3g14460, At1g62630, At4g27190 and RPM1, RPP13L2, RPP8L3, and RPPL1) and a TIR-NBS-LRR class disease resistance protein (SNC1) (Bakker et al., 2006;Xu et al., 2014). In addition, mdm-miR482b could also regulate the expression of two Rhotype GTPase-activating proteins, which were critical regulator of calcium-dependent signal transduction, H 2 O 2 -mediated cell death, and so on (Wu et al., 2000). On the other hand, atlastin GTPase 2 (ATL2) was found to could be cleaved by nine novel candidate miRNAs. ATL2 was rapidly activated after incubation with elicitors of the pathogen response and was associated with the plant defense signaling pathways by ubiquitin ligases (Serrano and Guzmán, 2004). In 2016, miR863-3p was found to silence two negative regulators (ARLPK1 and ARLPK2) of Arabidopsis defense at the early stage of infection of Pseudomonas FIGURE 8 | Relative transcription levels of miRNAs and corresponding target genes involved in the feedback regulation of sRNAs pathway in apple twig-Vm interaction. One DCL1, one AGO2, and two AGO1 transcripts were identified to be the target genes of mdm-miR162, mdm-miR403a, mdm-miR159c, and mdm-miR168b. Their transcript accumulation in apple twig bark tissues after challenge with V. mali was analyzed (12 and 48 hpi). The data were normalized to the inner reference gene apple translation elongation factor 1 alpha-subunit (EF-1a, Gene ID/Accession no. DQ341381). The mean value was calculated from three independent biological replicates. The relative expression level was analyzed as the fold-change of the mock-inoculated plants at that time point using the comparative 2 − CT method. syringae through mRNA degradation to promote immunity, and it could also silence SERRATE to positively regulate defense through translational inhibition at late infection stage (Niu et al., 2016). This result indicated that apple miRNAs might also play distinct roles by regulating the expression of the various target genes, possibly during the different developmental stages or the different stages of biological processes, such as the disease response. The expression of some selected miRNAs and target genes in apple twig bark in response to Vm were detected, and most were consistent with the results of sequencingbased estimations with up-/down-regulated expression trends. Combining the expression and the related reports, we speculated these regulated miRNAs/genes might contribute positively or negatively to the host-pathogen interaction.
Although many disease resistance-related miRNAs/genes were predicted and detected in this study, the host was still susceptible to Vm in reality. This indicated Vm might have evolved more "weapons" for attacking the host. Except for effectors, toxins, pectinases, and sRNAs from the pathogen could also enter the host cells to attack the host immunity by inhibiting the expression of critical genes in the disease response pathway, especially the sRNA pathway (Weiberg et al., 2013;Weiberg and Jin, 2015). For the sRNA pathway, short RNA duplexes are cleaved from long double-stranded RNAs by Dicer protein, and these sRNAs can guide gene silencing only when they are bound to the Argonaute proteins (AGOs) through a similar basepairing mechanism with their target RNA transcripts (Zamore et al., 2000;Bernstein et al., 2001;Chapman and Carrington, 2007;Kim et al., 2010). In this study, we found that the feedback regulation of the sRNA pathway was more complex than that in the other plants. The induced expression of Dicer genes may result in a change of the quantity of different kind sRNAs. The results of differences in the sRNAs classes observed in mock and inoculated samples offered a good evidence to support it. Meanwhile, AGO1 seem to play critical role in this interaction. It was found that mdm-MIR168b-p5, ppe-MIR168-p5, mdm-miR159c, and ppe-miR530_L-2R+2 could all target AGO1 in apple. However, AGO1 in Arabidopsis was regulated by miR168 (Vaucheret et al., 2004). Therefore, the role of apple miRNAs in the orderly regulation of the sRNA pathway and the entry of the Vm sRNAs into the apple twig cells to attack the host sRNA pathway (by hijacking Dicer or AGO from the host sRNAs) to disturb the normal order of regulation of disease resistance are extremely interesting subjects for research.
In conclusion, candidate miRNAs and corresponding target genes involved in the apple-Vm interaction were detected using Illumina and degradome sequencing technology. Many stress responses related miRNAs/genes were isolated and the expression of them were further detected, particularly the cross point in signal network of disease resistance. Moreover, the feedback regulation of the key component proteins in sRNA pathway regulated by miRNAs was also found to be associated with the disease response. The sRNA pathway may have a mysterious function in the apple twig-Vm interaction, which should be further explored.

AUTHOR CONTRIBUTIONS
HF and LH: contributed to the design of the work; HF and TZ: analyzed the sequencing data; MX and XG: drafted the work or revised it critically for important intellectual content; XZ: detected the expression of miRNAs and the corresponding target genes; LH: finally approved the version to be published, and was accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved.

ACKNOWLEDGMENTS
The study was funded by the National Natural Science Foundation of China (No. 31501591), and the Doctoral Scientific Research Foundation of Northwest A&F University.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2017. 02081/full#supplementary-material        Table S1 | Overview of reads from raw data to cleaned sequences.        Table S10 | Number statistics of target for each miRNA through predicting and degradome sequencing in both IVm and BMd.
Table S11 | Number statistics of mRNA with miRNA targeting sites through predicting and degradome sequencing in both IVm and BMd.
Table S12 | The primer sequences for quantitative real time PCR of candidate miRNAs and target genes.