Identification and Expression of miRNAs Related to Female Flower Induction in Walnut (Juglans regia L.)

Flower induction is an essential stage in walnut (Juglans regia L.) trees, directly affecting yield, yield stability, fruit quality and commodity value. The objective of this study was to identify miRNAs related to female flower induction via high-throughput sequencing and bioinformatics analysis. A total of 123 miRNAs were identified including 51 known miRNAs and 72 novel miRNAs. Differential expression was observed in 19 of the known miRNAs and 34 of the novel miRNAs. Twelve miRNAs were confirmed by RT-qPCR. A total of 1339 target genes were predicted for the differentially expressed miRNAs. The functions of 616 of those target genes had been previously annotated. The target genes of the differentially expressed miRNAs included: (i) floral homeotic protein APETALA 2 (AP2) and ethylene-responsive transcription factor RAP2-7 which were targeted by jre-miRn69; (ii) squamosa promoter-binding protein 1 (SPB1) and various SPLs (squamosa promoter-binding-like protein) which were targeted by jre-miR157a-5p; (iii) various hormone response factors which were targeted by jre-miR160a-5p (ARF18) and jre-miR167a-5p (ARF8) and (iv) transcription factor SCL6 which was targeted by jre-miR171b-3p, jre-miRn46 and jre-miRn49. The KEGG pathway analysis of the target genes indicated that the differentially expressed miRNAs were mainly enriched to ubiquitin mediated proteolysis, RNA degradation and various carbohydrate metabolism pathways. Many miRNAs were detected in J. regia during female flower induction. Some miRNAs (jre-miR157a-5p, jre-miR160a-5p, jre-miR167a-5p, miR171b-3p jre-miRn69 and jre-miRn49) were involved in female flower induction. The results of this experiment will contribute valuable information for further research about the function of miRNAs in flower induction of J. regia and other fruit trees.


Introduction
Walnut (Juglans regia L.) is an economically important nut tree, with more than 7000 years of evolutionary and domestication history in China [1,2]. The development of J. regia flowers begins with flower induction, flower initiation, and flower differentiation during the first growing season and ends with blooming in the second growing season. Flower induction involves the anatomical and morphological transition from vegetative meristems to reproductive meristems [3]. The date, intensity and quality of flower induction and differentiation directly influence the fruit tree's early fruiting, fruiting stability, fruit quality and commodity value. Flower induction in woody trees is closely related to the time when new shoots stop growing [4]. In J. regia, the induction of female flowers occurs 3 to 7 weeks after medium and short branches cease vegetative growth. Initiation and differentiation Table 1. Summary of small RNA sequencing and annotation in the four J. regia libraries. F_1, female flower buds before flower induction; F_2, female flower buds during flower induction; F_3, female flower buds after flower induction; JRL, leaf buds during flower induction. RNA categories, including known miRNA, novel miRNA, rRNA, tRNA, small nuclear RNAs (snRNA), small nucleolar RNAs (snoRNA), and ta-siRNA (Table 1). An average of 837,129 sRNAs (18.79% of the mapped sRNAs) were identified as known miRNAs, and an average of 92,751 reads (2.06%) were identified as novel miRNAs.

Identification of Known miRNAs in J. regia
To identify the known miRNAs in J. regia, mapped sRNAs in the transcriptome were compared to other plant miRNAs in the miRBase database (release 21). The results showed that 3,348,517 total reads representing 1013 unique sRNAs matched known miRNAs. Of these unique sRNAs, 51 mature miRNAs and 60 precursors were identified (Table S1). All of the precursors were able to adopt hairpin structures resembling the fold-back structure of miRNA. The analysis of nucleotide bias of known miRNAs showed that uracil appeared with high frequency at the 5′-end ( Figure S1). Among the known miRNAs, 38 miRNAs were detected in all four libraries. Two miRNAs (jre-miR170-3p and jre-miR166a-3p) were only identified in F_1. Jre-miR167c-5p was specific to F_2 and jre-miR393a-3p was specific to F_3 (Table S1). The expression levels of some miRNA families, such as MIR159 and MIR319, were much higher than those of other families (Figure 2A). The most abundant miRNA was jre-

