Rfam (http://www.sanger.ac.uk/software/Rfam)and the and B.rapa actin gene was used as the reference gene. NCBI GenBank noncoding RNA (http://www.ncbi.nlm The primers used are showed in Table S9. .nih.gov/)database,the small RNA sequences exactly matching rRNA,tRNA,small nuclear RNA (snRNA), Measurements of hormone contents small nucleolar RNA (snoRNA)together with sequences containing poly (A)tails were removed(Gardner et al. Approximately 500 mg frozen samples of leaves before 2009).We then aligned the remainder of the unique and after vernalization were used for endogenous phyto- small RNA sequences against miRBase 21.0(http://www hormone extraction.IAA and GA3 contents were mea- .mirbase.org/index.shtml)to identify B.rapa known miR- sured by the previously described methods (Takei et al. NAs.The sequences that match miRBase with up to two 2001,Dobrev and Vankova 2012).To link the RNA mismatches were considered as known miRNAs.Novel sequencing results to hormone contents,the examined miRNAs were predicted using the program Mireap leaves were to the same ones than for sequencing.The (http://sourceforge.net/projects/mireap)with default hormone was isolated based on the previously published parameters.The secondary structures of candidate novel protocol (Agar et al.2006).High-performance liquid miRNAs were confirmed by Mfold(Zuker 2003). chromatography-mass spectrometry (AB 5500,Beijing, China)was used to detect and quantify hormone follow- ing the previously reported protocol(Pan et al.2010). Differential expression analysis of miRNAs Both standard IAA and GA3 sample were purchased To investigate significantly differentially expressed from Sigma-Aldrich (St.Louis,MO).The results were miRNAs between the two libraries(control and vernal- analyzed using six replicates. ization),the read counts of miRNAs were normalized as transcripts per million according to Bayesian meth- Results ods (Audic and Claverie 1997).If the normalized read count is zero,the value was modified to 0.001.The Global analysis of sequences from small RNA libraries fold-change and P-value were calculated by normalized expression according to previously reported methods To analysis the roles of miRNAs in the vernalization (Li et al.2011).When P<0.05 and fold change (log2 process,we constructed two sRNA libraries from leaves (vernalization/control)>2 or <0.5,it was considered up- of control and vernalized plants.30 677 485 clean reads or down-regulated during the process of vernalization. in control and 20504 643 clean reads in vernalization library were obtained from Solexa sequencing (Table S1). Prediction and annotation of potential targets Among them,17%total sRNAs were non-vernalization library-specific with 50.8%unique sRNAs,10%total We used miRNA candidate targets as query sequences sRNAs were vernalization library-specific with 33% to blast against Chiifu Chinese cabbage database(http:// unique sRNAs and 73%total sRNAs were present in brassicadb.org/brad/)by a custom Perl script,based on both with 16.2%unique sRNAs (Fig.1;Table S2).Then, the published rules (Allen et al.2005,Schwab et al. we annotated and filtered out undesired sRNAs,such as 2005).A BlastX search with the Arabidopsis genome was rRNAs,tRNAs and mRNA.The proportions and numbers then applied to the predicted targets,considering they for different categories of small RNA were summarized both belong to Brassicaceae.Moreover,to further study in Table S3. the biological functions of miRNA targets,we obtained The size distribution of reads displayed similar trends GO annotations through Blast2GO program from Uni- between the two libraries (Fig.2).The majority of the Gene database based on the BLAST searches against the sRNAs were 24 nt long,which accounted for 47.2 and Nr database in NCBI under an E-value threshold of less 47.1%of total reads in control and vernalization library, than 10-5(Conesa et al.2005) respectively,followed by the 21 nt(23.5 and 19.2%). This was consistent with previous studies on A.thaliana QPCR analysis (Rajagopalan et al.2006),Oryza sativa(Zhu et al.2008), Medicago truncatula(Szittya et al.2008)and Populus We performed qPCR to confirm the quality of small trichocarpa(Puzey et al.2012).However,a difference in RNA sequencing and expression patterns of miRNAs the size distribution was noticed.The relative abundance and their targets.Total RNA was extracted from the of 21 nt sRNAs in vernalization library was significantly previously frozen plant tissues.The gPCR assay and lower than those in control library,indicating that the relative expression calculation were carried out followed 21 nt sRNAs might be repressed in B.rapa vernalization our previous report(Huang et al.2016).The 5S rRNA process. Physiol.Plant.2018Rfam (http://www.sanger.ac.uk/software/Rfam) and the NCBI GenBank noncoding RNA (http://www.ncbi.nlm .nih.gov/) database, the small RNA sequences exactly matching rRNA, tRNA, small nuclear RNA (snRNA), small nucleolar RNA (snoRNA) together with sequences containing poly (A) tails were removed (Gardner et al. 2009). We then aligned the remainder of the unique small RNA sequences against miRBase 21.0 (http://www .mirbase.org/index.shtml) to identify B. rapa known miR￾NAs. The sequences that match miRBase with up to two mismatches were considered as known miRNAs. Novel miRNAs were predicted using the program Mireap (http://sourceforge.net/projects/mireap) with default parameters. The secondary structures of candidate novel miRNAs were confirmed by Mfold (Zuker 2003). Differential expression analysis of miRNAs To investigate significantly differentially expressed miRNAs between the two libraries (control and vernal￾ization), the read counts of miRNAs were normalized as transcripts per million according to Bayesian meth￾ods (Audic and Claverie 1997). If the normalized read count is zero, the value was modified to 0.001. The fold-change and P-value were calculated by normalized expression according to previously reported methods (Li et al. 2011). When P ≤0.05 and fold change (log2 (vernalization/control) >2 or <0.5, it was considered up￾or down-regulated during the process of vernalization. Prediction and annotation of potential targets We used miRNA candidate targets as query sequences to blast against Chiifu Chinese cabbage database (http:// brassicadb.org/brad/) by a custom Perl script, based on the published rules (Allen et al. 2005, Schwab et al. 2005). A BlastX search with the Arabidopsis genome was then applied to the predicted targets, considering they both belong to Brassicaceae. Moreover, to further study the biological functions of miRNA targets, we obtained GO annotations through Blast2GO program from Uni￾Gene database based on the BLAST searches against the Nr database in NCBI under an E-value threshold of less than 10−5 (Conesa et al. 2005). QPCR analysis We performed qPCR to confirm the quality of small RNA sequencing and expression patterns of miRNAs and their targets. Total RNA was extracted from the previously frozen plant tissues. The qPCR assay and relative expression calculation were carried out followed our previous report (Huang et al. 2016). The 5S rRNA and B. rapa actin gene was used as the reference gene. The primers used are showed in Table S9. Measurements of hormone contents Approximately 500 mg frozen samples of leaves before and after vernalization were used for endogenous phyto￾hormone extraction. IAA and GA3 contents were mea￾sured by the previously described methods (Takei et al. 2001, Dobrev and Vankova 2012). To link the RNA sequencing results to hormone contents, the examined leaves were to the same ones than for sequencing. The hormone was isolated based on the previously published protocol (Agar et al. 2006). High-performance liquid chromatography-mass spectrometry (AB 5500, Beijing, China) was used to detect and quantify hormone follow￾ing the previously reported protocol (Pan et al. 2010). Both standard IAA and GA3 sample were purchased from Sigma-Aldrich (St. Louis, MO). The results were analyzed using six replicates. Results Global analysis of sequences from small RNA libraries To analysis the roles of miRNAs in the vernalization process, we constructed two sRNA libraries from leaves of control and vernalized plants. 30 677 485 clean reads in control and 20 504 643 clean reads in vernalization library were obtained from Solexa sequencing (Table S1). Among them, 17% total sRNAs were non-vernalization library-specific with 50.8% unique sRNAs, 10% total sRNAs were vernalization library-specific with 33% unique sRNAs and 73% total sRNAs were present in both with 16.2% unique sRNAs (Fig. 1; Table S2). Then, we annotated and filtered out undesired sRNAs, such as rRNAs, tRNAs and mRNA. The proportions and numbers for different categories of small RNA were summarized in Table S3. The size distribution of reads displayed similar trends between the two libraries (Fig. 2). The majority of the sRNAs were 24 nt long, which accounted for 47.2 and 47.1% of total reads in control and vernalization library, respectively, followed by the 21 nt (23.5 and 19.2%). This was consistent with previous studies on A. thaliana (Rajagopalan et al. 2006), Oryza sativa (Zhu et al. 2008), Medicago truncatula (Szittya et al. 2008) and Populus trichocarpa (Puzey et al. 2012). However, a difference in the size distribution was noticed. The relative abundance of 21 nt sRNAs in vernalization library was significantly lower than those in control library, indicating that the 21 nt sRNAs might be repressed in B. rapa vernalization process. Physiol. Plant. 2018