Small RNA and Degradome Deep Sequencing Reveals the Roles of microRNAs in Seed Expansion in Peanut (Arachis hypogaea L.)

Seed expansion in peanut is a complex biological process involving many gene regulatory pathways. MicroRNAs (miRNAs) play important regulatory roles in plant growth and development, but little is known about their functions during seed expansion, or how they contribute to seed expansion in different peanut lines. We examined seed miRNA expression patterns at 15 and 35 days after flowering (DAF) in two peanut eighth-generation recombinant inbred lines (RIL8); 8106, a medium-pod variety, and 8107, a super-pod variety. Using high-throughput sequencing, we identified 1,082 miRNAs in developing peanut seeds including 434 novel miRNAs. We identified 316 differentially expressed miRNAs by comparing expression levels between the two peanut lines. Interestingly, 24 miRNAs showed contrasting patterns of expression in the two RILs, and 149 miRNAs were expressed predominantly in only one RIL at 35 DAF. Also, potential target genes for some conserved and novel miRNAs were identified by degradome sequencing; target genes were predicted to be involved in auxin mediated signaling pathways and cell division. We validated the expression patterns of some representative miRNAs and 12 target genes by qPCR, and found negative correlations between the expression level of miRNAs and their targets. miR156e, miR159b, miR160a, miR164a, miR166b, miR168a, miR171n, miR172c-5p, and miR319d and their corresponding target genes may play key roles in seed expansion in peanut. The results of our study also provide novel insights into the dynamic changes in miRNAs that occur during peanut seed development, and increase our understanding of miRNA function in seed expansion.


INTRODUCTION
As one of the most important oil crops in the world, peanut (Arachis hypogaea L.) is grown on six continents but mainly in Asia, Africa, and America (Zhang et al., 2012;Zhao et al., 2015). Yield of peanut is determined primarily by two vital factors; (1) seed expansion, and (2) subsequent seed filling (Stanciel et al., 2000;Ramakrishna et al., 2006). However, compared to other food crops, peanuts show a distinct pattern of seed development. After fertilization, a new organ called the peg differentiates from the ovary (Brennan, 1969). Little mitotic division occurs in the embryo or endosperm during active geotropic peg growth until the peg penetrates the soil.
Rapid embryo cell division then begins a few days later. This occurs approximately 10-12 days following fertilization (Smith, 1956). Seed developmental changes in water content, nucleic acid levels, enzyme activities, and storage-protein deposition patterns were studied from the earliest stage when kernels can be removed (Aldana et al., 1972;Basha and Cherry, 1976;Yamada et al., 1980). It has been shown that developmentally expressed genes in seeds are regulated at the level of transcription initiation (Falvey and Schibler, 1991). In addition, some important genes that contribute highly to seed development have also been identified in peanut Zhang et al., 2012). Although some seed development genes have been studied in peanut, the roles of the upstream regulators of these coding genes, such as miRNAs, have not been well studied.
MicroRNAs (miRNAs) are a class of small (19-24 nt), noncoding RNAs that can regulate gene expression by cleavage or translational repression of the target gene mRNAs (Van Ex et al., 2011). miRNAs have been shown to play important roles in plant growth, development, and the response to environmental stresses (Sunkar and Zhu, 2004;Ding et al., 2012;Feng et al., 2014;Zhao et al., 2016). For example, in Arabidopsis, 33 different miRNA families have been detected in mature pollen, and several (miR156, miR2939, miR158, and miR845) showed elevated expression levels in pollen compared with leaves (Grant-Downton et al., 2009). In rice, miR167 may play a role in rice grain filling through the auxin-miR167-ARF8-OsGH3.2 regulatory pathway (Xue et al., 2009); miR397 is highly expressed in young panicles and grains (Zhang et al., 2013). In peanut, studies on miRNAs also have been reported Chi et al., 2011), and several miRNAs (miR160, miR164, miR393, miR396, miR397, miR482, and miR2118) have possible roles in the response to pathogen infection (Zhao et al., 2015). Although some miRNAs have been identified in peanut, to date, the mechanism of seed expansion involving miRNAs is unclear.
To understand the role of miRNAs in peanut seed expansion, we used two representative peanut lines in this study; 8106 and 8107 are eighth-generation recombinant inbred lines (RIL8) from a cross between cultivars of Huayou7 (female) and Huayou4 (male), both of which are erect Virginia-types with 8-10 branches. The main difference between the two RILs is the pod size: line 8106 has medium-sized pods (3.2 cm long × 1.3 cm wide), and a 100-seed weight of 100 g, while line 8107 has superlarge pods (5.5 cm × 2.07 cm) with a corresponding 100-seed weight of 182 g. These two typical peanut lines were investigated to help us understand peanut seed expansion by analyzing seed miRNA expression patterns and their target genes.

