Unidirectional movement of small RNAs from shoots to roots in interspecific heterografts

Long-distance RNA movement is important for plant growth and environmental responses; however, the extent to which RNAs move between distant tissues, their relative magnitude and functional importance remain to be elucidated on a genomic scale. Using a soybean (Glycine max)–common bean (Phaseolus vulgaris) grafting system, we identified 100 shoot–root mobile microRNAs and 32 shoot–root mobile phased secondary small interfering RNAs (phasiRNAs), which were predominantly produced in shoots and transported to roots, and, in most cases, accumulated to a level similar to that observed in shoots. Many of these microRNAs or phasiRNAs enabled cleavage of their messenger RNA targets or phasiRNA precursors in roots. In contrast, most mobile-capable mRNAs were transcribed in both shoots and roots, with only small proportions transported to recipient tissues. These findings suggest that the regulatory mechanisms for small RNA movement are different from those for mRNA movement, and that the former is more strictly regulated and, probably, more functionally important than the latter. A study using soybean–common bean grafting systems reveals that mobile sRNAs are predominantly produced in shoots and accumulate in roots by long-distance transport, while mobile mRNAs are produced in both shoots and roots, with a small proportion moved.

H igher plants possess two specialized vascular tissues, the xylem and phloem, which conduct water and nutrients from the roots to the shoots and organic components produced by photosynthesis in the leaves throughout the plants as their respective basic functions. Unlike the xylem, which is composed of primarily dead cells, the phloem consists of columns of living cells, including sieve elements and companion cells connected by plasmodesmata, allowing molecules such as carbohydrates, proteins/ amino acids and hormones to travel locally or over long distances [1][2][3][4][5] . The phloem also harbours a variety of RNAs such as mRNAs, microRNAs (miRNAs) and small interfering RNAs (siRNAs) [1][2][3]6,7 . Early studies demonstrated that several mRNAs produced in leaves were transported to distal tissues to exert physiological functions underlying specific traits [8][9][10][11][12] . Recent genomic analyses of transcripts in heterografted plants identified hundreds to thousands of mRNAs involved in shoot-to-root or vice versa trafficking [1][2][3]13 . In the past decade, several miRNAs have been identified as systemic signals mediating long-distance communications [14][15][16][17] . For example, miR399 in Arabidopsis thaliana and miR395 in Brassica rapa move from shoots to roots in response to nutrient deficiencies, while miR2111 in Lotus japonicus travels from shoots to roots to regulate the susceptibility of uninfected roots to the soil bacteria rhizobia as a mechanism to balance rhizobial infection and nodulation events [18][19][20][21] . More recently, a genomic analysis of grafted tissues identified numerous transposon-derived or heterochromatic siRNAs (hcsiRNAs) that are capable of trafficking from shoots to roots to modulate genome-wide DNA methylation in the recipient root cells 22,23 .
Despite such progress, surprisingly, genome-scale identification and characterization of mobile miRNAs-important regulators of gene expression-have not, to our knowledge, been conducted in any organisms. PhasiRNAs, a subclass of siRNAs that require a trigger miRNA for their biogenesis, were known to be capable of cell-to-cell trafficking [24][25][26] , but whether they move over long distances is yet to be determined. Shoots-to-roots mobile hcsiRNAs have been identified at the whole-genome level 22,23 , but whether hcsiRNAs also move from roots to shoots to regulate DNA methylation has not been investigated. Long-distance mobile mRNAs have been identified on a genome scale in several plants 3 , but the abundance of mobile mRNAs in recipient tissues relative to those produced locally is unclear. Here, we integrate small RNA (sRNA)-sequencing (sRNA-seq), Illumina short-read RNA-sequencing (RNA-seq), Nanopore long-read RNA-seq and degradome-sequencing (degradome-seq) data from heterografted and homografted soybean (G. max) and common bean (P. vulgaris)-two economically important leguminous crops diverged from a common ancestor ~17 million years ago-to address these outstanding questions.
Mobile miRNAs are predominately synthesized in shoots and accumulate in roots via long-distance movement. A total of 622 miRNAs, including miRNA-5p and miRNA-3p, were identified in the eight samples. These include 161 Gm-miRNAs, 72 Pv-miRNAs and 389 miRNAs sharing identical sequences between soybean and common bean (Supplementary Table 4). Of the 161 Gm-miRNAs, 67 were detected in the Gm/Pv shoots but not in the Pv/Gm shoots ( Fig. 2a and Supplementary Table 4), and thus defined to be mobile only from shoots to roots; one was detected in both the Gm/Pv roots and Pv/Gm shoots and defined to be mobile bidirectionally; and the remaining 93 were absent in either the Gm/Pv roots or Pv/Gm shoots and thus considered to be immobile. Of the 72 Pv-miRNAs, 33 were detected in the Pv/Gm roots but not in the 10 Gm-hcsiRNAs  showing sequence variation and srNAs showing 'same' sequences between the two crops. c, Numbers of mobile soybean (blue) and common bean (orange) srNAs and directionality of movement. d, Abundance of mobile soybean (Gm-) and common bean (Pv-) srNAs ranging from 16 to 26 nucleotides (nt). e, Comparison of Gm-hcsirNA abundances between homografted soybean roots and heterografted common bean roots. f, Comparison of Pv-hcsirNA abundances between homografted common bean roots and heterografted soybean roots. g, Comparison of Gm-hcsirNA abundances between homografted and heterografted soybean roots. h, Comparison of Pv-hcsirNA abundances between homografted and heterografted common bean roots. The specific tissue in two grafted samples is underlined below using the shoot/root notation. Boxplots display the distribution of the log 2 -transformed abundances of hcsirNAs in individual samples. The minima, maxima, centre, bounds of box and whiskers, and percentiles are marked in boxplots.
Gm/Pv shoots (Fig. 2b and Supplementary Table 4), and thus defined to be mobile from shoots to roots, whereas 39 were absent in either the Pv/Gm roots or the Gm/Pv shoots and thus considered to be immobile (Supplementary Table 1). Of the 67 shoot-to-root mobile Gm-miRNAs and 33 shoot-to-root mobile Pv-miRNAs, only 4 (6.0%) and 14 (42.4%) were detected in the Pv/Gm and Gm/Pv roots, respectively, and overall their abundances were much lower than those in the Gm/Gm and Pv/Pv roots, respectively (Fig. 2c,d  and Supplementary Tables 3 and 4). When all the soybean-common bean distinguishable (mobile and putative immobile) miRNAs were included, the overall abundance of Gm-miRNAs in the Pv/Gm roots was much lower than that in the Gm/Gm roots, and, similarly, the overall abundance of Pv-miRNAs in the Gm/Pv roots was much lower than that in Pv/Pv roots (Extended Data Fig. 1c,d and Supplementary Table 3). Together, these observations indicate that the mobile miRNAs were predominantly produced in shoots, and that the miRNA biogenesis machinery was largely suppressed in roots in general, but not fully blocked. The mobility of several miR-NAs was validated by stem-loop PCR (Extended Data Fig. 3). Of the mobile miRNAs detected in this study, six (miR166i, miR1509a, miR1510a, miR1510b, miR5770a and miR5770b) were detected to have moved from shoots to roots in both soybean (as Gm-miRNAs) and common bean (as Pv-miRNAs) (Supplementary Table 5).
With a few exceptions that may reflect the specificity of miRNA biogenesis and mobility, the relative abundances of the mobile Gm-miRNAs in the Gm/Pv roots and Pv-miRNAs in the Pv/Gm roots were similar to those observed in the Gm/Gm and Pv/Pv roots, respectively (Fig. 2a, Table 3). The relative abundances of the mobile Gm-miRNAs in the Gm/Pv roots and Pv-miRNAs in the Pv/Gm roots were also correlated with those in the Gm/Pv and Pv/Gm shoots, respectively (Extended Data Fig. 4a,b and Supplementary Table 3). In addition, the relative abundances of the mobile Gm-miRNAs in the Gm/Gm shoots and Pv-miRNAs in the Pv/Pv shoots were correlated with those in the Gm/Pv and Pv/ Gm shoots, respectively (Extended Data Fig. 4c,d and Supplementary  Table 3). Together, these observations suggest that the type of grafting, whether heterografting or homografting, did not affect the biogenesis and mobility of the miRNAs in shoots, and that the mobile miRNAs accumulated in roots substantially via long-distance movement.

