Comparative Pro�ling of Small RNAs Originated From Two Varieties of Chinese Cabbage (Brassica Rapa) With Distinct Leaf Traits

The species Brassica rapa includes many important vegetable crops. For a genome-wide survey of small RNAs in B. rapa, we performed massively parallel small RNA deep sequencing by using two representative varieties of non-heading Chinese cabbage (B. rapa ssp. chinensis) that are distinct in leaf shape, size, color and curvature. In total, 13.30 million small RNA reads were generated from Huaq, a variety with up-curved blades and wide petioles, and 14.69 million for Wut, a variety with down-curved blades and narrow petioles. After normalization, about half of small RNA reads in each variety were mapped to the published reference genome of B. rapa. In Huaq seedlings, unique small RNAs were much more than in Wut seedlings. Among them, 45 miRNAs were up-regulated or down-regulated in one variety, compared with those in another variety. Numbers of ta-siRNAs (trans-acting siRNAs) and ra-siRNAs (repeat-associated siRNAs) in Huaq seedlings were more than those detected in Wut seedlings, while gene-derived siRNAs (siRNAs derived from the sense and antisense strands of annotated genes regeion) in Huaq seedlings were less than in Wut seedlings. Especially, bra-miR156, bra-miR166 and bra-miR319 and one of ta-siRNAs that play important role in leaf shape, leaf size, leaf curvature and phase transition were differentially expressed between two varieties. In addition, total number of small RNAs derived from chloroplast genome of Wut was more than that of Huaq. The differentially-expressed small RNAs on genome-wide level provides the clue for the study of small-RNA-mediated morphology of leaves and is useful for developing of small RNA markers for leaf traits desirable for high yield and quality.


Introduction
Brassic rapa (B.rapa) includes many kinds of important vegetable crops.The edible parts of these crops are leafy heads, leaf clusters, modi ed blades, and modi ed petioles.The distinct morphologies exhibited by subspecies of B. rapa represent one of the most spectacular evolutionary changes that illustrate the structural evolution of plants under domestication.There are multiple genetic and epigenetic variances between different genotypes in B. rapa and A. thaliana, both of which diverged from a common ancestor separated approximately 13-17 million years ago (Mun et al., 2009).The B. rapa genome contains triplicated homologous counterparts of corresponding segments comparing to the genome of Arabidopsis.thaliana due to the events of entire genome triplication (Wang et al., 2011b).
Small RNAs of Chinese cabbage are involved in heat responses (Wang et al., 2011a;Yu et al., 2011).It is not known how small RNAs contributes to diversity of leaf shape during evolution for different genotypes.
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 RNAs is a big family, and there are four major classes of endogenous small RNAs in plants: microRNAs(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 Argonaute-protein-containing effector complexes (Jamalkandi and Masoudi-Nejad, 2009).miRNA are produced from loci which encoded a fairly long transcript(pri-miRNA), which is processed by DICERLIKE1(DCL1), HYPONASTIC LEAVES 1(HYL1) and SERRATE(SE) to a shorter hairpin structure, the pre-miRNA.The mature miRNA is recruited by AGO1 into the RNAi silencing complex(RISC).RNAase activity of AGO1 is directed to the target of the miRNA (Voinnet, 2009).ta-siRNAs are a second class of siRNA produced under direction of miRNA.Like miRNAs, ta-siRNAs are taken up by AGO1 into a silencing complex.miRNA are found in plants and animals but ta-siRNA have so far been identi ed only in plants (Montgomery et al., 2008).
Leaf develops from the shoot apical meristem (SAM), a group of stem cells at the apex of the shoot.On the anks of the SAM a group of cells loses indeterminacy and becomes a leaf primordium.These cells take cues from their position in the SAM to begin to divide, initially outwards to form a peg (a group of cell of leaf primordium), and then divisions become localized to the margin between the upper and lower side to form the leaf blade or lamina (Steeves, 1989).Many factors can affect the development of leaf from leaf primordia to a mature leaf with adaxial side and abaxial side (Chitwood et al., 2009;Juarez et al., 2004;Nagasaki et al., 2007;Wenkel et al., 2007).
Huaq and Wut are two varieties of non-heading Chinese cabbage.They are distinct in shape, size, color and curvature of leaves although the shortened shoots and roots are similar.The yield and quality of Chinese cabbage are largely dependent on shape, size, color and curvature of leaves.Recently, several microRNAs (miRNAs) and siRNAs have been identi ed to regulate leaf shape, curvature and size in Arabidopsis since they are involved in leaf polarity, cell division, auxin response and phase transition (Liu et al.;Wu et al., 2007;Yu et al., 2005).A question remains whether and how small RNAs regulate leaf phenotypes of Chinese cabbage, a closely relatives of Arabidopsis.At now, Small RNA deep sequencing is widely used for genomics research and is an important tool for analysis on the functional complexity of transcriptome (Hollister et al., 2011).It provides an approach to understand why close relative species or varieties present diverse phenotypes, and types and expression levels of some small RNAs could be regarded as molecular markers for leaf traits as well.Identi cation of differently expression of small RNAs between Huaq and Wut offers the opportunity to understand the molecular mechanisms involved in leaf variance at transcriptional and post-transcriptional level.

Materials And Methods
Plant material and small RNA sequencing Huaq-1 and Wut-1 are the two inbred lines of Huaq (B.rapa.ssp.chinensis var.huaq) and Wut (B.rapa.ssp.chinensis var.wut).All plants were grown under a 16/8 h light/dark photoperiod at 22℃ 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 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).