Plant Materials
Seeds of the two peanut lines were sown on 15 May 2016 in the field at Henan Agricultural University (Zhengzhou, China; E113 • 41 , N34 • 49 , 94 m altitude), where the average temperature is 14.4 • C and the average rainfall is 632 mm per year. Three replicated plots of both peanut lines were planted in the field; row length was 6 m, the inter-row spacing was 40 cm, and the within-row spacing was 20 cm. Peanut pod development can be divided into two periods: pod expansion and pod filling. Generally, the first sign of pod development was seen at 15 days after flowering (DAF), and the pods enlarge to their maximum size at about 35 DAF, which was named the stereotyped fruit. The peanut pods mature at about 60 DAF (Figure 1). In our study, we selected seeds at 15 and 35 DAF in which to study seed expansion in peanut. Fresh seeds were harvested from plants of both peanut RILs at 15 DAF and 35 DAF. Four samples were frozen immediately in liquid nitrogen and stored at −80 • C until use. Three replicates of samples were collected in this study.

Small RNA Library Construction and DNA Sequencing
Total RNA was extracted using TRK-1001 following the manufacturer's instructions. RNA quantity and purity were assessed using a Bioanalyzer 2100 and the RNA 6000 Nano LabChip Kit (Agilent Technologies, Santa Clara, CA, United States) with RIN number > 7.0. Total RNA was ligated to the RNA 3 and RNA 5 adapters and reverse transcription, followed by PCR, was performed to make cDNA constructs of the small RNAs. The small cDNA fractions that ranged from 22 to 30 nt in length were then isolated via 6% denaturing polyacrylamide gel electrophoresis. Finally, the cDNA constructs were purified, and the library was validated. We then performed single-end sequencing (50 bp) on an Illumina Hiseq2500 at the LC-BIO (Hangzhou, China) following the vendor's recommended protocol.

Identification of Known and Potential Novel miRNAs
The raw sequencing reads were subjected to the Illumina pipeline filter (Solexa 0.3), and the dataset was then further processed using an in-house program, ACGT101-miR (LC Sciences, Houston, TX, United States) to remove adapter dimers, low complexity reads, common RNA families (rRNA, tRNA, snRNA, snoRNA), repeats in the NCBI GenBank 1 , Repbase 2 , and Rfam 3 ) databases. Subsequently, the remaining non-annotated sequences were mapped to specific species precursors in miRBase 21.0 4 by BLAST searches to identify known miRNAs. On the basis of sequence similarity to peanut genome, these miRNA sequences were classified into subgroups 1a, 1b, 2a, 2b, and 3 with high-to-mid confidence in order (Ambros et al., 2003;Li et al., 2017). The classification criteria for each subgroup were as follows: in gp1a, reads were mapped to miRNAs/pre-miRNAs of specific species in miRbase and the pre-miRNAs were further mapped to genome (known miRNAs); in gp1b, reads were mapped to (except for specific species) miRNAs/pre-miRNAs of selected species in miRbase and the pre-miRNAs were further mapped to genome (conserved miRNAs); in gp2a, reads were mapped to miRNAs/pre-miRNAs of selected species in miRbase and the mapped pre-miRNAs were not further mapped to a FIGURE 1 | Morphological changes in two peanut lines during seed development. 8106 and 8107 lines indicate the two peanut recombinant inbred lines (RILs). 15 DAF,25 DAF,35 DAF,45 DAF,and 60 DAF indicate 15,25,35,45, and 60 days after flowering, respectively.
genome, but the reads (and of course the miRNAs of the pre-miRNAs) were mapped to a genome. The extended genome sequences from the genomic loci may form hairpins (conserved miRNAs); in gp2b, reads were mapped to miRNAs/pre-miRNAs of selected species in miRbase and the mapped pre-miRNAs were not further mapped to a genome, but the reads (and of course the miRNAs of the pre-miRNAs) were mapped to a genome. The extended genome sequences from the genomic loci may not form hairpins (conserved miRNAs); in gp3, reads were mapped to miRNAs/pre-miRNAs of selected species in miRbase and the mapped pre-miRNAs were not further mapped to a genome, and the reads were also not mapped to a genome (conserved miRNAs). All non-annotated reads with lengths of 18-25 nt were mapped to the reference genome sequences of two diploid ancestors of cultivated peanut, A. ipaensis 5 and A. duranensis 6 using the Bowtie package; only perfectly matched sRNAs were used in further analyses. Novel miRNAs were identified using the MIREAP software 7 based on their precursors, and the hairpin RNA structures containing sequences were predicated from the flanking 120 nt sequences using RNAfold software 8 . The key criteria for miRNA prediction were based on those that had been previously reported in the literature (Meyers et al., 2008).