b and Supplementary
In contrast to the mobile Gm-miRNAs, most of which were not produced in roots (Fig. 2c), their precursors, Gm-MiRNAs, were expressed at similar levels between the Pv/Gm and Gm/Gm roots (Fig. 2e and Supplementary Table 3). In addition, the expression levels of the Gm-MiRNA in the Pv/Gm roots were highly corelated with those in the Gm/Pv shoots (Extended Data Fig. 4e and Supplementary Table 3). This contrasting pattern was also observed between the abundance of Pv-miRNAs and the levels of Pv-MiRNA expression in the common bean tissues (Fig. 2d,f, Extended Data Fig. 4f and Supplementary Table 3), reinforcing our deduction that the miRNA biogenesis machinery was largely suppressed in the root tissues (Fig. 2g).
Mobile miRNA-mediated systemic gene regulation through mRNA cleavage in recipient cells. In an attempt to understand whether the mobile miRNAs execute a regulatory role in recipient cells, we generated and analysed degradome data and RNA-seq data from the grafted samples. Of the 67 mobile Gm-miRNAs and 33 mobile Pv-miRNAs, 32 and 10 were detected to enable the cleavage of their putative mRNAs in soybean, or common bean, or both (Supplementary Table  6). Three soybean genes, whose mRNAs were cleaved by corresponding mobile Pv-miRNAs in the Pv/Gm roots, showed reduced levels of expression in the Pv/Gm roots compared with the Gm/Gm roots, in which such cleavages were not detected. Four common bean genes, whose mRNAs were cleaved by corresponding mobile Gm-miRNAs in Gm/Pv, exhibited reduced levels of expression in the Gm/Pv roots compared with the Pv/Pv roots, in which such cleavages were not detected ( Fig. 2h-j). In addition, seven soybean genes, whose mRNAs were cleaved by corresponding Gm-miRNAs, exhibited reduced expression in the Gm/Gm roots compared with the Pv/Gm roots, in which no cleavages of those mRNAs were detected ( Fig. 2h-j). One of the mobile Gm-miRNAs, gma-miR4415a/b-3p, which was produced in soybean shoots, moved to the Gm/Gm and Gm/Pv roots, and was able to target soybean gene Glyma.20G051900 and common bean gene Phvul.006g011700, respectively, resulting in reduced expression of respective targets in the Gm/Gm and Gm/Pv roots compared with the Pv/Gm and Pv/Pv roots, respectively (Fig. 2i).

