Comparison of miRNAs and Their Targets in Seed Development between Two Maize Inbred Lines by High-Throughput Sequencing and Degradome Analysis

MicroRNAs (miRNAs) play an important role in plant growth, development, and response to environment. For identifying and comparing miRNAs and their targets in seed development between two maize inbred lines (i.e. PH6WC and PH4CV), two sRNAs and two degradome libraries were constructed. Through high-throughput sequencing and miRNA identification, 55 conserved and 24 novel unique miRNA sequences were identified in two sRNA libraries; moreover, through degradome sequencing and analysis, 137 target transcripts corresponding to 38 unique miRNA sequences were identified in two degradome libraries. Subsequently, 16 significantly differentially expressed miRNA sequences were verified by qRT-PCR, in which 9 verified sequences obviously target 30 transcripts mainly involved with regulation in flowering and development in embryo. Therefore, the results suggested that some miRNAs (e.g. miR156, miR171, miR396 and miR444) related reproductive development might differentially express in seed development between the PH6WC and PH4CV maize inbred lines in this present study.

Maize (Zea mays L.), one of the most important crops in the world, is widely used as a model plant for biological research [22]. Over recent decades, several published reports about miRNAs in maize have focused on many biological processes, including leaf development, root development, seed development and response to stresses [23][24][25][26][27][28][29][30][31]. For instance, Juarez et al. showed miR166 constituted a highly conserved signal in maize leaf development [23]. Zhang et al. revealed that submergence-responsive miRNAs were involved in the regulation of metabolic, physiological and morphological adaptations of maize roots [24]. In addition, Kang et al. identified 125 and 127 known miRNAs from seeds and leaves in maize [25]. Ding et al. reported that 34 miRNAs belonging to 20 miRNA families were obtained in germinating maize seeds by high-throughput sequencing [26]. Furthermore, Li et al. used high-throughput sequencing to find that diverse and complex miRNAs were involved in the seed imbibition process [27]. Through degradome sequencing, Liu et al. detected that 72 genes targeted by 62 differentially expressed miRNAs might be attributed to the development of maize ears [28]. Jin et al. researched the dynamic expression patterns of miRNAs at 4 distinct developmental grain filling stages in maize [29]. Moreover, Wu et al. reported that miR811 and miR829 confer a high degree of resistance to Exserohilum turcicum [30]. Sheng et al. described the identification and characterization of novel miRNAs that are differentially expressed in drought-tolerant and drought-sensitive maize inbred lines [31].
At the early stage of the hybrid maize breeding in America, there were two major races, namely Southern Dent and Northern Flint [32]. After 1947, the maize races were divided into two major germplasm groups, namely Reid group (from Southern Dent) and Lancaster group (from Northern Flint) [32]. These two groups represent the main genetic diversity that is available for maize breeding in China and America [32][33][34]. In this study, miRNAs and their target transcripts in the PH6WC and PH4CV maize inbred lines, which were respectively from Reid and Lancaster groups [33][34], were investigated by using high-throughput sequencing and degradome analysis. The results indicated that some miRNAs (e.g. miR156, miR171, miR396 and miR444) differentially expressed in the seed development between PH6WC and PH4CV maize inbred lines under different genetic backgrounds.

Plant materials
Two maize inbred lines (i.e. PH6WC and PH4CV) were separately grown in experimental fields (JiangPu, Nanjing, China), without any artificial cultivation (for instance, fertilization and deinsectization). According to previous studies, the maximum value of the grain filling rate in maize seed development occurs between 21-25 days after pollination (DAP) [29,35]. Therefore, samples form three stages (15,25 and 45 DAP), which respectively represent early, medium and late stages in seed development, were selected and combined together for representing the whole process of the maize seed development in this study. In addition, all samples were frozen in liquid nitrogen and stored at -80°C. To minimize inter-individual differences, three samples from the same sampling location were mixed together.

RNA extraction
Samples from PH6WC and PH4CV lines were individually subjected to RNA extraction using the TRIzol 1 reagent (Invitrogen, Carlsbad, CA, USA). The quality of the extracted total RNAs was verified by using an Agilent 2100 Bioanalyzer (for concentration, 28S/18S and RIN detection; Agilent Technologies) and a NanoDrop 2000 spectrophotometer (for OD 260/280 and OD 260/230 detection; Thermo Fisher Scientific). In addition, the total RNAs were used for highthroughput sequencing, degradome analysis and qRT-PCR verification, identically.