Reverse Transcription Reactions
Reverse transcription reactions were performed using the One Step PrimeScript miRNA cDNA Synthesis Kit (TaKaRa Co., Tokyo, Japan) following the manufacturer's instructions. Firststrand cDNA was synthesized in 20-µl reaction volumes containing 1 µg total RNA, 10 µl 2X miRNA reaction buffer mix, 2 µl 0.1% BSA, 2 µl miRNA PrimeScript RT enzyme mix, and RNase-free dH 2 O (to 20 µl). The reactions were incubated at 37 • C for 60 min, followed by 85 • C for 5 min, and then stored at −20 • C until use.

Validation of Differentially Expressed miRNAs
The qPCR reactions were performed with a SYBR PrimeScript miRNA RT-PCR Kit on a Fluorescence detection system (TianGen Biotech, Beijing, China). Each 20 µl reaction contained 1 µl cDNA template (∼100 ng), 1 µl 10 µM PCR forward primer, 1 µl 10 µM Uni-miR qPCR primer, 10 µl 2× SYBR premix EX TaqII, and 7 µl ddH 2 O. The amplification reactions were performed by first incubating at 95 • C for 5 min, followed by 45 cycles of 95 • C for 15 s, 60 • C for 30 s, and 72 • C for 45 s. Following amplification, a threshold was set and the threshold cycle (CT) was recorded automatically. All reactions were performed in triplicate for each sample. The relative expression levels of the miRNAs were calculated using the 2 − CT method (Livak and Schmittgen, 2001) and the data were normalized to the CT values for the Actin gene. The primer sequences for 16 differentially expressed miRNAs are given in Supplementary Table S13-1.

Degradome Library Construction and Target Identification
Total RNA was extracted using TRK1001 following the manufacturer's instructions. The total RNA quantity and purity were determined using a Bioanalyzer 2100 and RNA 6000 Nano LabChip Kit (Agilent, Santa Clara, CA, United States) with RIN number > 7.0. Approximately 20 µg of total RNA was used to prepare the Degradome library. The method differed considerably from previous efforts (Addo-Quaye et al., 2008, 2009b and followed the method of (Ma et al., 2010) with some modifications: (1) the poly(A)+ RNA was bound to the mRNA Capture Beads, (2) the 5 adapter was ligated to only those RNAs containing 5 -monophosphates, (3) reverse transcription was performed using 3 -random primers, (4) the library was PCR amplified, and (5) the libraries were sequenced using the 5 adapter only.
The purified cDNA library was used for cluster generation on Illumina's Cluster Station and then sequenced on the Illumina Hiseq 2500 instrument following the manufacturer's instructions. The extracted sequencing reads were stored in file SampleA_RawData.txt and were then used in e standard data analysis, which is described in Mao et al. (2012). A publicly available software package, CleaveLand3.0, was used to analyze the sequencing data generated. The key functions performed by this software and the relevant analysis results are described by Addo-Quaye et al. (2009a).

Gene Ontology (GO) and KEGG Analysis of Target Genes
Gene ontology (GO) annotations of target genes corresponding to the differentially expressed miRNAs were downloaded from the Gene Ontology 10 , NCBI 11 , and UniProt 12 . The KEGG 13 pathway was analyzed through the ClueGO plug-in 14 and Cytoscape software V2.8.2 15 to identify the significant pathways of the differential genes. GO terms and KEGG pathways were regarded to be significantly enriched with the corrected P-value ≤0.05, which was calculated using hypergeometric test (Hou et al., 2017).