Mobile phasiRNAs are predominately synthesized in shoots and accumulate in roots substantially via long-distance movement.
To determine whether phasiRNAs are long-distance mobile, we first identified all PHAS loci using the degradome-seq and sRNA-seq data generated in this study. A total of 29 soybean PHAS loci and 13 common bean PHAS loci were detected to have generated pha-siRNAs. Of all the phasiRNAs generated from the 42 loci, 61 were Gm-phasiRNAs, 18 were Pv-phasiRNAs and 25 were not distinguishable between soybean and common bean. Of the 61 Gm-phasiRNAs, 23 were detected in the Gm/Pv roots but not in the Pv/Gm shoots, and thus defined to be mobile only from shoots to roots, and 38 were absent in either the Gm/Pv roots or the Pv/Gm shoots and thus considered to be immobile ( Fig. 3a and Supplementary Table 7). Of the 18 Pv-phasiRNAs, 14 were detected in the Pv/Gm roots but not in the Gm/Pv shoots, and thus defined to be mobile only from shoots to roots, and four were absent in either the Pv/Gm roots or the Gm/Pv shoots and thus considered to be immobile ( Fig. 3b and Supplementary Table 7). Of the 23 shoot-to-root mobile Gm-phasiRNAs and 14 shoot-to-root mobile Pv-phasiRNAs, only one (4.3%) and four (28.6%) were detected in the Pv/Gm and Gm/Pv roots, respectively, and overall their abundances were lower than those in the Gm/Gm and Pv/Pv roots, respectively (Fig. 3c,d and Supplementary Table 7). When all the soybean-common bean distinguishable (mobile and putative immobile) phasiRNAs were included, the overall abundance of Gm-phasiRNAs in the Pv/Gm roots was much lower than that in the Gm/Gm roots, and, similarly, abundances between homografted soybean roots and heterografted common bean roots. b, Comparison of Pv-mirNA abundances between homografted common bean roots and heterografted soybean roots. c, Comparison of Gm-mirNA abundances between homografted and heterografted soybean roots. d, Comparison of Pv-mirNA abundances between homografted and heterografted common bean roots. e, Comparison of the expression levels of Gm-MiRNAs between homografted and heterografted soybean roots. f, Comparison of the expression levels of Pv-MiRNAs between homografted and heterografted common bean roots. The specific tissue in two grafted samples is underlined below using the shoot/root notation. Boxplots display the distribution of the log 2 -transformed abundances of mirNAs or MiRNAs in individual samples. The minima, maxima, centre, bounds of box and whiskers, and percentiles are marked in boxplots. g, Modes of mobile mirNA biogenesis and systemic regulation of target genes. h, Graphic demonstration of mobile mirNA-mediated downregulation of target genes in recipient tissues revealed by comparison of mrNA abundances between homografted and heterografted roots of a same crop. i, Exemplification of mobile mirNA-mediated downregulation of target genes in recipient tissues; for example, downregulation of Pv-mrNAs (or Gm-mrNAs) by mobile Gm-mirNAs (or Pv-mirNAs) in the Gm/Pv (or Pv/Gm) roots (type 1), and downregulation of Gm-mrNAs by mobile Gm-mirNAs in the Gm/Gm roots (type 2). The curved lines with arrowheads connect mirNAs and respective target genes. j, Exemplification of confirmed cleavages (sites and frequencies indicated by arrows and ratios, respectively) of mrNAs by mobile mirNAs in recipient tissues of heterografted and homografted plants.
the overall abundance of Pv-phasiRNAs in the Gm/Pv roots was much lower than that in Pv/Pv roots (Extended Data Fig. 1e,f and  Supplementary Table 3). Together, these observations indicate that the mobile phasiRNAs were predominantly produced in shoots, and that the phasiRNA biogenesis machinery was largely suppressed in roots in general, but not fully blocked.  Phvul.006G011700

gma-miR4996
Glyma.13G183500 Overall, the relative abundances of the mobile Gm-phasiRNAs in the Gm/Pv roots and Pv-phasiRNAs in the Pv/Gm roots were similar to those observed in the Gm/Gm and Pv/Pv roots, respectively ( Fig. 3a,b and Supplementary Table 3). The relative abundances of the mobile Gm-phasiRNAs in the Gm/Pv roots and Pv-phasiRNAs in the Pv/Gm roots were also correlated with those in the Gm/ Pv and Pv/Gm shoots, respectively (Extended Data Fig. 5a,b and Supplementary Table 3). In addition, the relative abundances of the mobile Gm-phasiRNAs in the Gm/Gm shoots and Pv-phasiRNAs in the Pv/Pv shoots were correlated with those in the Gm/Pv and Pv/Gm shoots (Extended Data Fig. 5c,d and Supplementary Table 3), respectively. Together, these observations suggest that the type of grafting, whether heterografting or homografting, did not affect the biogenesis and mobility of the phasiRNAs in shoots, and that the mobile phasiR-NAs accumulated in roots substantially via long-distance movement.

