MiRNA-seq-based profiles of miRNAs in mulberry phloem sap provide insight into the pathogenic mechanisms of mulberry yellow dwarf disease

A wide range of miRNAs have been identified as phloem-mobile molecules that play important roles in coordinating plant development and physiology. Phytoplasmas are associated with hundreds of plant diseases, and the pathogenesis involved in the interactions between phytoplasmas and plants is still poorly understood. To analyse the molecular mechanisms of phytoplasma pathogenicity, the miRNAs profiles in mulberry phloem saps were examined in response to phytoplasma infection. A total of 86 conserved miRNAs and 19 novel miRNAs were identified, and 30 conserved miRNAs and 13 novel miRNAs were differentially expressed upon infection with phytoplasmas. The target genes of the differentially expressed miRNAs are involved in diverse signalling pathways showing the complex interactions between mulberry and phytoplasma. Interestingly, we found that mul-miR482a-5p was up-regulated in the infected phloem saps, and grafting experiments showed that it can be transported from scions to rootstock. Based on the results, the complexity and roles of the miRNAs in phloem sap and the potential molecular mechanisms of their changes were discussed. It is likely that the phytoplasma-responsive miRNAs in the phloem sap modulate multiple pathways and work cooperatively in response to phytoplasma infection, and their expression changes may be responsible for some symptoms in the infected plants.

Mulberry trees that have long been cultivated for sericulture are susceptible to many diseases, among which yellow dwarf disease caused by phytoplasma is one of the most devastating 1 . Phytoplasmas are wall-less, obligate intracellular plant pathogens in the class Mollicutes 2 that infect several hundred economically important plants and cause devastating losses in agriculture and forestry 3 . The inability to culture phytoplasmas in vitro makes it difficult to characterize the plant pathogens at the molecular level, and the underlying molecular mechanisms of their pathogenicity are still poorly understood 4 .
When subjected to pathogen infection, the host plant activates sophisticated response mechanisms to reprogramme the expression of genes, proteins and metabolites 5 . The gene, protein and metabolite profiles in some host plants challenged with phytoplasmas have been investigated by differential methods [6][7][8][9][10][11][12][13][14][15][16] . MiRNAs functioning as negative regulators of gene expression are involved in the control of plant development and immunity 17 . Increasing evidence showed that miRNAs serve as an important mechanism for mediating gene expression during plant-pathogen interactions, and many miRNAs have been linked to resistance responses in plants [18][19][20][21] . A group of bacteria-responsive miRNAs and their target genes have been identified, and their regulatory functions have been extensively characterized in model plant species [21][22][23][24][25][26][27][28][29] . However, to our knowledge, only three studies have explored phytoplasma-responsive miRNAs in Mexican lime (Citrus aurantifolia L.), mulberry (Morus multicaulis Perr.) and Paulownia fortunei [30][31][32] . Although some miRNA families are conserved among various plant species, every species will have their own specific miRNAs, and the functions of these miRNAs may also be species-specific 33,34 . Moreover, different miRNAs, in addition to individual miRNAs in the same family, may be expressed differentially in various tissues and have different functions in response to the same pathogen infection 35 . Therefore, phytoplasma-responsive miRNAs have not been fully explored, and their mediating mechanisms for gene expression in response to phytoplasma are largely unknown.
Phloem is not only the major route for the translocation and distribution of organic metabolites but also an important mediator of whole-plant communication involved in whole plant events, including stress responses and long-distance signalling 36 . Recently, many miRNAs have been identified from the phloem exudates from Brassica napus 37 , Malus domestica (apple) 38 , and Lupinus albus 39 . Although the role of phloem miRNAs is not yet clear, some miRNAs in the phloem have been demonstrated to use long-distance signalling and have a role in mediating plant developmental patterning and stress responses 38,[40][41][42][43][44][45][46][47][48] . In plants, phytoplasmas are restricted to the sieve elements of phloem tissues 49 . Therefore, it is reasonable to assume that the phloem is the site where the host immediate defence response against phytoplasma occurs and is involved in the coordination of the defence response at the whole plant level. Identification and characterization of phytoplasma-responsive miRNAs in the phloem sap promises to enhance our understanding of the molecular mechanisms involved in yellow dwarf disease symptom development.
In the present study, based on transcriptome information for mulberry, we employed high-throughput sequencing to profile the miRNAs in the phloem sap during the response of mulberry to phytoplasma infection. Phytoplasma-responsive miRNAs were identified, and their potential target genes were predicted. In addition, the translocation and functions of the mul-miR482a-5p involved in the response of mulberry to phytoplasma infection were discussed. The results reported here may facilitate our understanding of phytoplasma pathogenicity.