Validation of the Target Genes
Expression analysis of selected target genes was performed by qPCR. First-strand cDNA was synthesized from 1 µg of RNA using TransScript First-Strand cDNA Synthesis SuperMix (TransGen Co., Beijing, China) following the manufacturer's instructions. Each 20 µl qPCR reaction contained 1 µl cDNA, 0.5 µl 10 µM forward primer, 0.5 µl 10 µM reverse primer, 10 µl 2× TransStart Top Green qPCR SuperMix, and 8 µl double-distilled water. The amplification conditions were an initial denaturation step of 95 • C for 5 min, followed by 45 cycles of 95 • C for 15 s, 56 • C for 30 s, and 72 • C for 45 s. All reactions were performed in triplicate for each sample. Relative expression levels of the targets were quantified using the 2 − CT method (Livak and Schmittgen, 2001) normalized to Actin CT values. The sequences of the primer pairs for the target genes are given in Supplementary Table S13-2.

Sequencing and Annotation of Peanut miRNAs
Lines 8106 and 8107 showed obvious difference in seed size (Figure 1). To investigate the dynamic variation of miRNAs during peanut seed expansion, two seed developmental stages, 15 DAF and 35 DAF, were selected to sequence the small RNAs using Solexa high-throughput sequencing technology. From RIL 8106 at the 15 DAF (C1) and 35 DAF (T1), we obtained 17,065,477 and 15,320,325 unfiltered sequence reads, respectively; RIL 8107 at 15 DAF (C2) and 35 DAF (T2) yielded 13,887,551 and 13,804,333 unfiltered sequence reads, respectively, from the miRNA libraries. We obtained 7,145,081, 6,008,924, 5,339,568, and 8,578,303 unique sequences from the C1, T1, C2, and T2 libraries, respectively (Supplementary Table S1). After discarding low-quality reads and further filtering the Rfam (rRNA, tRNA, snoRNA, snRNA, and other Rfam RNAs) and Repbase sequences, a total of 6,466,666, 5,121,789, 4,548,767, and 7,786,168 unique valid small RNA sequences remained, respectively (Supplementary Table S1). The lengths of the unique valid reads ranged from 18 to 25 nucleotides (nt), and the 21-24 nt sequences were predominant in the C1, T1, C2, and T2 libraries, with the 24 (nt) sequences being the most common (Supplementary Figure S1).

Identification of Known, Conserved, and Novel miRNAs
To identify the known miRNAs at different stages of seed expansion, the valid reads from each dataset were mapped to the precursor sequences in the peanut miRNA database available in miRBase (release 21). Based on the screening criteria of the miRNAs, we identified a total of 72 known miRNAs and 576 conserved miRNAs from the four small RNA libraries (Supplementary Table S2). Among these, 35 known miRNAs and 320 conserved miRNAs were shared in common (Supplementary Figures S2A,B). A total of 89 known and conserved miRNAs were confirmed in miRbase. In addition, 337 miRNAs were also confirmed, but their sequences were different from those reported in miRbase, including miR156b-5p_R+1, miR167-3p_L-1, miR394_L+1_2, miR408-5p_L+1R-1, and miR3508_2ss20CT21AT, among others. Furthermore, 222 new miRNAs were identified for newly reported 5p or 3p sequences, such as miR3510-p5_1ss22GT, miR3513-p5, miR3515-p3, miR3516-p5, and miR3517-p5 (Supplementary Table S2).
To identify the novel or predicted candidate miRNAs (PC miRNAs), the criteria for annotation of plant miRNAs (Meyers et al., 2008) were used in our study. The classification criteria for the PC miRNAs are as follows: reads were not mapped to pre-miRNAs of selected species in miRbase, but the reads were mapped to the genome and the extended genome sequences from the genome loci may form hairpins. We predicted 434 PC miRNAs from the four small libraries (332 miRNAs in C1, 332 miRNAs in T1, 241 miRNAs in C2, and 321 miRNAs in T2) (Supplementary Table S3). Among these, 172 PC miRNAs were found to be common to the four miRNA databases (Supplementary Figure S2C).

