Comparative Profiling of miRNAs and Target Gene Identification in Distant-Grafting between Tomato and Lycium (Goji Berry)

Local translocation of small RNAs between cells is proved. Long distance translocation between rootstock and scion is also well documented in the homo-grafting system, but the process in distant-grafting is widely unexplored where rootstock and scion belonging to different genera. Micro RNAs are a class of small, endogenous, noncoding, gene silencing RNAs that regulate target genes of a wide range of important biological pathways in plants. In this study, tomato was grafted onto goji (Lycium chinense Mill.) to reveal the insight of miRNAs regulation and expression patterns within a distant-grafting system. Goji is an important traditional Chinese medicinal plant with enriched phytochemicals. Illumina sequencing technology has identified 68 evolutionary known miRNAs of 37 miRNA families. Moreover, 168 putative novel miRNAs were also identified. Compared with control tomato, 43 (11 known and 32 novels) and 163 (33 known and 130 novels) miRNAs were expressed significantly different in shoot and fruit of grafted tomato, respectively. The fruiting stage was identified as the most responsive in the distant-grafting approach and 123 miRNAs were found as up-regulating in the grafted fruit which is remarkably higher compare to the grafted shoot tip (28). Potential targets of differentially expressed miRNAs were found to be involved in diverse metabolic and regulatory pathways. ADP binding activities, molybdopterin synthase complex and RNA helicase activity were found as enriched terms in GO (Gene Ontology) analysis. Additionally, “metabolic pathways” was revealed as the most significant pathway in KEGG (Kyoto Encyclopedia of Genes and Genomes) analysis. The information of the small RNA transcriptomes that are obtained from this study might be the first miRNAs elucidation for a distant-grafting system, particularly between goji and tomato. The results from this study will provide the insights into the molecular aspects of miRNA-mediated regulation in the medicinal plant goji, and in grafted tomato. Noteworthy, it would provide a basis how miRNA signals could exchange between rootstock and scion, and the relevance to diverse biological processes.

Introduction

Grafting is considered as one of the most successful techniques for plant and crop production, and protection. Literally, grafting is the union of a bud or stem or even a branch onto the root or stem or branch of another plant, while it is termed as a “distant-grafting” when it is carried out among different families, genus or species.

In recent years, distant-graft mutagenesis technology developed instantaneously and proved its ability to bypass the natural obstacles of incompatibility between distantly-related species. Notably, to integrate the rootstock of different sources and scion organically thus allow the scion to grow, develop and fruit normally with desired traits (Pan et al., 2012). Though it was thought initially that the genetic materials of rootstock integrate into the scion's genome creates the genetic variation induced by grafting (Taller et al., 1998), however, recent molecular biology research has identified that the stress-related retro-transposons have undergone transposition due to distant-grafting mechanism (Wei-Min, 2005), which is recognized as one of the effective mechanisms of genome rearrangement and gene mutation. Similarly, some other studies confirmed that the nucleic acid materials of rootstock can transmit through the graft-junction to the scion, where small RNA might silence, some special genes function or regulate expression of genes in scion (Stegemann and Bock, 2009). Undoubtedly, grafting confers a new physiological event to the involved plants though the molecular processes are relatively unknown. The gene expression changes occurred at transcriptional, post-transcriptional, and post-translational levels in the grafted plants. However, the mechanism of gene or loci expression regulation and the affiliation of small RNAs are largely unknown.

MicroRNAs (miRNAs) is a recently discovered small RNA class, which is short (20–24 nt), endogenous and non-coding, that regulates gene expression at post transcriptional level by cleaving or inhibiting translation of target gene transcripts (Llave et al., 2002; Joshi et al., 2010). Recent evidence indicates that they play a critical and pivotal role in the most significant biological processes and molecular pathways such as plant growth and development, including leaf morphology and polarity, lateral root formation, hormone signaling, transition between vegetative stage to flowering stage, regulation at flowering time, floral organ identity, and reproduction (Mallory and Vaucheret, 2006; Sunkar et al., 2007). Additionally, they have regulatory effects on nutrition homeostasis, abiotic stresses, and pathogen responses (Carrington and Ambros, 2003; Song et al., 2010). In plant, miRNAs are transcribed by RNA polymerase II to produce pri-miRNA (primary miRNA) transcripts and been cleaved to pre-miRNAs (miRNA precursors) followed by the catalytic activity of DCL1 (Dicer-like enzyme) (Subramanian et al., 2008; Naqvi et al., 2012). The pre-miRNA is additionally cleaved to a miRNA duplex (miRNA:miRNA*) and then one of these strand is incorporated into RISC (RNA induced silencing complex) (Szittya et al., 2008). Usually miRNA* is degraded, but also can be accumulated at a lower level in some exceptions (Wang et al., 2011). Then the remaining mature miRNA sequence guides AGO (Argonaute protein) toward the complementary target mRNA sequence resulting the silencing of the corresponding genes (Bartel, 2004; Zhu, 2008).

Recent studies of long-distance signal movement and functional analysis have proven that the macromolecules like proteins, mRNAs and miRNAs can translocate through the phloem. Moreover, these entities act as signals for the stress-related pathways, and for plant development (Kehr and Buhtz, 2008; Atkins et al., 2011). FLOWERING TIME protein is a well-studied molecule that moves from leaves to shoot apex in Arabidopsis (Corbesier et al., 2007). In few cases, movements of siRNAs (small interfering RNAs) were reported where induction of posttranscriptional gene silencing acted against viruses (Waterhouse et al., 2001). Previous report of biochemical analysis in pumpkin (Cucurbita maxima), has characterized a unique component of protein machinery CmPSRP1 (Phloem SMALL RNA BINDING PROTEIN1), which selectively binds to the small single-stranded RNA (ssRNA) species and mediates cell-to-cell trafficking, but not with dsRNA molecules (Hamilton et al., 2002). Further molecular study revealed that the endogenous CmPP16 moved from CC (companion cell) into SE (sieve element) through the pore-plasmodesmata when grafting was performed between pumpkin and cucumber (Tiedemann and Carstensbehrens, 1994).

