Gynophore miRNA analysis at different developmental stages in Arachis duranensis

Plant materials

Peanuts (A. duranensis) were grown in the field at Jiangsu Academy of Agricultural Sciences. Three different gynophore stages (A1: 5 days, A2: 10 days, and A3: 20 days) were collected and immediately stored in liquid nitrogen before RNA extraction. This study was approved by the Institutional Plant Care and Use Committee of Jiangsu Academy of Agricultural Sciences.

RNA extraction and high-throughput sequencing

Total RNA was obtained from the different gynophore stages (A1, A2, and A3) using Trizol agent (TaKaRa, Dalian, China), following the manufacturer instructions. Total RNA of peanut gynophores was isolated by an improved CTAB method with isopropanol instead of lithium chloride for RNA precipitation (Mi et al., 2008). Briefly, 1 g gynophore tissue was added to liquid nitrogen to obtain a fine powder. Subsequently, the powder was evenly mixed using 5 mL preheated (65°C) extraction buffer (2% CTAB, 2% PVP, 0.1 M Tris-HCl, 2.0 M NaCl, 25 mM EDTA, 2% beta-mercaptoethanol, pH 8.0). The mixture was incubated for 5 min at 65°C and shaken three times during the incubation time. After a short-time cooling, 2.5 mL isopropanol was added to the cooling mixture. After vortexing for 1 min, the mixture was centrifuged at 13,800 g for 15 min at 4°C. DNase was used to treat the extraction, and RNA was precipitated at 25°C for 10 min using 2.5 mL isopropanol. The extracted RNA was first resuspended with an equal volume of mixture (phenol:chloroform:isopropanol = 25:24:1), followed by resuspension using an equal mixture volume (chloroform:isopropanol = 24:1). Both 0.25mL 3 M NaOAC, pH 5.2, and 7.5 mL cooling ethanol were added to the mixture to precipitate the RNA overnight at -20°C. The quality of the prepared RNA from each sample was analyzed using Agilent 2100 (Agilent Technologies Co., Ltd., Palo Alto, CA, USA). The added fragmentation buffer (Ambion) was used to cut mRNA into short fragments (200-700 nt). First-strand cDNA was synthesized based on the short-fragment templates. DNA polymerase I (New England Biolabs, Beijing, China), dNTPs, RNase H (Invitrogen), and buffer were used to synthesize the second-strand cDNA. The fragments were purified using QiaQuick PCR kit and washed in EB buffer for end repair before adding polyA tails and adaptors. Fragments of suitable size were collected by agarose gel electrophoresis and amplified using polymerase chain reaction (PCR). The resulting products were sequenced by Illumina HiSeq™ 2000.

Prediction of novel miRNA

In order to predict novel gynophore miRNAs, the small RNA sequences were aligned to ESTs of peanut to collect precursor sequences. The structural features of the miRNA precursors were analyzed using the Mireap program constructed by the Beijing Genome Institute to find novel miRNA candidates. The collected structures were considered to be novel miRNA candidates when they conformed to the standard described by Allen et al. (2005). Mfold (Zuker, 2003) was used to examine the secondary structures of the collected pre-miRNA sequences. Furthermore, structures with the lowest energy were filtered for manual inspection.

Target gene prediction