Results
Purity assessment of phloem sap. To identify phloem-enriched sRNAs and ensure that the sRNAs identified in phloem sap did not result from contamination during sampling, the frequency of the ribulose bisphosphate carboxylase oxygenase (RuBisCo) large subunit gene that would be expected to in leaves, but not the sieve element-companion cell complex, was determined to assess the purity of the phloem sap sampled. The results detected no RuBisCo mRNA in the phloem sap samples, but this mRNA was clearly present in leaf tissue. Meanwhile, phloem-specific MmPP16 mRNA was detected in phloem sap samples (Fig. 1). This indicates that contamination from surrounding tissues in the collected phloem sap samples was very low.
Overview of small RNA in phloem saps. To examine the phytoplasma-responsive miRNAs in phloem sap, the sRNA libraries were constructed from phytoplasma-infected and healthy mulberry phloem saps and subjected to Solexa deep sequencing. After removing the rRNA, tRNA, and degradation products from the matching sequences, the remaining clean small RNA sequences were aligned with the miRNA precursors/mature miR-NAs in the miRBase database v21.0, and 86 known miRNAs members belonging to 78 families were identified ( Table 1). All non-annotated sRNA sequences were mapped to our mulberry transcriptome database, and 19 sequences were found to perfectly match the transcriptome sequences. The mapped RNAs were able to fold into hairpin structures and had negative folding free energies (from −21.0 to −95.65 kcal mol −1 with an average of about −49.05 kcal mol −1 ) according to Mfold; these values were lower than the folding free energies of rRNA (−33 kcal mol −1 ) and tRNA (−27.5 kcal mol −1 ) 50 . Therefore, these mapped RNAs were identified as candidate novel miRNAs (Table 2). Interestingly, some miRNA-3p sequences, such as mul-miR160b-3p, mul-miR166h-3p and mul-miR169p-3p, were identified, but the corresponding miRNA-5p sequences were not detected. This may indicate that these miRNA-3ps may be the authentic miRNAs. The size distribution of the small RNA sequences   Table 3. Expression profiling of novel miRNAs in mulberry phloem saps. IPS, phloem sap sampled from infected trees. HPS, phloem sap sampled from healthy trees.   identified in the two libraries showed that more than 80% of the mapped small RNAs were 20-24 nt long, with 24 nt and 21 nt as the major size groups (Fig. 2).

Expression profiling of miRNAs in response to phytoplasma-infection. The profiles of some miR-
NAs differed between the healthy and infected phloem sap libraries. A total of 30 known miRNAs and 13 novel miRNAs were found to be differentially expressed between phytoplasma-infected and healthy mulberry phloem sap libraries (Tables 1, 3). These differentially expressed miRNAs were considered phytoplasma-responsive miR-NAs, among which 15 miRNAs decreased and 28 miRNAs increased significantly in the infected phloem saps (P < 0.05, fold 2.0). These phytoplasma-response miRNAs consist of not only highly expressed miRNAs such as mul-miR166a, mul-miR166h-3p, mul-miR2199, and mul-miR5813 but also low-abundance miRNAs such as mul-miR160a, mul-miR167d-5p, mul-miR3630-3p, mul-miR3630-5p, and mul-miR396a. Although both the miRNA-5p and miRNA-3p of some miRNAs were detected, only miRNA-5p or miRNA-3p changed during phytoplasma-infection. This strongly suggested that single-strand miRNAs, and not miRNA-5p/miRNA-3p duplexes, are the phytoplasma infection-relevant molecular species. To validate the miRNA expression differences revealed by the sequencing experiments, we performed RT-qPCR analysis for 12 known miRNAs and 10 novel miRNAs covering different expression patterns. The results obtained by RT-qPCR showed a very strong correlation with read frequencies, demonstrating that our sequencing data are quantitative and reliable (Fig. 3).
Phytoplasma-responsive miRNAs related to diverse biologic processes. To understand the biological functions of phytoplasma-responsive miRNAs reported here, the target genes of these miRNAs were predicted and subjected to Gene Ontology (GO) analysis. Together, 131 target genes of 36 phytoplasma-responsive miRNAs were predicted and classified into nine categories according to their ontologies in Arabidopsis, based on KEGG functional annotations (Tables 4, 5; Fig. 4). The categories included metabolic process, transcription regulation, signalling pathway, stress and environmental response, development, response to hormones and hormone metabolism. Interestingly, 16% of target genes were found to be associated with signalling pathways. This indicated that many miRNAs in the phloem sap were associated with signal transduction and that diverse signalling pathways were involved in phytoplasma infection in the infected plants. In addition, many genes homologous to the sequences of unknown functions were predicted as targets of phytoplasma-responsive miR-NAs. Further analyses of these genes and miRNAs may reveal new biological functions for phloem. Thus, these phytoplasma-responsive miRNAs in the phloem saps might be related to diverse biologic processes, and the regulatory networks involved in the response to phytoplasma-infection in mulberry were intricate.

