《园艺作物育种学》课程教学资源（学术研究）Vernalization can regulate flowering time through microRNA mechanism in Brassica rapa

Physiologia Plantarum An International Journal for Plant Biology Physiologia Plantarum 2018 2018 Scandinavian Plant Physiology Society,ISSN 0031-9317 Vernalization can regulate flowering time through microRNA mechanism in Brassica rapa Feiyi Huang',Xiaoting Wut,Xilin Hou,Shuaixu Shao and Tongkun Liu* State Key Laboratory of Crop Genetics and Germplasm Enhancement,College of Horticulture,Nanjing Agricultural University,Nanjing 210095,China Correspondence Vernalization is an important process that regulates the floral transition *Corresponding author, in plants.MicroRNAs (miRNAs)are endogenous non-coding small RNA e-mail:liutk@njau.edu.cn (sRNA)molecules that function in plant growth and development.Despite Received 21 September 2017; that miRNAs related to flowering have previously been characterized,their revised 8 January 2018 roles in response to vernalization in pak-choi (Brassica rapa ssp.chinensis) has never been studied.Here,two sRNA libraries from B.rapa leaves doi:10.1111/ppl.12692 (vernalized and non-vernalized plants)were constructed and sequenced Two hundred eight known and 535 novel miRNAs were obtained,of which 20 known and 66 new miRNAs were significantly differentially expressed and considered as vernalization-related miRNAs.The corresponding targets were predicted on the basic of sequence homology search.In addition,11 miRNAs and eight targets were selected for real-time quantitative PCR to confirm their expression profiles.Functional annotation of targets using gene ontology and Kyoto encyclopedia of genes and genomes results suggested that most targets were significantly enriched in the hormone signaling pathway. Moreover,a decreased indole-3-acetic acid (IAA)and an increased GA3 hormone were detected after vernalization,indicating that the IAA and GA3 might response to vernalization.These results indicated that vernalization regulates flowering through microRNA mechanism by affecting endogenous hormone level in B.rapa.This study provides useful insights of promising miRNAs candidates involved in vernalization in B.rapa,and facilitates further investigation of the miRNA-mediated molecular mechanisms of vernalization in B.rapa. Introduction which can be stimulated by various environmental stimuli and endogenous signals.This process is con- Plant development consists of two phases,vegetative trolled by numerous genes,as well as some evolutionary and reproductive growth.Flowering is an important conserved microRNAs(miRNAs)(Irish 2010,Pose et al. biological process that ensures reproductive success, 2012).It has been determined that there are five major Abbreviation -AFB,AUXIN SIGNALING F-BOX;AP2,APETELA2;COLDAIR,cold assisted intronic noncoding RNA;COOLAIR, cold induced long antisense intragenic RNA;FT,FLOWERING LOCUS T;FLC,FLOWERING LOCUS C;GA,gibberellin;GO,gene ontology;IAA,indole-3-acetic acid;KEGG,Kyoto encyclopedia of genes and genomes;LFY,LEAFY;IncRNA,long non-coding RNA;MYB,MYB DOMAIN PROTEIN;miRNA,microRNA;nt,nucleotides;pre-miRNA,miRNA precursor;pri-miRNA,primary miRNA;gPCR,real-time quantitative PCR;RISC,the RNA-induced silencing complex;SAM,shoot apical meristem;SOC1 SUPPRESSOR OF OVEREXPRESSION OF CO 1;SPL,SQUAMOSA PROMOTER BINDING PROTEIN-LIKE;snRNA,small nuclear RNA;snoRNA,small nucleolar RNA;sRNA,small RNA:TIR1,Tolllinterleukin receptor 1. tThese authors equally contributed to this work. Physiol.Plant.2018

Physiologia Plantarum 2018 © 2018 Scandinavian Plant Physiology Society, ISSN 0031-9317 Vernalization can regulate flowering time through microRNA mechanism in Brassica rapa Feiyi Huang†, Xiaoting Wu†, Xilin Hou , Shuaixu Shao and Tongkun Liu* State Key Laboratory of Crop Genetics and Germplasm Enhancement, College of Horticulture, Nanjing Agricultural University, Nanjing 210095, China Correspondence *Corresponding author, e-mail: liutk@njau.edu.cn Received 21 September 2017; revised 8 January 2018 doi:10.1111/ppl.12692 Vernalization is an important process that regulates the floral transition in plants. MicroRNAs (miRNAs) are endogenous non-coding small RNA (sRNA) molecules that function in plant growth and development. Despite that miRNAs related to flowering have previously been characterized, their roles in response to vernalization in pak-choi (Brassica rapa ssp. chinensis) has never been studied. Here, two sRNA libraries from B. rapa leaves (vernalized and non-vernalized plants) were constructed and sequenced. Two hundred eight known and 535 novel miRNAs were obtained, of which 20 known and 66 new miRNAs were significantly differentially expressed and considered as vernalization-related miRNAs. The corresponding targets were predicted on the basic of sequence homology search. In addition, 11 miRNAs and eight targets were selected for real-time quantitative PCR to confirm their expression profiles. Functional annotation of targets using gene ontology and Kyoto encyclopedia of genes and genomes results suggested that most targets were significantly enriched in the hormone signaling pathway. Moreover, a decreased indole-3-acetic acid (IAA) and an increased GA3 hormone were detected after vernalization, indicating that the IAA and GA3 might response to vernalization. These results indicated that vernalization regulates flowering through microRNA mechanism by affecting endogenous hormone level in B. rapa. This study provides useful insights of promising miRNAs candidates involved in vernalization in B. rapa, and facilitates further investigation of the miRNA-mediated molecular mechanisms of vernalization in B. rapa. Introduction Plant development consists of two phases, vegetative and reproductive growth. Flowering is an important biological process that ensures reproductive success, Abbreviation – AFB, AUXIN SIGNALING F-BOX; AP2, APETELA2; COLDAIR, cold assisted intronic noncoding RNA; COOLAIR, cold induced long antisense intragenic RNA; FT, FLOWERING LOCUS T; FLC, FLOWERING LOCUS C; GA, gibberellin; GO, gene ontology; IAA, indole-3-acetic acid; KEGG, Kyoto encyclopedia of genes and genomes; LFY, LEAFY; lncRNA, long non-coding RNA; MYB, MYB DOMAIN PROTEIN; miRNA, microRNA; nt, nucleotides; pre-miRNA, miRNA precursor; pri-miRNA, primary miRNA; qPCR, real-time quantitative PCR; RISC, the RNA-induced silencing complex; SAM, shoot apical meristem; SOC1, SUPPRESSOR OF OVEREXPRESSION OF CO 1; SPL, SQUAMOSA PROMOTER BINDING PROTEIN-LIKE; snRNA, small nuclear RNA; snoRNA, small nucleolar RNA; sRNA, small RNA; TIR1, Toll/interleukin receptor 1. †These authors equally contributed to this work. which can be stimulated by various environmental stimuli and endogenous signals. This process is con￾trolled by numerous genes, as well as some evolutionary conserved microRNAs (miRNAs) (Irish 2010, Pose et al. 2012). It has been determined that there are five major Physiol. Plant. 2018

