Identiﬁcation of Grafting-Responsive MicroRNAs Associated with Growth Regulation in Pecan [ Carya illinoinensis (Wangenh.) K. Koch]

: Pecan [ Carya illinoinensis (Wangenh.) K. Koch] is an economically important nut tree and grafting is often used for clonal propagation of cultivars. However, there is a lack of research on the e ﬀ ects of rootstocks on scions, which are meaningful targets for directed breeding of pecan grafts. MicroRNAs (miRNAs) play an important role in many biological processes, but the mechanism underlying the involvement of miRNAs in grafting-conferred physiological changes is unclear. To identify the grafting-responsive miRNAs that may be involved in the regulation of growth in grafted pecan, six small RNA libraries were constructed from the phloem of two groups of grafts with signiﬁcantly di ﬀ erent growth performance on short and tall rootstocks. A total of 441 conserved miRNAs belonging to 42 miRNA families and 603 novel miRNAs were identiﬁed. Among the identiﬁed miRNAs, 24 (seven conserved and 17 novel) were signiﬁcantly di ﬀ erentially expressed by the di ﬀ erent grafts, implying that they might be responsive to grafting and potentially involved in the regulation of graft growth. Ninety-ﬁve target genes were predicted for the di ﬀ erentially expressed miRNAs; gene annotation was available for 33 of these. Analysis of their targets suggested that the miRNAs may regulate auxin transport, cell activity, and inorganic phosphate (Pi) acquisition, and thereby, mediate pecan graft growth. Use of the recently-published pecan genome enabled identiﬁcation of a substantial population of miRNAs, which are now available for further research. We also identiﬁed the grafting-responsive miRNAs and their potential roles in pecan graft growth, providing a basis for research on long-distance regulation in grafted pecan.


Introduction
As a traditional clonal propagation technique, grafting is widely used in horticultural crops. It combines materials from two different plants: the bottom part of one plant, or rootstock, which contributes roots and support; and the upper part, or scion, from another plant, contributing stems, leaves, flowers, and fruits [1]. Many phenotypic features of selected rootstocks profoundly influence pecan [29]; this is the only currently existing report concerning small miRNAs in pecan. Notably, the pecan genome was revealed recently [30]. The use of this genome for miRNA analysis enables us to achieve results with improved accuracy and reliability. In the present work, grafts with different growth vigor on two types of rootstocks (short and tall) were subjected to Illumina-based deep sequencing and analysis of miRNA. The conserved and novel miRNAs were identified by blast against the miRBase and pecan genomic sequences. Significantly differentially-expressed miRNAs in the two group of grafted pecans were analyzed. Their target genes were predicted, and these genes were subjected to Gene Ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis. Lastly, miRNAs which were significantly differentially expressed by our two growth variants were detected using the Quantitative Reverse Transcription-Polymerase Chain Reaction (qRT-PCR) technique. The main objective of our study was to enrich our knowledge of the miRNA population in pecan and explore grafting-responsive miRNAs associated with different growth vigor. This study will provide a basis for further investigation of the mechanisms involved in grafting-induced changes and the functions of miRNAs in the regulation of plant growth and development.