Mul-miR482a-5p accumulated in phloem sap under phytoplasma-infection is mobile. It was
reported that some miRNAs accumulated in the phloem sap are translocatable 51 . Our data showed that the mul-miR482a-5p accumulated strongly in the mulberry phloem sap under phytoplasma infection, so it is reasonable to suspect that it was mobile in the phloem sap. To investigate whether mul-miR482a-5p was mobile in the phloem, we performed grafting experiments using mul-miR482a-5p-overexpression and hen1-1 mutant Arabidopsis thaliana. After the establishment of graft unions, different parts of the successful grafts were used to analyse the mul-miR482a-5p abundance by RT-qPCR. As expected, the translocation of mul-miR482a-5p from overexpressing scions to hen1-1 rootstocks was observed in various independently grafted plants both with and without scions suffering Pseudomonas syringae pv. tomato DC3000 (Pst. DC3000) infection. However, little mul-miR482a-5p was detected in the hen1-1 scions grafted with mul-miR482a-5p overexpressing rootstock for scions with or without Pst. DC3000 infection (Fig. 5). This result indicates that mul-miR482a-5p can be transported efficiently across the graft junction from scions to rootstock under both infective and uninfective conditions. However, mul-miR482a-5p can scarcely be transported in the opposite direction.
To provide further evidence that the mul-miR482a-5p is mobile, the promoter of MUL-MIR482A-5p was cloned and fused to the reporter gene encoding β-glucuronidase (GUS) to analyse the expression pattern of mul-miR482a-5p in various tissues. Staining results showed that GUS activity was predominantly observed in stems and flowers, and the reporter signal in roots was very low, with no signal detected in the leaves and siliques (Fig. 6A). Corresponding to the staining results, when the pri-mul-miR482a transcript was examined in mulberry, more pri-mul-miR482a transcript was detected in the stem than in the root (Fig. 6B). However, this did not result in higher accumulation of mature mul-miR482a-5p in the stem than the root. In contrast, more mature mul-miR482a-5p accumulated in the root than the stem (Fig. 6C). The opposite accumulation of pri-mul-miR482a and mature mul-miR482a-5p in stems and roots implies that MUL-MIR482A was highly expressed in stem and mature mul-miR482a-5p was transported to roots. However, the abundance of mature mul-miR482a-5p increased in the infected roots compared to healthy roots (Fig. 6C), and the abundance of mul-miR482a-5p primary transcripts did not increase in the infected roots (Fig. 6B). Therefore, mul-miR482a-5p, but not its precursors, is mobile within the phloem. Although the level of mul-miR482a-5p increased in the infected stem barks, the abundance of mul-miR482a-5p did not increase in the infected leaves (Fig. 6C). This may confirm that mul-miR482a-5p can be transported from upper to lower parts but can scarcely be transported in the opposite direction.
Mul-miR482a-5p accumulated in phloem sap has physiological functions. The regulator of chromosome condensation family protein (RCC1) gene, trehalose 6-phosphate synthase gene and inositol 1,3,4-trisphosphate 5/6-kinase family protein gene were predicted as targets of mul-miR482a-5p and were experimentally verified by 5′-RLM RACE analyses (Fig. 7). To examine whether the translocation of mul-miR482a-5p in the phloem had physiological functions, the expression levels of its target genes in leaves, phloem saps, and roots were analysed by RT-qPCR (Fig. 8). The results showed that the expression levels of the three target genes have SCientifiC REPORtS | (2018) 8:812 | DOI:10.1038/s41598-018-19210-7 no significant change between infected and healthy leaves. This was consistent with the level of mul-miR482a-5p, which did not differ between infected and healthy leaves. The expression levels of the three target genes showed no significant change between infected and healthy phloem sap, although the level of mul-miR482a-5p was significantly increased in the infected phloem sap. This may be because there was no RNase, which was necessary for the cleavage of target mRNAs. Therefore, mul-miR482a-5p may have no physiological function that directs cleavage of its target mRNAs in the phloem sap. However, the expression of all three target genes was down regulated in the infected roots compared to healthy roots, and this coincided with the changes in mul-miR482a-5p in the roots. Thus, the translocation of mul-miR482a-5p in the phloem might have physiological cleavage functions for its target genes in the roots.
To further explore the physiological functions of mul-miR482a-5p, one of the target genes, RCC1, was also cloned and transformed into Arabidopsis thaliana, and RCC1-overexpressing and wild type Arabidopsis thaliana were inoculated with Pst. DC3000. The results showed that the RCC1 transgenic Arabidopsis plants showed stronger resistance to Pst. DC3000 than wild-type plants, suggesting that the RCC1 gene may be a positive regulator of defence responses (Fig. 9). Therefore, mul-miR482a-5p may repress the expression of RCC1 in infected mulberry and reduce host resistance to biotic stress.