pathways regulating the floral transition in Arabidop- Pak-choi (Brassica rapa ssp.chinensis),a subspecies of sis thaliana,namely the photoperiod,vernalization B.rapa,is a major vegetable crop and widely cultivated autonomous,gibberellins (GAs)and ageing pathways in Asia (Tian et al.2004).The cultivar Wuyueman',a (Srikanth and Schmid 2011).The vernalization pathway cold-requiring B.rapa which flowers later than other B. mainly mediates the response to the appropriate environ- rapa cultivars,was used in this study.Despite that numer- mental conditions.In temperate climates,flowering is ous miRNAs have been identified in B.rapa,no sys- usually coordinated with temperature and photoperiod. tematic studies were available for the roles of miRNAs In Arabidopsis,the models of flowering-time regulation in the vernalization pathway.Herein,we constructed suggest that vernalization and photoperiod influence two sRNA libraries(vernalization and non-vernalization) different signal genes.A series of'integrator'genes,such and performed the high-throughput Solexa sequencing as FT,SOC1,causes interactions between vernalization to analysis miRNAs and targets in B.rapa.The expres- and photoperiod response signals(Samach et al.2000). sion profiles of vernalization-related miRNAs and their Plant hormones and internal regulatory factors also func- targets were validated and analyzed by real-time quan- tion in the floral transition process(Li et al.2015).Recent titative PCR(qPCR).Then,we further analyzed B.rapa studies reported that GAs function to induce floral tran- miRNAs targets functions on the basic of the information sition and regulate floral development(Mutasa-Gottgens provided by a Blastx(Beijing,China)search against the and Hedden 2009).In addition to GAs,indole-3-acetic Arabidopsis genome and gene ontology(GO)and Kyoto acid (IAA)also play a vital role in the reproductive encyclopedia of genes and genomes(KEGG)analysis. phase transition.GAs signaling triggers early flowering, This work offers information to understand the regula- while IAA signaling leads to late-flowering slightly in tory roles and molecular mechanisms of miRNAs in B. Agapanthus praecox ssp.orientalis(Zhang et al.2014). rapa vernalization process. Increasing studies have been reporting that miRNAs play vital roles in controlling various plant develop- Materials and methods mental processes,included flower development (Zhao et al.2007,Yant et al.2010).MiRNAs are endogenous Plant materials 18-25 nt single-stranded non-coding small RNAs Young emerging leaves of B.rapa were used for (sRNAs),that functioned as post-transcriptional and tran- sequencing.Seeds were sowed in pots in a green- scriptional negative regulators(Voinnet 2009).They are house (22/18C,16/8 h light/dark).For vernalization first transcribed to long primary miRNAs (pri-miRNAs) treatment,1-month-old seedlings were transferred to a by RNA polymerase ll in the nucleus that fold to a growth chamber exposed to 4C for 6 weeks and then secondary structure(stem-loop)with partly complemen- transferred to the greenhouse(22/18C,16/8 h light/dark tarity.The pri-miRNAs are processed into pre-miRNAs daily)for 1 week.Seedlings grown in the greenhouse via Dicer-like one enzyme and further cleaved to small without vernalization treatment were used as control double-stranded RNA duplexes.The duplexes are then Leaves of 10-week-old vernalized and non-vernalized transported to the cytoplasm,and cleaved by Dicer, forming a miRNA(Bartel 2004,Kurihara and Watanabe (control)seedlings grown in the greenhouse were har- vested and frozen immediately in liquid nitrogen and 2004).The miRNA strand is loaded into RNA-induced stored at-80°C. silencing complex(RISC,Baumberger and Baulcombe 2005).The miRNA acts as a guide for RISC to degrade targets or inhibit their translation (Vazquez et al.2010). Construction of two small RNA library To date,nine conserved miRNA families have a ver and bioinformatics analysis of sequencing data ified function in flowering,such as miR156,miR159 Total RNA isolation and construction of two sRNA and miR172 (Luo et al.2013).MiR156 is required for libraries (vernalized and control seedlings)were com- flower transition by targeting SQUAMOSA PROMOTER pleted following reported procedures (Sunkar et al. BINDING PROTEIN-LIKE (SPLs)in Arabidopsis,maize 2008).Solexa technology was used to sequence sRNAs and rice (Chuck et al.2007,Wang et al.2009,Jiao et al. at BGI(Shenzhen,China).Small RNA reads were gen- 2010).MiR172 targets AP2 to regulate organ identity and erated by an Illumina Genome Analyzer at BGl.After flowering time (Varkonyi-Gasic et al.2012).MiR159 is removing the undesired raw reads,the extracted clean essential for normal anther development through con- reads of 18-30 nt were obtained for subsequent analy- trolling MYBs expression (Tsuji et al.2006).Gener- sis.The retained miRNA reads were mapped to B.rapa ally,these reports revealed that miRNA-mediated gene genome using SOAP software under default parameters regulation probably play crucial roles in plant flower (Li et al.2008).Only perfectly matched sequences were development. kept.The mappable sequences were search against Physiol.Plant.2018

pathways regulating the floral transition in Arabidop￾sis thaliana, namely the photoperiod, vernalization, autonomous, gibberellins (GAs) and ageing pathways (Srikanth and Schmid 2011). The vernalization pathway mainly mediates the response to the appropriate environ￾mental conditions. In temperate climates, flowering is usually coordinated with temperature and photoperiod. In Arabidopsis, the models of flowering-time regulation suggest that vernalization and photoperiod influence different signal genes. A series of ‘integrator’ genes, such as FT, SOC1, causes interactions between vernalization and photoperiod response signals (Samach et al. 2000). Plant hormones and internal regulatory factors also func￾tion in the floral transition process (Li et al. 2015). Recent studies reported that GAs function to induce floral tran￾sition and regulate floral development (Mutasa-Gottgens and Hedden 2009). In addition to GAs, indole-3-acetic acid (IAA) also play a vital role in the reproductive phase transition. GAs signaling triggers early flowering, while IAA signaling leads to late-flowering slightly in Agapanthus praecox ssp. orientalis (Zhang et al. 2014). Increasing studies have been reporting that miRNAs play vital roles in controlling various plant develop￾mental processes, included flower development (Zhao et al. 2007, Yant et al. 2010). MiRNAs are endogenous 18–25 nt single-stranded non-coding small RNAs (sRNAs), that functioned as post-transcriptional and tran￾scriptional negative regulators (Voinnet 2009). They are first transcribed to long primary miRNAs (pri-miRNAs) by RNA polymerase II in the nucleus that fold to a secondary structure (stem-loop) with partly complemen￾tarity. The pri-miRNAs are processed into pre-miRNAs via Dicer-like one enzyme and further cleaved to small double-stranded RNA duplexes. The duplexes are then transported to the cytoplasm, and cleaved by Dicer, forming a miRNA (Bartel 2004, Kurihara and Watanabe 2004). The miRNA strand is loaded into RNA-induced silencing complex (RISC, Baumberger and Baulcombe 2005). The miRNA acts as a guide for RISC to degrade targets or inhibit their translation (Vazquez et al. 2010). To date, nine conserved miRNA families have a ver￾ified function in flowering, such as miR156, miR159 and miR172 (Luo et al. 2013). MiR156 is required for flower transition by targeting SQUAMOSA PROMOTER BINDING PROTEIN-LIKE (SPLs) in Arabidopsis, maize and rice (Chuck et al. 2007, Wang et al. 2009, Jiao et al. 2010). MiR172 targets AP2 to regulate organ identity and flowering time (Varkonyi-Gasic et al. 2012). MiR159 is essential for normal anther development through con￾trolling MYBs expression (Tsuji et al. 2006). Gener￾ally, these reports revealed that miRNA-mediated gene regulation probably play crucial roles in plant flower development. Pak-choi (Brassica rapa ssp. chinensis), a subspecies of B. rapa, is a major vegetable crop and widely cultivated in Asia (Tian et al. 2004). The cultivar ‘Wuyueman’, a cold-requiring B. rapa which flowers later than other B. rapa cultivars, was used in this study. Despite that numer￾ous miRNAs have been identified in B. rapa, no sys￾tematic studies were available for the roles of miRNAs in the vernalization pathway. Herein, we constructed two sRNA libraries (vernalization and non-vernalization) and performed the high-throughput Solexa sequencing to analysis miRNAs and targets in B. rapa. The expres￾sion profiles of vernalization-related miRNAs and their targets were validated and analyzed by real-time quan￾titative PCR (qPCR). Then, we further analyzed B. rapa miRNAs targets functions on the basic of the information provided by a Blastx (Beijing, China) search against the Arabidopsis genome and gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) analysis. This work offers information to understand the regula￾tory roles and molecular mechanisms of miRNAs in B. rapa vernalization process. Materials and methods Plant materials Young emerging leaves of B. rapa were used for sequencing. Seeds were sowed in pots in a green￾house (22/18∘C, 16/8 h light/dark). For vernalization treatment, 1-month-old seedlings were transferred to a growth chamber exposed to 4∘C for 6 weeks and then transferred to the greenhouse (22/18∘C, 16/8 h light/dark daily) for 1 week. Seedlings grown in the greenhouse without vernalization treatment were used as control. Leaves of 10-week-old vernalized and non-vernalized (control) seedlings grown in the greenhouse were har￾vested and frozen immediately in liquid nitrogen and stored at −80∘C. Construction of two small RNA library and bioinformatics analysis of sequencing data Total RNA isolation and construction of two sRNA libraries (vernalized and control seedlings) were com￾pleted following reported procedures (Sunkar et al. 2008). Solexa technology was used to sequence sRNAs at BGI (Shenzhen, China). Small RNA reads were gen￾erated by an Illumina Genome Analyzer at BGI. After removing the undesired raw reads, the extracted clean reads of 18–30 nt were obtained for subsequent analy￾sis. The retained miRNA reads were mapped to B. rapa genome using SOAP software under default parameters (Li et al. 2008). Only perfectly matched sequences were kept. The mappable sequences were search against Physiol. Plant. 2018

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 P2 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.2018