Plant Materials
Experimental plants were selected from the approximately 200,000 1-year-old seedlings available in the nursery of Nanjing Green Universe Pecan Science & Technology Co. in December 2015. We obtained 180 short seedlings (average height ± standard deviation, 16.7 ± 2.8 cm) and 180 tall seedlings (68.3 ± 8.5 cm). They were transplanted as the split plot design to the experimental farm of Nanjing Forestry University (both of the row and plant spacing are 50 cm). Three blocks were set up, each split into four plots of seedlings as follows: 30 short seedlings for later grafting with the cultivar 'Pawnee'; 30 short seedlings for grafting with 'Shaoxing'; 30 tall seedlings for grafting with 'Pawnee'; and 30 tall seedlings for grafting with 'Shaoxing'. The areas containing all 60 seedlings of each size within a plot were designated as the main plots for each block; the groups of 30 seedlings were designated as secondary plots. 'Pawnee' originated in 1963 from a controlled cross of 'Mohawk' and 'Starking Hardy Giant' [31], and 'Shaoxing' originated from a seedling planted in Zhejiang Province, China [32].
In September 2016, seedling heights were measured and a second selection was made based on this data. In each block, those seedlings shorter than the average for short seedlings within that block were selected as short rootstocks (SR) for grafting; conversely, those seedlings taller than the average for tall seedlings within that block were selected as tall rootstocks (TR). The two types of selected rootstocks had significantly different heights (30.0 ± 4.5 cm for SR, 94.6 ± 8.2 cm for TR; p < 0.05; Figure S1). Patch budding was used to graft scion wood of the appropriate cultivar onto rootstocks, as specified above, at approximately 10 cm above ground level. All patches for each scion cultivar were collected from a single tree.
At the end of the 2017 growing season, the three shortest grafts in each SR secondary plot were marked, as were the three tallest grafts in each TR secondary plot. In September 2018, phloem samples 10-20 cm above the graft union were collected from the marked trees. The samples were immediately frozen in liquid nitrogen and then stored at −80 • C for use. At the end of the 2018 growing season, the shortest graft was selected among the three marked trees in each SR secondary plot, as were the tallest graft among the marked trees in each TR secondary plot. The graft height (from the graft union to the top) and diameter (1 cm above the graft union) of the selected grafts, both of which were 'Pawnee' grafts, are shown in Figure S2. The growth indexes of the SG grafts were significantly lower than those of the LG grafts in in both 2017 and 2018 (p < 0.05). Phloem from the two groups of 'Pawnee' grafts with significant differences in growth (SG: grafts with small growth increment on short rootstocks; LG: grafts with large growth increment on tall rootstocks) were used as the materials for miRNA sequencing and analysis, and the three grafts in each selected group were treated as three biological replicates.

Small RNA Library Construction and Deep Sequencing
Total RNA was extracted from the phloem samples using Trizol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's instructions. Small RNA libraries were constructed by TruSeq Small RNA Sample Prep Kits (Illumina, San Diego, CA, USA). In accordance with the manufacturer's protocol, 3 adapters and 5 adapters were first ligated to the RNA with T4 RNA ligase. The ligation products were then reverse-transcribed and PCR amplification was performed on the cDNAs produced by reverse transcription. Subsequently, the PCR fragments (138-146 bp) were recovered using PAGE (polyacrylamide gel electrophoresis) and were sequenced by an Illumina Hiseq2000/2500 (Illumina, San Diego, CA, USA).

Analysis of Differential Expressed MiRNAs
To analyze the differentially expressed miRNAs between the two groups of pecan grafts with different increment, the deep-sequencing counts were normalized as the procedures described in previous study [35]. The normalized value was subjected to Student's t-tests to analyze the significance of difference in miRNA expression according to the criteria of p-value < 0.05 [36]. To better display the data, Z-value was used to plot. It was calculated as the following formula [37]: Z-value = [log 2 (norm of sample i) − mean (log 2 (norm) of all samples)]/standard deviation (log 2 (norm) of all samples), norm representing the normalized value of miRNA counts.