Discussion
Component complexity of miRNAs in phloem sap. Although a number of miRNAs have been previously identified in the phloem sap of several herbaceous plants, to the best of our knowledge, no data on phloem miRNAs is available for woody perennials except for apple. Our study demonstrates that mulberry phloem sap contains small RNAs that are major contributors to the phloem sap RNA population. Some highly expressed miR-NAs such as miR166, miR167, and miR172 identified in mulberry phloem sap were also detected in the phloem sap of B. napus 37 , apple (M. domestica "Royal Gala") 39 , and pumpkin (Cucurbita maxima) 40 , suggesting that some miRNAs were conserved across plant species. However, some miRNAs, such as miR171, which were detected in mulberry, B. napus 37 39 , and pumpkin 40 phloem sap but were not present in mulberry. In addition, some novel miRNA candidates were identified in our data. Therefore, the compositions of phloem sap miRNAs differed between herbaceous and woody plants and even among woody species. Furthermore, our data showed that 43 miRNAs were differentially expressed in mulberry phloem sap in response to phytoplasma-infection. It was also reported that the composition of phloem sap miRNAs differed under different nutrient stresses. Thus, the miRNA composition of phloem sap is complex, and identification and characterization of the phloem miRNAs of mulberry may enhance current knowledge of the miRNA composition of phloem sap and help to discover new candidates that have significant action on phloem functions.

