Comparative profiling of genome-wide small RNAs in non-heading Chinese cabbage [Brassica rapa subsp. chinensis (L.) Hanelt]

Non-heading Chinese cabbage is characterized with diversity of leaf shape and of leaf architecture. Small RNAs are broadly involved in leaf development. To understand how small RNAs regulated leaf shape and leaf architecture, we choose two representative genotypes of non-heading Chinese cabbage (Brassica. rapa ssp. chinensis cv Wut) with distinct leaf shapes, and performed small RNA deep sequencing. Huaq, a genotype with long blades and wide petioles, generated 13.30 million small RNA reads while Wut, a genotype with round blades and long petioles, produced 14.69 million. After normalizing to RP10M, they share 0.5 million small RNAs with about 4 million reads which are mapping to the genome. Four types of small RNAs were selected from the libraries: miRNAs, ta-siRNAs, mRNA-derived siRNAs and rasiRNAs. 45 miRNA members among 75 display different expressional level such as miR166 and miR319, which play important role in leaf development. It indicated the different leaf phenotype between Huaq and Wut that possibly regulated by miRNA dosage. One of the small RNAs originated from TAS3 is Huaq dominant. About 2 thousand small RNAs derived from genes were expressionally different between genotypes, and some of them derived from antisense chains of genes. Chloroplast derived small RNAs were found that Wut were more than Huaq, which is consistent with its greener leaves. Comparative analysis of small RNAs on genome-wide level will broaden our view of molecular events underling diversity of leaf shape and architecture.