Target Gene Prediction and Enrichment Analysis (GO and KEGG) of Targets
To predict the genes targeted by the significantly differently expressed miRNAs, Target Finder (https://github.com/carringtonlab/TargetFinder) was used to identify miRNA binding sites against the genomic sequences of pecan. GO enrichment analysis was conducted to explore the functions of the target genes (http://www.geneontology.org/). KEGG analysis was performed for understanding of regulation pathway that the targets are potentially involved in (https://www.genome.jp/kegg/). p-value < 0.05 was used as a threshold to indicate the significantly enriched KEGG pathways and this p-value according to Fisher's exact test was calculated as the formula reported by Zhang et al. [38].

Validation of MiRNAs Using qRT-PCR
qRT-PCR was carried out to detect the miRNA expression and validate the results from deep sequencing. Total RNA was extracted from the six grafts as the methods described above and the first-strand cDNA was synthesized by TUREscript 1st Stand cDNA SYNTHESIS Kit (Aidlab, Beijing, China). U6 was used as the reference gene. The PCR amplification reactions were conducted in a volume of 10 µL, including 5 µL 2× SYBR ® Green Supermix (Bio-Rad, Hercules, CA, USA), 0.5 µL of each forward primer (4 µM), 1 µL cDNA and 3 µL ddH 2 O. They were performed on analytikjena-qTOWER2.2 Fluorescent (Analytik Jena AG, Jena, Germany) as the following procedures: 95 • C for 3 min, and 40 cycles of 95 • C for 10 s and 60 • C for 30 s. The qRT-PCR analysis of each sample was performed in three technical replicates. All of the primers were listed in Table S1. The relative expression levels of selected miRNAs were calculated using the 2 −∆∆Ct method [39]. for unique valid reads in the six grafts, respectively. The unique valid reads were mainly distributed between 21 and 24 nt, and 24-nt reads were the most abundant class ( Figure 1).

Overview of Sequencing Data
Forests 2020, 11,196 5 of 16 min, and 40 cycles of 95 °C for 10 s and 60 °C for 30 s. The qRT-PCR analysis of each sample was performed in three technical replicates. All of the primers were listed in Table S1. The relative expression levels of selected miRNAs were calculated using the 2 −ΔΔCt method [39].

Overview of Sequencing Data
A total of 21,056,396, 13,652,899, 10,785,140, 11,708,837, 15,829,590, and 18,327,912 raw reads were produced from the six grafts (Table 1)  LG: grafts with large growth increment on tall rootstocks. The statistics of length distribution were based on unique valid reads. LG: grafts with large growth increment on tall rootstocks. The statistics of length distribution were based on unique valid reads. Table 1. Analysis of small RNA sequences from the libraries of two groups of pecan grafts with significantly different growth increment. SG: grafts with small growth increment on short rootstocks; LG: grafts with large growth increment on tall rootstocks.

Identification of Conserved and Novel MiRNAs
The valid reads were compared to the reported precursor sequences in miRBase 22.0 and a total of 441 conserved miRNAs were identified from the six pecan graft libraries (Table S2), of which 357 miRNAs were identified as belonging to 42 miRNA families (Table S3). The number of members in each family ranged from 1 to 59 members. The MIR159 family had the most members, followed by MIR396 (31 members), MIR171_1 (27 members), MIR166 (26 members), and MIR156 (21 members). Four MIR166 members (mes-miR166a, cpa-miR166a, ppe-miR166a, and ptc-miR166a), two MIR482 members (ptc-miR482c-3p_2ss8GT20AT and ptc-miR482c-3p_2ss10GC20AT), mtr-miR159a, csi-miR482d-3p_1ss11AG, and gma-miR2118a-3p_R+1_2ss6GA21TA were among the top 10 conserved miRNAs by expression level in each of the six samples (Table S2). The conservation of the identified miRNAs was also analyzed using miRNA precursor statistics for other species. The results showed that there were over 100 miRNA precursors conserved between pecan and any of Glycine max L., M. domestica, Populus trichocarpa Torr. & A. Gray, Manihot esculenta Crantz, and Vitis vinifera L., while only one precursor could be matched in each of Helianthus paradoxus Heiser, Acacia mangium Willd., Helianthus ciliaris DC., Gossypium herbaceum L., Helianthus argophyllus Torr. & A. Gray, and Populus euphratica Oliv. (Table S4).
To predict novel miRNAs, the non-conserved valid reads were tested for matches to the genomic sequences of pecan. A total of 603 novel miRNAs corresponding to 579 unique precursor sequences was identified (Table S5). The lengths of these novel miRNAs ranged from 19 nt to 25 nt, with 24 nt as the most abundant. The precursors of all novel miRNAs presented a typical stem-loop hairpin secondary structure. These predicted precursors were 58-233 nt in length and their free energies were in range of −21.6 kcal/mol to −205.4 kcal/mol. The minimal folding free energy index (MFEI) of the precursors ranged from 0.9 to 2.5. Our sequencing results also showed that the maximum expression level of novel miRNAs was 3120.63 and only 35 miRNAs had expression levels of >100 reads in at least one sample, of which 12 miRNAs were expressed with reads of >100 in all of the six samples, including PC-3p-678_2589, PC-3p-396_3886, PC-5p-541_3068, PC-3p-1196_1708, and PC-5p-266_5602.

Differentially Expressed MiRNAs in Pecan Grafts with Different Growth Performance
To identify grafting-responsive miRNAs that were potentially involved in pecan graft growth regulation, we compared miRNA expression in the SG and LG groups. In total, 24 miRNAs were significantly differentially expressed (p < 0.05), which included seven conserved miRNAs and 17 novel miRNAs ( Figure 2 and Table S6). Among the conserved miRNAs, five were significantly down-regulated and two were up-regulated in SG compared to LG. There were four conserved miRNAs corresponding to miRNA families: cme-MIR319b-p5_2ss8TC21GC (MIR159 family), cme-MIR160c-p5_2ss13AG17AG (MIR160 family), cme-MIR160c-p3_2ss13AG17AG (MIR160 family), and ptc-miR399e_1ss21CT (MIR399 family). Among the novel miRNAs that were significantly differentially expressed, 16 miRNAs were significantly up-regulated and only one was down-regulated in SG compared to LG.

Go and KEGG Analysis of the Target Genes of Differentially Expressed MiRNAs
The target genes of significantly differentially expressed miRNAs were predicted by matching the miRNA sequences with pecan genomic data. A total of 95 genes were predicted to be targeted by eight of these miRNAs, among which 33 genes were annotated (Table S7). We found that cme-MIR160c-p3_2ss13AG17AG (cme-MIR160c-p5_2ss13AG17AG) targeted the highest number of genes (53), while each of two novel miRNAs (PC-3p-66417_64 and PC-3p-10652_328) matched only one target. Among the target genes, six targets (CIL1078S0020, CIL0962S0040, CIL1063S0032, CIL0937S0069, CIL0311S0011, and CIL0173S0004) were complementary to three novel miRNAs; other genes were targeted by the conserved miRNAs. A total of 33 targets had gene annotation.
To better understand the biological functions of the differentially expressed miRNAs, GO enrichment and KEGG pathway analysis were conducted. The results showed that a total of 91 targets were annotated with GO terms and one target can possess multiple GO annotations (Table S7). The GO terms can be classified into three categories: biological process, cellular component, and molecular function. The GO analysis showed that the miRNA targets were assigned to 11 biological processes, with the highest number of genes associated with cellular processes (31 genes), metabolic processes (26 genes) and response to stimulus (18 genes). The targets were involved in regulation activities at six cellular components with the most genes enriched at cell (74 genes) and membrane (29 genes). The targets also fell into six molecular functions with most genes annotated with catalytic activity (44 genes) and binding (42 genes) (Figure 3 and Table S7).

Go and KEGG Analysis of the Target Genes of Differentially Expressed MiRNAs
The target genes of significantly differentially expressed miRNAs were predicted by matching the miRNA sequences with pecan genomic data. A total of 95 genes were predicted to be targeted by eight of these miRNAs, among which 33 genes were annotated (Table S7). We found that cme-MIR160c-p3_2ss13AG17AG (cme-MIR160c-p5_2ss13AG17AG) targeted the highest number of genes (53), while each of two novel miRNAs (PC-3p-66417_64 and PC-3p-10652_328) matched only one target. Among the target genes, six targets (CIL1078S0020, CIL0962S0040, CIL1063S0032, CIL0937S0069, CIL0311S0011, and CIL0173S0004) were complementary to three novel miRNAs; other genes were targeted by the conserved miRNAs. A total of 33 targets had gene annotation.
To better understand the biological functions of the differentially expressed miRNAs, GO enrichment and KEGG pathway analysis were conducted. The results showed that a total of 91 targets were annotated with GO terms and one target can possess multiple GO annotations (Table S7). The GO terms can be classified into three categories: biological process, cellular component, and molecular function. The GO analysis showed that the miRNA targets were assigned to 11 biological processes, with the highest number of genes associated with cellular processes (31 genes), metabolic processes (26 genes) and response to stimulus (18 genes). The targets were involved in regulation activities at six cellular components with the most genes enriched at cell (74 genes) and membrane (29 genes). The targets also fell into six molecular functions with most genes annotated with catalytic activity (44 genes) and binding (42 genes) (Figure 3 and Table S7). Whereas GO analysis can identify enriched gene categories, KEGG analysis can be used for elucidation of biological pathways. The results of KEGG pathway enrichment analysis showed that a total of 34 target genes for seven differentially expressed miRNAs were classified under 28 KEGG pathways ( Figure 4 and Table S8). Those targets were mainly involved in four pathways: ATPbinding cassette (ABC) transporters, steroid biosynthesis, sesquiterpenoid and triterpenoid biosynthesis, and pyrimidine metabolism (p < 0.05).  Whereas GO analysis can identify enriched gene categories, KEGG analysis can be used for elucidation of biological pathways. The results of KEGG pathway enrichment analysis showed that a total of 34 target genes for seven differentially expressed miRNAs were classified under 28 KEGG pathways ( Figure 4 and Table S8). Those targets were mainly involved in four pathways: ATP-binding cassette (ABC) transporters, steroid biosynthesis, sesquiterpenoid and triterpenoid biosynthesis, and pyrimidine metabolism (p < 0.05). Whereas GO analysis can identify enriched gene categories, KEGG analysis can be used for elucidation of biological pathways. The results of KEGG pathway enrichment analysis showed that a total of 34 target genes for seven differentially expressed miRNAs were classified under 28 KEGG pathways ( Figure 4 and Table S8). Those targets were mainly involved in four pathways: ATPbinding cassette (ABC) transporters, steroid biosynthesis, sesquiterpenoid and triterpenoid biosynthesis, and pyrimidine metabolism (p < 0.05).

Expression Validation of Significantly Differentially Expressed MiRNAs Using qRT-PCR
We performed qRT-PCR to validate the sequencing data in the LG and SG groups. Five significantly differentially expressed miRNAs were randomly selected, consisting of two conserved (cas-miR5139_L-1 and eun-MIR482c-p3_1ss8GA) and three novel miRNAs (PC-5p-124408_31, PC-5p-42706_98, and PC-5p-96955_42). The results confirmed the existence of the novel miRNAs and that the miRNAs with low read number were detectable. It also showed that the five miRNAs were generally upregulated in SG samples compared with LG samples, which was consistent with the sequencing data ( Figure 5).

Expression Validation of Significantly Differentially Expressed MiRNAs Using qRT-PCR
We performed qRT-PCR to validate the sequencing data in the LG and SG groups. Five significantly differentially expressed miRNAs were randomly selected, consisting of two conserved (cas-miR5139_L-1 and eun-MIR482c-p3_1ss8GA) and three novel miRNAs (PC-5p-124408_31, PC-5p-42706_98, and PC-5p-96955_42). The results confirmed the existence of the novel miRNAs and that the miRNAs with low read number were detectable. It also showed that the five miRNAs were generally upregulated in SG samples compared with LG samples, which was consistent with the sequencing data ( Figure 5).

Discussion
Rootstocks are known to alter various scion behaviors in grafted plants. Rootstock-conferred effects on growth vigor has been largely elucidated in fruits trees on the dwarfing rootstocks [40]. Although gene expression patterns have been investigated in scions of different growth vigor on different rootstocks, there are no reports on the involvement of miRNAs in these differences. A single previous study has characterized pecan miRNAs, which identified 47 conserved miRNAs and 39 novel miRNAs present during graft union development [29]. The pecan genome was recently published [30] and in this study it was used in miRNA analysis to obtain results with more accuracy and reliability. In this present study, we identified the conserved and novel miRNAs in pecan and explored the differentially expressed miRNAs in the SG and LG groups. Then, we analyzed their targets to understand the functions of the miRNAs in regulation of graft growth.

miRNA Population Identified in Pecan
In the present study, 441 conserved miRNAs representing 42 families and 603 novel miRNAs were identified, indicating a larger number of miRNAs than that discovered in our earlier research [29]. This divergence may be due to use of the pecan genome to identify miRNAs, which permitted the discovery of a wider range of miRNAs than the trinity sequence in transcriptional data used as a reference in our previous work. The effects of different processes at work in pecan during graft union

Discussion
Rootstocks are known to alter various scion behaviors in grafted plants. Rootstock-conferred effects on growth vigor has been largely elucidated in fruits trees on the dwarfing rootstocks [40]. Although gene expression patterns have been investigated in scions of different growth vigor on different rootstocks, there are no reports on the involvement of miRNAs in these differences. A single previous study has characterized pecan miRNAs, which identified 47 conserved miRNAs and 39 novel miRNAs present during graft union development [29]. The pecan genome was recently published [30] and in this study it was used in miRNA analysis to obtain results with more accuracy and reliability. In this present study, we identified the conserved and novel miRNAs in pecan and explored the differentially expressed miRNAs in the SG and LG groups. Then, we analyzed their targets to understand the functions of the miRNAs in regulation of graft growth.

miRNA Population Identified in Pecan
In the present study, 441 conserved miRNAs representing 42 families and 603 novel miRNAs were identified, indicating a larger number of miRNAs than that discovered in our earlier research [29]. This divergence may be due to use of the pecan genome to identify miRNAs, which permitted the discovery of a wider range of miRNAs than the trinity sequence in transcriptional data used as a reference in our previous work. The effects of different processes at work in pecan during graft union development and graft growth on the expression of miRNAs should also be considered. According to the sequencing results, the conserved miRNAs, such as miR2118, miR395, miR168, and miR398, had not previously been detected in pecan; other conserved pecan miRNAs were not found in this experiment, but have been reported previously, such as miR860, miR818, miR862, and miR5998 [29]. The results also suggested that miR166, miR159, miR482, and miR2118 were expressed at a relatively high level in all of the six graft samples. Notably, miR166 and miR159 also had high expression levels during graft union development in pecan, and they are highly conserved across diverse plant species [41]. By analyzing the expression of the novel miRNAs, we observed that they were generally present at lower expression levels than those of conserved miRNAs, which was in agreement with previous reports [20,29,42]. Each plant species has its specific miRNAs that are not present in other, even closely related, species [43]. It can be speculated that these novel miRNAs may be involved in important regulatory pathways in pecan and further research is needed to reveal their regulatory roles.
In small RNA research, length distribution analysis is regarded as an effective assessment of the composition of sRNAs. Our results showed that the sRNAs were mainly distributed between 21 nt and 24 nt, and 24-nt sRNAs were the most abundant class, which was consistent with previous reports on pecan [28]. In some plants, such as Liriodendron chinense (Hemsl.) Sarg. [44], Gossypium hirsutum L. [45], and Camellia sinensis (L.) O. Kuntze [46], major peaks at 24-nt were also seen in sRNA length distribution patterns. Other studies, in different plants, identified 21-nt sRNAs as the most abundant class [47,48]. Additionally, differences in distribution patterns have also been observed in different lines of the same plant [49] and different developmental stages [50]. This discrepancy in the length distribution of sRNAs probably reflects variation in their responses in different species or individuals, and during different biological processes.

Grafting-Responsive MiRNAs Involved in Graft Growth and Development
In our study, a total of 24 miRNAs showed significant changes in expression levels between the LG and SG grafts (p < 0.05), which may be in response to grafting on rootstocks with different characteristics. Notably, 18 of the differentially expressed miRNAs (including 16 novel) were up-regulated in the SG grafts compared with LG grafts, with one novel miRNA down-regulated. In the study on watermelon, Liu et al. identified 45 novel miRNAs that were significantly differentially expressed in the hetero-grafting combinations [20]. These results suggests that the novel miRNAs play important roles in responding to grafting in pecan and may be involved in species-specific miRNA regulatory mechanisms. In addition, the generally opposing expression patterns of the novel and conserved miRNAs seen in this study imply different roles in response to grafting on different rootstocks. miRNAs bind to their target sequences with perfect or near-perfect complementarity and are known primarily as repressors of gene expression [51]. The analysis of targets for the differentially expressed miRNAs in two different groups of pecan grafts enables us to understand the regulatory network of growth in pecan.
Auxin is an essential regulator of plant growth and development. Its transport and changes in auxin metabolism are both involved in modulation of plant development [52]. By providing directional and positional signal information, polar transport of auxin has important effects on developmental processes, such as apical dominance, vascular differentiation, tropic growth, and organ development [53]. Polar transport of auxin between cells is active, requiring transporter molecules [54]. Auxin can be transported into cells via the H + -symport activity of the auxin permeases AUX1/LAX, in which LAX1, LAX2, and LAX3 are three proteins closely related to AUX1 [55]. Members of another family, the ATP-binding cassette (ABC) proteins such as ABCB1, ABCB4, and ABCB19 have also been shown to function as auxin transporters [56]. In the present study, it was found that LAX3 (auxin transporter-like protein 3 gene) and ABC transporter B family member 4 (ABCB4)-like gene were targeted by miR160 and cas-miR5139_L-1, respectively. Thus, interactions of these two miRNAs with auxin transporters might be involved in the processes of stem elongation and vascular differentiation by regulating the distribution of auxin in pecan grafts.
As components of SCF ubiquitin-ligase complexes (Skp I, Cullin, F-box proteins), F-box proteins function in recognizing and selectively recruiting target proteins into complexes for ubiquitin-mediated proteolysis [57]. Protein degradation is a commonly-employed mechanism for controlling protein abundance and is a particularly effective method for promoting unidirectional cell cycle transitions [58]. This study suggested that miR399 targeted the Pof1 (F-box/WD repeat-containing protein)-like gene. Pof1 is homologous to the F-box protein Met30 in budding yeast (Saccharomyces cerevisiae Hansen) [59], and SCF Met30 can trigger degradation of the transcriptional activator Met4 and regulate cell viability [60,61]. It was also reported that under normal growth conditions, SCF Pof1 is responsible for mediation of ubiquitylation and for degradation of Zip to ensure low levels. This strategy maintains normal cell division [59]. Given the existing evidence, Pof1 is thought to be closely related to the regulation of cell viability and growth through degradation of its target transcription factor, although its functions are not fully understood. Based on the previous report on G. hirsutum, cell growth related factors are involved in the regulation of plant height [45]. Therefore, it can be speculated that in our pecan grafts on different rootstocks, the altered expression of miR399 may play regulatory roles in controlling cellular activity via affecting Pof1-like gene expression and then mediated the graft growth.
Phosphorus (P) is an essential nutrient for plant growth. Although the total level of P is high in soil, the content of inorganic phosphate (Pi), an available P source, is often insufficient for plants. E2 conjugase, encoded by PHO2, is known to be a key component in maintaining Pi homeostasis [62]. Pi deprivation is reported to induce miR399, and overexpression of miR399 in Pi-replete conditions represses E2 conjugase expression and leads to high leaf Pi concentrations [63]. In this study, miR399 was not found to target PHO2, while a novel miRNA, PC-3p-10652_328, was observed to target PAH1 encoding phosphatidate phosphatase (PAP), which is associated with the availability of Pi during Pi limitation [64]. Membrane lipid remodeling is a major adaptative response to Pi deficiency, in which non-phosphorus galactolipids and sulfolipids substitute for the phospholipids so that the plant can utilize them as Pi reserves under conditions of shortage [65]. PAPs encoded by PAH1 and PAH2 are reported to be involved in the eukaryotic galactolipid biosynthesis pathway and the membrane lipid remodeling mediated by these two phosphatases is regarded as an adaptive strategy to cope with phosphate insufficiency [64]. Thus, PC-3p-10652_328 may mediate membrane lipid remodeling, which changes the available Pi levels and could regulate graft growth.

Long-Distance miRNA Transport in Grafts
The mobility of mRNA and small RNAs in phloem has attracted much attention from researchers. Several miRNAs, such as miR399, miR395, miR172 and miR156, have been identified as graft-transmissible signals sent through the phloem of plants. As Pant et al. reported, miR399 acts as a long-distance signal transmitted from shoot to root in A. thaliana to regulate phosphate homeostasis by targeting PHO2 [16]. In this study, miR399 was found to be expressed differentially by grafts on different rootstocks. One of its targets is the pof1-like gene, associated with cell activity; no targets were found to be related to Pi regulation. We thus hypothesize that miR399 might act as a graft-transmissible signal which mediates other biological activities. Other miRNAs with different expression patterns might also be long-distance signals transported between rootstocks and scions, and they may play important roles in the regulation of graft growth. In the future, we can conduct miRNA sequencing on the rootstocks and investigate miRNA mobility in grafts by comparison of miRNAs or precursor miRNAs expressed in scions and rootstocks. To better understand the long-distance regulation in grafts, further research is also needed to study their biological functions. The biological experimental methods, such as degradome sequencing and qRT-PCR, can be used to identify and validate the targets of miRNAs [66,67]. Furthermore, miRNAs overexpression or silencing can be induced to reveal the roles of the miRNAs in regulating the targets [68,69].

Conclusions
This study constructed six sRNA libraries from grafts with significantly different growth performance on short and tall rootstocks. By reference to the newly-available pecan genome, 441 conserved miRNAs and 603 novel miRNAs were identified, greatly exceeding the number found in our previous study. Of the identified miRNAs, seven conserved and 17 novel miRNAs exhibited significant differential expression patterns between growth categories, implying that the expression changes in these miRNAs may be in response to grafting and thus potentially involved in the regulation of graft growth. Analysis of the gene targets of these miRNAs suggested that miR160 and cas-miR5139_L-1 might be involved in the regulation of graft growth via auxin transport and distribution, that miR399 might play an important role in regulation of cell activity, and that PC-3p-10652_328, a novel miRNA, might participate in the regulation of available Pi.
Supplementary Materials: F The following are available online at http://www.mdpi.com/1999-4907/11/2/196/s1, Figure S1: Height of short rootstocks (SR) and tall rootstocks (HR) selected for grafting, Figure S2: Two groups (SG and LG; three grafts per group) of seedlings grafted on short and tall rootstocks (SR and TR) exhibited significantly different increment in two growth seasons (2017 and 2018), Table S1: Primers used in this study, Table S2: Conserved miRNAs identified in pecan, Table S3: MiRNA families in pecan, Table S4: Conservation of the identified miRNAs in pecan with other species, Table S5: Novel miRNAs identified in pecan, Table S6: Differentially expressed miRNAs in two groups of pecan grafts with significantly different increment, Table S7: Targets predicted for significantly differentially expressed miRNAs and GO enrichment analysis, Table S8: KEGG pathway analysis of the targets of the significantly differentially expressed miRNAs.
Author Contributions: F.P. conceived and designed the study. Z.L. collected experimental data, analyzed, and wrote the manuscript. F.L. and G.F. participated in collection of samples. P.T., K.Z., and Z.M. provided help in data analysis and improving manuscript. Y.L. provided the rootstocks. All authors have read and agreed to the published version of the manuscript.