Comparison of Differentially Expressed miRNAs Between the Two Peanut RILs
We compared the frequencies of occurrence of differentially expressed miRNAs at the 15 and 35 DAF stages between the two lines based on a Poisson distribution approach (Audic and Claverie, 1997). We identified 130, 177, 139, and 175 differentially expressed miRNAs in the T1 vs. C1, T2 vs. C2, C2 vs. C1, and T2 vs. T1 comparisons, respectively (Supplementary Tables S4-1-S4-4). Among these, 26 differentially expressed miRNAs were common to the four comparisons (Figure 2). In order to dissect the miRNAs at the pod enlargement stage, we mainly focused on differentially expressed miRNAs at the 15 DAF and 35 DAF stages between the two peanut RILs. We identified 86 known or conserved miRNAs from RIL 8106 and 119 known or conserved miRNAs from RIL 8107 that were differentially expressed at the 35 DAF stage compared to 15 DAF (Supplementary Tables S4-1,4-2). Of these, 48 miRNAs were common to the two peanut lines (Supplementary Figure S3A). Moreover, 38 miRNAs and 71 miRNAs were specifically expressed in RILs 8106 and 8107, respectively (Supplementary Figure S3A). We used the following criteria to identify differentially expressed miRNAs: the adjusted p-value was <0.05 in at least one of the two comparisons. Finally, 143 miRNAs were differentially expressed at 35 DAF in the two lines (Supplementary Table S5), and 11 miRNAs showed up-regulated expression ( Table 1), while another 21 miRNAs were down-regulated ( Table 2). There are 16 miRNAs showed contrasting expression patterns in the two lines (Table 3); among them, 10 miRNAs, such as miR159a-5p_2ss19TA20TA, miR167-3p_L-1, and miR172c-5p, were down-regulated in RIL 8106 but up-regulated in RIL 8107, while six miRNAs including miR166a_1ss5AN, miR166b_L+1_1ss22AC, and miR6300_R+4_1 were down-regulated in RIL 8107 but were up-regulated in RIL 8106 ( Figure 3A). Moreover, 95 miRNAs were expressed predominantly in only one of the RILs at the 35 DAF stage (Table 4 and Supplementary  Table S6).

Validation of Differentially Expressed miRNAs by qPCR
To confirm the high-throughput sequencing data and further comparative analyses, we verified the expression patterns of 16 randomly chosen miRNAs by qPCR. The qPCR results coincided with those of the high-throughput sequencing (Figure 4). For example, we confirmed that the levels of miR156e and miR530_1ss21CA were up-regulated in both peanut lines, whereas miR164a_L+1, miR319d_L+1R-2, PC-3p-12707_579, and PC-3p-55057_146 were down-regulated in both lines. Similarly, expression of miR172c-5p was shown by both methods to be down-regulated in RIL 8106 but up-regulated in RIL 8107. Conversely, miR166b_L+1_1ss22AC was shown by both methods to be down-regulated in 8107 but up-regulated in 8106. Moreover, miR160a, miR164a, miR168a, and miR171n_1ss20TC were expressed predominantly in RIL 8107. These results indicate that the frequency of occurrence as determined by highthroughput sequencing gave a reliable prediction of miRNA expression patterns.

Identification of Target Genes of miRNAs by Degradome Analysis
Target gene validation is important to further understand the biological functions of miRNAs. In this study, we constructed     and 2,979,274 annotated peanut genes in the DS1 and DS2 libraries, respectively. We identified cleaved targets for miRNAs based on a method in the Cleaveland pipeline (Addo-Quaye et al., 2009b), in which a host gene with an alignment score of 4 or less was considered to be a potential target. In total, 1,766 and 1,616 targets were identified from the DS1 and DS2 libraries, respectively (Supplementary Tables S9, S10). For these targets, 501, 5, 859, 46, and 355 belonged to categories 0, 1, 2, 3, or 4 in the DS1 library, and 495, 45, 646, 40, and 390 belonged to a category <4 in the DS2 library (Supplementary Tables S9, S10). As expected, most of the transcripts targeted by the highly conserved miRNAs were associated with conserved target genes. For example, miR156e targets the SPL6 and SPL18 genes; miR164a_L+1 targets the NAC021 gene; miR166b_L+1_1ss22AC targets ATHB-15 gene; miR319d_L+1R-2 targets the TCP3 and TCP4 genes; miR159b_R-1_1ss7GT targets the GAM1 gene, miR160a targets the ARF10, ARF16, and ARF18 genes; miR164a targets the NAC100 gene; miR168a targets the AGO1 gene; and miR171n_1ss20TC targets the SCL22 gene. Also, the expression analysis of the 12 representative target genes by qPCR showed negative correlations with the levels of their corresponding miRNAs (Figure 5). These results show that several miRNAs may be directly or indirectly involved in peanut seed expansion by regulating expression of their target gene(s).
There were 1,401 differentially expressed target genes between the DS1 and DS2 libraries (Supplementary Table S11). Examples of "target plots" (T-plots) for the targets for peanut miRNAs are shown in Figure 6. In this T-plot, NAC100 (XM_016332293.1) and SCL22 (XM_016328632.1) were cleaved at 613 nt and 1196 nt in the two peanut lines by miR164a and miR171n_1ss20TC, respectively (Figure 6). Interestingly, the read numbers of these miRNAs that cleaved the target genes in RIL 8107 showed a significant decrease compared to that in RIL 8106, which is consistent with the expression level of miRNAs identified by RNA-Seq and qPCR (Figure 4). Also, higher expression levels of these two targets, NAC100 and SCL22 were observed in RIL 8107 compared to the levels in RIL 8106 (Figure 5). The results showed that these two miRNAs may be involved in the seed expansion process in peanut by negatively regulating NAC100 and SCL22.