Pv
Although most of the mobile Gm-phasiRNAs and Pv-phasiRNAs were not produced in roots (Fig. 3c,d), their precursors, Gm-PHAS loci and Pv-PHAS loci, were expressed at similar levels between the Pv/Gm and Gm/Gm roots and between the Gm/Pv and Pv/Pv roots (Fig. 3e,f), respectively. The expression levels of the Gm-PHAS loci in the Pv/Gm roots and the Pv-PHAS loci in the Gm/Pv roots were also highly correlated with those in the Gm/Pv and Pv/Gm shoots (Extended Data Fig. 5e,f and Supplementary Table 3), respectively. Biogenesis of phasiRNAs from their PHAS loci requires trigger miRNAs. As exemplified in Fig. 3g, gma-miR1510b-3p was detected to have triggered the production of Gm-phasiRNAs from the transcripts of Glyms.04G219600 and Glyma.15G232600 in soybean shoots (Gm/Gm or Gm/Pv), but no such Gm-phasiRNAs were detected in the Pv/Gm roots, in which both gma-miR1510b-3p and the transcripts of two soybean genes were present (Fig. 3g,h). Together, these observations reinforce our deduction that pha-siRNA biogenesis was severely suppressed in the root tissues. Glyma.04G219600 Glyma.15G232600 Gmax. Gmax. Gmax. Gmax. Gmax. Gmax. Gmax. Gmax. Gmax.

Gmax.
Glyma.16G215400 Glyma.16G137200  Fig. 3 | Relative abundance and systemic gene regulation of mobile phasiRNAs. a, Comparison of Gm-phasirNA abundances between homografted soybean roots and heterografted common bean roots. b, Comparison of Pv-phasirNA abundances between homografted common bean roots and heterografted soybean roots. c, Comparison of Gm-phasirNA abundances between homografted and heterografted common bean roots. d, Comparison of Pv-phasirNA abundances between homografted and heterografted common bean roots. e, Comparison of the expression levels of Gm-PHAS loci between homografted and heterografted soybean roots. f, Comparison of the expression levels of Pv-PHAS loci between homografted and heterografted common bean roots. The specific tissue in two grafted samples is underlined below using the shoot/root notation. Boxplots display the distribution of the log 2 -transformed abundances of phasirNAs or PHAS transcripts in individual samples. The minima, maxima, centre, bounds of box and whiskers, and percentiles are marked in boxplots. g, Exemplification of regulatory cascades involving mobile phasirNAs. gma-mir1510b-3p triggered production of a cluster of phasirNAs from Glyma.04G219600 in soybean shoots, three of which enabled cleavage of their own precursor, with one validated to have implemented the cleavage in the Gm/Gm roots but not in the Pv/Gm roots. Cleavage site and frequency are indicated by arrows and ratios, respectively. h, relative abundances of gma-mir1510b-3p, all phasirNAs produced by the PHAS loci and the PHAS loci in the eight tissues. The curved lines with arrowheads connect the mirNA and its PHAS targets, PHAS loci and phasirNAs from the loci, or phasirNAs and their respective target genes.
targets of mobile phasiRNAs based on the degradome-seq data. Three of the mobile phasiRNAs produced from their precursor, the transcripts of Glyms.04G219600, enabled cis-directed cleavage of the precursor in the Gm/Gm roots. One of the three cleavage sites was also further confirmed by 5ʹ rapid amplification of complementary DNA ends (RACE)-PCR (Fig. 3g,h). By contrast, such cleavage sites were not detected in the Pv/Gm roots, exemplifying mobile phasiRNA-mediated systemic regulation of their own precursor genes in the recipient tissues.
According to the degradome-seq data, mRNAs from 19 genes were detected to be the targets of 17 phasiRNAs including four mobile Gm-phasiRNAs, one mobile Pv-phasiRNA, seven putative immobile Gm-phasiRNAs and five phasiRNAs that were not distinguishable between the two plants (Supplementary Table 8). Of the four mobile Gm-phasiRNAs, one was able to target two soybean genes (Glyma.16G050500, Glyma.19G100200) in the Gm/Gm roots and a common bean gene (Phvul.001G087000) in the Gm/Pv roots, and three were detected to target four soybean genes in the Gm/Gm roots. The mobile Pv-phasiRNA was detected to target a common bean gene (Phvul.001G087000) in the Pv/Pv roots and a soybean gene (Glyma.19G100200) in the Pv/Gm roots. As exemplified in Fig. 3h, some phasiRNAs were able to target not only their precursor PHAS loci, but also non-phasiRNA-producing genes.
Most mobile-capable mRNAs are transcribed in both shoots and roots, but mobile mRNAs make up only small proportions of the total mRNA abundance of a given gene in shoots or roots. In an attempt to assess the relative completeness of mobile mRNAs, we decoded the transcriptomes of the eight grafted samples by Illumina short-read RNA-seq and Oxford Nanopore Technology (ONT) long-read RNA-seq, which together detected expression of 38,241 soybean genes and 24,062 common bean genes. The Illumina RNA-seq data revealed the mobility of 1,322 mobile soybean mRNAs and 874 common bean mRNAs. Of the 1,322 soybean mRNAs, 1,167 moved from shoots to roots, 130 moved from roots to shoots and 25 moved bidirectionally. Of the 874 common bean mRNAs, 684 moved from shoots to roots, 153 moved from roots to shoots and 37 moved bidirectionally ( Fig. 4a and Supplementary  Tables 1 and 9). The ONT RNA-seq data revealed the mobility of 163 soybean mRNAs and 129 common bean mRNAs, of which ~56% were detected to cover the full length of respective coding sequences (Supplementary Table 9). Of the 163 soybean mRNAs, 117 moved from shoots to roots, 35 moved from roots to shoots and 11 moved bidirectionally. Of the 129 common bean mRNAs, 90 moved from shoots to roots, 29 moved from roots to shoots and 10 moved bidirectionally. In total, five pairs of orthologous genes in soybean and common bean were detected to have produced mobile mRNAs in both species. The trafficking of mRNAs from a number of genes including Glyma.11G114000, Glyma.11G228000, Phvul.001G259000, Phvul.001G113800, Phvul.001G229500, Phvul.010G023000 and Phvul.008G285100, which are predicted to be involved in various signalling pathways underlying plant growth, photosynthesis and/or vesicle-mediated transport, was further confirmed by quantitative PCR with reverse transcription (RT-qPCR) (Extended Data Fig. 6). Eight mobile soybean mRNAs and 14 mobile common bean mRNAs were detected by both Illumina RNA-seq and ONT RNA-seq.
It is proposed that soybean underwent a whole-genome duplication event ~13 million years ago, leaving ~16,500 duplicated gene pairs retained in the current genome 27 . We found that, of the 1,477 mobile mRNAs in soybean, 300 were from singletons, 1,075 were from one copy of duplicated gene pairs and 51 were from both copies of duplicated gene pairs (Supplementary Table 10). Thus, no apparent biases of the mobility against either duplicates or singletons were observed (3.6% in duplications versus 2.0% in singletons). Among all mobile mRNAs detected in soybean and common bean, only 72 were from orthologous genes in the two species, although enriched biological pathways of the mobile mRNAs in the two species are similar (Extended Data Fig. 7 and Supplementary Table 11).
To understand the relative magnitude of mobile mRNAs and the efficiency of mRNA movement, we examined and compared the relative abundances of mobile Gm-mRNAs and Pv-mRNAs in the eight grafted samples (Fig. 4b,c and Supplementary Table 9). The relative abundance of the shoot-to-root mobile Gm-mRNAs detected in the Gm/Pv roots was extremely low compared with that detected in their source tissue, the Gm/Pv shoots; and the relative abundance of the root-to-shoot mobile Gm-mRNAs detected in the Pv/Gm shoots was extremely low compared with their source tissue, the Pv/Gm roots. These observations indicate that only small proportions of the mobile-capable Gm-mRNAs produced in their source (shoot or root) tissues were transported to their respective recipient (root or shoot) tissues. Consistent with these observations, most of the shoot-to-root mobile-capable Gm-mRNAs showed higher or similar levels of abundance between the Gm/Gm and Pv/Gm roots, and most of the root-to-shoot mobile-capable Gm-mRNAs showed higher or similar levels of abundance between the Gm/Gm and Gm/Pv shoots. Similar observations were obtained for the mobile Pv-mRNAs (Fig. 4b,c and Supplementary Table 9). Together, these results suggest that the contribution of the mobile mRNAs from the source tissues to the recipient tissues was minimal. This scenario contrasts with that for the mobile sRNAs described earlier under three categories (that is, hcsiRNAs, miRNAs and pha-siRNAs), and further illustrated as a single sRNA category (Fig. 4b,c and Supplementary Tables 2, 4 and 7), which demonstrate extensive accumulation of mobile sRNAs in their recipient root tissue.