Identification of Known miRNAs in J. regia
To identify the known miRNAs in J. regia, mapped sRNAs in the transcriptome were compared to other plant miRNAs in the miRBase database (release 21). The results showed that 3,348,517 total reads representing 1013 unique sRNAs matched known miRNAs. Of these unique sRNAs, 51 mature miRNAs and 60 precursors were identified (Table S1). All of the precursors were able to adopt hairpin structures resembling the fold-back structure of miRNA. The analysis of nucleotide bias of known miRNAs showed that uracil appeared with high frequency at the 5 -end ( Figure S1). Among the known miRNAs, 38 miRNAs were detected in all four libraries. Two miRNAs (jre-miR170-3p and jre-miR166a-3p) were only identified in F_1. Jre-miR167c-5p was specific to F_2 and jre-miR393a-3p was specific to F_3 (Table S1). The expression levels of some miRNA families, such as MIR159 and MIR319, were much higher than those of other families (Figure 2A). The most abundant miRNA was jre-miR159a, accounting for 70.72, 71.09, 70.66 and 67.66% of the total sequences in F_1, F_2, F_3 and JRL, respectively. Jre-miR159b-3p, jre-miR166a-3p, jre-miR319a, jre-miR396b-5p, jre-miR319c and jre-miR170-5p had high expression levels, with more than 10,000 reads in each of the four libraries. Jre-miR159c, jre-miR162a-3p, and jre-miR396a-5p had moderate abundance, with more than 1000 reads. There were significant differences in expression levels among different members within the same miRNA family. For instance, in the MIR166 family, the mean expression of jre-miRNA166a-3p was 70,522 reads, whereas the mean expression of jre-miR166a-5p was only 0.25 reads. An average of 10,129 reads was obtained for jre-miR170-5p; however, only one read was identified for jre-miR170-3p. Furthermore, jre-miR166a-5p and jre-miR170-3p were only found in F_1. Jre-miR159b-3p, jre-miR166a-3p, jre-miR319a, jre-miR396b-5p, jre-miR319c and jre-miR170-5p had high expression levels, with more than 10,000 reads in each of the four libraries. Jre-miR159c, jre-miR162a-3p, and jre-miR396a-5p had moderate abundance, with more than 1000 reads. There were significant differences in expression levels among different members within the same miRNA family. For instance, in the MIR166 family, the mean expression of jre-miRNA166a-3p was 70,522 reads, whereas the mean expression of jre-miR166a-5p was only 0.25 reads. An average of 10,129 reads was obtained for jre-miR170-5p; however, only one read was identified for jre-miR170-3p. Furthermore, jre-miR166a-5p and jre-miR170-3p were only found in F_1.

Identification of Novel miRNAs in J. regia
Stem-loop structures, which are unique to miRNAs, were used to predict novel miRNAs in this study. A total of 371,006 novel miRNA reads were identified from the four libraries, representing 3489 unique sRNAs. A total of 75 miRNA precursors and 72 novel miRNAs were predicted from the four libraries (Table S2). The length of these potential novel miRNA precursors ranged from 65 to 295

Identification of Novel miRNAs in J. regia
Stem-loop structures, which are unique to miRNAs, were used to predict novel miRNAs in this study. A total of 371,006 novel miRNA reads were identified from the four libraries, representing 3489 unique sRNAs. A total of 75 miRNA precursors and 72 novel miRNAs were predicted from the four libraries (Table S2). The length of these potential novel miRNA precursors ranged from 65 to 295 nt. Most (47 of 72) of the novel miRNAs were located in the 5 -end of the stem-loop. The minimal folding free energy (MFE) of these precursors ranged from −14.91 to −202 kcal/mol, with an average of −54.69 kcal/mol. We also detected the complementary miRNA* sequences for each novel miRNA. The results showed that 45 of 72 novel miRNAs had complementary miRNA*s, providing more evidence that these miRNAs were real. In addition, most miRNA*s were present at lower copy numbers than those of their corresponding miRNAs. The length of the novel miRNAs ranged from 18-to 25-nt. The 21-and 24-nt miRNAs were most abundant, accounting for 37.5 and 34.72% of the total novel miRNAs, respectively. Uracil was the first nucleotide in nearly half of the novel miRNAs ( Figure S1). The expression levels of most novel miRNAs were generally less than those of the known miRNAs. The exceptions were jre-miRn3, jre-miRn16 and jre-miRn17, which all had more than 10,000 reads. Figure 2B shows the ten novel miRNAs with the highest expression. The most abundant miRNA was jre-miRn3 which accounted for 34.83% of the total novel reads in F_1, 48.52% in F_2, 43.71% in F_3, and 33.88% in JRL.
A total of 24 families were identified among these known and novel miRNA precursors. These families were conserved across 69 plant species (Table S3). In J. regia, the biggest family was MIR159/319 with nine members, followed by MIR156/157 and MIR171/170 which both had eight members. Moreover, MIR156/157 was found in 47 plant species, making it the most common family. MIR159/319, MIR166 and MIR396 were found across 42 plant species. The MIR171_1, MIR160 and MIR172 families were identified in 37, 36, and 33 different plant species, respectively.