MiRNAs in phloem sap and leaves are distinct in complement and expression pattern in response to phytoplasma infection.
Plant miRNAs often show differential expression among various tissues 39 , and many miRNAs present in the phloem sap of mulberry were identified in this study. When the identified miRNAs were compared to previously published collections of miRNAs in the leaves of mulberry in the response to phytoplasma infection 31 , there were 52 miRNAs identified in the phloem sap but not the leaves, and 134 miRNAs were identified in the leaves but not the phloem sap (Fig. 10). Among the 53 miRNAs common to the leaves and phloem sap, the relative levels of expression of some miRNAs, such as mul-miR2199, mul-miR2916, mul-miR5813 and mul-miR6300, were high in the phloem sap but low in the leaves. This demonstrates that  phloem sap contains a specific set of miRNAs distinct from leaves and that a set of phloem-enriched sRNAs exists. Moreover, among the 43 phytoplasma-responsive miRNAs identified in phloem sap, only 10 miRNAs were expressed differently in healthy and infected leaves. This may be because not all miRNAs in phloem sap can be translocated, and mobile miRNAs might be translocated in different directions. Therefore, the expression pattern of miRNAs in phloem sap was distinct from leaves in response to phytoplasma infection, and different miRNAs might have distinct localizations and functions. Interestingly, the 10 common phytoplasma-responsive miRNAs between phloem sap and leaves showed the same expression pattern. These miRNAs could potentially act as a long-distance information transmitters in response to phytoplasma infection in mulberry. Further experiments are required to uncover the translocatability and functions of these miRNAs in modulating the response of mulberry to phytoplasmas.
Role of miRNAs in phloem sap. Although many miRNAs have been detected in many plant phloem saps, only a few miRNAs have been shown to translocate between cells and over long distances [42][43][44]47 . It is not clear whether all differentially expressed miRNAs in the phloem sap are mobile. Since phytoplasmas are restricted to sieve elements of phloem tissues, the sieve elements may experience phytoplasma infection earlier than other tissues. Thus, the translocated phytoplasma-responsive miRNAs in phloem sap might have roles in long-distance signalling in response to phytoplasma infection and serve to coordinate physiological responses with other plant parts that are not yet infected. Even if some differentially expressed miRNAs in the phloem sap were not (B and C) Measurement of pri-mul-miR482a (B) and mature mul-miR482a-5p (C) in various tissues of mulberry, respectively. Pri-mul-miR482a and mul-miR482a-5p abundance were detected by RT-qPCR, and relative abundance was evaluated using comparative Ct method using actin (Accession No. DQ785808) and U6 as the reference, respectively. Values are given as the mean ± SD of three experiments in each group.
mobile, it was reported that the miRNAs in the phloem sap may act to prevent translation and movement of their target mRNAs 39 . Although the putative target genes were bioinformatically predicted, the role of differentially expressed miRNAs to prevent translation and movement of their target mRNAs, which may be involved in signalling in response to phytoplasma infection, remains to be studied. The elucidation of the roles of these phytoplasma-responsive miRNAs in the phloem sap may reveal the mechanisms underlying the interactions of phytoplasma and mulberry. The most typical symptoms of phytoplasma diseases indicate perturbations in plant hormonal balance [52][53][54][55] . In this study, we also found expression changes for several phloem sap miRNAs involved in auxin signalling and auxin metabolism, e.g., differentially expressed mul-miR319a was predicted to target the MYB transcription factor, which redirects auxin signal transduction by interacting with ARFs and plays a role in plant hormone responses 27 . Differential mul-miR1223e was predicted to target genes coding for O-fucosyltransferase family protein and SAUR-like auxin-responsive protein, which are associated with auxin metabolism (Table 4). In addition to the auxin signalling pathway, mul-miR1223e was found to target the genes involved in salicylic acid signalling, and mul-miR157a was predicted to target the gene coding galactose oxidase/kelch repeat superfamily protein associated with brassinosteroid biosynthesis. Meanwhile, mul-miR529b was found to target the gene involved in abscisic acid and jasmonic acid signalling pathways. In addition, mul-miR894 and mul-miRn25-5p were found to be associated with the ethylene signalling pathway, and mul-miRn23-5p was predicted to target the lateral root primordium (LRP) protein-related gene and SOB five-like 1 gene, which was related to the gibberellic acid-mediated signalling pathway and cytokinin metabolic processes. These data are consistent with  Symptoms induced in infected plants suggest that phytoplasma infection may modulate developmental processes within the plant host 56 . Our data showed that some phytoplasma-responsive miRNAs target the transcription factors involved in development. For example, mul-miR156a was predicted to target the Squamosa Promoter-Binding Protein-Like (SPL) family, which plays important roles in flower and fruit development, plant architecture and phase transition in plants 57,58 . Meanwhile, mul-miR319a was predicted to target the MYB transcription factor, which was crucial to the control of proliferation and differentiation in a number of cell types and key factors in regulatory networks controlling development 59 . In addition, our results showed differentially expressed miRNAs-such as mul-miR1223e, mul-miR165b-5p, mul-miR529, mul-miRn23-5p-that target the genes associated with modulating plant development through various pathways. The changes in these miRNAs in the phloem sap may disorder the expression of many genes involved in diverse development processes, causing symptoms of phytoplasma disease in the infected plants.
As intracellular parasites, phytoplasmas have lost many metabolic genes and must obtain essential metabolites from their hosts 60 , which has a great impact on the metabolome of infected plants 7,14,52 . Our results showed that phytoplasma infection alters the profiles of a number of miRNAs involved in metabolism in phloem sap. These differentially expressed miRNAs target genes associated with protein metabolism, CHO metabolism DC3000 in infected Arabidopsis leaves. Bacterial numbers were calculated at 3 days after inoculation and represented as CFU per gram leaf tissue, and CFU of Pst. DC3000 in infected leaves was counted in a 1/1000fold bacterium solution. Bioassays were performed three times, each with three replicates, and each value is mean ± SD of three experiments. Asterisks indicate significant difference based on Student's t-test (**P < 0.01). WT, Wild type Arabidopsis Col-0; OE, Transgenic RCC1 Arabidopsis plants. and lipid metabolism. Furthermore, some miRNAs targeting the genes associated with secondary metabolism were also differentially expressed. For example, mul-miR156a and mul-miR529-3p were predicted to target UDP-glucosyltransferase genes, which were involved in the flavonoid biosynthesis pathway 61 . Therefore, differentially expressed miRNAs in the phloem sap may disturb many metabolic processes and have a significant effect on the response against phytoplasmas.
Some phytoplasma-responsive miRNAs modulate overlapping signalling of biotic and abiotic stresses. Plants have evolved sophisticated mechanisms to sense and respond to diverse biotic stresses 62 .
Our miRNA expression analysis showed many differentially expressed miRNAs targeting genes associated with defence response genes. The changes in these miRNAs may down-regulate or up-regulate their target gene expression and alter plant resistance to phytoplasma. Interestingly, we also found that some differentially expressed miRNAs were responsive to abiotic stresses. For example, mul-miR397a were predicted to target the gene casein kinase II beta chain 3, which has roles in response to light stimulus response 63 , and mul-miR529b, which was predicted to target the 3-ketoacyl-CoA synthase 19 gene involved in response to cold stress 64 . In addition, mul-miR529-3p was predicted to target the ARM repeat superfamily protein associated with salt stress and shoot gravitropism 65 . Meanwhile, several phytoplasma-responsive miRNAs were reported to be differentially expressed under various abiotic stresses in other plant species. For instance, miR156, miR157, and miR390 were reported to be responsive to salt, drought, cold, and heat stress in many other plants 29,62 and were detected to be differentially expressed in phloem sap in response to phytoplasma infection in this study. This may be because phytoplasma infection had a great influence on the growth and development of mulberry, and these miRNAs may contribute to modulation of the necessary growth and developmental adjustments to adapt to phytoplasma-infected conditions. In conclusion, our results suggested that phytoplasma infection may cause both biotic and abiotic stress in the mulberry, and some phytoplasma-responsive miRNAs involved in both biotic and abiotic stress signalling converge upstream of phytoplasma infection. The infected plant can modulate protective responses by controlling the abundance of these miRNAs via overlapping signalling. However, further experiments are required to uncover signalling in the phloem sap that modulates the response of mulberry to phytoplasmas.