Discussion
Accumulating evidence has demonstrated the importance of plant mobile RNAs as signal molecules in shoot-root communications 3 , but no previous studies, to our knowledge, have investigated the mobility of multiple types of RNAs simultaneously at the whole-genome level, perhaps partly due to the limitations of the individual grafting systems used in those studies. For example, grafting of the wild-type Arabidopsis shoots with the roots of Arabidopsis mutants lacking the functional RNA polymerase IV required for 24-nucleotide sRNA biogenesis identified mobile 24-nucleotide hcsiRNAs 29,30 , but it would be incapable or less capable of detecting most of the other types of mobile RNAs, such as miRNAs and mRNAs, on a genomic scale. Similarly, grafting of different Arabidopsis ecotypes was able to identify a large number of mobile mRNAs 7 , but it would be less effective or ineffective in detecting mobile sRNAs showing identical sequences between the ecotypes. As plant mutants blocking miRNA biogenesis such as the Arabidopsis dcl1 mutants are generally lethal or infertile 31 , it would be impractical to graft such mutants with wild-type plants for genome-wide identification of mobile miRNAs. Although a number of heterografting systems with highly diverged plant species were used to identify genome-wide mobile mRNAs 3 , the effectiveness of those systems in identifying mobile sRNAs was not explored. In our study, the high degree of sequence divergence between the two leguminous crops enabled identification of 5,232 mobile hcsiRNAs, 100 mobile miRNAs and 32 mobile phasiRNAs, as well as 2,466 mobile mRNAs, throughout the two genomes. Given that less than a dozen mobile miRNAs and phasiRNAs have been previously reported 3 , the effectiveness of our experiments in identifying mobile sRNAs, particularly miRNAs and phasiRNAs, is laudable. Notably, for all types of mobile sRNAs investigated in this study, particularly the Gm-sRNAs, their relative abundances in the heterografted roots (Gm/Pv) were correlated with those in the homografted (Gm/Gm) roots, and their relative abundances in heterografted (Gm/Pv) shoots were also correlated with those in the homografted (Gm/Gm) shoots, suggesting that heterografting did not affect their biogenesis or their mobility any more than did homografting. In other words, the picture of RNA movement revealed by utilizing the heterografted plants reflects what occurs in homografted or natural plants.
We would like to point out that our approach also has its own limitations. Because only one sample of each tissue was sequenced, without biological replicates, there remains uncertainty as to whether those putative 'immobile' RNAs are truly immobile if replicates are included or more sequencing data are generated. In addition, it was only able to determine the mobility of 233 miR-NAs, including miRNA-5p and miRNA3-p, that are unambiguously assigned to the soybean or common bean genomes, which account for ~37% of all the miRNAs detected in the investigated shoot and root tissues. Whether this subset of miRNAs is representative of all miRNAs expressed in the two crops remains unclear. Nevertheless, several individual miRNAs showing sequence polymorphisms between soybean and common bean and shoot-to-root trafficking are also shared by other legumes, in which they have been demonstrated to play important roles in root and nodule development. For example, in Medicago truncatula, miR166 and miR1509 were found to regulate the root and nodule development by targeting mRNAs encoding HD-ZIP transcription factors and by triggering the production of phasiRNAs from transcripts of an APETALA2 homologue, respectively 32,33 . In addition, miR1510 was found to be a young miRNA specific to the Phaseoleae tribe of legumes and a predominant trigger for the production of phasiRNA from transcripts of NB-LRR genes that may underlie plant responses to biotic stresses 32,34,35 . Thus, although these miRNA sequences are diverged among species, they appear to be functionally conserved among these legumes and execute their functional roles systemically through long-distance trafficking.
It is noteworthy that our study revealed a few contrasting features between mobile sRNAs and mobile mRNAs. One of such features is the directionality of RNA movement. Previously, an elegant Arabidopsis shoot-root grafting experiment was designed to trace shoot-root movement of a green fluorescent protein (GFP)-derived hairpin RNA inverted repeat using GFP as a reporter 22 . This experiment demonstrated the movement of the GFP-derived silencing signal from shoot to root and vice versa, with the former more efficient than the latter. Nevertheless, the ≥15-nucleotide GFP-derived sRNAs transported to the roots accounted for only 0.06-0.08% of all of those detected in the shoots; the 21-22-nucleotide GFP-derived siRNAs transported to the roots were three orders of magnitude less abundant than those observed in the shoots 22 . The low efficiency of movement of this GFP-derived sRNA may not be representative of all mobile sRNAs present in the genome. Thus, whether sRNAs move from root to shoot remained largely elusive, although a few studies had shown movement of a silencing signal from root to shoot 22,36 . We demonstrate that the mobile sRNAs, including hcsiRNAs, miRNAs and phasiRNAs, were capable of moving only from shoots to roots, with few exceptions (5 of 11,291, <0.05%). Since root-shoot grafts carry a piece of stem from the rootstock, it is unclear whether such exceptions were contributed by that piece of stem rather than the root. In contrast to the nearly exclusive shoot-to-root movement of sRNAs, ~14% of the mobile mRNAs were detected to move from roots to shoots. Movement of mRNAs from shoots to roots and vice versa was commonly observed in other grafting experiments, but the numbers of mobile mRNAs as well as the proportions of mobile mRNAs trafficking in the two directions vary greatly among different experiments 37 . Such variation may be associated with growth conditions and developmental stages of the grafted plants, the time points after grafting for tissue collection and sequencing coverage 3 .
Another feature is the distinct efficiencies of mRNAs and siR-NAs. In most cases, the mobile sRNAs transported to and accumulated in the recipient tissues (that is, roots) were highly abundant and not produced locally, although the machinery for sRNA biogenesis was not fully blocked. By contrast, the proportions of mRNAs transported to and accumulated in the recipient tissues (either shoots or roots) were detected to be minimal (Fig. 4b,c). This may partly explain why only a small portion of mobile mRNAs was detected by both the Illumina sequencing and ONT sequencing of mRNAs, although such a low level of overlap may also be partly due to the uneven and biased coverages of species-specific short reads along individual genes-a pattern that was also revealed in a previous study 7 , as well as the relatively low coverages of the ONT long reads. In addition, the capability of the applied ONT RNA-seq approach to assess the integrity of mobile mRNAs remains limited, often with technical effects from RNA-seq sample preparation and library construction 38 as well as possible turnover of mRNAs in living cells 39 . Nevertheless, the majority (~56%) of the mobile mRNAs detected by the long reads cover the full length of their respective coding sequences. Given that some long-distance mobile mRNAs have been demonstrated to enable synthesis of functional proteins 3 , it would be reasonable to speculate that some of the mobile mRNAs are in full length and translatable. However, because mobile mRNAs transported to the recipient tissue of a homografted plant can be produced locally, the biological importance of mobile mRNAs remains to be elucidated. Together, these contrasting features suggest that the long-distance trafficking of sRNAs is more strictly regulated than that of mRNAs and that the mobile sRNAs may be more functionally important than mobile mRNAs.