High-throughput sequencing and miRNA identification
Two small RNA (sRNA) libraries were constructed using the Illumina TruSeq Small RNA Preparation Kit (LC Sciences, Hangzhou, China). The total sRNAs were ligated to 3p and 5p adapters (ADTs), and the corresponding cDNA was obtained by reverse-transcription PCR. Following purification, the cDNA from the two sRNA libraries was sequenced using an Illumina HiSeq 2000 (LC Sciences, Hangzhou, China). Removing low-quality data, the raw reads were obtained using the Illumina Pipeline v1.5 (LC Sciences, Hangzhou, China). After removing ADTs, sequences with lengths <18 and >25 nt, junk data, mRNA fragments, Rfam and Repeats, the clean reads were subjected to miRNA identification by using the selected Gramineae pre-miRNAs/miRNAs database in miRBase 21.0 and the maize genome database. Three mismatches were allowed between the reads and the known pre-miRNAs/miRNAs sequences. As results, the reads that mapped to known pre-miRNAs/miRNAs and also mapped to the maize genome were identified as conserved miRNAs. In addition, the reads that did not map to known pre-miRNAs/miRNAs but mapped to the maize genome were considered as novel miRNAs. Furthermore, the secondary structures of all identified and potential pre-miRNAs in the maize genome were predicted by using the UNAFold software [36]. The minimal folding energy indexes (MFEIs) of the novel miRNAs should be equal or greater than 0.9 [37][38][39].

Degradome sequencing and target identification
Two degradome libraries were constructed based on published methods [17,19,40]. Poly-A RNAs were obtained and ligated to a 5p adapter, and the cDNA was obtained by PCR. Following purification, the cDNA was sequenced through using an Illumina HiSeq 2000 (LC Sciences, Hangzhou, China). Removing low-quality data, the raw reads were obtained by using the Illumina Pipeline v1.5 (LC Sciences, Hangzhou, China). After removing ADTs and reads with lengths <15 nt, the remaining reads were compared with a cDNA library from the maize genome database. The mapped cDNA reads were then compared with the identified miRNAs to perform an alignment analysis by using CleaveLand 3.0 (LC Sciences, Hangzhou, China). The alignment scores 4 were considered. Furthermore, based on the number of degradome sequences and their abundance values, the miRNA targets were classified into 5 categories (0, 1, 2, 3 and 4, S5 Table) in accordance with reported method [17,19,40]. To further elucidate the potential functions, these miRNA targets were annotated through Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) [41].

qRT-PCR
Total RNA extraction was performed and subsequently used for quantitative real-time PCR (qRT-PCR). qRT-PCR was conducted by using SYBR Premix Ex Taq™ (Takara, Dalian, China) system on a BIOER Line-Gene K RT-PCR machine (BIOER, Hangzhou, China). The primers were listed in S6.1 Table, and U6 was used as internal reference [42][43]. In addition, reactions were performed in triplicate, and relative expression levels were quantified by using 2 -ΔΔCt method (S6.2 Table).

Statistical analysis
Log 2 (ratio) test and Chi-square 2×2 test were performed to identify differences in miRNA expression between the PH6WC and PH4CV libraries. Moreover, p values from Chi-square 2×2 test were adjusted to False Discovery Rate (FDR) as previous studies [44][45].

Data availability
The datum of the high-throughput sequencing and degradome sequencing were deposited in the short read archive (SRA) in National Center of Biotechnology Information (NCBI). Their numbers are SRX1686992, SRX1686966, SRX1684509 and SRX1684462.

Analysis of sRNA libraries
To identify miRNAs in seed develoment between two maize inbred lines (Fig 1a), two sRNA libraries were constructed. Following high-throughput sequencing, a total of 5,677,694 (from the PH6WC library) and 8,992,803 (from the PH4CV library) raw reads were generated (Fig 1b, S1 Table). After data filtering, 2,297,642 and 3,838,297 clean reads corresponding to 628,701 and 730,862 unique reads were obtained, respectively (Fig 1b, S1 Table). The distribution of the clean reads lengths was mainly between 21 and 24 nt (Fig 1c, S2 Table), which is consistent with previous reports [25][26][27][28][29][30][31]. Identification of conserved and novel miRNAs A total of 55 known unique miRNA sequences corresponding to 103 miRNAs were identified as conserved miRNAs, in which the miRNA sequences belonging to maize miRbase were identified as conserved known miRNAs (abbreviated as Con-K) and the sequences belonging to the selected Gramineae miRbase (except Zea mays) were identified as conserved novel miRNAs (abbreviated as Con-N) (S3 Table). Furthermore, these 55 known unique miRNA sequences were classified into 22 miRNA families (Fig 1d, S3 Table). On the other hand, a total of 24 novel unique miRNA sequences could form stem-loop structures with MFEI 0.9 were identified in two libraries (S3 Table). For expression comparison between the PH6WC and PH4CV libraries, unique miRNA sequences were analyzed through Log 2 (ratio) test and Chi-square 2×2 tests based on their normalized reads (Fig 2a, S3 Table). Following significant difference standard (p < 0.05 and |log 2 (PH4CV/PH6WC)| 1), 60 differentially expressed unique sequences were detected in two libraries (Fig 2, S3 Table). Comparing with the PH6WC library, 31 unique sequences were at up-expressed level and 29 unique sequences were at down-expressed level in the PH4CV library (Fig 2b, S3 Table).