Mul-miR482a-5p might negatively regulate mulberry resistance to phytoplasma-infection.
The miR482 superfamily is a group of plant-specific miRNAs targeting the NBS-LRR gene, and several reports have demonstrated that the miR482-NBS-LRR regulatory loop is part of the immune response induced by pathogens [66][67][68][69] . Moreover, miR482 was also found to participate in guidance for the biosynthesis of secondary pha-siRNAs, which are involved in controlling immune-response genes 66,70 . To date, miR482 has been confirmed to be distributed in more than 20 species, and approximately sixty primary transcripts of miR482 have been identified. Among the primary transcripts identified, there are approximately 20 primary transcripts that can be processed to generate both miR482-5p and miR482-3p (http://www.mirbase.org). To date, there are many reports of miR482-3p, but the function of miR482-5p is still not well understood. Plants infected with pathogens showed a reduced level of miR482 and an increased level of miR482 target mRNAs, suggesting that the miR482-mediated silencing cascade is suppressed by pathogen attack and may be a defence response of plants 71 . However, the level of mul-miR482a-5p, not mul-miR482a-3p, was changed significantly in the phloem sap infected by phytoplasma (Table 1). Mul-miR482a-5p was predicted to target the RCC1 gene, which is the guanine nucleotide exchange factor for the nuclear GTP binding protein Ran, and is probably involved in various biological processes, but the role of this gene under stress is currently not clear 72 . Our results showed that the RCC1 gene may be a positive regulator of defence responses (Fig. 9). Therefore, the up-regulation of mul-miR482a-5p may repress the expression of RCC1 in the infected plant and reduce host resistance to phytoplasma. This is consistent with the report that RCC1 family proteins were down-regulated in the incompatible interaction between soybean and Phytophthora sojae 73 . It was suggested that nucleocytoplasmic trafficking plays an essential role in the expression of disease resistance 74 , and nuclear localization of some disease resistance (R) proteins, such as members of CC-NB-LRR and TIR-NB-LRR proteins, is essential for their resistance function 75,76 . As RCC1 plays a major role in nucleocytoplasmic transport 72 , silencing of RCC1 may result in partial impairment of nucleocytoplasmic trafficking and loss of resistance function of some disease resistance (R) proteins. Thus, when mulberry plants were infected by phytoplasma, the mul-miR482a-3p did not decrease, suggesting that the miR482-NBS-LRR regulatory loop was not induced. Moreover, increased mul-miR482a-5p might repress the resistance to phytoplasma mediated by RCC1. Therefore, phytoplasma-derived suppression of RNA silencing may repress the whole host immune system during infection and potentially enhance phytoplasma colonization and amplification. However, further research is required to elucidate the regulatory mechanisms of the RCC1 gene.
In conclusion, the characterization of miRNA-Seq-based expression profiling of miRNAs allowed for the identification of many phytoplasma-responsive miRNAs in mulberry phloem sap. Future investigation will explore the functions and regulatory networks of these miRNAs. The information provided here will be particularly useful for a complete understanding of the function of miRNAs in phloem sap and will lay the foundation to reveal the mechanisms underlying phytoplasma pathogenicity.

Methods
Plant material. One-year-old cutting seedlings collected from Husang 32 (Morus multicaulis Perr.) were used as rootstock and grafted to scions from healthy or phytoplasma-infected mulberry trees (Husang 32). All establishment graft unions were incubated in a greenhouse, and plants showing Witches' broom disease symptoms were confirmed by PCR assay with an amplified fragment of the phytoplasma 16S rRNA gene (GenBank Accession No. EF532410) using the primers (P16mF: 5′-TAAAAG ACCTAGCAATAGG-3′ and P16mR: 5′-CAATCCGAACTGAGACTGT-3′) as previously described 53 . Phloem saps collection and purity assessing. Phloem sap was collected from infected and healthy mulberry plants using the shoot exudation method 77 . First, shoots were excised with a sterile razor blade, the first droplets were discarded, and the cut surface was blotted with sterile filter paper several times to avoid contamination. Exuding phloem sap was collected into an Eppendorf tube containing TRIzol reagent (Invitrogen, Carlsbad, CA, USA). Leaf and phloem sap RNAs were isolated using a TRIzol kit (Invitrogen) following the manufacturer's instructions. To detect RuBisCO and MmPP16 transcripts, cDNA was synthesized using oligo (dT)18 primer (GACTCTAGACGACATCGA(T)15) and ReverTra Ace M-MLV RTase (TaKaRa, Dalian, China). RT-PCR assays were performed in 25-µl reaction volumes containing 20 ng cDNA and 150 nM forward and reverse primers. The primers used for amplification of RuBisCO and MmPP16 genes are shown in Supplementary Table 1.
Small RNA library construction and high-throughput sequencing. Isolated phloem sap RNAs were used to prepare a small RNA library according to the protocol of the TruSeq Small RNA Sample Prep Kits (Illumina, San Diego, USA). Single-end sequencing (36 bp) was performed on an Illumina Hiseq2500 instrument following standard protocols. Three independent libraries each (biological replicates) were analysed for infected and healthy phloem saps.
MiRNAs identification. The raw sequences tags obtained from HiSeq sequencing were cleaned to remove adapter dimers, junk, low complexity, common RNA families (rRNA, tRNA, snRNA, snoRNA) and repeats, and the length distribution of the clean tags was summarized. The trimmed reads longer than 18 nt were annotated into different categories, and the sequences of 18-25 nucleotides were compared to a miRBase database v21.0 (http://www.mirbase.org/). The sequences with identical or related sequences from other plants were regarded as conserved miRNAs. All remaining unmapped sequences were BLASTed against our mulberry transcriptome database, and the hairpin RNA structures containing sequences were predicted using RNAfold software and used to predict novel miRNAs using Mireap (http://sourceforge.net/projects/mireap/). Differential expression analyses of miRNAs. The frequency of miRNA was normalized by the total number of miRNAs in every sample, where normalized expression = (Actual miRNAs sequencing reads count/ Total clean reads count) × 1,000,000. The fold change between infected (IPS) and healthy phloem sap (HPS) was calculated as follows: fold-change = log2(IPS/HPS). Statistical analysis was performed according to Poisson distribution, and the P value was calculated based on the formula: A fold-change ≥ 2 and P ≤ 0.05 were used as criteria to identify differentially expressed miRNAs, and an miRNA was designated as significantly differentially expressed if its expression value varied more than two-fold and P ≤ 0.05 between infected and healthy phloem saps. Target prediction of differential miRNAs. The target genes of the differentially expressed miRNAs were predicted using the software psRNATarget (http://plantgrn.noble.org/psRNATarget/) by submitting miRNA sequences to a search against our in-house mulberry transcriptome data following the criteria of (i) maximum expectation less than 3.0; (ii) multiplicity of target sites 2; (iii) range of central mismatch for translational inhibition 9-11 nt; and iv) SCientifiC REPORtS | (2018) 8:812 | DOI:10.1038/s41598-018-19210-7 maximum mismatches at the complementary site ≦ 4 without any gaps. All predicted target genes were aligned with the reference Arabidopsis thaliana database downloaded from TAIR (http://www.arabidopsis.org/; version TAIR10) to annotate their functions, and the GO terms of these targets were also annotated based on their TAIR GO categories.
Quantitative real-time PCR analysis for miRNAs and mRNAs. RNA was extracted using the TRIzol ® reagent following the manufacturer's recommendations (Invitrogen, Carlsbad, CA, USA) and digested with DNase I. Real-time PCR analyses for miRNAs and mRNAs were performed using the PrimeScriptTM miRNA qPCR Starter Kit Ver.2.0 (TaKaRa, Dalian, China) and the SYBR Premix Ex TaqTM kit (TaKaRa, Dalian, China) on the Rotor-Gene 3000A system (Bio-Rad, Munich, Germany), respectively, according to the manufacturer's protocol for the Rotor-Gene 3000A system. The U6 and actin genes were used as reference genes for miRNA and mRNA normalization, respectively. The U6 gene was amplified using the primer (5′-ATGGCCCCTGCGTAAGGATG-3′), and actin was amplified using primer pair (F: 5′-CAGTGCTTCTCACTGAGGCTC-3′ and R: 5′-GGAAGAGGACTTCTGGGCATC-3′). The primers (Supplementary Tables 1, 2) used to amplify the genes and miRNAs were designed based on our available mulberry transcriptome data. The relative expression levels of miRNA and mRNA were evaluated using the Comparative cycle threshold (Ct) method 78 . All samples were assayed in triplicate.
Target validation. For miRNA target validation, a modified gene-specific 5′ RNA ligase-mediated rapid amplification of cDNA ends (5′ RLM-RACE) was performed using the GeneRacer Kit (Invitrogen, Carlsbad, CA, USA). Total RNA was isolated from mulberry seedlings using the TRIzol kit (Invitrogen, Carlsbad, CA, USA) following the manufacturer's instructions, ligated to the RNA oligo adapter (5′-CGACUGGAGCACGAGGACACUGACAUGGACUGAAGGAGUAGAAA-3′) and reverse transcribed with SuperScript III reverse transcriptase using oligo(dT) primer. The resulting cDNA was PCR-amplified with GeneRacer 5′ primer (5′-CGACTGGAGCACGAGGACACTGA-3′) and each respective gene-specific outer primer (shown in Supplementary Table 3). The PCR product was further amplified by nested PCR using GeneRacer 5′ nested primer (5′-GGACACTGACATGGACTGAAGGAGTA-3′) and each respective gene-specific inner primer (shown in Supplementary Table 3). The final PCR product was gel-purified and cloned into a pMD18-T vector (Invitrogen, Carlsbad, CA, USA) for sequencing.
Micrografting experiments. Four-day-old seedlings of hen1-1 mutant and transgenic Arabidopsis thaliana were used for micrografting experiments. The seedlings were cut transversely using a sterile razor blade and combined inside silicon tubing (0.3 mm internal diameter). The graft unions were cultured on 1.5% (w/v) agar plates with half-strength MS medium for 9 days and were hydroponically cultured for 10 d. Grafted plants without adventitious roots were selected, and the scions and stocks of the selected grafted plants were harvested for RNA isolation. Then, the RNAs were used for RT-qPCR of the miRNA. To investigate whether the miRNA moves during infective conditions, the leaves of the scions of graft unions were spray-inoculated with Pst. DC3000 suspensions at 10 8 CFU mL −1 in 10 mM MgCl 2 with 0.04% (v/v) Silwet L-77. Three days after inoculation, the RNAs of the scions and stocks of the grafted plants were isolated and used for RT-qPCR of the miRNA.
Promoter activity analysis. The promoter was cloned using a Tail-PCR method and ligated into the vector pBI121 to replace the cauliflower mosaic virus (CaMV) 35S promoter and fused to the GUS (β-glucuronidase) reporter gene to create the promoter expression vector pMIR482::GUS. The derived construct vector was introduced into Agrobacterium tumefaciens strain GV3101, and the WT Arabidopsis plants were transformed with a floral dipping method. Histochemical staining for GUS activity was performed referring to the previously described method 79 .
Detection of resistance against Pst. DC3000. The RCC1 gene coding sequence was cloned and ligated into the vector pBI121, and the derived construct vector was introduced into A. tumefaciens strain GV3101 under the control of 35S. The WT Arabidopsis plants were transformed with floral dipping method. After transformation, the transformed plants were selected. Four-week-old transgenic and wild-type Arabidopsis seedlings were spray-inoculated with Pst. DC3000 suspensions at 10 8 CFU mL −1 in 10 mM MgCl 2 with 0.04% (v/v) Silwet L-77 or vacuum-infiltrated with bacterial suspensions at 10 7 CFU mL −1 with a syringe. Three days after inoculation, disease symptoms were recorded using a camera, and bacterial numbers were calculated and represented as colony-forming units (CFU) per gram leaf tissue in a 1/1000-fold bacterium solution. Bioassays were performed three times, with three replicates each. Data Availability. The datasets analyzed during the current study are available from the corresponding author on reasonable request.