Target Prediction and Function Analysis of miRNAs in J. regia
To understand the function of the differentially expressed miRNAs, their target genes were predicted by psRobot (http://omicslab.genetics.ac.cn/psRobot/target_prediction_1.php). A total of 1339 target genes were predicted for the 19 known miRNAs and 33 novel miRNAs. No target gene was predicted for jre-miRn44 (Table S5).
Most of the miRNAs targeted more than one gene, suggesting that these miRNAs had diverse regulatory roles. Jre-miRn3 targeted only one gene, indicating that it played a unique regulatory function. The predicted target genes were searched against the non-redundant database. A total of 616 of these target genes had an annotated function. These predicted target genes belonged to numerous transcription factors and genes with different biological functions. The target genes were involved in various functions, including floral development and flowering time. For example, floral homeotic protein APETALA 2 (AP2) and ethylene-responsive transcription factor RAP2-7 were predicted to be modulated by jre-miRn69 (Table 2 and Table S5). Squamosa promoter-binding protein 1 (SPB1) and squamosa promoter-binding-like protein (SPL4/6/9/13A/16/18) were predicted to be targeted by jre-miR157a-5p. Additionally, 16 target unigenes were involved in the hormone signaling pathway. Specifically, auxin response factor 18 (ARF18) was predicted to be targeted by jre-miR160a-5p, and auxin response factor 8 (ARF8) was predicted to be targeted by jre-miR167a-5p. Scarecrowlike protein 6 (SCL6) was predicted to be targeted by jre-miR171b-3p, jre-miRn46 and jre-miRn49. Jre-miRn49 was also predicted to target indole-3-acetic acid-amido synthetase GH3.1.

Target Prediction and Function Analysis of miRNAs in J. regia
To understand the function of the differentially expressed miRNAs, their target genes were predicted by psRobot (http://omicslab.genetics.ac.cn/psRobot/target_prediction_1.php). A total of 1339 target genes were predicted for the 19 known miRNAs and 33 novel miRNAs. No target gene was predicted for jre-miRn44 (Table S5).
Most of the miRNAs targeted more than one gene, suggesting that these miRNAs had diverse regulatory roles. Jre-miRn3 targeted only one gene, indicating that it played a unique regulatory function. The predicted target genes were searched against the non-redundant database. A total of 616 of these target genes had an annotated function. These predicted target genes belonged to numerous transcription factors and genes with different biological functions. The target genes were involved in various functions, including floral development and flowering time. For example, floral homeotic protein APETALA 2 (AP2) and ethylene-responsive transcription factor RAP2-7 were predicted to be modulated by jre-miRn69 (Table 2 and Table S5). Squamosa promoter-binding protein 1 (SPB1) and squamosa promoter-binding-like protein (SPL4/6/9/13A/16/18) were predicted to be targeted by jre-miR157a-5p. Additionally, 16 target unigenes were involved in the hormone signaling pathway. Specifically, auxin response factor 18 (ARF18) was predicted to be targeted by jre-miR160a-5p, and auxin response factor 8 (ARF8) was predicted to be targeted by jre-miR167a-5p. Scarecrow-like protein Molecules 2018, 23, 1202 7 of 15 6 (SCL6) was predicted to be targeted by jre-miR171b-3p, jre-miRn46 and jre-miRn49. Jre-miRn49 was also predicted to target indole-3-acetic acid-amido synthetase GH3.1. The Gene Ontology (GO) enrichment of the target genes showed that the differentially expressed miRNAs were involved in various biological processes and pathways. A total of 1612 GO terms were enriched (Table S6) and 302 significant GO terms (collected p-value < 0.05) were identified, including 185 "biological process" terms, 54 "cellular component" terms and 63 "molecular function" terms. Among the biological processes, most of the enriched GO terms were involved in biological regulation, including regulation of cellular process, macromolecule metabolic process, primary metabolic process, cellular metabolic process, and nitrogen compound metabolic process. For the cellular component category, the most represented GO terms were related to "cell", "cell part", "intracellular" and "intracellular part". Based on the molecular function category, the represented GO terms were related to "nucleic acid binding transcription factor activity", "transcription factor activity, sequence-specific DNA binding", "ADP binding", "nucleoside-triphosphatase activity" and "pyrophosphatase activity". GO terms were related to "nucleic acid binding transcription factor activity", "transcription factor activity, sequence-specific DNA binding", "ADP binding", "nucleoside-triphosphatase activity" and "pyrophosphatase activity". The Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis revealed 65 biochemical pathways. With the threshold of corrected p-value at <0.05, six pathways were significantly enriched involving 50 unigenes (Table S7). "Ubiquitin mediated proteolysis" (ko04120) was the most significantly enriched term with 21 unigenes, followed by "RNA degradation" (ko03018) with 10 unigenes, and "Other glycan degradation" (ko00511) and "Inositol phosphate metabolism" (ko00562) both with 7 unigenes (Figure 5). In addition, some pathways regulated by differentially expressed miRNAs were not significantly enriched; however they may still be important for flower induction in J. regia. Some examples include "plant circadian rhythm" (ko04712), "photosynthesis" (ko00195 and ko00196), "plant hormone signal transduction" (ko04075) and various carbohydrate metabolism pathways (ko00500, ko00040, ko00520, ko00030, ko00051, ko00053, ko00630, ko00650 and ko00010). The Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis revealed 65 biochemical pathways. With the threshold of corrected p-value at <0.05, six pathways were significantly enriched involving 50 unigenes (Table S7). "Ubiquitin mediated proteolysis" (ko04120) was the most significantly enriched term with 21 unigenes, followed by "RNA degradation" (ko03018) with 10 unigenes, and "Other glycan degradation" (ko00511) and "Inositol phosphate metabolism" (ko00562) both with 7 unigenes (Figure 5). In addition, some pathways regulated by differentially expressed miRNAs were not significantly enriched; however they may still be important for flower induction in J. regia. Some examples include "plant circadian rhythm" (ko04712), "photosynthesis" (ko00195 and ko00196), "plant hormone signal transduction" (ko04075) and various carbohydrate metabolism pathways (ko00500, ko00040, ko00520, ko00030, ko00051, ko00053, ko00630, ko00650 and ko00010).

Discussion
The sRNAs ranging from 21 to 24-nt were common in this study, which was typical of mature miRNAs in plants [20,21]. The 24-nt sRNA had the highest abundance in all four libraries, accounting for an average of 50.32% of the total reads. Similar results were obtained in other species, such as Arabidopsis [22], Citrus trifoliate [23], C. sativus [24], Carya cathayensis [19], Malus hupehensis [25], Pyrus bretschneideri [26] and Triticum turgidum [27]. There was a difference in the number of 21-nt and 24-nt sRNAs between leaf buds (JRL) and female flower buds (F_1, F_2 and F_3) (Figure 1), suggesting that the 21-nt and 24-nt sRNAs had different regulatory functions during flower induction. The analysis of nucleotide bias showed that uracil dominated the first three positions in the known miRNAs. In addition, uracil was the first nucleotide in nearly half of the novel miRNAs (Table S1). This result was similar to previous studies about Populus [28], Glycine max [29] and C. cathayensis [19].
The expression of miRNA families differs across species (Table S3). In J. regia, jre-miR159a was the most abundant miRNA, followed by jre-miR159b-3p, jre-miR166a-3p, jre-miR319a, jre-miR396b-5p, jre-miR319c and jre-miR170-5p. All of these miRNAs had >10,000 reads, suggesting that they had important roles in flower induction. MIR159/319 was the most abundant family with nine members, followed by MIR156/157 and MIR171/170 which both had eight members. A previous report indicated that MIR166 was the largest family in C. cathayensis (13 members) [19] and Cymbidium ensifolium (eight members) [21]. In this study, MIR166 had two members. Interestingly, there were significant differences in the expression of miRNAs within the same family. For example, jre-miR166a-3p had an average of 70,522 reads, whereas jre-miR166a-5p only had one read (the read was

Discussion
The sRNAs ranging from 21 to 24-nt were common in this study, which was typical of mature miRNAs in plants [20,21]. The 24-nt sRNA had the highest abundance in all four libraries, accounting for an average of 50.32% of the total reads. Similar results were obtained in other species, such as Arabidopsis [22], Citrus trifoliate [23], C. sativus [24], Carya cathayensis [19], Malus hupehensis [25], Pyrus bretschneideri [26] and Triticum turgidum [27]. There was a difference in the number of 21-nt and 24-nt sRNAs between leaf buds (JRL) and female flower buds (F_1, F_2 and F_3) (Figure 1), suggesting that the 21-nt and 24-nt sRNAs had different regulatory functions during flower induction. The analysis of nucleotide bias showed that uracil dominated the first three positions in the known miRNAs. In addition, uracil was the first nucleotide in nearly half of the novel miRNAs (Table S1). This result was similar to previous studies about Populus [28], Glycine max [29] and C. cathayensis [19].
The expression of miRNA families differs across species (Table S3). In J. regia, jre-miR159a was the most abundant miRNA, followed by jre-miR159b-3p, jre-miR166a-3p, jre-miR319a, jre-miR396b-5p, jre-miR319c and jre-miR170-5p. All of these miRNAs had >10,000 reads, suggesting that they had important roles in flower induction. MIR159/319 was the most abundant family with nine members, followed by MIR156/157 and MIR171/170 which both had eight members. A previous report indicated that MIR166 was the largest family in C. cathayensis (13 members) [19] and Cymbidium ensifolium (eight members) [21]. In this study, MIR166 had two members. Interestingly, there were significant differences in the expression of miRNAs within the same family. For example, jre-miR166a-3p had an average of 70,522 reads, whereas jre-miR166a-5p only had one read (the read was in F_1). Jre-miR170-5p had an average of 10,129 reads, and jre-miR170-3p was only found in F_1 with one read. The large differences among miRNA family members reflect the diversity of their potential physiological roles during flower induction.
We identified 72 novel miRNAs and 45 novel miRNA*s which corresponded to 75 precursors. The lengths of these potential novel miRNA precursors ranged from 65-to 295-nt, which was close to the lengths of miRNA precursors in Musa acuminata [30]. The average MFE of these miRNA precursors was −54.69 kcal/mol. This value was consistent with the low MFE of miRNA precursors in P. bretschneideri (−51.48 kcal/mol) [26] and Poncirus trifoliate (−52.41kcal/mol) [18] but less than that in C. cathayensis (−36.9 kcal/mol) [19]. According to the criteria for annotation of novel plant miRNA, the detection of miRNA* gave strong biogenesis evidence for the novel miRNAs [31]. In this study, we detected 45 miRNA*s. The expressions of these miRNA*s were significantly lower than their corresponding miRNAs, which agreed with the previous report that miRNA* sequences were degraded and therefore usually occurred at significantly lower levels than mature miRNAs [19,32,33].
Analysis of the target genes for the differentially expressed miRNAs revealed 1339 potential unigenes. The functions of 616 of the target genes had been previously annotated, and many genes were involved in flower development and flowering time. As reported in previous studies, miR156 regulated squamosa promoter-binding (SPB) genes and miR172 regulated AP2 genes [10,11]. The over-expression of miR156 prolonged the expression of juvenile vegetative traits and delayed flowering both in Arabidopsis [10,11,13] and in Zea mays [12]. In contrast, the over-expression of miR172 accelerated flowering in Arabidopsis. Four MIR156/157 members and five MIR172 members were identified in our study. Only jre-miR157a-5p was significantly up-regulated in F_1 vs. F_2. Jre-miRn69, which was up-regulated in F_1, F_2 and F_3 compared with JRL, was predicted to regulate AP2 and RAP2-7. The result, which was validated by RT-qPCR (Figure 4), suggested that jre-miR157a-5p and jre-miRn69 both play important roles in J. regia female flower induction.
In Arabidopsis, miR167 targeted ARF6 and ARF8, which regulated gynoecium and stamen development in immature flowers [34,35]. In C. cathayensis, cca-miR167 was up-regulated during the late stage of flower development, suggesting that the ARF family is important for flower development [19]. As shown in Table 2, 16 unigenes were involved in various hormone signaling pathways. The results of RT-qPCR showed that (i) jre-miR160a-5p was significantly down-regulated both in F_2 vs. F_1 and in F_3 vs. F_1 and (ii) jre-miR167a-5p was significantly up-regulated in F_1 vs. JRL. This suggests that jre-miR160a-5p and jre-miR167a-5p both play important roles during female flower induction in J. regia.
Transcription factor SCL, which belongs to the GRAS protein family, is targeted by miR171 and plays important roles in flowering control and apical meristem development in Arabidopsis [36]. In our study, SCL6 was targeted by jre-miR171b-3p, jre-miRn46, and jre-miRn49. The result of RT-qPCR confirmed that jre-miR171b-3p and jre-miRn49 were both down-regulated in F_2 vs. F_1, suggesting that jre-miR171b-3p and jre-miRn49 may play important roles in flower induction.

Plant Materials
The J. regia (cv. Xinxin2) trees in this study were growing at the Luntai Plant Germplasm Resource Garden, Xinjiang Academy of Agricultural Science. Five trees were selected for sampling. Female flower buds were collected before (F_1), during (F_2) and after (F_3) the female flower induction. Leaf buds (JRL) were collected during female flower induction. The samples from the five trees were pooled, immediately frozen in liquid N, and then stored at −70 • C.

Small RNA Library Construction and Deep Sequencing
Total RNA was extracted from the female flower buds (F_1, F_2 and F_3) and leaf buds (JRL) using a Trizol kit (Invitrogen, Carlsbad, CA, USA). The plant samples for extraction each contained three flower or leaf buds. Total RNA was extracted from three samples (replications) for each flower development stage and then pooled. A 3 µg aliquot of the pooled sample of total RNA was used for sequencing. The four sequencing libraries were generated using NEBNext ® Multiplex Small RNA Library Prep Set for Illumina ® (NEB, Acton, MA, USA) following the manufacturer's recommendations. Index codes were added to attribute sequences to each sample. Briefly, NEB 3 SR Adaptor was directly and specifically ligated to the 3 end of the small RNAs. After the 3 ligation reaction, the SR RT Primer was hybridized to the excess 3 SR Adaptor, transforming the single-stranded DNA adaptor into a double-stranded DNA molecule. Next, a 5 end adaptor was ligated to the 5 end of the small RNAs. The first strand cDNA was synthesized using M-MuLV Reverse Transcriptase (RNase H − ). Next, PCR amplification was performed and the products were purified on 8% polyacrylamide gel (100 V, 80 min). The cDNA fragments with length of 140-160 bp were recovered and dissolved in 8 µL elution buffer. Finally, the quality of the cDNA libraries was assessed on the Bioanalyzer 2100 system (Agilent, Palo Alto, CA, USA) using DNA High Sensitivity Chips. The clustering of the index-coded samples was performed on a cBot Cluter Generation System using a Truseq SR Cluster Kit v3-cBot-Hs (Illumina, San Diego, CA, USA) according to the manufacturer's instructions. The prepared libraries were sequenced on an Illumina Hiseq 2500/2000 platform. All of the processes were carried out by Novogene (Beijing, China).

Bioinformatics Analysis of Small RNAs
After discarding adapter contaminants and low-quality reads from the raw data, clean reads with lengths between 18 and 30-nt were chosen for analysis. The sRNAs were mapped to the J. regia transctriptome using the Bowtie method [37] without mismatch. Mapped sRNAs were used to search for known miRNA in the miRBase database (Release 21) (http://www.mirbase.org/). Modifications of mirdeep2 [38] and srna-tools-cli [39] software were used to obtain the potential miRNAs and to draw the secondary structures. We removed the sequences matching non-coding RNAs including rRNAs, tRNAs, snRNAs, and snoRNAs deposited in the Rfam database (http://rfam.xfam.org/). MiREvo [40] and mirdeep2 were integrated to predict novel miRNA by exploring the secondary structure, the Dicer cleavage sites, and the minimum free energy of sRNAs that were not annotated in the previous steps.

Expression Analysis of miRNA
The miRNA expression levels were estimated by TPM (transcript per million) using the following criteria [41]: Normalized expression = mapped read count/total reads × 1,000,000. Differential expression analysis of two samples was performed using the DEGseq (2010) R package [42]. The p-value was adjusted using the q-value. A q-value < 0.005 and |log2 (fold change)| > 1 were set as thresholds for significant differential expression.

Validation and Expression of miRNAs by RT-qPCR
The expression levels of 12 of the most promising miRNAs were determined using RT-qPCR. Total RNA samples were the same as the sRNA-seq samples. The addition of poly (A + ) tails to sRNAs by poly (A + ) polymerase and the synthesis of cDNA were both performed using a miRcute Plus miRNA First-strand cDNA Synthesis Kit (Tiangen, Beijing, China). The qPCRs were performed using a miRcute Plus miRNA qPCR Detection Kit (Tiangen) in a 20 µL reaction mixture consisting of 2 µL of cDNA (ten-fold dilution), 0.4 µL of 10 µM forward primer and 0.4 µL of 10 µM reverse primer and 10 µL of 2× miRcute Plus miRNA Premix (with SYBR&ROX). The reactions were carried out as follows: 95 • C for 15 min, and 40 cycles of 94 • C for 20 s and 60 • C for 34 s. Each sample was processed three times. The reference gene was chosen based on the stability analysis of eleven putative reference genes: jre-miR394a, jre-miR159a, jre-miR159c, jre-miR167d, jre-miR159b-3p, jre-miRn3, jre-miRn7, jre-miRn59, jre-miR142, U6 and 5.8s rRNA (data not shown). Jre-miR394a was the most stable gene in both the leaf buds and the female flower buds and was therefore used as the reference gene in this study. The relative expression levels of the miRNAs were quantified using the 2 −∆∆Ct method [43]. The expression level of F_1 was set to 1. The primer sequences are listed in Table S8.

Target Gene Prediction and Function Annotation
The putative target genes of differentially expressed miRNAs were predicted by psRobot_tar in psRobot [44]. To determine the function of the target genes, GO enrichment analysis was performed using the GOseq R package [45]. The GO terms with corrected p-values < 0.05 were considered significantly enriched. Pathway analysis of the target genes was performed using the KEGG database (http://www.genome.ad.jp/kegg/) [46]. KOBAS software [47] was used to test the statistical enrichment of the target genes in the KEGG pathways at a corrected p-value < 0.05.

Conclusions
In summary, high-throughput sequencing was used to identify 51 known miRNAs and 72 novel miRNAs from four small RNA libraries from female flower buds and leaf buds of J. regia. The miRNAs belonged to 24 miRNA families and were conserved across 69 plant species. Pairwise analysis indicated that 19 of the known miRNAs and 34 of the novel miRNAs were differentially expressed. Twelve miRNAs were validated by RT-qPCR and the results were similar to those detected by sRNA-seq. A total of 1339 target genes were predicted for the differentially expressed miRNAs. The functions of 616 of the target genes have been annotated. The target genes included (i) AP2 and ethylene-responsive transcription factor RAP2-7 targeted by jre-miRn69; (ii) SPB1 and various SPLs targeted by jre-miR157a-5p; and (iii) various hormone responding factors targeted by jre-miR160a-5p and jre-miR167a-5p. Transcription factor SCL6 was targeted by jre-miR171b-3p, jre-miRn46 and jre-miRn49. The KEGG pathway analysis of the target genes indicated that the differentially expressed miRNAs were mainly enriched to ubiquitin mediated proteolysis, RNA degradation and various carbohydrate metabolism pathways. These research results contribute valuable information for further research about the function of miRNAs in the flower induction of J. regia and other fruit trees.
Supplementary Materials: The following are available online. Figure S1: Examination of nucleotide bias within known miRNAs; Table S1: Known miRNAs in four Juglans regia libraries; Table S2: Novel miRNA in four Juglans regia libraries; Table S3: The conserved miRNA families among species; Table S4: Differentially expressed miRNAs identified by pairwise analysis; Table S5: Target gene prediction and function annotation for differentially expressed miRNAs; Table S6: GO enrichment for target genes of differentially expressed miRNAs; Table S7: KEGG pathway analysis for target genes of differentially expressed miRNAs; Table S8: Reverse transcription quantitative real-time PCR (RT-qPCR) primer sequences for candidate miRNAs.
Author Contributions: J.X. and L.Z. conceived and designed the experiments; L.Z. performed the experiments, analyzed the data and wrote the paper; S.Q., H.X., and L.M. contributed materials and data analysis.