Function of the Potential miRNA Targets
Gene ontology categories were assigned to 488 target genes for the 217 differentially expressed miRNAs  Table S12). Twelve classes of biological processes were identified, with the two most frequent being "regulation of transcription" and "cell division and differentiation" (Figure 3B). Moreover, other growth and development-related genes were also identified as miRNA targets, including "auxin mediated signaling pathway" and "cell cycle" (Figure 3B). These results imply the possible function of miRNAs in the regulation of biological processes involved in peanut seed expansion. In addition, based on the KEGG analysis, 208 target genes were significantly enriched in 15 pathways including "plant hormone signal transduction, " "spliceosome, " and "basal transcription factors" (Supplementary Figure S4). Examples of the plant hormone signal transduction pathway and the corresponding miRNAs are shown in Figure 3C. In this pathway, miR160a targets the ARF gene to respond to the IAA signal. Another miRNA, miR171n, is involved in ubiquitinmediated proteolysis, and controls the accumulation level of the DELLA protein in the GA pathway. These findings indicate that miRNAs have a significant effect on the regulation of peanut seed expansion by effecting hormone signaling transduction pathways.

Functions of miRNAs During Peanut Seed Expansion
Peanut seed expansion, including embryogenesis, cell division, and differentiation, and seed enlargement, is a complex biological process regulated by the expression of many genes. Previous studies have indicated that miRNAs and their predicted target mRNAs play important regulatory roles in plant growth and development (Jiang et al., 2014). Nevertheless, the function of miRNAs and their targets during peanut seed expansion has not been reported. In our study, we identified 72 known miRNAs, 576 conserved miRNAs, and 434 PC miRNAs from the four libraries (Supplementary Tables S2, S3). Interestingly, 222 miRNAs (including miR3513-p5, miR3515-p3, and miR3516-p5) were first identified in peanut, indicating that they are preferentially expressed and are specific to peanut seed development. In addition, several PC miRNAs accumulated only in the 15 DAF or 35 DAF stages, for example, PC-5p-94368_80, PC-5p-69905_114, PC-3p-63255_126, and PC-5p-118138_60, indicating that these miRNAs may play key roles in peanut seed expansion.

