“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.

According to the instruction of Trizol™ Reagent (Invitrogen, Carlsbad, CA), total RNA was extracted from each sample, and RNA purity was evaluated using NanoPhotometer® spectrophotometer (IMPLEN, CA, USA). RNA Nano 6000 Assay Kit and Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA) were used to assess RNA concentration and integrity.

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).

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.

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).

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).

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).

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.

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).

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.

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 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).

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²⁺-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.

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).

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).

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.

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 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.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.