Many studies reported on grafting, but very limited information is available on plant miRNAs mobility. Some other reviews have summarized presence of miRNAs in phloem exudates of different species, but negligible information reported on their mobility (Chuck and O'Connor, 2010; Kehr and Buhtz, 2013). However, RNA silencing signal movement was first convincingly reported in tobacco experiments (Palauqui et al., 1997). In a report of root patterning, miRNA165/166 appeared to move to adjacent cell layers from its biogenesis site. This is an example of a short-distance movement of miRNAs (Carlsbecker et al., 2010; Miyashima et al., 2011). However, recent evidence has indicated that only four miRNAs (miR399, miR395, miR172, and miR156) have been demonstrated as long-distance mobile signals in plants (Bhogale et al., 2014). Pant et al. (2008) reported that miR399 acted as a movable signal in Arabidopsis which regulated the phosphate homeostasis (Pant et al., 2008), whereas miR395 was reported to move from wild-type scions to rootstocks miRNA processing mutant hen1-1 under sulfate stress (Buhtz et al., 2010). In another study, miR172 was proposed as a long-distance mobile signal responsible for tuberization in potato (Solanum tuberosum) (Martin et al., 2009). A recent report by Bhogale et al. confirmed miR156 as a graft-transmissible miRNA which modulated potato plant architecture and tuberization in a potato micro-grafting method (Bhogale et al., 2014).

Our study was on grafting between tomato and goji. Both tomato (Solanum lycopersicum L.) and goji (Lycium chinense Mill.), both belonging in the night shade family (Solanaceae) of different genera. Goji, commercially called wolfberry/goji berry, is the common name to the fruit of two very closely related species: Lycium barbarum L. (Ningxia gouqi) and Lycium chinense Mill. (gouqi) (Flint, 1997). Tomato and goji are well recognized and high valued species for their unique and diversified uses and qualities. Goji fruit has a significant position at the glorious Chinese traditional medicine history.

According to the Chinese traditional medicine study, goji fruits are effective in nourishing liver and kidney, enhancing eyesight, enriching blood, invigorating sex, reducing rheumatism and so on. Some other important functions have been confirmed in modern clinical researches like immunity improvement, anti-oxidation, anti-aging, growth stimulation, hemopoiesis enhancing, incretion regulating, blood sugar reducing and bearing improvement. Goji is also widely used in brewing, beverage and many other food products (Dong et al., 2008; Zhao et al., 2013).

On the other hand, tomato is the world's top canned and the largest vegetable crop after potato and sweet potato (Fawusi, 1978). It is a very cheap source of vitamins and important condiment in the most diets. Along with tremendous invaluable nutritional functions of tomato, recent years, “lycopene” (a phytochemical), has drawn centric attention for its disease-fighting abilities, which is abundantly found in tomato. Studies in human have resulted that lycopene is effective at a wide range of cancers, including prostate of course, but also colorectal, breast, lung, endometrial, pancreatic, bladder, cervical, and skin cancers. Lycopene furthermore helpful to prevent heart disease and may slow the development of cataracts and muscular degeneration, an age-related vision problem that can lead to blindness (Agarwal and Rao, 2000).

By far the most conserved plant miRNAs have been identified through traditional Sanger sequencing method, but high-throughput sequencing technology is replacing Sanger and computational prediction methods vigorously. In computation-based approach, miRNAs can be identified with high conservation, but it cannot be characterized and likely have a probability to produce false-positive results (Xu et al., 2010). Recent studies reported that species-specific miRNAs of non-model plants are also very important. Moreover, non-conserved miRNAs are found at a low level of expression (Moxon et al., 2008). The deep sequencing technology can detect miRNAs with a high precision even less abundantly expressed miRNAs with a lower cost. Thus, miRNAs detection through small RNA library construction of high-throughput sequencing is considered as a feasible way in small RNA research (Song et al., 2010).

MiRNAs have been studying extensively, but this study might be the first miRNA elucidation of distant-grafting, particularly between medicinally important goji berry and tomato. However, miRNAs profiling from grafted tomato never been reported. Our present study focused on the miRNA-regulated gene expressions based on high-throughput small RNA libraries of grafted and control tomato samples. Obtained results would provide a unique basis for further unraveling the mechanism of miRNA transshipment between rootstock and scion and to build up the relevance to diverse biological processes between these plants.

Materials and Methods

Plant Material

Lycium chinense Mill. (cv. Large leave goji), was selected as a rootstock for this study and a commercial tomato cultivar, TA209 was selected as the scion. Healthy, fresh, and strong stem of goji was collected and grown in soil under control condition to develop mother stock plants. After several months, the stock plants sprouted a good number of branches and healthy leaves. Approximately, 2 months old, non-lignified and soft woody branch was chosen as rootstock. On the other hand, tomato seeds were germinated and grown under control conditions and the healthy shoot of approximately 6–7 weeks old was selected as scion.

Grafting Procedure

The stem of goji was cut about 15–20 cm above the soil level and made a vertical cut of 1–2 cm at the center. 3–5 cm stature scions were prepared by cleaning leaves and trimmed the base as a wedge followed by cleft grafting method (Sonoda and Nishiguchi, 2000). Graft junction was wrapped and covered immediate after grafting. The junction was wrapped by a polythene stripe and water was spraying three times in a day to avoid the desiccation of tender tomato shoot. Three weeks later, the bag has been removed permanently. Since the expression level of plant miRNAs changes on the spatio-temporal basis (Wang et al., 2012), we, followed the homo-grafting method in tomato to avoid the cutting/mechanical effect upon the miRNAs expression pattern. So, finally, we have produced two types of tomato plants for our small RNA study; (a) grafted tomato plant (where rootstock is goji) and (b) control tomato plant (where rootstock is a tomato itself).

Sample Tissue Selection, Total RNA Extraction, High-Throughput Sequencing and Small RNA Library Construction

For the high-throughput RNA sequencing and the small RNA library construction; we have studied four samples. Two types of samples (shoot and fruit) were studied from grafted and control tomato plants, thus four libraries were constructed such as: SFC (solanaceae fruit control); SFG (solanaceae fruit graft); SSC (solanaceae shoot control) and SSG (solanaceae shoot graft). For the shoot tissues, we have considered a shoot tip of 2 months old with the adjacent tender leaves and for fruit tissues, three fruits developmental stages (mature green, color breaking and ripening red) were considered. Target tissues were collected into sealer bags followed by snap freezing immediately in liquid N₂ and stored into liquid N₂ until crashing. The total RNAs of shoot were isolated from the shoot tip along with the adjacent leaves and for fruit, the total RNAs were isolated from the three fruit developing stages separately followed by the equal bulk mixing for the RNA sequencing and subsequent validation experiments. Target tissues were stored at −80°C ultra-refrigerator for later or further uses.

The total RNAs were isolated from target tissues by the TRIzol method, with a slight modification using RNAiso Plus (TAKARA BIO INC.) extraction kit, and the methods that were described previously (Hafner et al., 2008). RNA quantification and qualification were ensured before constructing small RNA libraries. One percent agarose gel electrophoresis was performed to monitor RNA degradation and contamination. Nano photometer (Implen, CA, USA) was used to check RNA purity. RNA concentration was measured using Qubit 2.0 flurometer (Life Technologies, CA, USA) and the integrity was assessed by Agilent bioanalyzer 2100 (Agilent Technologies, CA, USA). The Illumina sequencing was used at Novogene Company, Beijing, China. In short, the total RNAs were separated on 8% PAGE for small RNAs selection followed by ligation with Solexa 5′ and 3′ adapters sequentially. Adapter-ligated small RNAs were reverse transcribed after ligation and purification, finally PCR products were sequenced on an Illumina Hiseq 2500/2000 platform. The 50 bp single-end reads were generated.

Bioinformatics Analyses of Sequencing Data

From the sequenced, raw reads (raw data); only clean reads were used for the downstream analyses (a flowchart of the bioinformatics is included in the Figure S1). Raw reads of the fastq format were processed through Perl and python scripts. Then the reads were screened out to separate clean reads by removing contaminating reads, adapters containing sequences, insert tags, reads with poly-A tails and the sequences beyond the 18–40 nt length. Clean reads of these sRNA tags were mapped to the tomato genome (http://solgenomics.net/organism/Solanum_lycopersicum/genome) using Bowtie (Langmead et al., 2009) with no mismatches. Non-coding RNAs (rRNAs, tRNAs, snRNAs, snoRNAs) matching sequences was BLASTN searched with the sequences of Rfam (http://www.sanger.ac.uk/software/Rfam) and NCBI GenBank databases (http://www.ncbi.nlm.nih.gov/blast/Blast.cgi) and were excluded to perform the downstream analyses (Benson et al., 2011).

The remaining sRNA tags were used to search the evolutionary conserved or known miRNAs. Reads were searched against miRBase 20.0 considering maximum two mismatch (Griffiths-Jones et al., 2008). To obtain the potential miRNAs and drawing the secondary structures, a modified software, miRDeep2 (Friedlander et al., 2012) and srna-tools-cli were used. The characteristic of the hairpin structure of miRNA precursors was used to predict the potential novel miRNAs. By the integrating operation of miREvo (Wen et al., 2012) and miRDeep2 software, novel miRNAs were predicted through exploring the secondary structures, DCL1 cleavage sites and the minimum free energy (MFE) of the sRNA tags.

To predict miRNA families, we have used known miRNAs in miFam.dat (http://www.mirbase.org/ftp.shtml) and the novel miRNA precursor were aligned to Rfam (http://rfam.xfam.org/). For the prediction of potential target genes of miRNAs, psRobot software was used (Wu et al., 2012).

Quantification is very important to estimate the expression of miRNAs, so after prediction of the reads corresponding to the target genes was normalized by TPM (transcripts per million) according to Zhou et al. (2010). The normalization formula is: normalized expression = mapped read count/total reads*1,000,000.

For the detection of differential expression of the miRNAs between two samples, we have used DEGseq (2010) R package (to minimize the positive false rate for the sequencing data for without biological replicates). The corresponding P-values were adjusted using q-values (Storey, 2003). The threshold level of q-value was set as q < 0.01 and |log₂(foldchange)|>1 to estimate the significantly differentially expressed miRNAs.

The dataset studied in this article is available in the NCBI (SRA) public repository under the accession number SRP051797.

Validation of miRNA Expression Using qRT-PCR

Poly (A)-tailed RT-qPCR was used to measure the gene expression variation and validate the deep sequencing results of miRNAs. Reverse-transcription reaction and real-time PCR primer design were conducted according to previous report (Shi and Chiang, 2005) and the manufacturer's instructions (Takara-Primescript miRNA qPCR starter kit, 2.0), where each PCR reaction was performed in a volume of 25 μl containing 12.5 μl of SYBR Premix Ex Taq II (2×), 1 μl of each forward primer, 1 of μl universal reverse primer and reverse-transcribed cDNA from ~100 pg of total RNA. The PCR protocol was 5 s at 95°C, 40 cycles of 95°C for 5 s, 60°C for 20 s and 72°C for 1 min. Real- time PCR was performed on the Applied Biosystem 7500 Fast detection system using the SYBR Green I method, and all the reactions were run in biological triplicates. The melting curve was used to determine the specificity of PCR products (primer amplicons).

The delta-C_T (corresponding cycle threshold) method was used to calculate the relative expressional levels of miRNAs (Livak and Schmittgen, 2001). In this method, the C_T values of the wanted miRNAs and reference genes were first transformed for measurement using delta-C_T followed by dividing the quantities of wanted miRNAs by the geometric mean of the reference genes. The standard deviation and mean values were calculated using to triplicates RT-qPCR assays. The C_T was calculated using the machine accessory software and converted into relative copy numbers using a standard curve as previously described reports (Jones-Rhoades and Bartel, 2004; Shi and Chiang, 2005). The gene 5.8 S rRNA was used as reference gene in the qPCR detection of miRNAs. Student's t-test was used for the statistical analysis of the RT-qPCR.

Functional Assignments of Potential Target Genes

GO enrichment analysis was used for the target gene candidates of differentially expressed miRNAs. GOseq based Wallenius non-central hyper-geometric distribution (Young et al., 2010) was used in the GO enrichment analysis. KEGG database was used for further understanding of the functional diversity of “target gene candidates” (Kanehisa et al., 2008). Finally, we have used KOBAS (Mao et al., 2005) software to test the statistical enrichment of the target gene candidates in KEGG pathways.

Validation of Predicted Targets Using RLM 5′-Race

For the validation of candidate targets, modified RNA Ligase-Mediated 5′-RACE (RLM 5′-RACE) was performed. 5′-Full RACE Kit (Takara) was used according to the manufacturer's instructions, with slight modifications. Briefly, total RNA was directly ligated to the 5′ adaptor followed by reverse transcription with the oligo (dT) primers. PCR was performed with 5′ primers and 3′ gene-specific primers using the cDNA as the template (Table S1). The 5′-RACE PCR products were purified using Takara PCR product recovery kit, cloned, and sequenced.

Results

Summary of Small RNAs Generated through Deep Sequencing in Grafted and Control Tomato Plants

For the investigation of small RNAs, especially miRNAs' regulatory network in grafted tomato plant, four small RNA libraries were constructed from the shoot and fruit tissues of the control and grafted tomato plant followed by prediction of involved targets and important regulatory biological pathways. More than 35.22 million raw reads were obtained from four libraries and approximately 34.38 million (97.63%) found with high quality (Table 1).

TABLE 1

Table 1. Statistical summary of data generated by high-throughput small RNA sequencing in tomato.

After removing different adapters, poly-(A/T/G/C) sequences and reads beyond the 18–40 nt range, 31,034,648 (88.10% of total raw reads) total small RNA reads with 18,497,982 mapped sRNA sequences were obtained from the four libraries (Table 1). These sequences were comprised of known miRNAs, putative novel miRNAs, rRNAs, tRNAs, snRNAs, snoRNAs, repeat associated RNAs, TAS (Transacting small interfering RNAs or TASi) and unannotated fragments (Figure 1). These total clean sRNA tags were then mapped to the tomato genome without any mismatches by Bowtie software to analyze their expression and distribution on the reference (Langmead et al., 2009). The highest proportion of reads mapped to the tomato genome found from SSG (Solanaceae shoot graft) library 4,577,167 (70.48%), and the lowest proportion yielded from SFG library 3,669,417 (40.15%). Important observation was found from the unique sRNA reads in the grafted samples. Grafted tomato fruit has shown notably higher accumulation of unique sRNA reads over control (39.32 and 17.63% of clean sRNA reads, respectively). Similarly, grafted shoot also has demonstrated exceptionally higher frequency than control (83.72 and 48.01% of clean sRNA reads, respectively; Table 1). This phenomenon strongly supports that a tremendous changes have occurred in transcriptomic level of tomato tissues from the rootstock of goji. In this study, 76.67% unique reads were found common between SFC and SFG libraries where 20.11% were SFC specific and 3.22% were SFG specific. In SSC and SSG libraries only 10.91% were common and 52.19 and 36.89% were found to be SSC and SSG specific, respectively (Figure 2).

FIGURE 1

Figure 1. Small RNA classification. (A) Small RNAs in SFC library (Solanaceae fruit control); (B) Small RNAs in SFG (Solanaceae fruit graft); (C) Small RNAs in SSC (Solanaceae shoot control); (D) Small RNAs in SSG (Solanaceae shoot graft). Known miRNA, perfect matching mature miRNA sequences in miRBase miRNAs repository; Novel miRNA, new mature miRNA sequences detected by characteristic stem-loop hairpins and minimum free energy (MFE); rRNA, tRNA, snRNA, snoRNA are non-coding sRNAs category; TAS, transacting small interfering RNAs; Unannotated, reads that had no known information.

FIGURE 2

Figure 2. Common and specific unique sRNA sequence in SFC vs. SFG (A) and SSC vs. SSG (B). SFC, Solanaceae fruit control; SFG, Solanaceae fruit graft; SSC, Solanaceae shoot control; SSG, Solanaceae shoot graft.

For the length distribution analysis, 21–24 nt long RNA reads found significantly abundant throughout the sRNA libraries which are a typical sRNA distribution pattern of high-throughput RNA sequencing. However, 24 nt long sequences showed a remarkably higher abundance only in the shoot tissues (between control and grafted shoot) and comprised approximately 40 and 37% of total, which is four times and two times higher than their corresponding fruit tissues. However, grafted fruit has represented higher abundance (18%) in 24 nt long sequences than the control one (10%). This is an opposite picture of shoot tissues where the control tomato shoot exhibited 3% higher than the grafted shoot. This scenario suggested that 24 nt long sRNAs has a remarkable abundance in shoot tissues while fruit sRNAs are steadily distributed from 21 to 24 nt long sequences. Even though, 24 nt long sRNAs were found the most abundant class across the libraries (Figure 3).

FIGURE 3

Figure 3. Length distribution of small RNAs in tomato (A), X-axis, length of sRNAs distribution; Y-axis, corresponding percentage of raw reads; The distribution of member number for each evolutionary known miRNA families in tomato (B), X-axis, name of the miRNAs family; Y-axis, representing miRNA numbers.

Identification of Evolutionary Known miRNAs

To analyze and predict the evolutionary known miRNAs, our sRNA libraries were BLASTN searched against the plant known mature miRNAs deposited in miRBase (v20, June 2013). A total of 68 evolutionary known miRNAs have been identified that have orthologs in other plant species and widely distributed over 37 diversified miRNA families (Figure 3). Between 37 miRNA families, SSG represented all the miRNA families followed by SSC and SFG, each of these represented 34 evolutionary known families. The number of miRNA family found less dramatically in SFC library, which represented only 25. The highest number of known miRNAs was expressed in SSG (68) followed by SSC (65) and SFG (62) and the lowest (44) was expressed in SFC. The results have identified four known miRNAs (miR171c, miR171d, miR530, and miR5302) that showed only shoot specific expression among the sRNA libraries and no known sequences identified as fruit tissues specific (Table 2 and details in Table S2).

TABLE 2

Table 2. Summarized information of miRNAs in sRNA libraries of tomato.

The member of evolutionary known miRNA families was further analyzed and distributed in Figure 3, where the majority of families (23 of 37) contained single members. Micro RNA families MIR156, MIR172, and MIR5303 contained highest five members while 11 families viz, MIR162, MIR166, MIR167, MIR168, MIR171, MIR1919, MIR319, MIR398, MIR482, MIR6024, and MIR7997 contained several members (2–4).

Abundance of and expressions of plant miRNAs always vary on spatio-temporal differences (Wang et al., 2012). Results have shown that the abundance of miRNA families like MIR159, MIR166, MIR6022, MIR162, MIR482, and some other were consistently higher throughout the libraries. But SFC has shown comparatively lower abundance for all the miRNA families. At the same time, MIR397, MIR5302, MIR8007 has infrequently sequenced. The possible reason is these miRNAs may be expressed at the lower level under certain condition and in specific tissues. This inconsistent sequencing frequency might be suggested that the conserved miRNAs are well expressed than non-conserved miRNAs in plant tissues; moreover, they may have some distinct physiological roles with diversified regulatory effects in grafted tomato. For example, miR159 is a highly conserved miRNA which controls the expression of some MYB transcription factors and regulates germination, vegetative growth, flowering time, anther development and seed shape (Palatnik et al., 2007). Conversely, some other miRNAs, such as, miR530, miR8007, miR1916, and miR8021 were found at a low-level pattern of expression. This scenario indicates that they have a narrow range of specificity to plants. These miRNAs have less conserved homologs, and their functions are also not clear yet (Liu et al., 2013).

Prediction of Potentially Novel miRNAs and Nucleotide Bias

Another most important feature of deep sequencing is the ability to detect novel miRNAs from sRNA transcriptome. Based on the available tomato genome sequence, we have detected flanking sequences of the grafted tomato candidate miRNAs. Potential novel miRNAs were identified through the exploration of secondary precursor structures, Dicer cleavage sites and by measuring the (MFE). Custom scripts were used to obtain the counts of predicted miRNAs as well as to detect the first nucleotide bias on the first position with certain length along with each position. According to these criteria, a total of 168 novel mature miRNA sequences were identified. Most of the novel miRNA candidates (166) were expressed in shoot tissues (SSC and SSG library) followed by SFG (162). Surprisingly, SFC expressed the lowest (101), which almost half of its corresponding grafted fruit tissues (Table 2). However, the novel miRNAs were found to be consistent between shoot tissues, but it was a tremendous enlargement in grafted fruit compared to its control. This big change at transcription level suggests that the regulatory reaction of miRNA sequences occurred in tomato fruits after grafting on goji and this inflation of transcriptions and miRNAs might be because of the pulling up of miRNAs from root stock to scion.

Predicted novel miRNAs found with a low abundance which supports the notion that the non-conserved miRNAs are often expressed at a lower level than the conserved ones (Figure 4). The lowly expressed miRNAs might suggest their specific roles under various habitats, developmental phases or in specific tissues, and their regulation nature remains to be investigated (Table 3).

FIGURE 4

Figure 4. Expression level (number of raw sequence reads) of evolutionary known miRNAs in each library in tomato. X-axis, name of the miRNAs family; Y-axis, representing sRNA reads. (A) SFC and SFG libraries. SFC, Solanaceae fruit control; SFG, Solanaceae fruit graft. (B) SSC and SSG libraries. SSC, Solanaceae shoot control; SSG, Solanaceae shoot graft.

TABLE 3

Table 3. Summary of differentially expressed miRNAs among the sRNA libraries in tomato.

Data of the novel candidate miRNAs were summarized with the nucleotide bias at the first position of all sRNA sequences of certain lengths. This summary could often be helpful to assess the prediction accuracy according to the base bias of known miRNAs. From the libraries, it was observed that the majority candidates started with a 5′ U (uracil) and with a length of 21 nt (Figure S2).

Differentially Expressed Known and Potentially Novel miRNAs between Grafted and Control Tomato

MiRNAs are expressed in a spatio-temporally manner or in the response of environmental stimuli. So their biogenesis is highly regulated. The miRNAs are differentially expressed between libraries and provide a clue at the molecular events related to grafting. Thus, they make a correlation between rootstock and scion. Firstly, the read density measurement (transcripts per million-TPM) was normalized followed by P < 0.01 and the absolute value of log2 ratio fold-change > 1.0 were maintained as the threshold for the statistical significance level of miRNAs.

From the Table 4, 33 known miRNAs found significantly differentially expressed between grafted and control fruit sRNA libraries and 11 in shoot libraries. A total of 123 miRNAs (26 known and 97 novels) were up-regulated in the grafted fruit and only 28 miRNAs (8 known and 20 novels) were up-regulated in grafted shoot. Moreover, 133 miRNAs (33 known and 100 novel) found with ≥ 3 log2 fold-changes between fruit samples, while only 8 miRNAs (2 known and 6 novel) in shoot samples. This picture suggests that the fruit of tomato is highly responsive by the grafting technique on goji rootstock and there might be a possible exchange of biologically important phytochemicals from root stock to scion. These grafting results imply miRNAomes, in the form of sRNAs transcripts, a very complex issue (Lu C. et al., 2005), though miRNA sequences thought to be explicitly conserved through plant species (Voinnet, 2009).

TABLE 4

Table 4. Significantly differentially expressed novel fruit miRNAs with 3 or more fold change.

The comparative analysis of miRNAs expression level, a great variation of transcriptional accumulation was observed between grafted and control tomato fruit. Among the evolutionary known miRNAs, 33 were found highly significant with more than 3 fold change with a large range of abundance from 7.09 to 191,474 TPM. Highest fold-change (10.4) observed in miR172a with up-regulation in SFG. Moreover, miR6022, miR166a-5p, and miR6024 were found as the most abundant families with 191,474.2, 4,180.5, and 2,990.12 transcripts, respectively. Conversely, miR166a-5p, miR6022, and miR7997c found as the most representing miRNA families in SFC library with 41,753.1, 5,917.10, and 3,805.84 transcripts, respectively. However, miR166a-5p found consistently higher representing miRNA in both the grafted and control fruit libraries. On the other hand, in the shoot sRNA libraries (SSG vs. SSC) only two known miRNAs (miR398a-3p and miR8021) were found in a fold change more than 3. In case of novel miRNAs, three (tono23, tono38, and tono39) were found with more than 3 fold-change. The highest fold change was identified in tono23, and it was also down-regulated in SSG and other two were up-regulated in control shoot (SSC) (Table S3).

Target Genes Prediction, GO, and KEGG Analysis of Grafted Tomato Fruits

Micro RNAs and their corresponding target genes exhibit a high sequence complementarity in plant genomes and thus explore target genes directly and successfully (Jones-Rhoades and Bartel, 2004). In the expression regulation of a target gene mechanism, miRNAs usually get hybridized to the mRNA transcripts of a target gene to promote RNA degradation, inhibit translation or both. For the better understanding of the biological functions of the differentially expressed miRNAs in grafted tomato, GO enrichment and KEGG database were analyzed. A total of 6255 genes were found to be associated in the functional annotation that was targeted by 70 significantly differentially expressed miRNAs in this study and three targets (Solyc07g062840.2, Solyc11g027650.1, and Solyc03g121000.2) were complementary to three novel miRNAs (tono50, tono95, and tono97) (Table S4). The rest of the targets were found as the complementary sequence to the known miRNAs. The potential targets of the most miRNAs were varied with a span from few to several hundreds. However, three miRNAs named tono97, miR167a and miR298b-3p have targeted the least number, one in each, while the highest number of genes (981) were targeted by miR5303. Noteworthy, many miRNAs did not match with any target at all; suggesting that this might be due to the absence of targets in the reference genome and/or low expression of these miRNAs in the sRNA libraries, nonetheless, indicates the scope of their functional characterization in tomato grafting.

GO and KEGG databases are very important for the functional annotation of gene enrichment analysis of the potential targets of significantly differentially expressed miRNAS. A total of 163 miRNAs (33 known and 130 novels) in fruit sRNA libraries (SFG vs. SFC) and 43 (11 known and 32 novels) miRNAs shoot (SSG vs. SSC) libraries were considered for the functional annotation of the targets following some previous reports (Kanehisa and Goto, 2000; Mao et al., 2005; Young et al., 2010). There are three components of GO, such as BP (biological process), CC (cellular component), and MF (molecular function). Over-represented P-values followed by corrected P-values were employed to identify the involved enriched gene categories in GO enrichment analysis. KEGG analysis was assigned for the elucidation and better understanding of different biological pathways. However, GO annotation was assigned to the all putative targets as the basis of GO components. In SFG vs. SFC libraries, only the “molecular function” was found as the enriched GO term. Among these “ADP binding regulators,” targets were identified as highly significant and over-representing, which is 2.28% of the total MF related targets and regulated by 40 genes (Figure 5 and Table S5). However, SSG vs. SSC libraries showed diversified regulatory activities in cellular component and also in molecular functions. A total of 45 genes were involved in regulation at the nucleus followed by ADP binding activities and 23 genes are involved, moreover, molybdopterin synthase complex and RNA helicase activity were as well found significantly enriched, associated with 3 and 8 genes, respectively.

FIGURE 5

Figure 5. Gene ontology terms and numbers of predicted target genes for differentially expressed miRNAs from SFG vs. SFC (A) and SSG vs. SSC (B). SFG, Solanaceae fruit graft; SFC, Solanaceae fruit control; SSG, Solanaceae shoot graft; SSC, Solanaceae shoot control. CC, Cellular Component; and MF, Molecular Function. Right hand side scale, targeted gene numbers corresponding to the GO terms; left hand side scale, percent of targeted gene numbers corresponding to the GO terms.

For the elucidation of important biological pathways' interaction and the target genes of differentially expressed miRNAs, we further have followed KEGG pathway enrichment analysis. After analyzing the target genes of miRNAs and the corrected P-values, a highly diversified 113 biochemical pathways were found to be enriched in fruit tissues. From these, 20 pathways have been identified on the basis of rich factor analysis; corrected q-value and number of involved genes in SFG vs. SFC libraries. A total of 298 target genes of the miRNAs were associated. However, “metabolic pathways” was found as the most significant pathway term in respect of rich factor and was regulated by the highest number of genes 172, followed by the ubiquitin mediated proteolysis (16 genes) and ribosome biogenesis in eukaryotes (11 genes). In SSG vs. SSC libraries, only 29 genes were found to be involved in the top 20 enriched biological pathways. Moreover, the number of involved genes in SSG vs. SSC was not more than two in any pathways. However, the enriched pathways were related to selenocompound metabolism, sulfur metabolism, and some other (Figure 6). Details of the pathway terms, id name, gene name and hyperlink of pathways are enlisted in Table S5.

FIGURE 6

Figure 6. KEGG analysis with 20 most enriched pathways from SFG vs. SFC (A), and SSG vs. SSC (B). Coloring of q-values indicates the significance of rich factor; circle indicates target genes involved and size is proportional to the gene number.

Analysis of Graft-Transmissible miRNAs As a Whole miRNA Movable Molecule

In our previous study on sRNAs characterization of goji control plants (Khaldun et al., 2015), 62 miRNAs found to be significantly differentially expressed in shoot and fruit tissues. We have further analyzed that information with present study (where goji is the rootstock) to identify common and specific significantly differentially expressed miRNAs in both shoot and fruit tissues of grafted tomato plant and their corresponding control plants. Figure 7 representing that only one significantly differentially expressed miRNA (tono154) was specific to SSG vs. SSC combination, indicating that there was no remarkable difference in miRNAs number between control and grafted tomato shoot, whereas 10 miRNAs found that were shown expression only in grafted and control tomato fruit (SFG vs. SFC) and was not found in other combinations. Interesting observation is that SSC vs. SFC showed specific difference in 28 miRNAs and LSC vs. LFC showed in 22, suggesting that numbers of miRNAs and regulating patterns widely variy between shoot and fruit tissues. However, six significantly differentially expressed miRNAs were found common in both shoot and fruit tissues of tomato and Lycium (LSC vs. LFC and SSC vs. SFC). Two miRNAs (goji novel28 and miR8007a-5p) were common between shoot and fruit in tomato and Lycium (LSC vs. LFC and SFG vs. SFC) whereas three (miR395a, miR6023, and miR8021) were common in both shoot and fruit of tomato (SFG vs. SFC and SSG vs. SSC). However, no miRNAs found common between grafted fruit and shoot of tomato and the fruit and shoot of Lycium. After analyzing the RNA reads count, three miRNAs (goji novel28, goji novel9, and miR8007a-5p) found as the transported from rootstock of goji to grafted tomato plant (Table 5). In the qRT-PCR experiment these miRNAs were expressed in the control tomato samples in a very low read counts. So, the sequencing data for these three miRNAs somehow showed expression deviation in the qRT-PCR results. But, one thing is clear these miRNAs shown strong expression in goji samples whereas they have shown a very weak level of expression. However, further experiment can identify additional information for these miRNAs as a whole miRNA movement from rootstock to the scion.

FIGURE 7

Figure 7. Common and specific response differential miRNAs in grafted tomato plants and their controls from compared libraries. SFC, Solanaceae fruit control; SFG, Solanaceae fruit graft; SSC, Solanaceae shoot control; SSG, Solanaceae shoot graft; LSC, goji shoot control; LFC, goji fruit control.

TABLE 5

Table 5. Expression pattern of the graft-transmissible miRNAs in the small RNA libraries.

Expression Validation of Sequencing Data with qRT-PCR and RLM 5′-Race

Quantitative RT-PCR was used to validate the sequencing data experimentally. For this, 19 miRNAs were selected, 14 of which were selected from the differentially expressed between grafted and control fruits and 5 between grafted and control shoot sRNA libraries after determining the appropriate reference genes (Figure 8). The C_T values and RQ values of these miRNAs were obtained, and 17 miRNAs had a similar expression pattern in the grafted and control tomato samples compared with the sequencing results. However, miR1919-5p and tono2 showed acceptable inconsistency in the expression profile between the sequencing data and RT–qPCR. The possible reasons for this inconsistency might be the difference to the amplification performance of primers or other unknown reasons. The primers that were used for this study are listed in supporting Table S6. Additionally, six important targets were validated by modified RLM 5′-RACE with their miRNA sequencing numbers and cleavage sites (Figure S3).

FIGURE 8

Figure 8. Expression ratios (grafted fruit/control fruit) of miRNAs in (A) qRT-PCR and (B) sequencing results; and (grafted shoot/control shoot) ratios in (C) qRT-PCR and (D) sequencing results. X-axis, name of the miRNAs that were selected for qRT-PCR; column above the X-axis, miRNAs that were up-regulated in the grafted fruit and shoot; column below the X-axis, miRNAs that were up-regulated in the grafted and control shoot.

Discussion

An accurate mechanism of the plant development and response to the environmental stimuli could be achieved by the gene regulation through sequence-specific interaction between miRNAs and their targets. Therefore, a holistic effort is required for the unraveling of the genetic mechanisms, especially addressing the roles of miRNAs in the grafting system. Discovery of whole sets of miRNAs and their targets may play a centric role in elucidation of complicated miRNA-mediated regulatory network for the graft system that controls the plant growth and development, and other physiological and biological processes. Undoubtedly, the advent of deep sequencing platforms have significantly widened the ability of miRNA exploration and in return provided a new horizon to identify the conserved, non-conserved, lowly expressed as well as species-specific miRNAs in a large scale. By deep sequencing, not only the conserved miRNAs but also the novel miRNAs were identified. To the best of authors' knowledge, this is the first report on high-throughput RNA sequencing of grafting between two distantly related medicinally important herb and shrub. Moreover, the differential expression was observed for the long-distance translocation of miRNAs over a long period of time (from young shoot to mature fruit). Quality control data and the size distribution have confirmed the reliability of throughput sequencing data, which was further validated by the qRT-PCR and the modified 5′RLM-RACE experiments.

miRNAome in Shoot and Fruit Tissues of Grafted and Control Tomato

High-throughput sequencing is rapidly replacing the conventional Sanger sequencing platforms though the most conserved plant miRNAs have been identified through traditional approaches. By now, the most plant miRNAs have been discovered in Arabidopsis, rice and poplar through the traditional methods (Lu S. F. et al., 2005). Some recent studies have been reported that non-conserved, and species-specific miRNAs are normally found at a lower abundance compare to conserved miRNAs, which is literally hard to reveal efficiently and cost-effectively (Moxon et al., 2008). However, deep sequencing has been successfully applied both in model (Lu et al., 2006; Rajagopalan et al., 2006; Fahlgren et al., 2007) and non-model plants (Morin et al., 2008; Qiu et al., 2009) addressed as blessings for miRNAs identification study (Song et al., 2010).

Populations and the distribution patterns are distinguishing features of sRNA libraries. Composition of sRNAs often indicates the roles of different type of sRNAs in specific tissue types, or even particular species. Most of the miRNAs with known functions are 20–24 nt long. In our four libraries, both of the shoot libraries (SSC and SSG) comprised of the 24 nt sRNAs with a remarkably higher frequency (Figure 2). On the other hand, both the fruit libraries (SFC and SFG) are comprised of 21 nt long sRNAs. This is a good indication that 24 nt long sRNAs play more roles in shoot or vegetative tissues of tomato, on the contrary; 21 nt sRNAs have important roles in fruit tissues. In some plant species like grapevine, Pinus cordata and Populus balsamifera, 21 nt has substantially more than 24 nt sRNAs (Barakat et al., 2007; Morin et al., 2008; Pantaleo et al., 2010). The opposite scenario also reported in some other studies where 24 nt class shown higher-frequency (Rajagopalan et al., 2006; Fahlgren et al., 2007; Szittya et al., 2008). The high frequency of 24 nt sequence may reflect the complexity of tissue samples since 24 nt sequences are considered as siRNAs. siRNAs are known to be involved in heterochromatin modification, specially genomes with content of repeat sequences (Herr, 2005). Moreover, a previous study was reported that 24 nt sRNAs was predominantly comprised of repeat and transposons (Lippman and Martienssen, 2004).

Until date, nearly 22 miRNA families were reported as frequently distributed at least in 20 plant species thus considered as highly conserved and found in dicots and monocots (Sunkar and Jagadeeswaran, 2008). Well-conserved miRNA families have retained homologous target interactions and performed analogous molecular functions in the plant kingdom in evolutionary timescales (Axtell et al., 2007). Gene duplication event is thought to generate plant MIR genes, which then evolved by random mutations into short, imperfectly paired hairpins (Allen et al., 2004; Axtell, 2008). Those miRNAs are mostly non-conserved and are believed to be evolutionarily young and mainly represented by a single copy MIR gene (Jones-Rhoades et al., 2006). Recent studies reported on many non-conserved miRNAs in several species. These non-conserved miRNAs are found among distinct species and in different phylogenetic families. We also have found several non-conserved miRNAs in our study, although they were found in a very low level of expression. Nevertheless, this situation does not rule out that they also could be expressed highly in other tissues. Apart from conserved miRNAs, the abundance of novel miRNAs in our study found at a lower level suggesting their involvement more specific processes to grafted tomato. This also supports the results of previous reports (Liu et al., 2013) where species-specific miRNAs were believed to be newly evolved and lowly expressed than the conserved ones. Therefore, further functional characterization of these novel miRNAs may reveal interesting and useful information about their role in signal transduction and development in grafted tomato on goji. Micro RNAs usually undergo some minimum annotating criteria regarding biogenesis and expression and convincingly detected by either of the methods like northern blotting, qRT-PCR or sequencing analysis (Song et al., 2010). Additionally sequence detection or cloning data of miRNA* are required as a biogenesis proof since they are complementary to mature miRNA sequences (Jones-Rhoades et al., 2006; Meyers et al., 2008). Among the 168 novel miRNAs, 98 were found with their complementary miRNA* sequences (Table 2).

Target Annotation of Differentially Expressed miRNAs

Some miRNAs such as, miR156, miR162, miR164, miR166, miR172, miR397, and miR398 were reported to be highly conserved in tomato fruit and developmental stages (Zuo et al., 2012; Karlova et al., 2013). In our present study, homologs of miR162, miR164, miR166, and miR397 were identified as significantly differentially expressed between grafted and control tomato fruits but not in the shoot tissues (Table S7). All these miRNAs found to be up-regulated in grafted fruit. Surprisingly, these miRNAs were detected as significantly differentially expressed between goji fruit and shoot, where they were up-regulated in fruit (Khaldun et al., 2015). This accumulation difference of these miRNAs in grafted tomato fruit may be due to the long-distance movement after grafting. MiR398 exhibited an interesting expression pattern where it was found as significantly expressed in grafted tomato fruit though it was down regulated in control tomato samples (Table S8). Moreover, it has found significantly differentially expressed in goji fruit, suggested this induction probably due to grafting method. This phenomenon could be a proof of long distance movement of miRNAs from rootstock to scion, though their expression was remarkably low in grafted tomato fruit. These results convincingly implicated that grafting could indeed affect miRNAs expression in scion.

Most of the cases in plants, targets of miRNA family mostly belong to the same gene family (Karlova et al., 2013). Additionally, the perfect sequence complementarity between miRNA and targets conferred the target prediction process straightforward and efficient (Jones-Rhoades and Bartel, 2004). Number of target genes of a single miRNA has counted from one to more than several hundred in our study. Mostly, a single miRNA can work on many targets; similarly, a single gene could be regulated by many miRNAs, making the functional annotation extremely complex. Thus, from high-throughput sequencing method enriched GO annotation and KEGG pathway analyses are efficiently used for the functional annotation of miRNAs with their targets since it is not possible “single gene functional characterization” from such huge data. “Metabolic pathways” was identified as the most enriched pathway between grafted and control fruit. Most of the known miRNAs were also involved in this regulatory mechanism. “Amino sugar and nucleotide sugar metabolism” has identified as the most enriched pathway between grafted shoot and control shoot, where only two genes and two miRNAs involved. So, these genes and miRNAs could be a good source of biological pathway regulation process and should be paid special attention in the future study. The targeted genes of graft-mediated miRNAs encoded proteins of diverse functions such as transcription factors, mostly, and others were involved in metabolism, signal transduction, growth, development and other biological processes (Table S5). Hormones, nucleic acids, metabolics, and proteins can act as systemic signals that could provide useful exchange information between rootstock and scion. It is expected for the conserved miRNAs that their targets should also be well-conserved. Most of the conserved miRNAs (such as miR156, miR159, miR160, miR164, miR167, miR171, miR172, miR319, and some others) usually target a range of transcription factors like MYBs, ARFs, SBPs, NACs, AP2-like factors, GRFs, and GRASs, and their miRNAs-mediated regulations are important for plant growth and development and may act in the core gene expression networks (Liu et al., 2013). Molecular information can move bidirectional in grafting system; however, signaling between rootstock-scion and scion-rootstock requires further studies. In comparison with grafted tomato, and control samples of tomato and goji we have found many miRNAs were differentially expressed.

Conclusions

Present study reported the miRNAs profiling and target annotation in grafted tissues between two important plants L. chinense and tomato. This report, to our knowledge, is the first elucidation of miRNAs in grafting between distantly related two plants. High-throughput sequencing has employed, and 68 evolutionary known miRNAs identified belonging to 37 diversified families. Additionally, 168 putative novel miRNAs were identified from shoot and fruit tissues of grafted and control tomato. A total of 163 (33 known and 130 novels) miRNAs were explored as significantly differentially expressed between fruit sRNA libraries. Moreover, GO analysis has detected “ADP binding activities,” “molybdopterin synthase complex,” and “RNA helicase activity” as the most enriched terms, and KEGG analysis has identified “metabolic pathways” as the most enriched pathway involving 172 genes in grafted and control fruit tissues. Therefore, this study provided a unique insight of graft-oriented miRNAs and their complex regulatory functions in Lycium and tomato, which would be helpful for studying the various biological pathways as well as the miRNA translocation in a distant-grafting system.

Author Contributions

AK and YW conceived and designed the experiments. AK, WH, HL, and SL performed the experiments. AK, WH, and SZ analyzed the data. YW, HL, and SL contributed chemicals/reagents/materials/analysis tools. AK, YW, and WH wrote the paper.

Conflict of Interest Statement

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.

Acknowledgments

This work was supported by grants from the National Natural Science Foundation of China (31200242 and 31470391) and the Scientific Project of Ningxia Agriculture Comprehensive Development (znnfkj2015). Youth Innovation Promotion Association CAS (2015286). We would like to thank Mr. Chengyu Wang for teaching the grafting techniques. We also thank Dr. Adeniyi Charles Adeola extending his helping hand for revising the language of the manuscript.

Supplementary Material

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fpls.2016.01475

Figure S1. Flowchart of the small RNAs indentification bioinformatics.

Figure S2. First nucleotide bias of novel miRNA candidates in tomato. (A) miRNA first nucleotide bias in SFC, (B) miRNA first nucleotide bias in SFG, (C) miRNA first nucleotide bias in SSC, (D) miRNA first nucleotide bias in SSG. Number on top of the bars indicates the number of sequences corresponding to the miRNA length; X-axis, length of miRNAs; Y-axis, representing percent of each nucleotide.

Figure S3. Target validation of known tomato miRNAs. Arrows showing the 5′-ends of the cleavage products.

Table S1. The designed primers for modified RLM 5′-RACE and targets annotation in tomato.

Table S2. Statistics of evolutionary known and novel miRNAs in grafted and control tomato plants.

Table S3. Details of differentially expressed evolutionary known and novel miRNAs in grafted and control tomato plants.

Table S4. Differentially expressed miRNAs and their corresponding genes in grafted tomato.

Table S5. GO annotation, KEGG enriched pathways and corresponding genes in tomato.

Table S6. The primer designed and candidate list of miRNAs for qRT-PCR with differential expression through deep sequencing.

Table S7. Significantly up-regulated evolutionary known miRNAs with 2.4-, or more fold-change in grafted tomato.

Table S8. Significantly down-regulated evolutionary known miRNAs in grafted tomato.

References

Agarwal, S., and Rao, A. V. (2000). Tomato lycopene and its role in human health and chronic diseases. Can. Med. Assoc. J. 163, 739–744. Available online at: http://www.cmaj.ca/content/163/6/739.full

PubMed Abstract

Allen, E., Xie, Z., Gustafson, A. M., Sung, G.-H., Spatafora, J. W., and Carrington, J. C. (2004). Evolution of microRNA genes by inverted duplication of target gene sequences in Arabidopsis thaliana. Nat. Genet. 36, 1282–1290. doi: 10.1038/ng1478

PubMed Abstract | CrossRef Full Text | Google Scholar

Atkins, C. A., Smith, P. M., and Rodriguez-Medina, C. (2011). Macromolecules in phloem exudates—a review. Protoplasma 248, 165–172. doi: 10.1007/s00709-010-0236-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Axtell, M. J. (2008). Evolution of microRNAs and their targets: are all microRNAs biologically relevant? Biochim. Biophys. Acta 1779, 725–734. doi: 10.1016/j.bbagrm.2008.02.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Axtell, M. J., Snyder, J. A., and Bartel, D. P. (2007). Common functions for diverse small RNAs of land plants. Plant Cell 19, 1750–1769. doi: 10.1105/tpc.107.051706

PubMed Abstract | CrossRef Full Text

Barakat, A., Wall, P., and Diloreto, S. (2007). Conservation and divergence of microRNAs in Populus. BMC Genomics 8:481. doi: 10.1186/1471-2164-8-481

PubMed Abstract | CrossRef Full Text

Bartel, D. (2004). MicroRNAs: genomics, biogenesis, mechanism, and function. Cell 116, 281–297. doi: 10.1016/S0092-8674(04)00045-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Benson, D. A., Karsch-Mizrachi, I., Lipman, D. J., Ostell, J., and Sayers, E. W. (2011). GenBank. Nucleic Acids Res. 39, D32–D37. doi: 10.1093/nar/gkq1079

PubMed Abstract | CrossRef Full Text | Google Scholar

Bhogale, S., Mahajan, A. S., Natarajan, B., Rajabhoj, M., Thulasiram, H. V., and Banerjee, A. K. (2014). MicroRNA156: a potential graft-transmissible microRNA that modulates plant architecture and tuberization in Solanum tuberosum ssp. andigena. Plant Physiol. 164, 1011–1027. doi: 10.1104/pp.113.230714

PubMed Abstract | CrossRef Full Text | Google Scholar

Buhtz, A., Pieritz, J., Springer, F., and Kehr, J. (2010). Phloem small RNAs, nutrient stress responses, and systemic mobility. BMC Plant Biol. 10:64. doi: 10.1186/1471-2229-10-64

PubMed Abstract | CrossRef Full Text | Google Scholar

Carlsbecker, A., Lee, J.-Y., Roberts, C. J., Dettmer, J., Lehesranta, S., Zhou, J., et al. (2010). Cell signalling by microRNA165/6 directs gene dose-dependent root cell fate. Nature 465, 316–321. doi: 10.1038/nature08977

PubMed Abstract | CrossRef Full Text | Google Scholar

Carrington, J. C., and Ambros, V. (2003). Role of microRNAs in plant and animal development. Science 301, 336–338. doi: 10.1126/science.1085242

PubMed Abstract | CrossRef Full Text | Google Scholar

Chuck, G., and O'Connor, D. (2010). Small RNAs going the distance during plant development. Curr. Opin. Plant Biol. 13, 40–45. doi: 10.1016/j.pbi.2009.08.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Corbesier, L., Vincent, C., Jang, S., Fornara, F., Fan, Q., Searle, I., et al. (2007). FT protein movement contributes to long-distance signaling in floral induction of Arabidopsis. Science 316, 1030–1033. doi: 10.1126/science.1141752

PubMed Abstract | CrossRef Full Text | Google Scholar

Dong, J.-Z., Yang, J.-J., and Wang, Y. (2008). Resources of Lycium species and related research progress. Zhongguo Zhong Yao Za Zhi 33, 2020–2027. doi: 10.3321/j.issn:1001-5302.2008.18.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Fahlgren, N., Howell, M. D., Kasschau, K. D., Chapman, E. J., Sullivan, C. M., Cumbie, J. S., et al. (2007). High-Throughput sequencing of Arabidopsis microRNAs: evidence for frequent birth and death of MIRNA genes. PLoS ONE 2:e219. doi: 10.1371/journal.pone.0000219

PubMed Abstract | CrossRef Full Text | Google Scholar

Fawusi, M. O. A. (1978). Emergence and seedling growth of pepper as influenced by soil compaction, nutrient status and moisture regime. Sci. Hortic. 9, 329–335. doi: 10.1016/0304-4238(78)90042-0

CrossRef Full Text | Google Scholar

Flint, H. L. (1997). “Lycium barbarum” Landscape Plants for Eastern North America: Exclusive of Florida and the Immediate Gulf Coast. Chichester: John Wiley & Sons.

Friedlander, M. R., Mackowiak, S. D., Li, N., Chen, W., and Rajewsky, N. (2012). miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Res. 40, 37–52. doi: 10.1093/nar/gkr688

PubMed Abstract | CrossRef Full Text | Google Scholar

Griffiths-Jones, S., Saini, H. K., Van Dongen, S., and Enright, A. J. (2008). miRBase: tools for microRNA genomics. Nucleic Acids Res. 36, D154–D158. doi: 10.1093/nar/gkm952

PubMed Abstract | CrossRef Full Text | Google Scholar

Hafner, M., Landgraf, P., Ludwig, J., Rice, A., Ojo, T., Lin, C., et al. (2008). Identification of microRNAs and other small regulatory RNAs using cDNA library sequencing. Methods 44, 3–12. doi: 10.1016/j.ymeth.2007.09.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Hamilton, A., Voinnet, O., Chappell, L., and Baulcombe, D. (2002). Two classes of short interfering RNA in RNA silencing. Embo J. 21, 4671–4679. doi: 10.1093/emboj/cdf464

PubMed Abstract | CrossRef Full Text | Google Scholar

Herr, A. J. (2005). Pathways through the small RNA world of plants. FEBS Lett. 579, 5879–5888. doi: 10.1016/j.febslet.2005.08.040

PubMed Abstract | CrossRef Full Text | Google Scholar

Jones-Rhoades, M., and Bartel, D. (2004). Computational identification of plant microRNAs and their targets, including a stress-induced miRNA. Mol. Cell 14, 787–799. doi: 10.1016/j.molcel.2004.05.027

PubMed Abstract | CrossRef Full Text | Google Scholar

Jones-Rhoades, M., Bartel, D., and Bartel, B. (2006). MicroRNAs and their regulatory roles in plants. Annu. Rev. Plant Biol. 57, 19–53. doi: 10.1146/annurev.arplant.57.032905.105218

PubMed Abstract | CrossRef Full Text | Google Scholar

Joshi, T., Yan, Z., Libault, M., Jeong, D.-H., Park, S., Green, P. J., et al. (2010). Prediction of novel miRNAs and associated target genes in Glycine max. BMC Bioinformat. 11:S14. doi: 10.1186/1471-2105-11-s1-s14

PubMed Abstract | CrossRef Full Text | Google Scholar

Kanehisa, M., Araki, M., Goto, S., Hattori, M., Hirakawa, M., Itoh, M., et al. (2008). KEGG for linking genomes to life and the environment. Nucleic Acids Res. 36, D480–D484. doi: 10.1093/nar/gkm882

PubMed Abstract | CrossRef Full Text | Google Scholar

Kanehisa, M., and Goto, S. (2000). KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 28, 27–30. doi: 10.1093/nar/28.1.27

PubMed Abstract | CrossRef Full Text | Google Scholar

Karlova, R., Van Haarst, J. C., Maliepaard, C., Van De Geest, H., Bovy, A. G., Lammers, M., et al. (2013). Identification of microRNA targets in tomato fruit development using high-throughput sequencing and degradome analysis. J. Exp. Bot. 64, 1863–1878. doi: 10.1093/jxb/ert049

PubMed Abstract | CrossRef Full Text | Google Scholar

Kehr, J., and Buhtz, A. (2008). Long distance transport and movement of RNA through the phloem. J. Exp. Bot. 59, 85–92. doi: 10.1093/jxb/erm176

PubMed Abstract | CrossRef Full Text | Google Scholar

Kehr, J., and Buhtz, A. (2013). “Endogenous RNA constituents of the phloem and their possible roles in long-distance signalling,” in Phloem. Molecular Cell Biology, Systemic Communication, Biotic Interactions, eds G. A Thompson and A. J. E. Van Bel (Oxford: Wiley-Blackwell Publishing), 186–208.

Khaldun, A., Huang, W., Liao, S., Lv, H., and Wang, Y. (2015). Identification of microRNAs and target genes in the fruit and shoot tip of Lycium chinense: a traditional chinese medicinal plant. PLoS ONE 10:e0116334. doi: 10.1371/journal.pone.0116334

PubMed Abstract | CrossRef Full Text | Google Scholar

Langmead, B., Trapnell, C., Pop, M., and Salzberg, S. L. (2009). Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 10:R25. doi: 10.1186/gb-2009-10-3-r25. Available online at: http://link.springer.com/article/10.1186/gb-2009-10-3-r25?view=classic

PubMed Abstract | CrossRef Full Text

Lippman, Z., and Martienssen, R. (2004). The role of RNA interference in heterochromatic silencing. Nature 431, 364–370. doi: 10.1038/nature02875

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, N., Yang, J., Guo, S., Xu, Y., and Zhang, M. (2013). Genome-wide identification and comparative analysis of conserved and novel microRNAs in grafted watermelon by high-throughput sequencing. PLoS ONE 8:e57359. doi: 10.1371/journal.pone.0057359

PubMed Abstract | CrossRef Full Text | Google Scholar

Livak, K. J., and Schmittgen, T. D. (2001). Analysis of relative gene expression data using real-time quantitative PCR and the 2(−ΔΔCT) method. Methods 25, 402–408. doi: 10.1006/meth.2001.1262

PubMed Abstract | CrossRef Full Text | Google Scholar

Llave, C., Xie, Z., Kasschau, K., and Carrington, J. (2002). Cleavage of scarecrow-like mRNA targets directed by a class of Arabidopsis miRNA. Science 297, 2053. doi: 10.1126/science.1076311

PubMed Abstract | CrossRef Full Text | Google Scholar

Lu, C., Kulkarni, K., Souret, F. F., Muthuvalliappan, R., Tej, S. S., Poethig, R. S., et al. (2006). MicroRNAs and other small RNAs enriched in the Arabidopsis RNA-dependent RNA polymerase-2 mutant. Genome Res. 16, 1276–1288. doi: 10.1101/gr.5530106

PubMed Abstract | CrossRef Full Text | Google Scholar

Lu, C., Tej, S. S., Luo, S. J., Haudenschild, C. D., Meyers, B. C., and Green, P. J. (2005). Elucidation of the small RNA component of the transcriptome. Science 309, 1567–1569. doi: 10.1126/science.1114112

PubMed Abstract | CrossRef Full Text | Google Scholar

Lu, S. F., Sun, Y. H., Shi, R., Clark, C., Li, L. G., and Chiang, V. L. (2005). Novel and mechanical stress-responsive microRNAs in Populus trichocarpa that are absent from Arabidopsis. Plant Cell 17, 2186–2203. doi: 10.1105/tpc.105.033456

PubMed Abstract | CrossRef Full Text | Google Scholar

Mallory, A., and Vaucheret, H. (2006). Functions of microRNAs and related small RNAs in plants. Nat. Genet. 38, S31–S36. doi: 10.1038/ng1791

PubMed Abstract | CrossRef Full Text | Google Scholar

Mao, X. Z., Cai, T., Olyarchuk, J. G., and Wei, L. P. (2005). Automated genome annotation and pathway identification using the KEGG Orthology (KO) as a controlled vocabulary. Bioinformatics 21, 3787–3793. doi: 10.1093/bioinformatics/bti430

PubMed Abstract | CrossRef Full Text | Google Scholar

Martin, A., Adam, H., Díaz-Mendoza, M., Żurczak, M., González-Schain, N. D., and Suárez-López, P. (2009). Graft-transmissible induction of potato tuberization by the microRNA miR172. Development 136, 2873–2881. doi: 10.1242/dev.031658

PubMed Abstract | CrossRef Full Text | Google Scholar

Meyers, B., Axtell, M., Bartel, B., Bartel, D., Baulcombe, D., Bowman, J., et al. (2008). Criteria for annotation of plant microRNAs. Plant Cell 20, 3186–3190. doi: 10.1105/tpc.108.064311

PubMed Abstract | CrossRef Full Text | Google Scholar

Miyashima, S., Koi, S., Hashimoto, T., and Nakajima, K. (2011). Non-cell-autonomous microRNA165 acts in a dose-dependent manner to regulate multiple differentiation status in the Arabidopsis root. Development 138, 2303–2313. doi: 10.1242/dev.060491

PubMed Abstract | CrossRef Full Text | Google Scholar

Morin, R. D., Aksay, G., Dolgosheina, E., Ebhardt, H. A., Magrini, V., Mardis, E. R., et al. (2008). Comparative analysis of the small RNA transcriptomes of Pinus contorta and Oryza sativa. Genome Res. 18, 571–584. doi: 10.1101/gr.6897308

PubMed Abstract | CrossRef Full Text | Google Scholar

Moxon, S., Jing, R., Szittya, G., Schwach, F., Pilcher, R. L. R., Moulton, V., et al. (2008). Deep sequencing of tomato short RNAs identifies microRNAs targeting genes involved in fruit ripening. Genome Res. 18, 1602–1609. doi: 10.1101/gr.080127.108

PubMed Abstract | CrossRef Full Text | Google Scholar

Naqvi, A. R., Sarwat, M., Hasan, S., and Roychodhury, N. (2012). Biogenesis, functions and fate of plant microRNAs. J. Cell. Physiol. 227, 3163–3168. doi: 10.1002/jcp.24052

PubMed Abstract | CrossRef Full Text | Google Scholar

Palatnik, J., Wollmann, H., Schommer, C., Schwab, R., Boisbouvier, J., Rodriguez, R., et al. (2007). Sequence and expression differences underlie functional specialization of Arabidopsis microRNAs miR159 and miR319. Dev. Cell 13, 115–125. doi: 10.1016/j.devcel.2007.04.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Palauqui, J. C., Elmayan, T., Pollien, J. M., and Vaucheret, H. (1997). Systemic acquired silencing: transgene-specific post-transcriptional silencing is transmitted by grafting from silenced stocks to non-silenced scions. Embo J. 16, 4738–4745. doi: 10.1093/emboj/16.15.4738

PubMed Abstract | CrossRef Full Text | Google Scholar

Pan, X., Sun, X., Zhang, F., Zhao, C., Zhang, X., and Du, W. (2012). Optimization of distant grafting mutagenesis technology in soybean. Soybean Science 31, 237–241. doi: 10.3969/j.issn.1000-9841.2012.02.015

CrossRef Full Text | Google Scholar

Pant, B. D., Buhtz, A., Kehr, J., and Scheible, W. R. (2008). MicroRNA399 is a long-distance signal for the regulation of plant phosphate homeostasis. Plant J. 53, 731–738. doi: 10.1111/j.1365-313X.2007.03363.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Pantaleo, V., Szittya, G., Moxon, S., Miozzi, L., Moulton, V., Dalmay, T., et al. (2010). Identification of grapevine microRNAs and their targets using high-throughput sequencing and degradome analysis. Plant J. 62, 960–976. doi: 10.1111/j.1365-313X.2010.04208.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Qiu, D., Pan, X., Wilson, I. W., Li, F., Liu, M., Teng, W., et al. (2009). High throughput sequencing technology reveals that the taxoid elicitor methyl jasmonate regulates microRNA expression in Chinese yew (Taxus chinensis). Gene 436, 37–44. doi: 10.1016/j.gene.2009.01.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Rajagopalan, R., Vaucheret, H., Trejo, J., and Bartel, D. (2006). A diverse and evolutionarily fluid set of microRNAs in Arabidopsis thaliana. Genes Dev. 20, 3407–3425. doi: 10.1101/gad.1476406

PubMed Abstract | CrossRef Full Text

Shi, R., and Chiang, V. L. (2005). Facile means for quantifying microRNA expression by real-time PCR. Biotechniques 39, 519–525. doi: 10.2144/000112010

PubMed Abstract | CrossRef Full Text | Google Scholar

Song, C., Wang, C., Zhang, C., Korir, N., Yu, H., Ma, Z., et al. (2010). Deep sequencing discovery of novel and conserved microRNAs in trifoliate orange (Citrus trifoliata). BMC Genomics 11:431. doi: 10.1186/1471-2164-11-431

PubMed Abstract | CrossRef Full Text | Google Scholar

Sonoda, S., and Nishiguchi, M. (2000). Graft transmission of post-transcriptional gene silencing: target specificity for RNA degradation is transmissible between silenced and non-silenced plants, but not between silenced plants. Plant J. 21, 1–8. doi: 10.1046/j.1365-313x.2000.00645.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Stegemann, S., and Bock, R. (2009). Exchange of genetic material between cells in plant tissue grafts. Science 324, 649–651. doi: 10.1126/science.1170397

PubMed Abstract | CrossRef Full Text | Google Scholar

Storey, J. D. (2003). The positive false discovery rate: a Bayesian interpretation and the q-value. Ann. Stat. 31, 2013–2035. doi: 10.1214/aos/1074290335

CrossRef Full Text | Google Scholar

Subramanian, S., Fu, Y., Sunkar, R., Barbazuk, W. B., Zhu, J.-K., and Yu, O. (2008). Novel and nodulation-regulated microRNAs in soybean roots. BMC Genomics 9:160. doi: 10.1186/1471-2164-9-160

PubMed Abstract | CrossRef Full Text

Sunkar, R., Chinnusamy, V., Zhu, J., and Zhu, J.-K. (2007). Small RNAs as big players in plant abiotic stress responses and nutrient deprivation. Trends Plant Sci. 12, 301–309. doi: 10.1016/j.tplants.2007.05.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Sunkar, R., and Jagadeeswaran, G. (2008). In silico identification of conserved microRNAs in large number of diverse plant species. BMC Plant Biol. 8:37. doi: 10.1186/1471-2229-8-37

PubMed Abstract | CrossRef Full Text | Google Scholar

Szittya, G., Moxon, S., Santos, D., Jing, R., Fevereiro, M., Moulton, V., et al. (2008). High-throughput sequencing of Medicago truncatula short RNAs identifies eight new miRNA families. BMC Genomics 9:593. doi: 10.1186/1471-2164-9-593

PubMed Abstract | CrossRef Full Text | Google Scholar

Taller, J., Hirata, Y., Yagishita, N., Kita, M., and Ogata, S. (1998). Graft-induced genetic changes and the inheritance of several characteristics in pepper (Capsicum annuum L.). Theor. Appl. Genet. 97, 705–713. doi: 10.1007/s001220050946

CrossRef Full Text | Google Scholar

Tiedemann, R., and Carstensbehrens, U. (1994). Influence of grafting on the phloem protein patterns in Cucurbitaceae. I. additional phloem exudate proteins in Cucumis sativus grafted on two Cucurbita species. J. Plant Physiol. 143, 189–194. doi: 10.1016/S0176-1617(11)81685-8

CrossRef Full Text | Google Scholar

Voinnet, O. (2009). Origin, biogenesis, and activity of plant microRNAs. Cell 136, 669–687. doi: 10.1016/j.cell.2009.01.046

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, K., Li, M., Gao, F., Li, S., Zhu, Y., and Yang, P. (2012). Identification of conserved and novel microRNAs from Liriodendron chinense floral tissues. PLoS ONE 7:e44696. doi: 10.1371/journal.pone.0044696

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, L., Liu, H., Li, D., and Chen, H. (2011). Identification and characterization of maize microRNAs involved in the very early stage of seed germination. BMC Genomics 12:154. doi: 10.1186/1471-2164-12-154

PubMed Abstract | CrossRef Full Text | Google Scholar

Waterhouse, P. M., Wang, M.-B., and Lough, T. (2001). Gene silencing as an adaptive defence against viruses. Nature 411, 834–842. doi: 10.1038/35081168

PubMed Abstract | CrossRef Full Text | Google Scholar

Wei-Min, X. (2005). Isolation and Analysis of the Ty1/copia Retrotransposons of Mungbean and Its Possible Role on the Variation of Graft [XL]. Doctoral, Beijing University.

Wen, M., Shen, Y., Shi, S. H., and Tang, T. (2012). miREvo: an integrative microRNA evolutionary analysis platform for next-generation sequencing experiments. BMC Bioinformat. 13:140. doi: 10.1186/1471-2105-13-140

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, H.-J., Ma, Y.-K., Chen, T., Wang, M., and Wang, X.-J. (2012). PsRobot: a web-based plant small RNA meta-analysis toolbox. Nucleic Acids Res. 40, W22–W28. doi: 10.1093/nar/gks554

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, Q., Liu, Y., Zhu, A., Wu, X., Ye, J., Yu, K., et al. (2010). Discovery and comparative profiling of microRNAs in a sweet orange red-flesh mutant and its wild type. BMC Genomics 11:246. doi: 10.1186/1471-2164-11-246

PubMed Abstract | CrossRef Full Text | Google Scholar

Young, M. D., Wakefield, M. J., Smyth, G. K., and Oshlack, A. (2010). Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biol. 11:R14. doi: 10.1186/gb-2010-11-2-r14

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, S., Tuan, P., Li, X., Kim, Y., Kim, H., Park, C., et al. (2013). Identification of phenylpropanoid biosynthetic genes and phenylpropanoid accumulation by transcriptome analysis of Lycium chinense. BMC Genomics 14:802. doi: 10.1186/1471-2164-14-802

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, L., Chen, J., Li, Z., Li, X., Hu, X., Huang, Y., et al. (2010). Integrated profiling of microRNAs and mRNAs: microRNAs located on Xq27.3 associate with clear cell renal cell carcinoma. PLoS ONE 5:e15224. doi: 10.1371/journal.pone.0015224

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, J.-K. (2008). Reconstituting plant miRNA biogenesis. Proc. Natl. Acad. Sci. U.S.A. 105, 9851–9852. doi: 10.1073/pnas.0805207105

PubMed Abstract | CrossRef Full Text | Google Scholar

Zuo, J., Zhu, B., Fu, D., Zhu, Y., Ma, Y., Chi, L., et al. (2012). Sculpting the maturation, softening and ethylene pathway: the influences of microRNAs on tomato fruits. BMC Genomics 13:7. doi: 10.1186/1471-2164-13-7

PubMed Abstract | CrossRef Full Text | Google Scholar