Presumptive miRNA-mRNA Modules Involved in Peanut Seed Expansion
In this study, we identified 11 up-regulated miRNAs ( Table 1) at 35 DAF in the two peanut lines. For example, the target gene of miR156e encodes the Squamosa promoter binding protein-like (SPL) transcription factor. Previous studies have shown that SPL family members (SPL3, SPL4, and SPL5) make a secondary contribution to the regulation of flowering and appear to function mostly in the control of flowering time and developmental phase change (Wu and Poethig, 2006;Wang et al., 2009). Another study also showed that OsSPL16 encodes a protein that promotes cell division, also with positive consequences for grain width and yield in rice . In this work, miR156e was up-regulated at the 35 DAF stage, while the expression pattern of the SPL gene was significantly down-regulated. The reduced expression of SPL might prolong the developmental phase associated with seed expansion.
Expression of 33 miRNAs showed down-regulation in the two peanut RILs at 35 DAF ( Table 2 and Supplementary Table S8-1).
The target gene of miR164a_L+1encodes the NAC transcription factor, which plays a role in the early auxin response (Guo et al., 2005;Zhang et al., 2006). Another down-regulated miR319d_L+1R-2 is a member of the miR319 family and targets TCP4-like (TCP4) genes, which encode plant-specific transcription factors involved in jasmonic acid (JA) biosynthesis (Nag et al., 2009;Schommer et al., 2014). Previous studies have demonstrated that TCP positively regulates the JA level in plants (Schommer et al., 2008;Hao et al., 2012). JA has also been reported to affect grain filling in rice (Kim et al., 2009). In addition, miR319 might influence sweet corn seed vigor (Gong et al., 2015). In our present work, we found that the level of TCP4-specific mRNA increased compared with the decrease in miR319d_L+1R-2 (Figures 4, 5). The increased expression of TCP4 may enhance seed vigor by increasing the level of JA, which could then contribute to seed expansion in peanut.
There were 24 miRNAs that showed contrasting patterns of differential expression in the two peanut lines at 35 DAF ( Table 3 and Supplementary Table S8-2). For example, miR172c-5p was down-regulated in the RIL 8106, but was up-regulated in RIL 8107 (Table 3 and Figure 4). The target gene of miR172c-5p is the splicing factor U2AF large subunit B-like (U2AF), which was up-regulated in RIL 8106 and downregulated in RIL 8107 at the 35 DAF stage (Figure 5), indicating that U2AF may also be involved in peanut seed expansion. However, miR166b_L+1_1ss22AC was up-regulated in RIL 8106, but down-regulated in RIL 8107 (Table 3 and Figure 4). miR166b_L+1_1ss22AC is a member of the miR166 family and targets the homeobox-leucine zipper protein ATHB-15 gene. In this study, we observed that expression of ATHB-15 increased with the decrease in miR166b_L+1_1ss22AC expression in RIL 8107 at 35 DAF, indicating that ATHB-15 may also play a role in peanut seed expansion by regulating early embryo development.
We also found that 95 known or conserved miRNAs ( Table 4 and Supplementary Table S6) and 54 PC miRNAs (Supplementary Table S8-3) were expressed predominantly in only one of the two peanut RILs. For example, miR159b_R-1_1ss7GT and PC-3p-79977_98 were down-regulated in RIL 8106 at 35 DAF (Table 4 and Supplementary Table S8-3). miR159 negatively regulates the expression of GAMYB genes at the posttranscriptional level; GAMYB was first identified as a downstream GA signaling target in aleurone cells of barley (Hordeum vulgare L.) (Gubler et al., 1995). MYB family members were also found to play an important role in response to the presence of abscisic acid (ABA) during seed development in rice (Aldridge and Probert, 1993;Reyes and Chua, 2007). Our results indicate that the absence of miR159b_R-1_1ss7GT in RIL 8106 lead to increased expression of the GAMYB gene, which would be consistent with its role in peanut seed development. The target gene of PC-3p-79977_98, a novel miRNA, is the ubiquitin-conjugating enzyme E2-C (UCE2) which contributes to APC-dependent protein ubiquitylation in vivo during the cell cycle (Criqui et al., 2002) and cell differentiation (Wen et al., 2006). In this study, we found that UCE2 was up-regulated in concert with the down-regulation of PC-3p-79977_98 in RIL 8106 (Figures 4, 5), indicating that PC-3p-79977_98 may also participate in peanut seed expansion by regulating the cell cycle. Our findings suggest that the strong expression of these miRNAs in RIL 8106 could be related to a specific seed development mechanism in this line.
We also found that the expression of a number of miRNAs was down-regulated only in RIL 8107, including miR160a, miR168a, and miR171n_1ss20TC, etc. (Table 4 and Figure 4). The target gene of miR160a is a member of the auxin response factor (ARF) gene family. miR160 negatively regulates ARF10 (Liu et al., 2007), ARF16 (Wang et al., 2005), and ARF17 (Mallory et al., 2005) of the repressor ARF family, and plays a critical role in maintaining the process of seed germination and the normal developmental programs for embryos, roots, leaves, and floral organs. Furthermore, ARFs regulate the expression of early auxin responsive genes, including the AUX/IAA genes (Ulmasov et al., 1997;Tiwari et al., 2003), and frees ARFs from repression by AUX/IAA proteins. The accumulation of ARF16 resulting from the down-regulation of miR160a might enhance the auxin response and thus enhance seed development in peanut. Transcription of miR168a also was significantly decreased in RIL 8107, and its target is Argonaute 1 (AGO1) which is the core component of the RNA-induced silencing complex (RISC) (Vaucheret et al., 2006) and plays crucial roles in controlling cotyledon expansion and hormone signaling (Vaucheret et al., 2004;Kidner and Martienssen, 2005). Additionally, scarecrowlike protein 22 (SCL22) was predicted to be the target gene of miR171n_1ss20TC, and it can control seed development (Kim and Ahn, 2014). We found that miR171n_1ss20TC was downregulated in RIL 8107, indicating that SCL may also be involved in peanut seed expansion. These results suggest that genotypespecific regulation of miRNAs might explain why the seeds of RIL 8107 are larger than seeds of RIL 8106. Altogether, our results suggest that at least 10 regulatory modules of these miRNAtargets play an important role in seed expansion in peanut.