The psRNATarget program (http://bioinfo3.noble.org/psRNATarget) was used to analyze the potential targets of the peanut miRNAs. The novel peanut miRNA fragments were used as custom miRNA sequences. Arachis transcript/genomic libraries (GSS, EST, and nucleotide databases) were used as custom plant databases. The analyzed genes were scored, and the standard roles used in this study were as follows: every G:U wobble pairing received 0.5 points, 2.0 points were assigned for every indel, and all non-canonical Watson-Crick pairings received 1.0 point each. The complete score of an alignment was calculated based on 20 nt. Scores of all possible consecutive 20 nt subsequences were calculated when the fragment was longer than 20 nt, and the minimum score was collected as the total score of the query-subject alignment. Given that the targets’ complementary sequence at the miRNA 5'-end is critical, mismatches other than G:U wobbles at positions 2-7 at the 5'-end received a 0.5 point deduction of the final score (Griffiths-Jones et al., 2005). The sequences whose final scores were less than 3.0 points were considered to be the miRNA targets. Redundant sequences were deleted using the Vector NTI program, when the potential target mRNA sequences were collected. The target sequences were subjected to BLASTx analysis and the functions of the potential targets were analyzed in the NCBI database.

Analysis of GO pathway

BLASTx analysis was performed using the target sequences and the NCBI database to better understand the function of the miRNA targets and the metabolic regulatory networks related to the identified peanut miRNAs. A BLASTx of the Interpro and KEGG databases was used to collect the predicted target proteins with an E value of 1e-30. The target gene function and metabolic pathway of the miRNAs were verified using the best hits. Finally, we obtained molecular function, biological process, and cellular component of the target genes from the GO and Interpro database.

Data analysis

Clean reads were filtered from the raw reads obtained after sequencing, by removing sequences containing adapters. The collected reads were used for de novo assembly, which was performed using the SOAPdenovo program. Firstly, the filtered reads with suitable overlap length were combined to form contigs. Secondly, the filtered reads were re-mapped to contigs, and then scaffolds were made. Lastly, pair-end reads were used again for gap filling to construct unigenes. To show the alternation of each DEG during the whole development stages, their fold-change data were imported into MultiExperiment Viewer (MeV v4.8, http://www.tm4.org/mev/) for hierarchical clustering analysis with the average linkage method. The functional annotation, GO functional analysis, and KEGG pathway analysis of these unigenes were carried out using BLASTx (E value < 0.00001) against the Nr, Swiss-Prot, KEGG, and COG databases.

Sequence analysis

In the present study, high-throughput sequencing was used to identify low-abundance candidate miRNAs in gynophores at different developmental stages. After filtering out low-quality tags, excluding adaptor/acceptor sequences, and removing contamination caused by adaptor-adaptor ligation, multiple-clean reads were obtained. All sequences were found to be small RNAs, including miRNA, noncoding RNA, tRNA, rRNA, snRNA, and snoRNA. Of the small RNAs, 70.62% were common to stages A1 and A2, whereas 16.28 and 13.10% were unique to stages A1 and A2, respectively. Furthermore, 60.45% of the small RNAs were common to stages A1 and A3, whereas 15.72 and 23.83% were unique to stages A1 and A3, respectively. Finally, 64.17% of the small RNAs were common to stages A2 and A3, 12.33% were unique to stage A2 and 23.50% were unique to stage A3 (Table 1). The numbers and kinds of small RNAs identified in the gynophore at different stages are shown in Table 1. A large number of known miRNAs were detected in all three developmental stages of the peanut gynophore. More miRNAs were found in stage A2 than in stages A1 or A3. Interestingly, stage A1 contained the most different kinds of miRNAs, whereas stage A3 had the fewest kinds of miRNAs. The length distribution of small RNAs at different stages often determines their roles in the particular stage and biogenetic machines. As shown in Figure 1, the majority of the small RNAs in the stage A1, A2, and A3 libraries were 21 nt long and accounted for 29.9, 34.0, and 32.5% of the total sequence numbers, respectively. In addition, the number of small RNAs (24 nt) from stages A1, A2, and A3 were 26.4, 22.7, and 30.4%, respectively, of the total number of RNAs. The 21 and 24 nt small RNAs were the most prevalent in the gynophore, which was consistent with the results reported for other peanut tissues (Chi et al., 2011). There were more 21 nt small RNAs than 24 nt RNAs in gynophore. Furthermore, as found in other studies (Chi et al., 2011), the high percentage of 24 nt small RNAs in the peanut gynophore could be indicative of the complexity of the peanut genome.

Figure 1

Tag length distribution of the three different gynophore stages.

Known and novel miRNAs

There were 266 known miRNAs detected in the three different stages of the gynophore (Table S1). As shown in previous studies, most miRNAs from the gynophore were highly conserved in different plant species (Sunkar et al., 2008), which indicates that the ancient regulatory pathways of conserved miRNAs are present in peanut too. miRNA members from different stages of the gynophore were compared with each other and great differences were found among them. The miRNAs miR166d-5p, miR7785-3p, miR5713, miR1080, miR8578.1, and miR482c-5p were present in more than 1000 reads in stage A1 but were not found at all in stages A2 and A3. Likewise, three miRNAs (miR166a, miR482e, and miR9722) were identified in more than 1000 reads from stage A2, whereas they were never observed in stages A1 and A3. The known miRNA, miR166u, was sequenced in more than 300,000 reads in stage A3, which indicates that this miRNA might play a significant and essential role in forming the peanut gynophore. Furthermore, a cluster analysis was performed to visualize the read differences between the three stages of the gynophore. As shown in Figure 2A, many known miRNAs showed great differences between reads of any two stages, suggesting that some play specific roles in different gynophore stages. A cluster analysis was used to analyze the differences of reads among the three gynophore stages (Figure 2B). Many novel miRNAs showed great read differences between any two stages, which suggests that there were some specific roles of the miRNAs in the different gynophore stages. In order to compare the miRNA expression from the three gynophore stages, the number of reads mapped to a gene was analyzed. As shown in Figure 3, the stage A1 score was the highest and stage A3 was the lowest. The number of reads mapped to a gene from stage A1 was higher than the other two stages suggesting that the miRNA function in stage A1 may be more complex compared to the other two stages. These results indicate that different miRNAs had different expression patterns, which may be due to stage-specific expression.Figure 2

Cluster analysis of miRNAs from three different gynophore stages (A1, A2, and A3) of known (A) and novel (B) miRNAs. Phylogeny scale bar shows the different levels of clusterings of miRNAs by their expression patterns in three stages.

Figure 3

Number of reads to one gene from different gynophore stages (A1, A2, and A3). The contrast between stages A1 and A2 of known (A) and novel (D) miRNAs, between stages A1 and A3 (B: known, E: novel miRNAs), and between stages A2 and A3 (C: known, F: novel miRNAs).

Prediction of function

In order to fully investigate the biological functions of the known and novel miRNAs, GO analysis was used to investigate the gene ontology of all targets regulated by the miRNAs (Du et al., 2010). The terms from different stages are listed in Table S2. The terms were divided into three categories; biological process, cellular component, and molecular function. For known miRNAs, we found that many genes were involved in the lignin catabolic process at stage A1, while no such genes were detected in stages A2 and A3. This suggests that the lignin catabolic process is important to insert the gynophore into the ground, since lignin hardens the gynophore enough to pierce the soil surface. The cellular component demonstrated that 46 genes were involved in apoplast formation at stage A1 but not in stages A2 and A3. This indicates that the apoplast development process is necessary for stage A1 to develop into stage A2. A large number of genes were involved in different binding processes in all three developmental stages. Many genes participated in the copper ion-binding process of stage A1, while none of these genes was found in stages A2 and A3. This indicates that the copper ion-binding process might play an important role in preparation before developing the gynophore in stage A3. The genes involved with protein kinase activity were fully functioning in stage A2, suggesting that protein kinases at this stage were more active than at the other two stages. As listed in Table S2, novel miRNA terms from stages A1, A2, and A3 were collected to analyze the biological functions. Several terms relating to energy consumption (ATP metabolic process, proton transport, and ATP hydrolysis-coupled proton transport) were found only at stage A3, which suggests that stage A3 may be a rapid development phase of the gynophore. In addition, genes that negatively regulate biological process were specific to stages A2 and A3, whereas no such genes were filtered in stage A1 (Figure 4). This suggests that some biological processes of the reproductive process at stage A1 may be suppressed in stages A2 and A3. There were fewer genes related to nutrient reservoir activity in stage A2 than in stages A1 and A3, indicating that stage A2 might be the consuming phase of the gynophore. As shown in Figure 5, the biological functions of novel miRNAs were similar to known miRNAs. The genes involved in growth processes were found only in stage A3, indicating that the rapid development of the gynophore occurs during stage A3.Figure 4

Function of known miRNAs from different gynophore stages (A1, A2, and A3).

Figure 5

Function of novel miRNAs from different gynophore stages (A1, A2, and A3).