Rfam (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

■Control specific in B.rapa were also identified and named as bra-miRn. ☐Common Most of the novel miRNAs were control-specific or vernalization-specific,which accounted for 47.1 (252) Vernalization specific and 28.4%(152),respectively.Only 131 novel miRNAs appeared in both control and vernalization libraries. Information about the novel miRNAs was listed in Table S5.We also found that conserved miRNAs level was higher than non-conserved and novel miRNAs level (Tables S4 and S5),which is in agreement with the pre- vious findings that evolutionary conservation may have a correlation with the expression level(Chi et al.2011). 50.89 16.2% 33.0% Identification of vernalization-related miRNAs in B.rapa After annotating known and novel miRNAs,we exam- ined vernalization-responsive miRNAs based on the differential expression analysis.We used a threshold Control Vernalization and absolute value of logzRatio >1 and P-value 10.These results indicated that these sRNAs.The two libraries are shown by different colors. identified vernalization-responsive miRNAs might play vital roles in B.rapa vernalization ldentification of known and novel miRNAs By sequencing,we obtained 43389 (control)and Prediction and annotation of vernalization-related 35 545 (vernalization)miRNA candidates.To isolate miRNA targets known miRNAs in B.rapa,the sequences of sRNAs To investigate the biological roles of the predicted targets were blasted against the known miRNAs of various plant in B.rapa,we searched against EST and cDNA sequences species deposited in miRBase 21.0.After BLASTn and in B.rapa genome annotation database(http://brassicadb sequence analysis,208 expressed known miRNAs were .org/brad/index.php)and Arabidopsis genome database identified in both libraries.Detailed sequences and (https://www.arabidopsis.org/index.jsp).A total of 714 reads of the known miRNAs are shown in Table S4.Five and 666 potential targets for 92 known and 213 new hundred thirty-five potential novel miRNAs candidates miRNAs were identified,with an average of 7.76 and Physiol.Plant.2018

Control Vernalization 50.8% 16.2% 33.0% Control specific Common Vernalization specific Fig. 1. Summary of the common and specific sequences of the unique sRNAs from control and vernalization library. Frequence (%) 50 45 40 35 30 25 20 15 10 5 0 Control Vernalization 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 Length (nt) Fig. 2. Length distribution and frequency percent of sRNA sequences in control and vernalization libraries in B. rapa. Y-axis represents percent￾ages of sRNAs identified in this study; x-axis represents the length of sRNAs. The two libraries are shown by different colors. Identification of known and novel miRNAs By sequencing, we obtained 43 389 (control) and 35 545 (vernalization) miRNA candidates. To isolate known miRNAs in B. rapa, the sequences of sRNAs were blasted against the known miRNAs of various plant species deposited in miRBase 21.0. After BLASTn and sequence analysis, 208 expressed known miRNAs were identified in both libraries. Detailed sequences and reads of the known miRNAs are shown in Table S4. Five hundred thirty-five potential novel miRNAs candidates in B. rapa were also identified and named as bra-miRn. Most of the novel miRNAs were control-specific or vernalization-specific, which accounted for 47.1 (252) and 28.4% (152), respectively. Only 131 novel miRNAs appeared in both control and vernalization libraries. Information about the novel miRNAs was listed in Table S5. We also found that conserved miRNAs level was higher than non-conserved and novel miRNAs level (Tables S4 and S5), which is in agreement with the pre￾vious findings that evolutionary conservation may have a correlation with the expression level (Chi et al. 2011). Identification of vernalization-related miRNAs in B. rapa After annotating known and novel miRNAs, we exam￾ined vernalization-responsive miRNAs based on the differential expression analysis. We used a threshold and absolute value of log2Ratio≥1 and P-value ≤0.05 to search for miRNAs up- or down-regulated during the process of vernalization. Twenty known and 66 novel miRNAs were observed to be significantly differentially expressed (Table S6). Among these, 43 miRNAs (14 known and 29 novel ones) were upregulated, and 43 miRNAs (6 known and 37 novel ones) were down￾regulated. The changes in expression levels of known miRNAs were much lower than the novel miRNAs. Of these known miRNAs, bra-miR1063a was remarkably downregulated with more than fourfold change, while bra-miR5714, bra-miR157a and bra-miR5077 were sharply upregulated with more than twofold change. The majority of the novel vernalization-responsive miRNAs were only detected in the control or vernal￾ization library, indicating that these novel miRNAs may have stage-specific functions. Of them, seven novel miRNAs (bra-miRn93, bra-miRn25 and bra-miRn345 were downregulated; bra-miRn426, bra-miRn458, bra-miRn503 and bra-miRn395 were upregulated) were significantly differentially expressed with an absolute value of log2ratio>10. These results indicated that these identified vernalization-responsive miRNAs might play vital roles in B. rapa vernalization. Prediction and annotation of vernalization-related miRNA targets To investigate the biological roles of the predicted targets in B. rapa, we searched against EST and cDNA sequences in B. rapa genome annotation database (http://brassicadb .org/brad/index.php) and Arabidopsis genome database (https://www.arabidopsis.org/index.jsp). A total of 714 and 666 potential targets for 92 known and 213 new miRNAs were identified, with an average of 7.76 and Physiol. Plant. 2018

3.13 target genes per miRNA,respectively (Table S7). miRNAs and their corresponding targets at different We viewed that the known miRNAs targeted more genes stages of vernalization (4C for 0,2,4 and 6 weeks), on average compared to novel miRNAs,and this find- 11 vernalization-related miRNAs (9 known and 2 ing suggested that the known miRNAs could function novel miRNAs)and 8 targets were randomly selected in a broader range of biological process than the novel for qPCR analysis.The expression profiles were dis- miRNAs.In addition,no targets for 116 known and played in Figs.4 and 5,and the primers were listed 322 novel miRNAs were found,possibly due to the tar- in Table S9.As Fig.4 shown,there exists a similar gets expression levels below the detection level dur- tendency between the qPCR and sRNA sequencing ing the process of vernalization.Additionally,114 tar- results of the selected miRNAs expressions,suggesting gets related to vernalization were isolated as the tar- that the results of sRNA sequencing were reliable.Nine gets of 33 vernalization-related miRNAs(11 known and miRNAs (bra-miR157,bra-miR159a,bra-miR393b 22 novel miRNAs).These targets were related to a bra-miR5020b, bra-miR5054, bra-miR5072 wide variety of biological processes,and several targets bra-miR5077,bra-miR5368 and bra-miRn429)expres- were transcription factors.For instance,MYBs(targeted sions were increased and two miRNAs (bra-miR397a by bra-miR159a)and SPLs (targeted by bra-miR157a), and bra-miRn93)expressions were reduced.The results which were in accord with the previous research in rice suggested that the levels of the tested miRNAs vary and Chinese cabbage (Sunkar et al.2008,Wang et al. significantly during the process of vernalization.Further- 2011).Additionally,many hypothetical and unknown more,some miRNAs showed stage-specific expression, function genes were targeted,indicating possible new probably involved in vernalization. functions for these miRNAs in B.rapa. Additionally,we also validated the expression pat- To functionally analyzed the annotation,the targets terns of eight corresponding targets.The qRCR results of vernalization-related miRNAs were subjected to showed that different targets of one miRNA expressed GO enrichment analysis by Blast2 GO program with varied during the process of vernalization and in differ- default parameters.GO analysis indicated that these ent tissues in B.rapa,and that most targets expression predicted targets could be classified into five molecular were relatively higher at flower or leaf tissues(Fig.5).For functions,16 biological processes and nine cellular example,NEDD1-Bra009155 exhibited highest expres- components (Fig.3).Most targets,which were clas- sion level in flower,while SAP130B-Bra034172 was sified as the binding category,encode transcription highly expressed in leaf,and both of them are targeted factors.Since transcription factors play an important by bra-miR5077.These results showed that different tar- role in gene expression,these targets are thought to gets of the same vernalization-related miRNA may have be the major targets of miRNAs (Eldem et al.2013). different function.In addition,there are inverse corre- Interestingly,we found that most targets were enriched lations between the expression levels of the selected in hormone-mediated signaling pathway,including miRNAs and those of their targets,suggesting that these GAs-mediated signaling pathway (Fig.S1).To further miRNAs might negatively regulate their targets.The rela- investigate the biological pathways influenced by vernal- tionships between miRNAs and their targets were shown ization in B.rapa,KEGG pathway enrichment analysis in Supplementary Table S10.These findings suggested was carried out.We found 106 pathways enriched with that miRNA-mediated silencing of their potential tar- targets of significantly differentially expressed miRNAs, get genes occur during the process of vernalization in and each target was assigned to at least one KEGG B.rapa. annotation.Among these pathways,plant hormone signal transduction involved the most miRNA targets (Bra0329.54,Bra003518,Bra007720.Bra026953and Bra016754,all targeted by bra-miR393b)(Table S8).The Endogenous hormone measurements functional annotations of these genes showed that they As mentioned above,GO and KEGG results suggested all acted in auxin signal pathway.The GO and KEGG that plant hormones may function in the process of B results suggested that plant hormone (GA and auxin) rapa vernalization.To further investigate the roles of signal transduction has vital influences on the process endogenous hormones in the vernalization process,the of B.rapa vernalization. IAA and GA contents were measured in B.rapa.Sam- ples were collected from the leaves before and after ver- qPCR analysis of miRNAs and their targets nalization.As Fig.6 shown,the IAA level decreased and the GA level increased after vernalization,indicating To confirm the deep sequencing data and quan- that the IAA and GA signaling might be involved in the tify the expression patterns of both the identified process of B.rapa vernalization. Physiol.Plant.2018

3.13 target genes per miRNA, respectively (Table S7). We viewed that the known miRNAs targeted more genes on average compared to novel miRNAs, and this find￾ing suggested that the known miRNAs could function in a broader range of biological process than the novel miRNAs. In addition, no targets for 116 known and 322 novel miRNAs were found, possibly due to the tar￾gets expression levels below the detection level dur￾ing the process of vernalization. Additionally, 114 tar￾gets related to vernalization were isolated as the tar￾gets of 33 vernalization-related miRNAs (11 known and 22 novel miRNAs). These targets were related to a wide variety of biological processes, and several targets were transcription factors. For instance, MYBs (targeted by bra-miR159a) and SPLs (targeted by bra-miR157a), which were in accord with the previous research in rice and Chinese cabbage (Sunkar et al. 2008, Wang et al. 2011). Additionally, many hypothetical and unknown function genes were targeted, indicating possible new functions for these miRNAs in B. rapa. To functionally analyzed the annotation, the targets of vernalization-related miRNAs were subjected to GO enrichment analysis by Blast2 GO program with default parameters. GO analysis indicated that these predicted targets could be classified into five molecular functions, 16 biological processes and nine cellular components (Fig. 3). Most targets, which were clas￾sified as the binding category, encode transcription factors. Since transcription factors play an important role in gene expression, these targets are thought to be the major targets of miRNAs (Eldem et al. 2013). Interestingly, we found that most targets were enriched in hormone-mediated signaling pathway, including GAs-mediated signaling pathway (Fig. S1). To further investigate the biological pathways influenced by vernal￾ization in B. rapa, KEGG pathway enrichment analysis was carried out. We found 106 pathways enriched with targets of significantly differentially expressed miRNAs, and each target was assigned to at least one KEGG annotation. Among these pathways, plant hormone signal transduction involved the most miRNA targets (Bra032954, Bra003518, Bra007720, Bra026953 and Bra016754, all targeted by bra-miR393b) (Table S8). The functional annotations of these genes showed that they all acted in auxin signal pathway. The GO and KEGG results suggested that plant hormone (GA and auxin) signal transduction has vital influences on the process of B. rapa vernalization. qPCR analysis of miRNAs and their targets To confirm the deep sequencing data and quan￾tify the expression patterns of both the identified miRNAs and their corresponding targets at different stages of vernalization (4∘C for 0, 2, 4 and 6 weeks), 11 vernalization-related miRNAs (9 known and 2 novel miRNAs) and 8 targets were randomly selected for qPCR analysis. The expression profiles were dis￾played in Figs. 4 and 5, and the primers were listed in Table S9. As Fig. 4 shown, there exists a similar tendency between the qPCR and sRNA sequencing results of the selected miRNAs expressions, suggesting that the results of sRNA sequencing were reliable. Nine miRNAs (bra-miR157, bra-miR159a, bra-miR393b, bra-miR5020b, bra-miR5054, bra-miR5072, bra-miR5077, bra-miR5368 and bra-miRn429) expres￾sions were increased and two miRNAs (bra-miR397a and bra-miRn93) expressions were reduced. The results suggested that the levels of the tested miRNAs vary significantly during the process of vernalization. Further￾more, some miRNAs showed stage-specific expression, probably involved in vernalization. Additionally, we also validated the expression pat￾terns of eight corresponding targets. The qRCR results showed that different targets of one miRNA expressed varied during the process of vernalization and in differ￾ent tissues in B. rapa, and that most targets expression were relatively higher at flower or leaf tissues (Fig. 5). For example, NEDD1-Bra009155 exhibited highest expres￾sion level in flower, while SAP130B-Bra034172 was highly expressed in leaf, and both of them are targeted by bra-miR5077. These results showed that different tar￾gets of the same vernalization-related miRNA may have different function. In addition, there are inverse corre￾lations between the expression levels of the selected miRNAs and those of their targets, suggesting that these miRNAs might negatively regulate their targets. The rela￾tionships between miRNAs and their targets were shown in Supplementary Table S10. These findings suggested that miRNA-mediated silencing of their potential tar￾get genes occur during the process of vernalization in B. rapa. Endogenous hormone measurements As mentioned above, GO and KEGG results suggested that plant hormones may function in the process of B. rapa vernalization. To further investigate the roles of endogenous hormones in the vernalization process, the IAA and GA3 contents were measured in B. rapa. Sam￾ples were collected from the leaves before and after ver￾nalization. As Fig. 6 shown, the IAA level decreased and the GA3 level increased after vernalization, indicating that the IAA and GA signaling might be involved in the process of B. rapa vernalization. Physiol. Plant. 2018

Number of genes 0 1020304050 60 cellular process metabolic process single-organism process response to stimulus developmental process biological regulation Biological process regulation of biological process signaling reproduction reproductive process multicellular organismal process cellular component organization or biogenesis localization establishment of localization growth multi-organism process cell cell part organelle Cellular component membrane organelle part macromolecular complex membrane part extracellular region membrane-enclosed lumen Molecular binding catalytic activity nucleic acid binding transcription factor activity enzyme regulator activity Fig.3.GO annotation of targets of vernalization-associated miR- function structural molecule activity NAs Physiol.Plant.2018

Fig. 3. GO annotation of targets of vernalization-associated miR￾NAs. cellular process metabolic process single-organism process response to stimulus developmental process biological regulation regulation of biological process signaling reproduction reproductive process localization multicellular organismal process cellular component organization or biogenesis establishment of localization Number of genes Biological process Cellular component Molecular function growth multi-organism process cell cell part organelle membrane organelle part macromolecular complex membrane part extracellular region membrane-enclosed lumen binding catalytic activity nucleic acid binding transcription factor activity enzyme regulator activity structural molecule activity 0 10 20 30 40 50 60 Physiol. Plant. 2018

5 1.2 9 bra-miR157a bra-miR159a bra-miR393b bra-miR397a 1.0 0.8 3 3 5 0.6 4 0.4 I 0.2 0 0 2 4 02 4 0 2 4 6 024 9 9 bra-miR5020b bra-miR5054 bra-miR5072 bra-miR5077 8 8 8 7 61 5 4 0 0 024 0 2. 1.2 bra-miR5368 bra-miRn429 bra-miRn93 1.0 41 0.8 0.6 0.4 0.2 0 0 2 4 6 02 46 Treatment time(week) Fig.4.The expression patterns of miRNAs during the process of vernalization by gRT-PCR.The transcript levels of miRNAs relative to the transcript levels of internal reference were quantified using 2-C.The expression level in the untreated sample(0h)was set to a value of 1.Each bar shows the mean +s of triplicate biological assays. Physiol.Plant.2018

5 4 3 2 1 0 0 2 46 bra-miR157a 5 4 3 2 1 0 0 2 46 bra-miR159a 5 4 3 2 1 0 0 2 46 bra-miR397a 9 8 7 6 3 2 1 0 0 2 46 bra-miR5020b 5 4 3 2 1 0 0 2 46 bra-miR5054 9 8 7 6 5 4 3 2 1 0 0 2 46 bra-miR5072 9 8 7 6 5 4 3 2 1 0 0 2 46 bra-miR5077 9 8 7 6 3 2 1 0 0 2 46 bra-miR5368 5 4 3 2 1 0 0 2 46 bra-miRn429 1.0 0.8 0.6 0.4 0.2 0 0 2 46 bra-miR393b 1.2 1.0 0.8 0.6 0.4 0.2 0 0 2 46 bra-miRn93 1.2 0 2 46 Relative expression Treatment time (week) Fig. 4. The expression patterns of miRNAs during the process of vernalization by qRT-PCR. The transcript levels of miRNAs relative to the transcript levels of internal reference were quantified using 2−ΔΔCT. The expression level in the untreated sample (0 h) was set to a value of 1. Each bar shows the mean±SE of triplicate biological assays. Physiol. Plant. 2018

complete small RNA libraries that will discover more 25 口0week ▣2week novel miRNAs,especially those specifically expressed 口4week ■6week in leaves after vernalization.In this work,we have obtained 208 known and 535 novel miRNAs involved in B.rapa vernalization and predicted their targets using bioinformatics approaches (Tables S5-S7).Most of the identified known miRNAs were evolutionarily conserved across diverse plant species,and 54 conserved miRNAs (belonging to 32 families)were detected.The expression patterns of several selected miRNAs and their targets were confirmed through qPCR analysis(Figs 4 and 5), validating the accuracy of the RNA-seg results in this 172 study. To date,nine flowering-related miRNA families were studied,which together with their target genes con- trolled the floral transition in Arabidopsis.Of these miRNAs,miR156/157 and miR159 were downregulated 14 during the process of floral transition,whereas miR172 was upregulated (Luo et al.2013).However,both sequencing and gPCR results showed that bra-miR157a was upregulated after vernalization,which targets 11 8 SPL genes (Table S10).Interestingly,previous studies reported that high levels of miR156 can prevent preco- cious flowering (Wang 2014).Vernalization just works as the first step to obtain the competence to flower.For most vernalization-requiring plants,like 'wuyueman', ra003305 91307 they are also LD plants,still needing long days after ver- nalization,which ensures that flowering does not begin during the decreasing day lengths of late fall and occurs Fig.5.The temporal and spatial expression patterns of targets during in the spring.The results indicated that vernalization the process of vernalization by qPCR.(A)The relative expression lev- could affect flowering through microRNA mechanism in els of targets during different stages of vernalization.(B)The relative B.rapa. expression levels of targets in various tissues(root,stem,leaf and flower). Each bar shows the mean +sE of triplicate assays.The targets detected GO and KEGG analyses revealed that the targets in qPCR were Bra034842 (MYB65),Bra032954 (AFB2),Bra026953 were significantly enriched in the hormone-mediated (AFB3),Bra009155(NEDD1),Bra034172(5AP130B).Bra016891(5PL9). Bra038324(5PL6),Bra003305(SPL15),Bra013077(CHR35. A B 16T 1200 GA3 OIAA Discussion With genomic and biocomputing tools widely devel- 800 oped,high-throughput sequencing becomes a straight- forward method to analyze small RNA profiling for 8 plants in most cases.To date,28645 miRNAs(miRBase database,Release June 21,2014)have been discovered 400 from various kinds of multicellular eukaryotes,including 4 humans,animals and plants.To analyze the roles of miRNAs during vernalization process in B.rapa,we constructed two sRNAs libraries from the leaves with or 0 0 without vernalization treatment.With Solexa sequenc- Control ing,we identified 25 million clean reads per library, Vernalized Control Vernalized more than in the previous study on B.rapa that generated about 9 million clean reads by deep sequencing (Wang Fig.6.Gibberellic acid 3(A)and auxine IAA(B)contents in leaves before et al.2015).This suggest that we generated a relatively and after vernalization in B.rapa. Physiol.Plant.2018

0 week 2 week 4 week 2.5 A B 2.0 1.5 6 week 1.0 0.5 0 Relative expression Bra034842 Bra032954 Bra009155 Bra034172 Bra016891 Bra038324 Bra003305 Bra013077 10 8 6 4 2 0 Relative expression Bra034842 Bra032954 Bra009155 Bra034172 Bra016891 Bra038324 Bra003305 Bra013077 root stem leaf 18 16 14 12 flower Fig. 5. The temporal and spatial expression patterns of targets during the process of vernalization by qPCR. (A) The relative expression lev￾els of targets during different stages of vernalization. (B) The relative expression levels of targets in various tissues (root, stem, leaf and flower). Each bar shows the mean±SE of triplicate assays. The targets detected in qPCR were Bra034842 (MYB65), Bra032954 (AFB2), Bra026953 (AFB3), Bra009155 (NEDD1), Bra034172 (SAP130B), Bra016891 (SPL9), Bra038324 (SPL6), Bra003305 (SPL15), Bra013077 (CHR35). Discussion With genomic and biocomputing tools widely devel￾oped, high-throughput sequencing becomes a straight￾forward method to analyze small RNA profiling for plants in most cases. To date, 28 645 miRNAs (miRBase database, Release June 21, 2014) have been discovered from various kinds of multicellular eukaryotes, including humans, animals and plants. To analyze the roles of miRNAs during vernalization process in B. rapa, we constructed two sRNAs libraries from the leaves with or without vernalization treatment. With Solexa sequenc￾ing, we identified 25 million clean reads per library, more than in the previous study on B. rapa that generated about 9 million clean reads by deep sequencing (Wang et al. 2015). This suggest that we generated a relatively complete small RNA libraries that will discover more novel miRNAs, especially those specifically expressed in leaves after vernalization. In this work, we have obtained 208 known and 535 novel miRNAs involved in B. rapa vernalization and predicted their targets using bioinformatics approaches (Tables S5–S7). Most of the identified known miRNAs were evolutionarily conserved across diverse plant species, and 54 conserved miRNAs (belonging to 32 families) were detected. The expression patterns of several selected miRNAs and their targets were confirmed through qPCR analysis (Figs 4 and 5), validating the accuracy of the RNA-seq results in this study. To date, nine flowering-related miRNA families were studied, which together with their target genes con￾trolled the floral transition in Arabidopsis. Of these miRNAs, miR156/157 and miR159 were downregulated during the process of floral transition, whereas miR172 was upregulated (Luo et al. 2013). However, both sequencing and qPCR results showed that bra-miR157a was upregulated after vernalization, which targets 11 SPL genes (Table S10). Interestingly, previous studies reported that high levels of miR156 can prevent preco￾cious flowering (Wang 2014). Vernalization just works as the first step to obtain the competence to flower. For most vernalization-requiring plants, like ‘wuyueman’, they are also LD plants, still needing long days after ver￾nalization, which ensures that flowering does not begin during the decreasing day lengths of late fall and occurs in the spring. The results indicated that vernalization could affect flowering through microRNA mechanism in B. rapa. GO and KEGG analyses revealed that the targets were significantly enriched in the hormone-mediated 16 12 8 4 0 GA3 GA3 content (ng g-1 FW) 1200 800 400 0 IAA A IAA content (ng g-1 FW) B Control Vernalized Control Vernalized Fig. 6. Gibberellic acid 3 (A) and auxine IAA (B) contents in leaves before and after vernalization in B. rapa. Physiol. Plant. 2018

signaling pathway,and most of which were involved in Vemalization GA and IAA mediated signaling pathway.We therefore measured the contents of IAA and GA in B.rapa dur- ing the process of vernalization.The results revealed that the content of IAA was reduced and the content of GA was increased after vernalization(Fig.6),suggesting that miR159 GA vernalization could induce GA3 and reduce IAA. Plant hormones,components of the flowering time syndrome,have been reported in model plants (Davis 2009).However,the regulatory roles of hormones GAMYBs in B.rapa vernalization process are still unclear. miR393 (MYB65 and MYB101) Auxin functions in the initiation of floral primordia and flower development (Alabadi et al.2009).In this work,we found five auxin-related genes were signifi- cantly changed,Bra032954,Bra026953,Bra016754, Bra003518 and Bra007720,which were homologous genes of AtAFB2,AtAFB3 and AtTIR1 respectively. AFB2.AFB3.TIR1 These five genes were all targeted by bra-miR393b. In Arabidopsis,miR393 is encoded by miR393a and miR393b (Jones-Rhoades and Bartel 2004).Exogenous IAA treatment could enhance miR393b transcription (Chen et al.2011).In our sequencing and qPCR results, bra-miR393b was reduced after vernalization,which was consistent with the decrease in the contents of IAA Prevent premature flowering after vernalization,indicating that the IAA signaling involved in flowering time regulation is regulated by Fig.7.Hypothetical model for flowering regulation by miRNAs asso- miR393 (Figs.4 and 6).MiR393 post-transcriptionally ciation with IAA and GA after vernalization in B.rapa.The up-and regulates AtTIR1,AtAFB2 and AtAFB3 to alter plant down-regulation are marked red and green,respectively. auxin responses(Xia et al.2012).Auxin response factors (ARFs)constituted a small subset of F-box-containing Brassica napus (Zanewich and Rood 1995).Interestingly, auxin receptors and played key roles in bud growth and GA3 level was also triggered by vernalization in B.rapa flower development (Nagpal et al.2005,Zhao 2010). (Fig.6).However,the increased GA3 did not lead to TIRI acted as auxin receptors and played a role in the early flowering,why? degradation of the Aux/IAAs repressors,which promotes In Arabidopsis,miR159 is involved in the floral AFB mediating auxin-related specific genes(Dharmasiri transition through targeting MYBs in response to GA et al.2005,Sato and Yamamoto 2008).Furthermore, signaling (Simpson and Dean 2002).GAs treatment osmiR393 overexpression resulted in early flowering directly degrades DELLA proteins and finally result in in rice (Xia et al.2012).These results suggested that an increase in LFY level (Yamaguchi and Abe 2012).In vernalization could reduce bra-miR393 expression by this study,we got seven targets of bra-miR159a,MYB65 downregulation of IAA,and finally prevent premature (Bra034842 and Bra002042),MYB81(Bra000531), flowering(Fig.7). MYB101 (Bra022888 and Bra021791)and MYB120 In addition to auxin,GAs also plays important roles (Bra035547 and Bra002938).Bra034842 is the homol- in flowering regulation.GAs could induce phase tran- ogous gene of AtMYB65,which is a GA-specific sitions of flowering and regulate floral development transcriptional regulator regulating the floral initia- (Mandaokar and Browse 2009).Vernalization could tion through activating LFY in the GA signaling pathway regulate GAs biosynthesis in vegetative plants,and in Arabidopsis (Gocal et al.2001).AtmiR159 overex- control regulatory signals which control the GAs respon- pression transgenic plants showed late flowering with sive genes expressions (Oka et al.2001).The ga1-3 decreasing GAMYBs(MYB33,MYB65 and MYB101)and mutant Arabidopsis,which contains a deletion in an LFY expression levels (Achard et al.2004),and a similar enzyme that acts on an early step of gibberellic acid result have been shown here (Figs.4 and 5;Table S10). biosynthesis,is responsive to vernalization (Michaels GAMYBs,a downstream component of GA response, and Amasino 1999).Additionally,the endogenous acts on binding to the promoter of LFY and then induces level of GAs was upregulated during vernalization in LFY expression (Gocal et al.2001).So,upregulated Physiol.Plant.2018

signaling pathway, and most of which were involved in GA and IAA mediated signaling pathway. We therefore measured the contents of IAA and GA3 in B. rapa dur￾ing the process of vernalization. The results revealed that the content of IAA was reduced and the content of GA3 was increased after vernalization (Fig. 6), suggesting that vernalization could induce GA3 and reduce IAA. Plant hormones, components of the flowering time syndrome, have been reported in model plants (Davis 2009). However, the regulatory roles of hormones in B. rapa vernalization process are still unclear. Auxin functions in the initiation of floral primordia and flower development (Alabadi et al. 2009). In this work, we found five auxin-related genes were signifi- cantly changed, Bra032954, Bra026953, Bra016754, Bra003518 and Bra007720, which were homologous genes of AtAFB2, AtAFB3 and AtTIR1 respectively. These five genes were all targeted by bra-miR393b. In Arabidopsis, miR393 is encoded by miR393a and miR393b (Jones-Rhoades and Bartel 2004). Exogenous IAA treatment could enhance miR393b transcription (Chen et al. 2011). In our sequencing and qPCR results, bra-miR393b was reduced after vernalization, which was consistent with the decrease in the contents of IAA after vernalization, indicating that the IAA signaling involved in flowering time regulation is regulated by miR393 (Figs. 4 and 6). MiR393 post-transcriptionally regulates AtTIR1, AtAFB2 and AtAFB3 to alter plant auxin responses (Xia et al. 2012). Auxin response factors (ARFs) constituted a small subset of F-box-containing auxin receptors and played key roles in bud growth and flower development (Nagpal et al. 2005, Zhao 2010). TIR1 acted as auxin receptors and played a role in the degradation of the Aux/IAAs repressors, which promotes AFB mediating auxin-related specific genes (Dharmasiri et al. 2005, Sato and Yamamoto 2008). Furthermore, osmiR393 overexpression resulted in early flowering in rice (Xia et al. 2012). These results suggested that vernalization could reduce bra-miR393 expression by downregulation of IAA, and finally prevent premature flowering (Fig. 7). In addition to auxin, GAs also plays important roles in flowering regulation. GAs could induce phase tran￾sitions of flowering and regulate floral development (Mandaokar and Browse 2009). Vernalization could regulate GAs biosynthesis in vegetative plants, and control regulatory signals which control the GAs respon￾sive genes expressions (Oka et al. 2001). The ga1-3 mutant Arabidopsis, which contains a deletion in an enzyme that acts on an early step of gibberellic acid biosynthesis, is responsive to vernalization (Michaels and Amasino 1999). Additionally, the endogenous level of GAs was upregulated during vernalization in Vernalization miR159 GA IAA GAMYBs miR393 (MYB65 and MYB101) LFY AFB2, AFB3,TIR1 Prevent premature flowering Fig. 7. Hypothetical model for flowering regulation by miRNAs asso￾ciation with IAA and GA after vernalization in B. rapa. The up- and down-regulation are marked red and green, respectively. Brassica napus (Zanewich and Rood 1995). Interestingly, GA3 level was also triggered by vernalization in B. rapa (Fig. 6). However, the increased GA3 did not lead to early flowering, why? In Arabidopsis, miR159 is involved in the floral transition through targeting MYBs in response to GA signaling (Simpson and Dean 2002). GAs treatment directly degrades DELLA proteins and finally result in an increase in LFY level (Yamaguchi and Abe 2012). In this study, we got seven targets of bra-miR159a, MYB65 (Bra034842 and Bra002042), MYB81 (Bra000531), MYB101 (Bra022888 and Bra021791) and MYB120 (Bra035547 and Bra002938). Bra034842 is the homol￾ogous gene of AtMYB65, which is a GA-specific transcriptional regulator regulating the floral initia￾tion through activating LFY in the GA signaling pathway in Arabidopsis (Gocal et al. 2001). AtmiR159 overex￾pression transgenic plants showed late flowering with decreasing GAMYBs (MYB33, MYB65 and MYB101) and LFY expression levels (Achard et al. 2004), and a similar result have been shown here (Figs. 4 and 5; Table S10). GAMYBs, a downstream component of GA response, acts on binding to the promoter of LFY and then induces LFY expression (Gocal et al. 2001). So, upregulated Physiol. Plant. 2018

expression of bra-miR159a and downregulated expres- microRNAs and short interfering RNAs.Proc Natl Acad sion of MYBs inhibited the promotion of LFY by GAs to Sci USA102:11928-11933 prevent premature flowering (Fig.7). Chen ZH,Bao ML,Sun YZ,Yang YJ,Xu XH,Wang JH,Han As we know,FLC,a key gene in vernalization pathway, N,Bian HW,Zhu MY(2011)Regulation of auxin is regulated by two IncRNAs,COLDAIR and COOLAIR response by miR393-targeted transport inhibitor (Swiezewski et al.2009,Heo and Sung 2011).In this response protein 1 is involved in normal development in work,we did not find any vernalization-related miRNA Arabidopsis.Plant Mol Biol 77:619-629 targeting FLC.Our findings showed that miRNAs mainly Chi X,Yang Q,Chen X,Wang J,Pan L,Chen M,Yang Z, act on interactions between vernalization and hormone He Y,Liang X,Yu S(2011)Identification and pathways,rather than directly act on vernalization.We characterization of microRNAs from peanut(Arachis therefore speculated that vernalization can regulate flow- hypogaea L.)by high-throughput sequencing.PloS One 6:e27530 ering by modulating hormones level via miRNAs and by modulating FLC at the epigenetic level.These results Chuck G,Cigan AM,Saeteurn K,Hake S(2007)The help us understand the roles of miRNAs in vernalization heterochronic maize mutant Corngrass1 results from overexpression of a tandem microRNA.Nat Genet 39: pathway in B.rapa.Whether these effects are produced 544-549 via miRNA-target modules requires further functional Conesa A,Gotz S,Garcia-Gomez JM,Terol J,Talon M, analysis. Robles M(2005)Blast2GO:a universal tool for annotation,visualization and analysis in functional Author contributions genomics research.Bioinformatics 21:3674-3676 Davis SJ(2009)Integrating hormones into the Conception and design of the study:F.H.and T.L.Per- floral-transition pathway of Arabidopsis thaliana.Plant formed the experiments:F.H.,S.S.and X.W.Wrote Cell Environ 32:1201-1210 the paper:F.H.Manuscript revision and approval:F.H., Dharmasiri N,Dharmasiri S,Weijers D,Lechner E, T.L.and X.H.All authors red and approved the final Yamada M,Hobbie L,Ehrismann JS,Jurgens G,Estelle M manuscript. (2005)Plant development is regulated by a family of auxin receptor F box proteins.Dev Cell 9:109-119 Dobrev Pl,Vankova R(2012)Quantification of abscisic Acknowledgements-This work was supported by grants acid,cytokinin,and auxin content in salt-stressed plant from the Fundamental Research Funds for the Central Uni- tissues.Methods Mol Biol 913:251-261 versities(Y0201700179),the Natural Science Foundation of Eldem V,Okay S,Unver T(2013)Plant microRNAs:new Jiangsu province (BK20171374)and National Natural Sci- players in functional genomics.Turk J Agric For 37:1-21 ence Foundation of China(Nos.31330067,31301782). Gardner PP,Daub J,Tate JG,Nawrocki EP,Kolbe DL, Lindgreen S,Wilkinson AC,Finn RD,Griffiths-Jones S, Eddy SR,Bateman A(2009)Rfam:updates to the RNA References families database.Nucleic Acids Res 37:D136-D140 Achard P,Herr A,Baulcombe DC,Harberd NP(2004) Gocal GF,Sheldon CC,Gubler F,Moritz T,Bagnall DJ, Modulation of floral development by a MacMillan CP,Li SF,Parish RW,Dennis ES,Weigel D gibberellin-regulated microRNA.Development 131: King RW(2001)GAMYB-like genes,flowering,and 3357-3365 gibberellin signaling in Arabidopsis.Plant Physiol 127: Agar G,Turker M,Battal P,Emre EM(2006)Phytohormone 1682-1693 levels in germinating seeds of Zea mays L.exposed to Heo JB,Sung S(2011)Vernalization-mediated epigenetic selenium and aflatoxines.Ecotoxicology 15:443-450 silencing by a long intronic noncoding RNA.Science Alabadi D,Blazquez MA,Carbonell J,Ferrandiz C, 331:76-79 Perez-Amador MA(2009)Instructive roles for hormones Huang F,Tang J,Hou X(2016)Molecular cloning and in plant development.Int J Dev Biol 53:1597-1608 characterization of BcCSP1,a Pak-choi (Brassica rapa Allen E,Xie Z,Gustafson AM,Carrington JC(2005) Ssp.chinensis)cold shock protein gene highly microRNA-directed phasing during trans-acting siRNA co-expressed under ABA and cold stimulation.Acta biogenesis in plants.Cell 121:207-221 Physiol Plant 38:1-8 Audic S,Claverie JM(1997)The significance of digital Irish VF(2010)The flowering of Arabidopsis flower gene expression profiles.Genom Res 7:986-995 development.Plant J 61:1014-1028 Bartel DP(2004)MicroRNAs:genomics,biogenesis, Jiao Y,Wang Y,Xue D,Wang J,Yan M,Liu G,Dong G, mechanism,and function.Cell 116:281-297 Zeng D,Lu Z,Zhu X,Qian Q,Li J (2010)Regulation of Baumberger N,Baulcombe D(2005)Arabidopsis OsSPL14 by OsmiR156 defines ideal plant architecture ARGONAUTE1 is an RNA Slicer that selectively recruits in rice.Nat Genet 42:541-544 Physiol.Plant.2018

expression of bra-miR159a and downregulated expres￾sion of MYBs inhibited the promotion of LFY by GAs to prevent premature flowering (Fig. 7). As we know, FLC, a key gene in vernalization pathway, is regulated by two lncRNAs, COLDAIR and COOLAIR (Swiezewski et al. 2009, Heo and Sung 2011). In this work, we did not find any vernalization-related miRNA targeting FLC. Our findings showed that miRNAs mainly act on interactions between vernalization and hormone pathways, rather than directly act on vernalization. We therefore speculated that vernalization can regulate flow￾ering by modulating hormones level via miRNAs and by modulating FLC at the epigenetic level. These results help us understand the roles of miRNAs in vernalization pathway in B. rapa. Whether these effects are produced via miRNA-target modules requires further functional analysis. Author contributions Conception and design of the study: F.H. and T.L. Per￾formed the experiments: F.H., S.S. and X.W. Wrote the paper: F.H. Manuscript revision and approval: F.H., T.L. and X.H. All authors red and approved the final manuscript. Acknowledgements – This work was supported by grants from the Fundamental Research Funds for the Central Uni￾versities (Y0201700179), the Natural Science Foundation of Jiangsu province (BK20171374) and National Natural Sci￾ence Foundation of China (Nos. 31330067, 31301782). References Achard P, Herr A, Baulcombe DC, Harberd NP (2004) Modulation of floral development by a gibberellin-regulated microRNA. Development 131: 3357–3365 Agar G, Turker M, Battal P, Emre EM (2006) Phytohormone levels in germinating seeds of Zea mays L. exposed to selenium and aflatoxines. Ecotoxicology 15: 443–450 Alabadi D, Blazquez MA, Carbonell J, Ferrandiz C, Perez-Amador MA (2009) Instructive roles for hormones in plant development. Int J Dev Biol 53: 1597–1608 Allen E, Xie Z, Gustafson AM, Carrington JC (2005) microRNA-directed phasing during trans-acting siRNA biogenesis in plants. Cell 121: 207–221 Audic S, Claverie JM (1997) The significance of digital gene expression profiles. Genom Res 7: 986–995 Bartel DP (2004) MicroRNAs: genomics, biogenesis, mechanism, and function. Cell 116: 281–297 Baumberger N, Baulcombe D (2005) Arabidopsis ARGONAUTE1 is an RNA Slicer that selectively recruits microRNAs and short interfering RNAs. Proc Natl Acad Sci USA 102: 11928–11933 Chen ZH, Bao ML, Sun YZ, Yang YJ, Xu XH, Wang JH, Han N, Bian HW, Zhu MY (2011) Regulation of auxin response by miR393-targeted transport inhibitor response protein 1 is involved in normal development in Arabidopsis. Plant Mol Biol 77: 619–629 Chi X, Yang Q, Chen X, Wang J, Pan L, Chen M, Yang Z, He Y, Liang X, Yu S (2011) Identification and characterization of microRNAs from peanut (Arachis hypogaea L.) by high-throughput sequencing. PloS One 6: e27530 Chuck G, Cigan AM, Saeteurn K, Hake S (2007) The heterochronic maize mutant Corngrass1 results from overexpression of a tandem microRNA. Nat Genet 39: 544–549 Conesa A, Götz S, García-Gómez JM, Terol J, Talón M, Robles M (2005) Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics 21: 3674–3676 Davis SJ (2009) Integrating hormones into the floral-transition pathway of Arabidopsis thaliana. Plant Cell Environ 32: 1201–1210 Dharmasiri N, Dharmasiri S, Weijers D, Lechner E, Yamada M, Hobbie L, Ehrismann JS, Jürgens G, Estelle M (2005) Plant development is regulated by a family of auxin receptor F box proteins. Dev Cell 9: 109–119 Dobrev PI, Vankova R (2012) Quantification of abscisic acid, cytokinin, and auxin content in salt-stressed plant tissues. Methods Mol Biol 913: 251–261 Eldem V, Okay S, Unver T (2013) Plant microRNAs: new players in functional genomics. Turk J Agric For 37: 1–21 Gardner PP, Daub J, Tate JG, Nawrocki EP, Kolbe DL, Lindgreen S, Wilkinson AC, Finn RD, Griffiths-Jones S, Eddy SR, Bateman A (2009) Rfam: updates to the RNA families database. Nucleic Acids Res 37: D136–D140 Gocal GF, Sheldon CC, Gubler F, Moritz T, Bagnall DJ, MacMillan CP, Li SF, Parish RW, Dennis ES, Weigel D, King RW (2001) GAMYB-like genes, flowering, and gibberellin signaling in Arabidopsis. Plant Physiol 127: 1682–1693 Heo JB, Sung S (2011) Vernalization-mediated epigenetic silencing by a long intronic noncoding RNA. Science 331: 76–79 Huang F, Tang J, Hou X (2016) Molecular cloning and characterization of BcCSP1, a Pak-choi (Brassica rapa Ssp. chinensis) cold shock protein gene highly co-expressed under ABA and cold stimulation. Acta Physiol Plant 38: 1–8 Irish VF (2010) The flowering of Arabidopsis flower development. Plant J 61: 1014–1028 Jiao Y, Wang Y, Xue D, Wang J, Yan M, Liu G, Dong G, Zeng D, Lu Z, Zhu X, Qian Q, Li J (2010) Regulation of OsSPL14 by OsmiR156 defines ideal plant architecture in rice. Nat Genet 42: 541–544 Physiol. Plant. 2018