Possible miRNA-Dependent Regulatory Pathways That Participate in Peanut Seed Expansion
Phytohormones can act as signaling compounds that promote and influence plant development. For example, auxin has been shown to be essential for plant growth and development by controlling cell division and cell elongation (Sauer et al., 2013). Gibberellic acid (GA) is involved in plant flowering and embryo development through its regulation of the expression of related genes (Blázquez and Weigel, 2000;White and Rivin, 2000). miR159 has been reported to regulate Arabidopsis floral development in the GA pathways; in addition, miRNA167 is involved in the auxin pathways during grain filling in rice (Achard et al., 2004;Xue et al., 2009). In this study, we also found that several miRNAs, including miR160a, miR171n, and miR156e, participate in auxin signal transduction, GA signal transduction, and BR signal transduction by regulating their target genes ( Figure 3C). These results suggest that these miRNAs may play key roles in peanut seed expansion.
A number of other conserved and novel miRNAs may also be involved in peanut seed expansion, and we further assessed the functions of their potential target genes in seed development. Examples are miR164a-mediated cleavage of NAC, miR319d-mediated cleavage of TCP, miR166b-mediated cleavage of ATHB, miR159b-mediated cleavage of GAMYB, PC-3p-79977_98-mediated cleavage of UCE2, miR168a-mediated cleavage of AGO1, and miR172c-5p-mediated cleavage of U2AF. These miRNA-mRNA modules might be involved in regulating seed expansion in peanut (Figure 7). In summary, miRNAtarget gene models, hormone transduction, and transcriptional regulation comprise a complex network that regulates seed expansion in peanut.

CONCLUSION
In this study, we found that 143 conserved miRNAs and 74 PC miRNAs were differentially expressed in the two peanut RILs at the 35 DAF stage. Of these, expression of 11 miRNAs was found to be up-regulated, while 33 miRNAs were down-regulated. Twenty-four miRNAs displayed contrasting expression patterns in the two peanut lines. Moreover, 149 miRNAs were expressed predominantly in only one of the two RILs. We also validated the expression of a number of representative miRNAs by qPCR. The potential miRNA target genes were verified by sequencing the two degradome libraries. Our results demonstrate that the expression patterns of some miRNAs can be very different even between two closely related inbred lines from the same species. The differentially expressed miRNAs and their target genes identified in this study could be important to the regulatory networks that control seed expansion in peanut.

AUTHOR CONTRIBUTIONS
DY designed and conceived the research. XM and ZX wrote the manuscript. XM and XZ analyzed the data. KZ, FL, KL, and LN performed the experiments. ZX, JH, and DY edited the manuscript. All authors read and approved the final manuscript.

FUNDING
This work was financially supported by grants from the National Natural Science Foundation of China (No. 31471525) and Key Scientific and Technological Project in Henan Province (Nos. 161100111000 and S2012-05-03).

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2018.00349/ full#supplementary-material FIGURE S1 | Size distribution of peanut small RNAs. C1 and C2 are small RNAs isolated from RILs 8106 and 8107 at 15 DAF, respectively. T1 and T2 are the small RNAs isolated from RILs 8106 and 8107 at 35 DAF, respectively.                S13 | All primer sequences in this study. (S13-1) Primer sequences used for qPCR analysis of 16 differentially expressed miRNAs. (S13-2) Primer pair sequences used for qPCR analysis of 12 target genes.