Identification of miRNA targets
The degradome analysis was performed in order to explore the potential miRNA targets. A total of 9,591,315 and 9,061,915 raw reads were generated from the PH6WC and PH4CV degradome libraries, respectively (S4 Table). After removing the ADTs and reads <15 nt, the remaining reads were compared with the Zea mays cDNA library. A total of 6,301,841 and 6,675,216 mapped cDNA reads were obtained from the two degradome libraries (S4 Table). The mapped cDNA reads were then compared with the identified miRNAs. Finally, a total of 137 miRNA targets corresponding to 38 unique miRNA sequences were discovered in the PH6WC and PH4CV degradome libraries (S5 Table). Among these, 74 targets were detected in both two libraries (S5 Table).

Discussion
High-throughput sequencing and degradome sequencing were used to identify differentially expressed miRNAs and their targets in seed development between PH6WC and PH4CV maize inbred lines. In total, 79 conserved and novel unique miRNA sequences corresponding to 127 miRNAs were identified in two inbred lines, in which 60 unique miRNA sequences existed different expression (p < 0.05 and |log 2 (PH4CV/PH6WC)| 1). Among these, 16 miRNA sequences showed significantly difference (normalized reads 50, p < 0.01 and |log2(PH4CV/ PH6WC)| 1), and were verified by qRT-PCR. Nine verified significantly differential miRNA sequences targeted 30 transcripts. By gene annotation, we found that most these targets were connected with maize reproductive development (Table 1). For instance, the targets of zma-miR156a,b,c,d,e,f,g,l-5p were predicted as SPL10 and SPL11 from Squamosa promoter-binding-like protein family, taking part in the regulation of timing of transition from vegetative to reproductive phase [46]. According to the previous studies, miR159 and miR319 probably from a homology family potentially regulate heterochronic development [47]. Similarly, zma-miR159c_L-1R-1-3p and zma-miR319a_R+1, b_R+1,c_R+1,d_R+1-3p in our study targeted to a same gene TCP2. Moreover, Gong et al. also reported that miR159 and miR319 might influence sweet corn seed vigor [48]. APUM1 gene, might regulate mRNA translation [49], was considered as the target of zma-miR167g_R+1,h_R +1,i_R+1-5p. Additionally, SCL6 gene was predicted as the target of zma-miR171c_1ss4CT,d,e, n_R+1_1ss12CT-3p. Previous studies declared SCL family is able to influence flowering time, which can control seed development [50][51]. In this present study, the major targets of zma-miR396a,b-5p and zma-miR396c_L-1 were GRFs (3, 5 and 6), members of growth-regulating factors, involved in regulation of cell expansion in leaf and cotyledons tissues [52]. This result agreed with previous study in miRNA identification in maize grain filling stages [29]. Besides, we found that QQT2 was another target of zma-miR396a,b-5p. QQT2 gene possibly is a regulator in embryo development [53]. EDA18 gene as the target of zma-miR444a_1ss12TC, b_1ss12TC might participate in pollen and embryo sac development [54][55]. zma-miR827_R-1-3p should be conjectured to take part in maize seed development although its target was unknown in our study, because miR827 had obvious differences in expression level in our study and another previous study [29]. We believed that zma-miR166b,d,e,g,h,i-3p probably relate to the maize seed development, because 1) it high expressed in our study and other previous study [29]; 2) its target, ATHB9 gene, is involved in the determination of adaxial-abaxial polarity in ovule primordium [56].

Conclusions
In this study, several differentially expressed miRNAs and their targets were identified in seed development in two maize inbred lines (i.e. PH4CV and PH6WC) by using high-throughput sequencing and degradome analysis. The results indicated that miR156, miR171, miR396 and miR444, which respectively targeted to SPL, SCL, QQT and EDA genes, might differentially expressed in the seed development in two maize inbred lines, especially involved in flowering  regulation and embryo development. Therefore, this preliminary results might improve our understanding on the regulatory roles of miRNAs in maize seed development under different genetic backgrounds.
Supporting Information S1 Table. Overview of sRNA sequences in the PH6WC and PH4CV libraries.  Table. qRT-PCR verification and primers of miRNAs in the PH6WC and PH4CV lines. (XLSX)