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 controlled 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 Arabidopsis 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 environmental 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 function in the floral transition process (Li et al. 2015). Recent studies reported that GAs function to induce floral transition 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 developmental 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 transcriptional 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 complementarity. 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 verified 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 controlling MYBs expression (Tsuji et al. 2006). Generally, 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 numerous miRNAs have been identified in B. rapa, no systematic 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 expression profiles of vernalization-related miRNAs and their targets were validated and analyzed by real-time quantitative 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 regulatory 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 greenhouse (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 harvested 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 completed following reported procedures (Sunkar et al. 2008). Solexa technology was used to sequence sRNAs at BGI (Shenzhen, China). Small RNA reads were generated 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 analysis. 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 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.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 miRNAs. 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 vernalization), the read counts of miRNAs were normalized as transcripts per million according to Bayesian methods (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 upor 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 UniGene 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 phytohormone extraction. IAA and GA3 contents were measured 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 following 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 <0.05 to search for miRNAs up-or down-regulated during the process of vernalization.Twenty known and 66 novel Fig.1.Summary of the common and specific sequences of the unique sRNAs from control and vernalization library. miRNAs were observed to be significantly differentially expressed (Table S6).Among these,43 miRNAs (14 known and 29 novel ones)were upregulated,and 43 01 ■Control miRNAs (6 known and 37 novel ones)were down- 5 ■Vernalization regulated.The changes in expression levels of known 40 miRNAs were much lower than the novel miRNAs.Of 50050 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. 15 The majority of the novel vernalization-responsive 10 miRNAs were only detected in the control or vernal- 5 ization library,indicating that these novel miRNAs may 0 have stage-specific functions.Of them,seven novel 15161718192021222324252627282930 miRNAs (bra-miRn93,bra-miRn25 and bra-miRn345 Length (nt) were downregulated;bra-miRn426,bra-miRn458, bra-miRn503 and bra-miRn395 were upregulated)were Fig.2.Length distribution and frequency percent of sRNA sequences in control and vernalization libraries in B.rapa.Y-axis represents percent- significantly differentially expressed with an absolute ages of sRNAs identified in this study:x-axis represents the length of value of logzratio>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 percentages 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 previous 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 examined 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 downregulated. 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 vernalization 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 finding 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 targets expression levels below the detection level during the process of vernalization. Additionally, 114 targets related to vernalization were isolated as the targets 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 classified 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 vernalization 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 quantify 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 displayed 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) expressions 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. Furthermore, some miRNAs showed stage-specific expression, probably involved in vernalization. Additionally, we also validated the expression patterns of eight corresponding targets. The qRCR results showed that different targets of one miRNA expressed varied during the process of vernalization and in different 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 expression 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 targets of the same vernalization-related miRNA may have different function. In addition, there are inverse correlations between the expression levels of the selected miRNAs and those of their targets, suggesting that these miRNAs might negatively regulate their targets. The relationships between miRNAs and their targets were shown in Supplementary Table S10. These findings suggested that miRNA-mediated silencing of their potential target 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. Samples were collected from the leaves before and after vernalization. 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