Methods
Plant materials and growth conditions. Seeds of soybean variety Williams 82 and common bean variety Jiulibai were soaked in water for 1 d and then sowed in soil to grow in the growth chamber under the cycles of 12 h light at 28 °C/12 h darkness at 24 °C, with humidity at 30%. The stems of the seedlings were cut 7 d (V0 stage) after sowing for the grafting experiments following a protocol described previously 40 . The grafting experiments include four distinct scion/rootstock combinations-soybean/common bean and common bean/soybean heterografted plants, and soybean/soybean and common bean/common bean homografted plants. As the connection of phloem generally takes place 3 d after initial grafting and the connection of xylem takes up to 7 d post-grafting in Arabidopsis 29 , we sampled the eight tissues homografted G. max (soybean) shoots (Gm/Gm), homografted G. max roots (Gm/Gm), heterografted G. max shoots (Gm/Pv), heterografted G. max roots (Pv/Gm), homografted P. vulgaris (common bean) shoots (Pv/Pv), homografted P. vulgaris roots (Pv/Pv heterografted P. vulgaris shoots (Pv/Gm) and heterografted P. vulgaris roots (Gm/Pv) 10 d after the grafting. Each of the eight samples processed for sRNA and mRNA sequencing was composed of the same tissue from four individual plants in each grafting experiment without replicates.
Sequencing of sRNA, mRNA and degradomes, and processing of raw data. Total RNAs were isolated from grafted shoot and root samples using TRIzol Reagent (Invitrogen/Life Technologies). The sRNA-seq libraries were constructed using the NEBNext Multiplex Small RNA Library Prep Set for Illumina (NEB) following the manufacturer's manual, and then sequenced using the Illumina HiSeq 2500 platform to generate the 50-base-pair (bp) single-end reads. The raw reads were processed with the fastx-toolkit (v.0.0.14, http://hannonlab.cshl.edu/fastx_toolkit) for removal of low-quality reads and adaptor sequences, and the processed reads were mapped to the soybean and common bean reference genomes (v.12.1, phytozome) using the Bowtie program 41 (v.1.1.1) with 0 mismatches (−v 0). The abundance matrix of sRNAs was normalized to c.p.m. (counts per million reads) based on the mapping results using in-house Perl scripts.
Total RNAs were processed using the Epicentre Ribo-zero rRNA Removal Kit (Epicentre) to deplete ribosomal RNAs, and the processed RNA samples were used to construct RNA-seq libraries using the NEBNext Ultra Directional RNA Library Prep Kit for Illumina (NEB). Then the RNA libraries were sequenced using the Illumina HiSeq 4000 platform to generate 150-bp paired-end reads. We developed a more conservative approach to determine RNA mobility which relied less on SNP detection than a previous study 7 . Briefly, the raw reads were processed with the fastx-toolkit program for removal of low-quality reads and adaptor sequences, and the processed reads from each library were subsequently mapped to the soybean or common bean reference genome 27,28 using STAR 42 (v.2.5.4b) with the default parameters. The reads uniquely mapped to individual genes were extracted using Samtools 43 (v.1.8), and then the abundance of mRNAs from each gene was determined and normalized to c.p.m. using the Bedtools program 44 (v.2.29.0), and differentially expressed genes (DEGs) between compared samples were defined using the 'EdgeR' packages in R 45 (v.3.28), based on the recommendation by its User's Guide for experiments without biological replicates. In an attempt to increase the accuracy in finding DEGs, we used a dispersion value of 0.2, which is much larger than the dispersion value previously estimated based on 40 highly confident soybean housekeeping genes 46 , and an adjusted false discovery rate-controlled Q value (P value) <0.05 to identify candidate DEGs. The DEG candidates targeted by a number of selected miRNAs and phasiRNAs were further validated by RT-qPCR in this method.
To detect long mobile RNA fragments, ONT sequencing of cDNAs with a strand-switching method was employed using the cDNA-PCR Sequencing Kit (SQK-PCS 109) and PCR Barcoding Kit (SQK-PBK004), following the manuals of the manufacturer (Oxford Nanopore Technologies). The base-called reads were collected and converted to FASTA format using the MinKNOW software (https://nanoporetech.com). The Nanopore reads were mapped to the soybean and common bean genomes using the minimap2 program 47 (v.2.11) with the parameters (-ax splice -uf -k14-secondary = no-splice-flank = no). The reads that mapped to both the genomes were excluded from further analyses, and only those mapped uniquely to one of the two genomes were processed for the subsequent analyses. Gene models were annotated de novo using the Stringtie program with the '-L' parameter specified for the long-read data 48 . The newly generated annotation files were merged and compared with the gene models defined in the soybean reference genome 27,28 using the Cuffmerge and Cuffcompare programs 49 , respectively. All the mobile mRNAs were further confirmed by alignment with the Illumina short reads, including both mobile and putative immobile ones generated in this study, using the Integrative Genomics Viewer (IGV) software 50 manually. The gene expression level was quantified based on the Illumina short reads in the same tissues using the Stringtie program.
Degradome-seq, also referred as to parallel analysis of RNA ends (PARE), is a modified 5ʹ RACE with high-throughput short-read sequencing method for target mRNA confirmation and cleavage site detection, and was conducted following the protocol described previously 51 with modifications made by LC-BIO (Hangzhou, China). The 47-bp single-end PARE sequences were sequenced using the Illumina HiSeq 2500 platform. A combination of the degradome-seq data, with the transcriptome data and the sRNA reads from soybean and common bean, was used to detect cleavage sites within mRNAs using the Paresnip2 program 52 (v.4.5). All kinds of sequencing information are summarized in Supplementary Table 12.
Sequence analyses for detection of mobile RNAs. sRNAs with sizes ranging from 16 nucleotides to 26 nucleotides and copy numbers ≥1 c.p.m. in at least one of the eight samples were kept for detection of mobile sRNAs. The 24-nucleotide sRNAs were defined as hcsiRNAs; the miRNAs were identified by searching against miRbase, the miRNA database 53 , with a focus on previously identified soybean and common bean miRNAs 33 . The newly identified miRNA-5p and miRNA-3p were considered to be different miRNAs. The phasiRNA precursor genes and phasiRNAs were identified using the PhaseTank program 54 (v.1.0) with two reference genomes (soybean and common bean) and other default parameters. Putative miRNA triggers for production of phasiRNAs were predicted using degradome data and psRNAtarget-a plant sRNA target analysis server 55 . mRNAs were identified by comparison with annotated genes in the soybean and common bean genomes 27,28 . The mobility of the sRNAs that are either soybean-or common bean-specific, or show sequence polymorphisms between the two crops, was determined by comparison among the heterografted samples. For example, a soybean sRNA was defined to be a shoot-to-root mobile sRNA when it was detected in Gm/Pv roots but not in the Pv/Gm shoots; a soybean sRNA was defined to be a root-to-shoot mobile sRNA when it was detected in the Pv/Gm shoots but not in the Gm/Pv roots; a soybean sRNA was defined to be mobile bidirectionally when it was detected in both the Gm/Pv roots and the Pv/Gm shoots; and a soybean sRNA was considered to be immobile when it was detected in neither the Gm/Pv roots nor the Pv/Gm shoots. The mobility of mobile mRNAs from annotated genes in the two reference genomes was identified using the same principle described above, but only individual Illumina reads that were uniquely mapped to only one of the two genomes were included in identification of mobile mRNAs and evaluation of their relative abundance. For the dataset from ONT long-read RNA-seq, only the reads uniquely mapped to one of the two genomes were aligned with the Illumina short reads generated from a same tissue for sequence confirmation and evaluation of relative expression of the mobile genes.
RT-qPCR, stem-loop RT-qPCR, and 5ʹ RACE-PCR and sequencing, and analysis of PCR fragments. RT-PCR, RT-qPCR, stem-loop RT-PCR, stem-loop RT-qPCR and 5ʹ RACE-PCR were performed as previously described 56,57 . In addition to RNA-seq, the relative abundance of mRNA from chosen genes was evaluated by RT-qPCR, in which the soybean gene GmCons4 (GenBank ID BU578186) 58 or common bean gene PvActin11 (Phvul.008G011000) 59 was used as a respective internal reference to quantify the relative expression levels of the soybean and common bean genes from three biological replicates. All the primers used in this study are listed in Supplementary Table 13.
Gene set enrichment and pathway analysis. Gene ontology enrichment analysis of a given gene group, such as DEGs and mobile mRNAs, was performed and visualized with the clusterProfiler R package 60 (v.3.10). The clusterProfiler's input files were formatted using in-house Perl scripts. Gene ontology term information of soybean and common bean was extracted from the gene annotation files in the Phytozome database. The gene ontology terms with false discovery rate-adjusted P < 0.05 were considered as significant enrichment.
Reporting Summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
The raw read sequences are deposited in the National Center for Biotechnology Information Sequence Read Achieve (http://www.ncbi.nlm.nih.gov/sra/) under the accession number PRJNA648759. The miRNAbase 53 database (Release 22.1, 2018) was used for identifying known soybean and common bean miRNAs and their precursors.

Code availability
Custom Perl codes used for data analysis in this paper can be found at https:// github.com/wang3283/sRNA-normalize-and-GO-anlaysis-input.

nature research | reporting summary
April 2020 Corresponding author(s): Jianxin Ma Last updated by author(s): Nov 30, 2020 Reporting Summary Nature Research wishes to improve the reproducibility of the work that we publish. This form provides structure for consistency and transparency in reporting. For further information on Nature Research policies, see our Editorial Policies and the Editorial Policy Checklist.

Statistics
For all statistical analyses, confirm that the following items are present in the figure legend, table legend, main text, or Methods section.

n/a Confirmed
The exact sample size (n) for each experimental group/condition, given as a discrete number and unit of measurement A statement on whether measurements were taken from distinct samples or whether the same sample was measured repeatedly The statistical test(s) used AND whether they are one-or two-sided Only common tests should be described solely by name; describe more complex techniques in the Methods section.
A description of all covariates tested A description of any assumptions or corrections, such as tests of normality and adjustment for multiple comparisons A full description of the statistical parameters including central tendency (e.g. means) or other basic estimates (e.g. regression coefficient) AND variation (e.g. standard deviation) or associated estimates of uncertainty (e.g. confidence intervals) For null hypothesis testing, the test statistic (e.g. F, t, r) with confidence intervals, effect sizes, degrees of freedom and P value noted

Software and code
Policy information about availability of computer code Data collection All sofeware including fastx-toolkit (version 0.0.14), MinKNOW, minimap2 program47 (version 2.11) used for data collection in the manuscript have been cited/referenced.

Data analysis
All sofeware including Bowtie program41 (version 1.1.1), STAR42 (version 2.5.4b), Bedtools program44 (version 2.29.0), "EdgeR" packages in R45 (version 3.28), Integrative Genomics Viewer (IGV), Paresnip2 program52 (version 4.5), PhaseTank program54 (version 1.0), clusterProfiler R package60 (version 3.10) used for the data analyses in the manuscript have been cited/referenced. For manuscripts utilizing custom algorithms or software that are central to the research but not yet described in published literature, software must be made available to editors and reviewers. We strongly encourage code deposition in a community repository (e.g. GitHub). See the Nature Research guidelines for submitting code & software for further information.

Data
Policy information about availability of data All manuscripts must include a data availability statement. This statement should provide the following information, where applicable: -Accession codes, unique identifiers, or web links for publicly available datasets -A list of figures that have associated raw data -A description of any restrictions on data availability All the raw read sequences generated in this study have been deposited in the National Center for Biotechnology Information Sequence Read Achieve (http:// www.ncbi.nlm.nih.gov/sra/) under the accession number PRJNA648759