between different genotypes. rapa and Arabidopsis thaliana, diverged (A. thaliana), diverged from a common ancestor separated approximately 13-17 million years ago (Mun et al. 2009), and the B. rapa genome contains triplicated homologous counterparts of corresponding segments comparing to the genome of A. thaliana due to the events of entire genome triplication .
Small RNAs are involved in heat response of Chinese cabbage . It is not known how small RNAs contributes to diversity of leaf shape. In general, small RNAs cause transcriptional gene silencing by guiding heterochromatin formation at homologous loci, whereas others lead to posttranscriptional gene silencing through mRNA degradation or translational inhibition. Some of small RNAs are derived from and target foreign nucleic acids such as viruses and transgenes, but most of them are derived from endogenous loci and regulate a multitude of developmental and physiological processes, and response to abiotic or biotic stress (Chen 2009).
Small RNA is a big family, and there are four major classes of endogenous small RNAs in plants: micro-RNAs (miRNA), transacting siRNAs (ta-siRNA), natural antisense transcript-derived siRNA (nat-siRNA) and heterochromatic siRNAs (hc-siRNA). Small RNAs can arise through distinct biogenesis routes and function through loading into Argonauteprotein-containing effector complexes (Jamalkandi 2009, Shukla et al.2008. miRNAs are produced from loci which encoded a fairly long transcript (pri-miRNA), this is processed by DICERLIKE1(DCL1), HYPONASTIC LEAVES 1(HYL1) and SERRATE (SE) to a shorter, hairpin structure, the pre-miRNA. The mature miRNA is produced by DCL1 processing of the pre-miRNA and is taken up by AGO1 into the RNAi silencing complex(RISC). AGO1's RNAase activity is directed to the target of the miRNA (Voinnet 2009). ta-siRNAs are a second class of siRNA produced under the direction of miRNA. ta-siRNAs are encoded by TAS loci. TAS loci transcripts are processed by miRNA directed cleavage into short single strand fragments. Second strand synthesis via RDR6 and DCL4 produces 21 nt long ta-siRNAs. Like miRNAs, ta-siRNAs are taken up by AGO1 into a silencing complex. miRNAs are found in plants and animals, but ta-siRNAs have so far been identified only in plants (Montgomery 2008).
Small RNAs have been implicated as factors in many aspects of plant development. They comprise several new layers of gene regulation, affecting transcription, stability, localization and translation. Such as miR156 control the vegetative phase and floral transition (Wang et al. 2008;Wut et al. 2006), miRNA 160 and miRNA 167 controls hormone biosynthesis and signaling (Mallory et al. 2005;Reyes et al. 2007), miRNA 165/166 effects leaf polarity formation and morphogenesis (Juarez et al. 2004, Sieber et al. 2007, and Others function in stress resistance such as drought (Liu et al. 2008), over-oxidation (Sunkar et al. 2006), and phosphate starvation (Fujii et al. 2005).
Leaves develop from the shoot apical meristem (SAM), a group of stem cells at the apex of the shoot. On the flanks of the SAM a group of cells loses indeterminacy and becomes a leaf primordia. These cells take cues from their position in the SAM to begin to divide, initially outwards to form a peg, and then divisions become localized to the margin between the upper and lower side to form the leaf blade or lamina (Steeves and Sussex 1989). Many factors can affect the development of leaf from leaf primordia to a mature leaf with adaxial side and abaxial side, several papers described the regulation of leaf polarity by the miR390/ARF pathway (Nagasaki et al. 2007;Chitwood et al. 2009). miR165/166 is conserved amongst seed plants and negatively regulates accumulation of members of the HD-ZIPIII family of transcription factors involved in meristem function, vascular patterning and adaxial cell fate (Juarez et al. 2004a, Wenkel et al. 2007a. miR164A is required for regulation of CUC gene expression in the developing leaf lamina (Nikovics et al. 2006).
Recently, small RNA deep sequencing is widely used for genomics research and is an important tool for analysis of the functional complexity of transcriptomes. It provides an approach to understand why close relative species or varieties present diverse phenotypes. In Arabidopsis, the two closely related species A. thaliana and A. lyrata differ greatly in Transposable Elements (TE) copy number, and every major TE family has more copies and TEs have been more active in A. lyrata since the divergence from A.thaliana (Hollister et al. 2011). It indicates that siRNA direct gene expression could depend on active ability of siRNAs which will differ among species. Small RNAs expression could be an index for plant development and response to environment changing, but there is little known about the difference of small RNAs between species with close genetic relation. Huaq and Wut were chosen as our research objects for their different morphotype and different abilities to adapt to environment. Huaq has long blades and wide petioles while Wut has round blades and long petioles. Wut are more tolerant to cold environment than Huaq. Identification of the small RNA expression differences between Huaq and Wut offers the opportunity to understand the epigenetic molecular mechanisms in transcriptional and post-transcriptional level.

Materials and methods
Plant material and small RNA sequencing Huaq and Wut, two inbred lines belonging to nonheading Chinese cabbage (Brassica. rapa ssp. chinensis cv Wut), was used for constructing small RNA libraries. All plants were grown under a 16/8 h light/dark photoperiod at 22°C for 3 weeks, the above-ground portions of the seedlings were harvested, quick-frozen immediately with liquid nitrogen and stored at -80°C. Total RNAs were extracted using TRIzol Reagent (Invitrogen). RNA concentration was measured by Eppendorf Biophotometer 6121. DNA was removed by digestion of using TURBO DNA-free kit (Ambion). RNA samples were sent to Keygene N.V., Wageningen in Netherlands where the samples were prepared by using the MIRVana MIRNA Isolation Kit (Ambion, Austin, TX, USA)and sequenced on Illumina GAII sequencer, according to the Alternative v1.5 Protocol (Illumina, San Diego, CA, USA).
Annotation of the origin of small RNAs Small RNAs which are miRNAs were identified by comparing to the former annotated miRNA of B.rapa . ta-siRNAs were aligned to the TAS genes in B.rapa homologous to Arabidopsis. Small RNAs which generated from genes were selected by comparing the location of small RNAs with that of Brassica gene models (http://brassicadb.org). Further, these genes were annotated with Arabidopsis gene name and Gene ontology (GO) analysis. Small RNAs which locus in genome are more than 50 were considered as ra-siRNA. For statistic analysis, small RNAs originated from other kinds of RNA, they were aligned to Arabidopsis chloroplast genome, chondriosome genome, transposon elements, ribosome RNA, tRNA and snoRNA.

Analysis of expressional difference
The abundance of all mapped reads was normalized to the number of reads per 10 million (RP10M). Small RNAs which only detected in one genotype with more than 10 reads was considered genotype enriched, the abundance of small RNAs which ratio is more than 2 folds change between genotypes were considered expressional different, and defined as genotype dominant.

The distinct phenotypes of Huaq and Wut
Huaq and Wut are two inbred lines of non-heading Chinese cabbage (Brassica rapa ssp. chinensis var. Huaq). To demonstrate difference in leaf shape and leaf architecture between two genotypes, we made careful observation and measurement by using the plants of 6-week old plants. The blades of Huaq were lightly green, relatively flat while the petioles are white, wide and thick (Fig. 1). The blades of Wut were deeply green, curved downward while the petioles are green, long and thick. Comparing to Huaq, Wut was more tolerant to low temperature. Meanwhile, Wut had much more rosette leaves that were clustered into many whorls.
Characterization and classification of small RNA reads Before analysis, all small RNA reads were filtered out with low-complexity, low quality and repetitive reads, and 13,295,284 reads in Huaq and 14,690,798 reads in Wut were obtained. For convenient comparison of abundance between two libraries with different total number, the reads of small RNAs in Huaq and Wut was normalized to the number of reads per 10 million (RP10M)(normalized read = raw read 9 10 million/total reads). the data of normalized reads was used for further analysis (Table 1).
For mapping to the genome, small proportion of small RNAs with length less than 17 nt were filtered out. The reads sequences were aligned against the Chinese cabbage genome (http://brassicadb.org) by using the SOAP2 software, with normalized number of reads which was set to threshold of less than one mismatch. In total reads, Huaq has 63.07% and Wut has 61.88% reads mapped to B. rapa genome (Table 1). About one third of total reads did not be mapped on genome which may because some gap in whole B. rapa genome or different sequence variance between the genotypes and reference genome.
6.30 million and 6.18 million small RNA reads in Huaq and Wut were mapped on B. rapa genome database. Among them, there were about 2.10 million unique reads for Huaq, and 1.73 million for Wut. 4,110,151 and 4,702,562 small RNA reads of Huaq and Wut shared 497,645 unique small RNAs, indicating the apparent difference on accumulation of common small RNAs between Huaq and Wu. Interestingly, 1.42 million reads were matched on sense or antisense strands of gene regions in Huaq whereas 1.50 million reads on gene regions in Wut. Between two genotypes, total numbers of miRNAs, ta-siRNAs, and rasiRNAs were different (Table 1).
In both genotypes, 75 miRNA families were detected. However, there were 0.45 million miRNA reads in Huaq and 0.30 million in Wut. These data revealed that more miRNAs were generated from Fig. 1 The picture of Huaq and Wut, a Huaq grows in field. b The display of leaves from Huaq seedling. c Wut grows in field. d The display of leaves from Wut seedling Huaq than Wut. More than 80 putative ta-siRNAs were generated from TAS3 precursor in both genotypes. Additionally, about 40-45 thousand small RNAs with about 0.2 million reads belonged to ra-siRNAs.
To distinguish small RNAs which exhibit different expressional level in Huaq and Wut, we found there are 15 thousand Huaq dominant small RNAs, and 18 thousand Wut dominant ( Although the total numbers of unique reads were approximately the same between Huaq and Wut, the abundance of small RNA reads categorized in different length were substantially different (Fig. 2). Among the small RNAs in the length from 17 to 32 nt, the distribution of total reads was consistent with that of unique reads. The peaks of length for unique small RNAs were 24 nt in both genotypes. The abundance 21 * 24 nt small RNAs in Huaq was more than that of Wut. However, from the abundance of 25 * 32 nt small RNAs in Wut was more than that of Huaq.
Still there was a large proportion is unknown except that some of small RNA fragments were derived from chloroplast small RNA (csRNA), chondriosome small RNA (chosRNA), rRNAs, tRNA, snRNAs and snoRNAs. For learning the proportion of these small RNAs, they were aligned to Arabidopsis Tair 9 chloroplast genome, chondriosome genome, and cDNA database (Fig. 3). At abundance distribution of different type RNAs, Wut has more csRNA and rRNA than Huaq.

Small RNA distribution in chromosome
We calculated the distribution of small RNA reads in 100-kilobase pair (100 kb) windows along ten chromosomes (Fig. 4). There were many hot spots on genome where small RNA reads are enriched. Interestingly, the peaks of the total small RNA abundance between Huaq and Wut were bearly noticeable

Difference in expression of individual miRNAs between Huaq and Wut
Recently, many species of higher plants have been identified to have species-specific miRNAs (Axtell and Bowman 2008). To verify whether miRNAs in Huaq was not present in Wut, we compared the sequences of miRNAs. Among 75 miRNAs, 3 were Wut-specifc. They were bra-miR403, bra-miR5711, bra-miR399f. Comparison of expressional difference miRNAs between Huaq and Wut showed that 76.9% miRNAs were expressed predominant for Huaq whereas 23.1% miRNAs predominant in Wut. When all miRNAs were grouped into 5 groups according to the range of abundance ([ 10,000, 1000-10,000, 250-1000, 100-250, \ 100 reads) (Fig. 5a-e), we recognized that in the group belonging to largest abundance, bra-miR156a were more abundant in Wut than in Huaq, while bra-miR159a and bra-miR166a accumulated more in Huaq than in Wut. miR156 is the most ancient miRNAs found in perhaps all land plant lineages (Axtell and Bowmen 2008). In other groups, bra-miR1885b.3 and bra-miRNA319a were upregulated 2.5 folds compared with those of Wut. and meawhile, bra-miRNA169b, bra-miR5711, bra-miR5725 and bra-miR395a were upregulated in Wut compared with those of Huaq.

ta-siRNA in Huaq and Wut
About ta-siRNAs in non-heading Chinese cabbage, only TAS3 transcripts were found by comparison with ta-siRNAs of Arabidopsis . We chose TAS3 precursors as reference, and aligned all small RNAs with them. About 80 small RNAs in each genotype were matched with TAS3 precursors sequence, and one of them was conserved with TAS3a(TTCTTGACCTTGTAAGACCCC), known ta-siRNA in Arabidopsis. The abundance of this ta-siRNA in Huaq and Wut was 10 and 3, respectively. It was aligned on chromosome 5th of the genome of B. rapa. It is not known whether the small RNAs derived from TAS3 precursors are ta-siRNAs. TAS3 ta-siRNAs are generated from miR390-mediated cleavage, and its target genes are two auxin response factor genes, ARF3 and ARF4, which promote abaxial identity and expression of adult traits in leaves (Adenot et al. 2006;Hunter et al. 2006).
Small RNA aligned to the annotated genes of B.rapa For analysis of function of small RNAs, we aligned small RNA reads to the B. rapa genes. More than 0.44 million small RNAs in Huaq and 0.41 million in Wut were matched on the sequences of the annotated genes ( Table 1). The half of them was located on antisense strands of the detected genes (Fig. 6). There were 44 Huaq-enriched small RNAs and 27 Wut-enriched small RNAs from antisense strand of genes. Using the threshold of 2 folds changing as standard, 238 small RNAs in Huaq were higher in abundance than in Wut, while 77 small RNAs in Wut was higher than in Huaq (Fig. 7a). To detect the putative target genes of these small RNAs, small RNAs were applied to match with their possible target gene. Potential targets (homologous to Arabidopsis genes) were displayed in supplementary data. These putative targets (filtered out unknown genes) included hormone-related genes such as CTR1, AREB3 and ATCHX4 (Fig. 7b). The abundance of these small RNA matched on CTR1, AREB3 were higher in Huaq, while that of ATCHX4 was higher in Wut. In addition, ATIPT3, a gene that regulates cytokinin biosynthesis, was more abundant in Huaq than in Wut. ATCOAA, MEE39, and ATCAP-E1 were also more abundant in Huaq. In contrast, ATCS-

Discussion
On the whole view of small RNAs from deep sequencing, the total numbers of 21-24nt small RNA were more abundant in Huaq than in Wut. This fact indicated that Huaq had strong ability to generate more classic siRNAs. Chloroplasts are important for photosynthesis and respiration. Abundance comparison of csRNA and rRNA revealed that small RNAs were more abundant in Wut than in Huaq. This was consistent with deep green leaves of Wut, which was indicating more chlorophylls synthesizing.
Ribosome biogenesis is a complex process that is comprised of transcription, modification, and processing of ribosomal RNA and production of ribosomal proteins (Deisenroth et al. 2010). Moreover, abundant chloroplast small RNAs in Wut may reflect strong activity of photosynsis, a process small RNAs are involved in. The difference in miRNA abundance between Huaq and Wut may cause variance of some important development processes such as adaxial/ abxial polarity, cell division and auxin response. Although the accurate roles of these miRNAs in formation of leaf shape and leaf architecture has not been identified, the difference of miRNA abundance provides a clue for exploring molecular reasons of leaf diversity in B. rapa.
There are many conserved miRNAs that effect the plant growth and development, Such as miRNA 156, which is first and one of the most ancient miRNA in land plant lineages (Axtell and Bowman 2008). Its target gene is 11 of 17 SPL genes, and SPL3, SPL4 and SPL5 promote vegetative phase change as well as floral transition. Overexpression of mi156 results in prolonged expression of juvenile characters and extremely delayed flowering (Chuck 2007.). Comparing to Huaq, more miR156 accumulated in Wut. That was coincidence with the different development period between Huaq and Wut, Wut had longer juvenile period and more rosette leaves than Huaq. For miR172, it displayed reversed expressional difference that more miR172 accumulated in Wut. In Arabidopsis, miR172 promotes flowering by repressing TOE1 and TOE2 (Aukerman and Sakai 2003), while the maize Cg1 mutant that overexpresses miR156 has reduced miR172 levels (Chuck et al. 2007). More evidence shows miR172 are directly regulated by miR156 so plays converse role to miR156 in flower transition. In the rosette leaves of Huaq and Wut, the absolute quantity of miR156 was largely more than miR172, and the expressional difference of miR172 was in opposite gradient to miR156 between Huaq and Wut. It was consistent that miR156 also played a critical role for phase change and repression of miR172 in no heading Chinese cabbage. b Fig. 4 The distribution and abundance show of small RNAs on the genome of B.rapa. Wut has a hot spot at end of chromosome 9 th , which disappeared at same place of Huaq The miR166 is conserved amongst seed plants and negatively regulates accumulation of members of the transcription factors HD-ZIPIII families involved in meristem function, vascular patterning and adaxial cell fate (Zhong and Ye 2007;Wenkel et al. 2007a, b). Accumulation of miR165/166 on the abaxial side of leaf primordia, delimits the HD-ZIPIII transcripts in adaxial side of leaf and effects the polarity of the leaf (Kidner et al. 2004, Juarez et al. 2004). Overexpression of miR165/166 will conduct adaxial character leaves. Huaq exhibited accumulated more mature miR166 than Wut (Fig. 5), whether lack of miR166 will lead Wut to present somewhat curled leaves still need more experiment investigation. In Arabidopsis, miR319a was initially characterized based on the dramatic leaf phenotype that results from overexpression of miR319a, as in the jaw-D allele that results in a jagged and wavy (jaw) leaf-phenotype (Palatnik et al. 2003). However, Wut seemed to accumulate less miR319. The expressional difference of target genes of miR319 needs to be testified. TCP genes as miR319 target sequences are found in a wide range of species, indicating that miRNA-mediated control of leaf morphogenesis is conserved in plants with very different leaf forms (Palatnik et al. 2003.). In this case, Huaq and Wut exhibit vivid leaf type.
Biogenesis of miRNAs must be rigorously regulated, and many miRNA express in a spatiotemporally regulated manner or in response to environmental stimuli (Chen 2009). For the purpose to investigate miRNA-mediated leaf development in Chinese cabbage, an attempt has been made in our lab to verify the morphological consequences of several miRNAs differentially expressed in Huaq and Wut.
We analyzed the miRNAs derived from antisense strand of Brassica genes, three genes Bar000375(-ATHMA1), Bra037319(ATHMA3) and Bra0003728(ATIPT3) (Fig. 7), with difference expresses. Arabidopsis HAIRY MERISTEM (HAM) family of transcription regulators act as conserved interacting cofactors with WUS/WOX proteins. HAM and WUS share common targets in vivo and their physical interaction is important in driving downstream transcriptional programs and in promoting shoot stem cell proliferation (Zhou Y et al. 2015). ATIPT3 Encodes cytokinin synthase involved in cytokinin biosynthesis. The accumulation of cytokinins in plants expressing AtIPT3 has a significant effect on root development, increased Cytokinin Accumulation Affects Plant Development and Cell Proliferation, over-expression of AtIPT3 increased leaf size was primarily the result of an increased cell number. (Galichet et al. 2008.), Huaq and Wut has differences in leaf shape could be different in expression of those genes.

Declarations
Conflict of interest We declare that we have no financial and personal relationships with other people or organizations that can inappropriately influence our work, there is no professional or other personal interest of any nature or kind in any product, Fig. 6 The distribution of small RNAs mapped to gene regions of Huaq and Wut. a The distribution of small RNAs derived from sense strand of gene regions. b The distribution of small RNAs derived from antisense strand of gene regions. The diagram was draw by logarithmic value of abundance service and/or company that could be construed as influencing the position presented in, or the review of, the manuscript entitled, ''Comparative profiling of genome-wide small RNAs in non-heading Chinese cabbage (Brassica rapa ssp. chinensis)''. Fig. 7 The abundance of small RNAs located in antisense strand of genes. a The number of the small RNAs located in antisense strand of genes. b The heat map for abundance of small RNAs located in antisense strand of Brassica genes with the name of orthologous genes of Arabidopsis