Alignments of raw reads with reference genome
Recently, genome sequence of Chinese cabbage (Chiifu-401-42) genotype has been reported (Wang et al., 2011b).This genome sequence has made us possible to match small RNA sequences on the good reference data.After removing adaptors, the small RNA with length of 17 ~ 36 nucleotides (nt) were mapped to genomic sequences (BRAD) 1 of Chinese cabbage  by SOAP2 (short oligonucleotide alignment program) (Li et al., 2009) with the 0 or 1 mismatch.For learning small RNAs derived from chloroplast, mitochondria, rRNAs, tRNA, TE (transposable elements) and snoRNAs (small nucleolar RNA), all reads were aligned with Arabidopsis TAIR9 cDNA database.
Annotation of the origin of small RNAs Small RNAs which are miRNAs were identi ed by comparing to the former annotated miRNA of B. rapa (Yu et al., 2011).The ta-siRNAs that are homologous to Arabidopsis homologs were aligned to the TAS genes in B. rapa.Small RNAs that were generated from the genes were selected by comparing the location of small RNAs with Brassica gene (http://brassicadb.org).Then, these genes were annotated on the basis of homologous genes of Arabidopsis and Gene Ontology (GO) analysis.Small RNAs that have more than 50 multiple locations on genome sequences were considered as ra-siRNA (repeat-associated siRNA).Meanwhile, small RNAs were aligned to Arabidopsis chloroplast genome, mitochondria genome, transposable 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 cultivar with more than 10 reads was considered cultivarenriched.The identical reads whose P-value cutoffs of 0.01 and 2-fold change of Huaq-to-Wut ratio were considered to be different in expression levels, and de ned as Huaq-predominant, and vice versa.The Pvalue was calculated based on the Bayes inference: P = 0.5 (N1+N2) where N1 and N2 are the reads of two samples, respectively.The small RNAs that were speci cally expressed in one variety rather than another were de ned as speci c reads.

Results
The distinct phenotypes of Huaq and Wut seedlings Huaq-1 and Wut-1 are two inbred lines of Huaq and Wut, respectively.To demonstrate difference in leaf shape and leaf architecture between two varieties, we carefully observed leaf shape, size, color and curvature of the two inbred lines at rosette stage (6-week old).For the variety Huaq, leaves were big; blade was oval in shape; and petioles were extremely wide and thick (Fig. 1A, B, Table 1).For the variety Wut, leaves were small; blades were round and dark-green; and petioles were long and narrow.Especially, Huaq petioles were characterized with inward curvature (Fig. 1C, D).
Table 1.Leaf parameters of Huaq and Wut.The data were recorded 6 weeks after sowing and were the mean of 20 individuals.
Characterization and classi cation of small RNA reads Before analysis, all small RNA reads were ltered out with low-complexity, low quality and repetitive reads.In total, 13,295,284 reads for Huaq and 14,690,798 reads for Wut were obtained 2 .For convenient comparison of abundance between two small RNA libraries, the raw reads of Huaq and Wut were normalized to the number of reads per 10 million (RP10M) (Table 2).For mapping to the genome, the sequence length of reads shorter than 17 nt were ltered out.The quali ed reads of small RNAs were nally aligned with the Chinese cabbage reference genome (http://brassicadb.org)by using the SOAP2 software, set one mismatch as threshold 2 .After analysis, 63.07% of Huaq and 61.88% of Wut within the raw reads were mapped to B. rapa genome (Wang et al., 2011b) (Table 2).As a result, 6.30 million and 6.18 million small RNA reads in Huaq and Wut were mapped on B. rapa reference genome, meanings that Huaq plants produce more small RNAs than Wut plants.Among them, there were about 2.10 million unique reads for Huaq, and 1.73 million for Wut.Interestingly, 1.42 million reads were matched on sense or antisense strands of annotated genes in Huaq, while 1.50 million reads on the genes in Wut.Between two varieties, total numbers of miRNAs, ta-siRNAs, and ra-siRNAs were different (Table 2).Huaq and Wut shared 497,645 unique reads, indicating the distinct difference on accumulation of common small RNAs between Huaq and Wut.
In total, 75 miRNA families were detected in Huaq and Wut 2 .However, 0.45 million miRNA reads were derived from Huaq and 0.30 million from Wut, revealing that more miRNAs were generated from Huaq than Wut.More than 80 putative ta-siRNAs were generated from TAS3 precursor in both varieties.Additionally, about 20000 small RNAs belonged to ra-siRNAs.
The expression of a number of small RNAs was different between Huaq and Wut.7725 reads were Huaqspeci c since they were not detected in Wut; and 4767 reads were Wut-speci c since they were not detected in Huaq (Table 3).Meanwhile, 15805 small RNAs were Huaq-predominant (P < 0.01) since they were more abundant in Huaq than in Wut; and 18048 small RNAs were Wut-predominant since they were more abundant in Wut than in Huaq.In addition, numbers of gene-derived small RNAs and ra-siRNAs were signi cantly different.Although the total numbers of unique reads were approximately the same between Huaq and Wut, the length of small RNA reads were substantially different (Fig. 2).The peaks of length for unique reads were the same at 24 nt in both varieties.However, the abundance of small RNAs with 21 ~ 24 nt in Huaq was more than that of Wut.In contrast, the abundance of 25 ~ 32 nt small RNAs in Wut was less than that of Huaq.
Many small RNA reads were perfectly matched with chloroplast genome.They were designated as chloroplast-related small RNAs.Among them, some were matched with nuclear sequences as well, and thus ltered out; and the rest were chloroplast-speci c small RNAs (csRNA).The situation was the same for mitochondria-speci c small RNAs (mitRNA).To detect the proportion of csRNA and mitRNA, small RNA reads were aligned to Arabidopsis Tair 9 chloroplast genome, mitochondria genome, and cDNA database (Fig. 3).The result showed that proportion of csRNAs in Wut seedlings was more than in Huaq seedlings in terms of number of unique reads and abundance of reads.

Small RNA distribution in chromosome
We calculated the distribution of small RNA reads in chromosome scale at 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 distributed similarly.In contrast, the hot sites of small RNA transcription between Arabidopsis species were very diverse, especially between two most related species Arabidopsis lytrata and Arabidopsis thaliana (Hollister et al., 2011).One exception was that at the end of chromosome 9th of Wut, there was a hot spot.

Differential expression of individual miRNAs between Huaq and Wut
Recently, many species of higher plants have been identi ed to have species-speci c miRNAs (Axtell and Bowman, 2008).To verify whether some miRNAs existed in one variety but not in another, we compared the sequences of miRNAs in Huaq and Wut seedlings.Among 78 miRNAs, 3 miRNAs were Wut-speci c (Supplementary Table 1).They were bra-miR403, bra-miR5711, bra-miR399f.
ta-siRNA in Huaq and Wut seedlings Alignment of small RNAs with ta-siRNAs of Arabidopsis indicated that only ta-siRNAs from B. rapa TAS3 (BrTAS3) transcripts were presented in Huaq and Wut datasets.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).We choose BrTAS3 precursors as reference to align all small RNAs.More than 80 small RNAs were matched with TAS3, but only one of them was conserved with TAS3a D7(+) (5'-TTCTTGACCTTGTAAGACCCC-3'), the known ta-siRNA in Arabidopsis.This ta-siRNA was located on the 5th chromosome of the genome of B. rapa.It was not known whether the other small RNAs derived from TAS3 precursors acted as ta-siRNAs in B. rapa.
Small RNAs aligned to the annotated genes of B. rapa For analysis of function of small RNAs, we aligned small RNA reads with the annotated genes of B. rapa.
More than 0.44 million small RNAs in Huaq and 0.41 million in Wut were matched on the sequences of the annotated genes 2 (Table 2).The half of them was located on antisense strands of the annotated genes (Fig. 6).They included 44 Huaq-speci c small RNAs and 27 Wut-speci c small RNAs.We set the threshold of 2-fold changes as standards of difference in abundance.As a result, 238 small RNAs in Huaq showed higher abundance than in Wut, while 77 small RNAs in Wut displayed higher abundance than in Huaq (Fig. 7A).
To detect the putative target genes of these small RNAs, the small RNAs were applied to match with their possible target genes of Arabidopsis and B. rapa.Potential targets (homologous to Arabidopsis genes) were displayed in Supplementary Table 2.The putative targets included hormone-related genes such as CTR1, AREB3 and ATCHX4 (Fig. 7B).The abundance of the small RNAs matched on CTR1, AREB3 genes were higher in Huaq, while that matched on ATCHX4 was higher in Wut.ATIPT3 regulates cytokinin biosynthesis.Small RNAs matched with ATIPT3 was more abundant in Huaq than in Wut.The small RNAs matched on ATCOAA, MEE39, and ATCAP-E1 were also more abundant in Huaq while those matched with ATCS-C (Wut-enriched), VHA-A, and ATNTH2 were more abundant in Wut than in Huaq.
These results showed that many small RNAs that may target at the important genes were different in abundance between Huaq and Wut.

Discussion
Comparative pro ling shows remarkable difference in types and abundance of small RNAs between different varieties of Chinese cabbage The whole view of small RNAs from small RNA deep sequencing provides information of non-coding RNAs in Chinese cabbage.These small RNAs are derived from nuclear, chloroplast and mitochondria genomes.In Huaq seedlings, the total numbers of 21-24 nt small RNAs are more abundant than in Wut.High amount of small RNAs in Huaq plants is mainly because of more types of small RNAs (unique small RNAs).Meanwhile, the gene-derived small RNAs in Huaq plants are more than those of Wut.The genederived small RNAs are mainly composed of small RNAs generated on antisense strands of the annotated genes.In the other words, the natural antisense transcripts (NATs) in Huaq plants generate more nat-siRNAs than in Wut plants.nat-siRNA plays an important role in plant development, stress response and disease resistance at transcriptional and post-transcriptional levels (Liu et al., 2008;Swiezewski et al., 2007).More nat-siRNAs in Wut plants indicate that development of Wut plants are vulnerable to stress and disease.Among 78 miRNAs, 3 miRNAs were Wut-speci c. 76.9% miRNAs were expressed predominantly in Huaq seedlings, and 23.1% miRNAs predominant in Wut seedlings.
Chloroplasts are important for photosynthesis and respiration.Abundance comparison of csRNA and rRNA reveals that small RNAs are more abundant in Wut seedlings than in Huaq seedlings, consistent with dark-green leaves of Wut seedlings.Ribosome biogenesis is a complex process comprising transcription, modi cation, and processing of ribosomal RNA and production of ribosomal proteins (Deisenroth and Zhang, 2010).Moreover, abundant rRNA-derived chloroplast small RNAs in Wut seedlings may re ect strong activity of photosynthesis, a process in which small RNAs are involved.
Genome-wide survey of small RNAs provides a basis for study of molecular mechanism governing leaf traits of Chinese cabbage 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.In Arabidopsis, several conserved miRNAs regulate shape, size, color and curvature of leaves.miR156 targets at 11 of 17 SPL genes (Axtell and Bowman, 2008).Overexpression of SPL3, SPL4 and SPL5 result in adult characters of leaf shape on juvenile leaves and thus promote vegetative phase change and oral transition.Overexpression of miR156 results in prolonged expression of juvenile characters due to extended period of juvenile leaf shape (Chuck et al., 2007).In Huaq seedlings, miR156 accumulated more than in Wut seedlings.Coincidently, Huaq plants maintains longer period of juvenile leaf shape than Wut plants.miR172 accumulates more in Wut seedlings than Huaq.In Arabidopsis, miR172 is antagonist of miR156 and promotes owering by repressing TOE1 and TOE2 (Glazinska et al., 2009), while the maize Cg1 mutant that overexpresses miR156 has reduced miR172 levels (Chuck et al., 2007).Higher level of miR172 in Wut may affect leaf traits by repression of miR156.The miR165/166 is conserved among seed plants and negatively regulates accumulation of members of the transcription factors HD-ZIP III families involved in meristem function, vascular patterning and adaxial cell fate (Wenkel et al., 2007;Zhong and Ye, 2007).Accumulation of miR165/166 on the abaxial side of leaf primordia delimits the HD-ZIP III transcripts in adaxial side of leaf and affects the polarity of the leaf (Juarez et al., 2004;Kidner and Martienssen, 2004).Overexpression of miR165/166 will conduct abaxial character of leaves.Huaq seedlings generate more miR166 than Wut seedlings.This difference could be re ected with leaf curvature.In Wut seedlings, miR319 accumulated less than in Huaq seedlings.It is interesting to determine whether this difference is linked to the phenotype of leaf margins of Wut.In Arabidopsis, miR319-mediated targets regulate leaf morphogenesis (Palatnik et al., 2003).
The correlation between miRNA expression and leaf traits in Chinese cabbage indicate their possible relation.It remains unknown whether these miRNAs are causally related to the certain leaf traits of Chinese cabbage.The functions of these small RNAs should be identi ed with the loss-of-function mutants and the transgenic plants, especially in the backgrounds of Chinese cabbage.Biogenesis of small RNAs must be rigorously regulated, and many small RNAs are expressed in a spatiotemporally regulated manner or in response to environmental stimuli (Chen, 2009).Total numbers of unique small RNAs in Huaq were much more than those of Wut, whereas total numbers of chloroplast-speci c small RNAs derived from chloroplast genome of Wut were more than those of Huaq.Possibly, the genes responsible for a subset of small siRNAs or miRNAs are different in sequences and function, and thus exhibit different ability to direct biogenesis of small RNAs.For the purpose to investigate miRNA-mediated leaf development in Chinese cabbage, an attempt has been made in our lab to verify the sequences and functions of the genes responsible for biogenesis of small RNAs to analyze the morphological consequences of several miRNAs differentially expressed in Huaq and Wut.These small RNAs are useful for understanding molecular mechanism of leaf development and stress resistance and for genetic improvement of major agronomic traits.
Small RNAs are useful for developing of small RNA makers for leaf traits desirable for high yield and quality The change in sequences and abundance of small RNAs may affect leaf development by altering expression of targeted genes.Among Huaq-or Wut-speci c small RNAs, there are a number of single nucleotide polymorphisms (SNPs).If some small RNAs are functional in leaf development, the SNPs may be linked with certain leaf traits.Then, the small RNAs with these SNPs may become the molecular makers for leaf shape, size, color or curvature.To nd out small RNAs markers, we have started to establish the relationship between SNPs in small RNAs and some quantitative trait loci (QTLs) responsible for leaf shape, size, color or curvature of Chinese cabbage.On the other hand, abundance of some small RNAs may be linked with certain leaf traits.Huaq has stronger ability to generate classic small RNAs than Wut.These small RNAs included those not only generated from the region corresponding to the known small RNAs but also from the externally and internally transcribed regions and space regions.Small RNAs are considered to be sequence-speci c, depending on nucleotides, such as DNA and histone methylation and posttranscriptional regulation at RNA levels (Chen, 2009).
Functional small RNAs are capable of affecting leaf phenotypes.It is possible for us to disclose the linkage relationship between abundance of small RNAs and leaf traits.Some functional small RNAs could be used as expression markers for leaf traits of Chinese cabbage.The yield and quality of Chinese cabbage are highly dependent on leaf shape, size, color and curvature.The vegetable breeders usually design these leaf traits according to the consumer's requirement but seldom success in breeding of the desirable cultivars.One explanation for this di culty is that the molecular makers for leaf shape, size, color and curvature are rare.The SNPs and abundance of small RNAs in Chinese cabbage should be a basis for exploring molecular markers of leaf variance.

Figures
Figure 1 The    The distribution and abundance of small RNAs on the genome of B.rapa.A hot spot appears at the end of 9th chromosome in Wut but not in Huaq.
Figure 6 The distribution of small RNAs mapped to the annotated genes of Huaq and Wut.(A) The distribution of small RNAs derived from sense strands of the annotated genes.(B) The distribution of small RNAs derived from the antisense strands of the annotated genes.The draft was draw by logarithmic value of abundance.
Figure 7 The abundance of small RNAs located in antisense strands of the annotated genes.(A) The number of the small RNAs located in the antisense strands.(B) The abundance of small RNAs located in antisense strands of Brassica genes (at the right), corresponding to orthologous genes of Arabidopsis (at the left).
seedlings of Huaq and Wut in the eld.(A) Huaq seedling.(B) The leaves of Huaq seedling in order (from inside to outside).(C) Wut seedling.(D) The leaves of Wut seedling in order (from inside to outside).Page 16/21

Table 2 .
Classi cation of small RNAs of unique reads and total reads in non-heading Chinese cabbage that mapped to the reference genome of B. rapa."shared" indicates common small RNAs between Huaq

Table 3 .
Classi cation of different expressional small RNAs in non-heading Chinese cabbage."speci c" indictates the small RNAs speci cally expressed in one variety."predominant" means that the abundance of the small RNAs in one variety are two times more than in another variety."Gene region derived" indicates the small RNAs matched on sense or antisense strand of the annotated coding genes.