BACKGROUND
the tumourous stem mustard, brassica juncea , is an important vegetable and the raw material for fuling mustard.
the genetics breeding, physiology, biochemistry and classification of mustards have been extensively studied, but little work has been done at the molecular level.
the three diploid  and three allotetraploid  species of brassica, related according to the u-triangle theory, are one of the best model systems to study polyploidy and genetic relationships  <cit> .
thus, investigating tumourous stem mustard at the molecular level may help to clarify the genetic relationships mentioned above.
development of the stem in tumourous stem mustard, which is pickled or eaten fresh, is directly related to the quality and yield of tumourous mustards.
stem swelling is a key characteristic of tumourous stem formation, and although many new cultivars of tumourous stem mustard have been bred, the molecular mechanism underlying stem swelling is unclear.
tumourous stem mustard is an annual plant, and the stem does not swell except in plants sown between mid-september and mid-october in chongqing and the other valleys of the yangtze river, china; thus, the production period of edible stems is limited  <cit> .
identification of the genes controlling stem swelling and the regulatory network would facilitate molecular breeding and increase the yield and quality of this crop.
tumourous stem mustard shares a close genetic relationship with the model plants arabidopsis thaliana and arabidopsis lyrata, whose genetic backgrounds are more clearly understood.
thus, a. thaliana may aid in our understanding of the mechanism of swelling in tumourous stem mustard.
in recent years, novel, high-throughput, deep sequencing technologies have allowed the efficient generation of large-scale ests and improved the speed of gene discovery  <cit> .
in addition, a. thaliana microarrays may be used to study stem swelling because of the close genetic relationship between that organism and tumourous stem mustard.
rna sequencing  generates absolute rather than relative gene expression measurements and provides greater insight and accuracy than microarrays  <cit> .
in this study, we generated over  <dig> billion bases of high-quality dna sequence using illumina technology and demonstrated the suitability of short-read sequencing for the de novo assembly and annotation of genes expressed without prior genome information.
in addition, we constructed five digital gene expression  libraries and compared the gene expression profiles of tumourous stem mustard at different developmental stages with those of different varieties.
the assembled annotated transcriptome sequences and gene expression profiles will facilitate the identification of genes involved in tumourous stem mustard swelling and be a useful reference for other cruciferae developmental studies.
methods
plant material
the tumourous stem mustard  strains yong'an  and a mutant variety dayejie  were collected from the chongqing fuling institute of agricultural sciences, chongqing, china.
after seeding on october 5th, dayejie and four developmental stages of yong'an were collected in february of the next year.
dayejie stems were collected  <dig> weeks after seeding .
the stems of yong'an were collected  <dig>   <dig>   <dig>  and  <dig> weeks after seeding .
after washing in physiological saline and  <dig> % depc-treated water, the fresh samples were stored in liquid nitrogen until the extraction of total rna.
cdna preparation and illumina sequencing
total rna was extracted using a plant rna kit according to the manufacturer's instructions .
the quantity of rna was verified by an ultraviolet spectrometer and electrophoresis on a denaturing formaldehyde agarose gel.
to obtain complete gene expression information, a pooled rna sample from tissues of different developmental stages  was used for transcriptome analysis.
oligo beads were used to isolate poly + mrna from the total rna.
fragmentation buffer was added to disrupt the mrna into short fragments.
taking these short fragments as templates, a random hexamer primer was used to synthesise first-strand cdna.
second-strand cdna was synthesised using buffer, dntps, rnase, and dna polymerase i. the resulting short fragments were purified with a qiaquick pcr extraction kit and resolved with eb buffer for end repair and addition of a poly tail.
afterwards, the short fragments were connected with sequencing adapters.
following agarose gel electrophoresis, suitable fragments were selected as templates for pcr.
finally, the library was sequenced using illumina hiseq™  <dig>
transcriptome analysis
the raw data from the images were collected using solexa ga pipeline  <dig>  by removing low-quality reads , adaptor sequence fragments, and empty reads.
next, de novo assembly of the transcriptome into unigenes was carried out with soapdenovo, a short-reads assembly program  <cit> .
briefly, soapdenovo first combines reads with a particular overlap to form longer fragments without n, which are called contigs.
then, the reads are mapped back to contigs.
the program is able to detect contigs from the same transcript as well as the sequences between these contigs using paired-end reads.
next, soapdenovo connects the contigs using n to represent unknown sequences between each two contigs, and scaffolds are made.
paired-end reads are used again for gap filling of the scaffolds to obtain sequences with the lowest ns and which cannot be extended on either end.
such sequences are defined as unigenes.
subsequently, a blastx alignment  was performed between the unigenes and protein databases, including the non-redundant , swiss-prot, kyoto encyclopaedia of genes and genomes , and cog databases, and the best alignments were used to decide the sequence direction of the unigenes.
if the results from the different databases conflicted, a priority order of nr, swiss-prot, kegg, and cog was followed to decide the sequence direction of the unigenes.
when a unigene happened to be unaligned in any of the above databases, estscan was used to predict its coding regions and decide its sequence direction  <cit> .
in the final step, using nr annotation, the blast2go program was used to obtain the gene ontology  and kegg annotations of the unigenes  <cit> .
after obtaining the go annotation for each unigene, wego software was used to classify the unigenes by function and to determine the distribution of gene functions in the species at the macro level  <cit> .
all raw transcriptome data have been deposited at the sequence read archive  of ncbi.
dge library preparation and sequencing
total rna was extracted from the stems of yongan1hao, yongan2hao, yongan3hao, yongan4hao, and daye3bianzhong using a plant rna kit  according to the manufacturer's instructions.
next, a dge library was prepared using an illumina gene expression sample prep kit.
briefly, total rna from the samples was used for mrna capture with magnetic oligo beads.
double-stranded cdnas were synthesised directly on the poly + rna-bound beads and then digested with nlaiii.
those cdna fragments with 3' ends were purified from the magnetic beads, and illumina adaptor  <dig> was added to their 5' ends.
after digestion with mmei, which recognises the junction between illumina adaptor  <dig> and the sequence catg, 21-bp tags containing adaptor  <dig> were ligated to the 3' ends of the tags to create a tag library.
the library was amplified by pcr over  <dig> cycles, and 85-bp strips were purified by page.
the single-stranded molecules were attached to the illumina chip for sequencing.
all raw tag data were deposited at the sra of ncbi.
dge analysis
the raw image data were transformed by base calling into sequence data.
to map the dge tags, the sequenced raw data were filtered to remove adaptor sequences, low-quality sequences , empty tags , and tags with only one copy number .
for tag annotation, clean tags containing catg and the 21-bp tag sequences were mapped to our transcriptome reference database, allowing no more than one nucleotide mismatch.
the clean tags were designated as unambiguous clean tags.
for gene expression analysis, the number of expressed tags was calculated and normalised to the number of transcripts per million tags.
to compare the differences in gene expression, the tag frequency in each dge library was statistically analysed according to the method of audic and claverie  <cit> .
we used a false discovery rate of <  <dig>  and an absolute value of the log  <dig> ratio of >  <dig> as the threshold for judging the significance of the gene expression differences.
next, the differentially expressed genes were subjected to go and kegg ontology  enrichment analysis.
enriched p-values were calculated according to the hypergeometric test:
p=1-∑i=0m-1min-mn-inn
in this equation, n represents the number of genes with a go/ko annotation, n represents the number of differentially expressed genes in n, m represents the number of genes in each go/ko term, and m represents the number of differentially expressed genes in each go/ko term.
for go enrichment analysis, all of the p-values were subjected to bonferroni correction.
a corrected p-value of <  <dig>  was selected as the threshold for determining significant enrichment of the gene sets.
quantitative real-time pcr  validation
total rna was extracted as described for dge library preparation and sequencing.
total rna  from each sample was reverse-transcribed in a 10-μl reaction using the amv rna pcr kit  <dig>  .
the sequences of the primers used are given in additional file 1: table s <dig>  the  <dig> s rrna gene of tumourous stem mustard was used as an internal control.
qrt-pcr was performed using sybr premix ex taq™ kit  according to the manufacturer's protocol.
the selected genes were verified using a bio-rad iq <dig> real-time pcr detection system with a cycling temperature of 57°c and with a single peak on the melting curve to ensure a single product.
at least three replicates were tested per sample.
RESULTS
illumina sequencing and sequence assembly
a total of  <dig> , <dig> reads  were generated through illumina sequencing and assembled into  <dig>  contigs.
using paired end joining and gap filling, the contigs were further assembled into  <dig>  scaffolds with a mean length of  <dig> bp.
after clustering the scaffolds together with the nucleotide sequences available from the ncbi, we obtained  <dig>  unigenes, including  <dig>  clusters and  <dig>  singletons, with a mean length of  <dig> bp.
the sequence clusters had two to seven scaffolds per cluster, and 92% of the clusters contained only one scaffold.
the size distribution indicated that the lengths of the  <dig>  unigenes exceeded  <dig>  bp .
to evaluate the quality of the dataset, the ratio of the gap length to the length of the assembled unigenes was analysed .
most of the unigenes showed gap lengths < 5% of the total length, which accounted for  <dig> % of the total number of unigenes .
annotation of predicted proteins
to annotate the unigenes, we first searched the reference sequences using blastx against the nr ncbi protein database with a cut-off e-value of 10- <dig>  a total of  <dig>  unigenes  returned a significant blast result .
the species distribution of the best match result for each sequence is shown in figure  <dig>  the sequences had a  <dig> % match with a. thaliana, followed by arabidopsis lyrata subsp.
lyrata , brassica , and oryza sativa .
cog classification
the assembled unigenes were compared against the cog database for the analysis of phylogenetically widespread domain families.
the results revealed  <dig>  unigenes with significant homology and assigned them to the appropriate cog clusters.
these cog classifications were grouped into  <dig> functional categories .
the five largest categories were 'general function' ; 'replication, recombination, and repair' ; 'transcription' ; 'translation, ribosomal structure, and biogenesis' ; and 'carbohydrate transport and metabolism' .
go classification
go assignments were used to classify the functions of the predicted tumourous stem mustard genes.
based on sequence homology,  <dig>  sequences were categorised into  <dig> functional groups .
in each of the three main categories  of the go classification, the major subcategories were as follows: six subcategories for biological process ; four subcategories for cellular component ; and two subcategories for molecular function .
only a few genes were clustered in terms of 'biological adhesion', 'cell killing', 'locomotion', viral reproduction', 'virion', 'virion part', 'auxiliary transport protein activity' and 'electron carrier activity'.
figure  <dig> shows the go classification of the tumourous stem mustard transcriptome.
functional classification using the kegg database
functional classification and pathway assignment was performed using the kegg database  <cit> .
first, the  <dig>  unigenes with an e-value ≤ 1e- <dig> were compared using blastx against the kegg database.
in total,  <dig>  unigenes were assigned to  <dig> kegg pathways .
the major pathways were 'metabolic pathways ', 'biosynthesis of secondary metabolites ', 'plant-pathogen interaction ', 'spliceosome ' and 'starch and sucrose metabolism '; the gene numbers and percentages assigned to these pathways were  <dig>  ,  <dig>  ,  <dig>  ,  <dig>   and  <dig>  , respectively.
dge library sequencing
five dge libraries for tumourous stem mustard were sequenced , which generated approximately  <dig> to  <dig> million high-quality reads for each library .
the percentage of clean reads among the raw reads in each library was > 97% .
among the clean reads, the number of sequences that could be mapped to unigenes ranged from  <dig>  to  <dig>  million, and the percentage of cleans reads ranged from  <dig>  to  <dig> % in the five libraries.
as table  <dig> shows, the vast majority of these mapped reads were uniquely matched to unigenes, and the percentage of multi-position matched reads was no more than  <dig> %.
note: the percentages of reads containing n, adaptors, low quality, clean reads.
the numbers in parentheses indicate the percentage of each type of read among the total raw reads
gene expression variations among the different samples
first, to evaluate the dge data, we analysed the distribution of unigene coverage in each sample, which is the number of clean reads that aligned to the reference unigenes.
as shown in figure  <dig>  most unigene coverage was > 50% .
second, the number of clean reads was calculated and the gene expression level was calculated using the reads per kb per million reads method for each unigene   <cit> .
third, the differentially expressed genes were identified using an algorithm developed by audic et al.
<cit>  at different developmental stages and in various samples.
variations in gene expression were identified based on comparisons of daye3bianzhong and yongan2hao to yongan4hao, and yongan1hao and yongan2hao to yongan4hao.
the results showed significantly differentially expressed genes in the above comparisons:  between daya3bianzhong and yongan2hao,  <dig>  and  <dig>  genes were up- and downregulated, respectively;  between daya3bianzhong and yongan3hao,  <dig>  and  <dig>  genes were up- and downregulated, respectively;  between daya3bianzhong and yongan4hao,  <dig>  and  <dig>  genes were up- and downregulated, respectively;  between yongan1hao and yongan2hao,  <dig>  and  <dig>  genes were up- and downregulated, respectively;  between yongan1hao and yongan3hao,  <dig>  and  <dig>  genes were up- and downregulated, respectively; and  between yongan1hao and yongan4hao,  <dig>  and  <dig>  genes were up- and downregulated, respectively .
the  <dig>  genes occurring simultaneously in the six comparisons above  were screened out, and the comparison ratio values of these genes were used to build a cluster tree.
as figure  <dig> shows, the genes were divided into four groups based on whether they showed a trend toward up-  or downregulation .
to investigate some of the most differentially up- and downregulated genes , eight up-regulated and five downregulated genes clustered and were screened out.
of these thirteen genes, only one downregulated gene had a defined function , wherea five upregulated genes had defined functions .
a total of seven genes among the thirteen have unknown functions .
next, the  <dig>  screened genes were analysed for function using the go annotation system.
the major subcategories were as follows: three subcategories  in the cellular component cluster; two subcategories  in the molecular function cluster; and three subcategories  in the biological process cluster .
the ko categories were further analysed, and the results indicated that the major pathways were 'metabolic pathways ', 'biosynthesis of secondary metabolites ', 'plant hormone signal transduction ' and 'microbial metabolism in diverse environments ' .
validation of rna-seq-based gene expression
to validate the expression profiles obtained by rna-seq, real-time rt-pcr was performed on eight genes selected at random with high or low expression levels.
expression comparisons were performed between yongan2hao and yongan1hao, yongan3hao and yongan1hao, and yongan4hao and yongan1hao by qrt-pcr.
for all of the genes, the trend in real-time rt-pcr expression was in agreement with the rna-seq data except for prh <dig> .
discussion
tumourous stem mustard, an important cash crop and the raw material for fuling pickles, is a world-famous vegetable crop.
presently, many varieties of tumourous stem mustard have been bred, but the regulatory pathway and molecular mechanism of mustard tumourous stem development are unclear, and molecular biological research on tumourous stem mustard is rare.
by transcriptome sequence analysis, we obtained  <dig> , <dig> reads corresponding to about  <dig> gb of raw sequence data.
the predicted  <dig>  unigenes were subjected to blast annotation, and 72% of the unigenes returned a significant blast result.
as expected, most unigenes shared the highest sequence similarity with crucifers .
the number of genes similar to brassica genes was lower than the number of similar genes between a. thaliana and a. lyrata subsp.
lyrata , possibly due to the characteristics of the b. juncea transcriptome in the ncbi protein database.
arabidopsis thaliana is an important model plant with a clear genetic background that is very useful for researching gene functions in tumourous stem mustard.
our transcriptome analysis is the first high-throughput sequencing of tumourous stem mustard and will serve as a basis for other studies.
to investigate the regulatory pathway and molecular mechanism of tumour swelling, we created five dge libraries from plants at different developmental stages and samples from a non-swollen mutant to analyse the gene expression patterns at various developmental stages.
the quality of the dge libraries was further confirmed by qrt-pcr analysis.
because the tumours continue to swell for a long period after initiation, different swelling stages were selected for each experimental group, and non-swollen stages of yong'an and the non-swollen mutant strain dayejie were used for comparison.
the dge profiles of the swelling stages were compared with the controls.
different genes may be involved in tumourous stem formation, and the common genes in the six comparison groups reduced the number of differentially expressed genes that might be related to tumourous stem development.
compared with the go annotation results of the dge screening of genes in the transcriptome data, we found no genes distributed in the three molecular function subcategories  in the dge group, indicating that these three subcategories are not related to tumour swelling.
a comparison of the results of additional file 3: table s <dig> and table s <dig> showed that the pathway order number of 'plant hormone signal transduction' was about  <dig> , suggesting that these genes were screened out and that the pathway 'plant hormone signal transduction' is related to tumour swelling.
mapping the dge data back to the transcriptome database revealed that about 30% of the reads were mapped and that > 60% remained transcribed sequences.
although large amounts of data were obtained by transcriptome sequencing, the reference sequences may still be insufficient and may have caused the lower mapped ratio, which could be resolved by increasing the sequencing depth and enhancing the accuracy of the assembly.
although  <dig>  differentially expressed genes were discovered using the above method, the key genes related to tumourous stem formation need to be analysed further.
the genes with the greatest changes in expression were selected for further study; a log  <dig> ratio value >  <dig> was used as a threshold to select thirteen genes for further analysis.
seven of the genes have unknown functions and six genes have a functional annotation based on sequence similarity.
of the six annotated genes, four gene functions require further clarification, and two genes, phytochrome kinase substrate  <dig>   and abr <dig>  , have functions whose details are relatively well known.
girdhar et al.
<cit>  determined that abr <dig> functions in the negative regulation of aba responses during seed germination and plays a role in the aba signalling pathway in arabidopsis.
in this study, the expression level of abr <dig> was significantly higher in the inflation stage of the stem tumour.
this suggests that aba signalling is related to stem inflation.
as additional file 6: table s <dig> shows, four genes may be involved in the 'plant hormone signal transduction ' pathway, suggesting that plant hormones are related to stem swelling; however, abr <dig> was not one of the four genes annotated in the 'plant hormone signal transduction' pathway in this study.
during the inflation stage of stem tumours, plant hormones play key roles in the number of cells in the stem, which increase rapidly accompanied by cell splintering.
pks <dig> expression was upregulated during the tumour inflation stage.
to our knowledge, pks <dig> is phosphorylated in a phytochrome-dependent manner and negatively regulates phytochrome light signalling in arabidopsis  <cit> .
other reports have shown that pks <dig> regulates root phototropism and gravitropism and leaf flattening and positioning, and that it affects the state of phytochrome a in etiolated arabidopsis seedlings  <cit> .
chen et al.
<cit>  found that illumination and temperature were the main factors affecting the formation of stem tumours, but they did not determine whether pks <dig> plays key roles in illumination.
is pks <dig> related to stem inflation in tumourous stem mustard?
whether or not the overexpression or knockout/down of abr <dig> and pks <dig> alters the trend of inflation requires further research.
CONCLUSIONS
although the molecular functions of individual tumourous mustard genes and the associated signal transduction pathways remain largely unknown, the present transcriptome analysis provides valuable information regarding tumourous stem development, which may facilitate future investigations of the detailed regulatory mechanisms and pathways.
additionally, we obtained about  <dig> gb of raw sequence data and predicted  <dig>  assembled unigenes, and annotated  <dig>  unigenes using illumina sequencing technology.
this is the first large-scale transcriptome and expression study of the tumourous stem mustard, brassica juncea var.
tumida tsen et lee for which the data has been deposited in genbank .
we believe that our data not only provide more molecular information for researchers, but will also help to accelerate gene expression and function research in brassica juncea as well as that on the evolutionary relationships of the cruciferae.
authors' contributions
cyf conceived the study.
cyf, sq, zxy participated in experiment materials preparation.
zxy participated in rna extraction.
sq analyzed data and performed qrt-pcr materials.
cyf, sq wrote the paper.
zgf provided the plant varieties and participated in experiment materials preparation, fyh, lyh, hxh, sjj, jhz, hdw, pz, xlx and yjp participated in part of experiment.
all authors read and approved the final manuscript.
supplementary material
additional file 1
table s <dig> primers used in qrt-pcr for validating differentially expressed genes.
click here for file
additional file 2
table s <dig> top hits obtained by blastx for the unigenes.
blastx against the nr, swiss-prot, kegg and cog protein database was used with a cutoff e-value < 1e- <dig> .
click here for file
additional file 3
table s <dig> ko's assigned of unigenes.
click here for file
additional file 4
table s <dig> expression level of unigenes in each library .
click here for file
additional file 5
table s <dig> discription of thirteen screened out genes.
expression ratios absolute values >  <dig> in six comparison groups were used for rule.
click here for file
additional file 6
table s <dig> ko's assigned of screening  <dig> different genes.
click here for file
additional file 7
table s <dig> the unigenes' id of screened out genes.
click here for file
acknowledgements
we acknowledge the beijing genomics institute at shenzhen for its assistance in original data processing and related bioinformatics analysis.
this work was supported by the natural science foundation of china  and chongqing natural science foundation.

 the tumourous stem mustard, brassica juncea, is an important vegetable and the raw material for fuling mustard. the genetics breeding, physiology, biochemistry and classification of mustards have been extensively studied, but little work has been done at the molecular level. the three diploid and three allotetraploid species of brassica, related according to the u-triangle theory, are one of the best model systems to study polyploidy and genetic relationships. thus, investigating tumourous stem mustard at the molecular level may help to clarify the genetic relationships mentioned above. development of the stem in tumourous stem mustard, which is pickled or eaten fresh, is directly related to the quality and yield of tumourous mustards. stem swelling is a key characteristic of tumourous stem formation, and although many new cultivars of tumourous stem mustard have been bred, the molecular mechanism underlying stem swelling is unclear. tumourous stem mustard is an annual plant, and the stem does not swell except in plants sown between mid-september and mid-october in chongqing and the other valleys of the yangtze river, china; thus, the production period of edible stems is limited. identification of the genes controlling stem swelling and the regulatory network would facilitate molecular breeding and increase the yield and quality of this crop. tumourous stem mustard shares a close genetic relationship with the model plants arabidopsis thaliana and arabidopsis lyrata, whose genetic backgrounds are more clearly understood. thus, a. thaliana may aid in our understanding of the mechanism of swelling in tumourous stem mustard. in recent years, novel, high-throughput, deep sequencing technologies have allowed the efficient generation of large-scale ests and improved the speed of gene discovery. in addition, a. thaliana microarrays may be used to study stem swelling because of the close genetic relationship between that organism and tumourous stem mustard. rna sequencing generates absolute rather than relative gene expression measurements and provides greater insight and accuracy than microarrays. in this study, we generated over billion bases of high-quality dna sequence using illumina technology and demonstrated the suitability of short-read sequencing for the de novo assembly and annotation of genes expressed without prior genome information. in addition, we constructed five digital gene expression libraries and compared the gene expression profiles of tumourous stem mustard at different developmental stages with those of different varieties. the assembled annotated transcriptome sequences and gene expression profiles will facilitate the identification of genes involved in tumourous stem mustard swelling and be a useful reference for other cruciferae developmental studies. methods plant material the tumourous stem mustard strains yong'an and a mutant variety dayejie were collected from the chongqing fuling institute of agricultural sciences, chongqing, china. after seeding on october 5th, dayejie and four developmental stages of yong'an were collected in february of the next year. dayejie stems were collected weeks after seeding. the stems of yong'an were collected and weeks after seeding. after washing in physiological saline and % depc-treated water, the fresh samples were stored in liquid nitrogen until the extraction of total rna. cdna preparation and illumina sequencing total rna was extracted using a plant rna kit according to the manufacturer's instructions. the quantity of rna was verified by an ultraviolet spectrometer and electrophoresis on a denaturing formaldehyde agarose gel. to obtain complete gene expression information, a pooled rna sample from tissues of different developmental stages was used for transcriptome analysis. oligo beads were used to isolate poly + mrna from the total rna. fragmentation buffer was added to disrupt the mrna into short fragments. taking these short fragments as templates, a random hexamer primer was used to synthesise first-strand cdna. second-strand cdna was synthesised using buffer, dntps, rnase, and dna polymerase i. the resulting short fragments were purified with a qiaquick pcr extraction kit and resolved with eb buffer for end repair and addition of a poly tail. afterwards, the short fragments were connected with sequencing adapters. following agarose gel electrophoresis, suitable fragments were selected as templates for pcr. finally, the library was sequenced using illumina hiseq™ transcriptome analysis the raw data from the images were collected using solexa ga pipeline by removing low-quality reads, adaptor sequence fragments, and empty reads. next, de novo assembly of the transcriptome into unigenes was carried out with soapdenovo, a short-reads assembly program. briefly, soapdenovo first combines reads with a particular overlap to form longer fragments without n, which are called contigs. then, the reads are mapped back to contigs. the program is able to detect contigs from the same transcript as well as the sequences between these contigs using paired-end reads. next, soapdenovo connects the contigs using n to represent unknown sequences between each two contigs, and scaffolds are made. paired-end reads are used again for gap filling of the scaffolds to obtain sequences with the lowest ns and which cannot be extended on either end. such sequences are defined as unigenes. subsequently, a blastx alignment was performed between the unigenes and protein databases, including the non-redundant, swiss-prot, kyoto encyclopaedia of genes and genomes, and cog databases, and the best alignments were used to decide the sequence direction of the unigenes. if the results from the different databases conflicted, a priority order of nr, swiss-prot, kegg, and cog was followed to decide the sequence direction of the unigenes. when a unigene happened to be unaligned in any of the above databases, estscan was used to predict its coding regions and decide its sequence direction. in the final step, using nr annotation, the blast2go program was used to obtain the gene ontology and kegg annotations of the unigenes. after obtaining the go annotation for each unigene, wego software was used to classify the unigenes by function and to determine the distribution of gene functions in the species at the macro level. all raw transcriptome data have been deposited at the sequence read archive of ncbi. dge library preparation and sequencing total rna was extracted from the stems of yongan1hao, yongan2hao, yongan3hao, yongan4hao, and daye3bianzhong using a plant rna kit according to the manufacturer's instructions. next, a dge library was prepared using an illumina gene expression sample prep kit. briefly, total rna from the samples was used for mrna capture with magnetic oligo beads. double-stranded cdnas were synthesised directly on the poly + rna-bound beads and then digested with nlaiii. those cdna fragments with 3' ends were purified from the magnetic beads, and illumina adaptor was added to their 5' ends. after digestion with mmei, which recognises the junction between illumina adaptor and the sequence catg, 21-bp tags containing adaptor were ligated to the 3' ends of the tags to create a tag library. the library was amplified by pcr over cycles, and 85-bp strips were purified by page. the single-stranded molecules were attached to the illumina chip for sequencing. all raw tag data were deposited at the sra of ncbi. dge analysis the raw image data were transformed by base calling into sequence data. to map the dge tags, the sequenced raw data were filtered to remove adaptor sequences, low-quality sequences, empty tags, and tags with only one copy number. for tag annotation, clean tags containing catg and the 21-bp tag sequences were mapped to our transcriptome reference database, allowing no more than one nucleotide mismatch. the clean tags were designated as unambiguous clean tags. for gene expression analysis, the number of expressed tags was calculated and normalised to the number of transcripts per million tags. to compare the differences in gene expression, the tag frequency in each dge library was statistically analysed according to the method of audic and claverie. we used a false discovery rate of < and an absolute value of the log ratio of > as the threshold for judging the significance of the gene expression differences. next, the differentially expressed genes were subjected to go and kegg ontology enrichment analysis. enriched p-values were calculated according to the hypergeometric test: p=1-∑i=0m-1min-mn-inn in this equation, n represents the number of genes with a go/ko annotation, n represents the number of differentially expressed genes in n, m represents the number of genes in each go/ko term, and m represents the number of differentially expressed genes in each go/ko term. for go enrichment analysis, all of the p-values were subjected to bonferroni correction. a corrected p-value of < was selected as the threshold for determining significant enrichment of the gene sets. quantitative real-time pcr validation total rna was extracted as described for dge library preparation and sequencing. total rna from each sample was reverse-transcribed in a 10-μl reaction using the amv rna pcr kit. the sequences of the primers used are given in additional file 1: table s the s rrna gene of tumourous stem mustard was used as an internal control. qrt-pcr was performed using sybr premix ex taq™ kit according to the manufacturer's protocol. the selected genes were verified using a bio-rad iq real-time pcr detection system with a cycling temperature of 57°c and with a single peak on the melting curve to ensure a single product. at least three replicates were tested per sample.  illumina sequencing and sequence assembly a total of, reads were generated through illumina sequencing and assembled into contigs. using paired end joining and gap filling, the contigs were further assembled into scaffolds with a mean length of bp. after clustering the scaffolds together with the nucleotide sequences available from the ncbi, we obtained unigenes, including clusters and singletons, with a mean length of bp. the sequence clusters had two to seven scaffolds per cluster, and 92% of the clusters contained only one scaffold. the size distribution indicated that the lengths of the unigenes exceeded bp. to evaluate the quality of the dataset, the ratio of the gap length to the length of the assembled unigenes was analysed. most of the unigenes showed gap lengths < 5% of the total length, which accounted for % of the total number of unigenes. annotation of predicted proteins to annotate the unigenes, we first searched the reference sequences using blastx against the nr ncbi protein database with a cut-off e-value of 10- a total of unigenes returned a significant blast result. the species distribution of the best match result for each sequence is shown in figure the sequences had a % match with a. thaliana, followed by arabidopsis lyrata subsp. lyrata, brassica, and oryza sativa. cog classification the assembled unigenes were compared against the cog database for the analysis of phylogenetically widespread domain families. the results revealed unigenes with significant homology and assigned them to the appropriate cog clusters. these cog classifications were grouped into functional categories. the five largest categories were'general function';'replication, recombination, and repair';'transcription';'translation, ribosomal structure, and biogenesis'; and'carbohydrate transport and metabolism'. go classification go assignments were used to classify the functions of the predicted tumourous stem mustard genes. based on sequence homology, sequences were categorised into functional groups. in each of the three main categories of the go classification, the major subcategories were as follows: six subcategories for biological process; four subcategories for cellular component; and two subcategories for molecular function. only a few genes were clustered in terms of'biological adhesion','cell killing','locomotion', viral reproduction','virion','virion part','auxiliary transport protein activity' and'electron carrier activity'. figure shows the go classification of the tumourous stem mustard transcriptome. functional classification using the kegg database functional classification and pathway assignment was performed using the kegg database. first, the unigenes with an e-value ≤ 1e- were compared using blastx against the kegg database. in total, unigenes were assigned to kegg pathways. the major pathways were'metabolic pathways','biosynthesis of secondary metabolites','plant-pathogen interaction','spliceosome' and'starch and sucrose metabolism'; the gene numbers and percentages assigned to these pathways were,,, and, respectively. dge library sequencing five dge libraries for tumourous stem mustard were sequenced, which generated approximately to million high-quality reads for each library. the percentage of clean reads among the raw reads in each library was > 97%. among the clean reads, the number of sequences that could be mapped to unigenes ranged from to million, and the percentage of cleans reads ranged from to % in the five libraries. as table shows, the vast majority of these mapped reads were uniquely matched to unigenes, and the percentage of multi-position matched reads was no more than %. note: the percentages of reads containing n, adaptors, low quality, clean reads. the numbers in parentheses indicate the percentage of each type of read among the total raw reads gene expression variations among the different samples first, to evaluate the dge data, we analysed the distribution of unigene coverage in each sample, which is the number of clean reads that aligned to the reference unigenes. as shown in figure most unigene coverage was > 50%. second, the number of clean reads was calculated and the gene expression level was calculated using the reads per kb per million reads method for each unigene. third, the differentially expressed genes were identified using an algorithm developed by audic et al. at different developmental stages and in various samples. variations in gene expression were identified based on comparisons of daye3bianzhong and yongan2hao to yongan4hao, and yongan1hao and yongan2hao to yongan4hao. the results showed significantly differentially expressed genes in the above comparisons: between daya3bianzhong and yongan2hao, and genes were up- and downregulated, respectively; between daya3bianzhong and yongan3hao, and genes were up- and downregulated, respectively; between daya3bianzhong and yongan4hao, and genes were up- and downregulated, respectively; between yongan1hao and yongan2hao, and genes were up- and downregulated, respectively; between yongan1hao and yongan3hao, and genes were up- and downregulated, respectively; and between yongan1hao and yongan4hao, and genes were up- and downregulated, respectively. the genes occurring simultaneously in the six comparisons above were screened out, and the comparison ratio values of these genes were used to build a cluster tree. as figure shows, the genes were divided into four groups based on whether they showed a trend toward up- or downregulation. to investigate some of the most differentially up- and downregulated genes, eight up-regulated and five downregulated genes clustered and were screened out. of these thirteen genes, only one downregulated gene had a defined function, wherea five upregulated genes had defined functions. a total of seven genes among the thirteen have unknown functions. next, the screened genes were analysed for function using the go annotation system. the major subcategories were as follows: three subcategories in the cellular component cluster; two subcategories in the molecular function cluster; and three subcategories in the biological process cluster. the ko categories were further analysed, and the results indicated that the major pathways were'metabolic pathways','biosynthesis of secondary metabolites','plant hormone signal transduction' and'microbial metabolism in diverse environments'. validation of rna-seq-based gene expression to validate the expression profiles obtained by rna-seq, real-time rt-pcr was performed on eight genes selected at random with high or low expression levels. expression comparisons were performed between yongan2hao and yongan1hao, yongan3hao and yongan1hao, and yongan4hao and yongan1hao by qrt-pcr. for all of the genes, the trend in real-time rt-pcr expression was in agreement with the rna-seq data except for prh. discussion tumourous stem mustard, an important cash crop and the raw material for fuling pickles, is a world-famous vegetable crop. presently, many varieties of tumourous stem mustard have been bred, but the regulatory pathway and molecular mechanism of mustard tumourous stem development are unclear, and molecular biological research on tumourous stem mustard is rare. by transcriptome sequence analysis, we obtained, reads corresponding to about gb of raw sequence data. the predicted unigenes were subjected to blast annotation, and 72% of the unigenes returned a significant blast result. as expected, most unigenes shared the highest sequence similarity with crucifers. the number of genes similar to brassica genes was lower than the number of similar genes between a. thaliana and a. lyrata subsp. lyrata, possibly due to the characteristics of the b. juncea transcriptome in the ncbi protein database. arabidopsis thaliana is an important model plant with a clear genetic background that is very useful for researching gene functions in tumourous stem mustard. our transcriptome analysis is the first high-throughput sequencing of tumourous stem mustard and will serve as a basis for other studies. to investigate the regulatory pathway and molecular mechanism of tumour swelling, we created five dge libraries from plants at different developmental stages and samples from a non-swollen mutant to analyse the gene expression patterns at various developmental stages. the quality of the dge libraries was further confirmed by qrt-pcr analysis. because the tumours continue to swell for a long period after initiation, different swelling stages were selected for each experimental group, and non-swollen stages of yong'an and the non-swollen mutant strain dayejie were used for comparison. the dge profiles of the swelling stages were compared with the controls. different genes may be involved in tumourous stem formation, and the common genes in the six comparison groups reduced the number of differentially expressed genes that might be related to tumourous stem development. compared with the go annotation results of the dge screening of genes in the transcriptome data, we found no genes distributed in the three molecular function subcategories in the dge group, indicating that these three subcategories are not related to tumour swelling. a comparison of the results of additional file 3: table s and table s showed that the pathway order number of'plant hormone signal transduction' was about, suggesting that these genes were screened out and that the pathway'plant hormone signal transduction' is related to tumour swelling. mapping the dge data back to the transcriptome database revealed that about 30% of the reads were mapped and that > 60% remained transcribed sequences. although large amounts of data were obtained by transcriptome sequencing, the reference sequences may still be insufficient and may have caused the lower mapped ratio, which could be resolved by increasing the sequencing depth and enhancing the accuracy of the assembly. although differentially expressed genes were discovered using the above method, the key genes related to tumourous stem formation need to be analysed further. the genes with the greatest changes in expression were selected for further study; a log ratio value > was used as a threshold to select thirteen genes for further analysis. seven of the genes have unknown functions and six genes have a functional annotation based on sequence similarity. of the six annotated genes, four gene functions require further clarification, and two genes, phytochrome kinase substrate and abr, have functions whose details are relatively well known. girdhar et al. determined that abr functions in the negative regulation of aba responses during seed germination and plays a role in the aba signalling pathway in arabidopsis. in this study, the expression level of abr was significantly higher in the inflation stage of the stem tumour. this suggests that aba signalling is related to stem inflation. as additional file 6: table s shows, four genes may be involved in the'plant hormone signal transduction' pathway, suggesting that plant hormones are related to stem swelling; however, abr was not one of the four genes annotated in the'plant hormone signal transduction' pathway in this study. during the inflation stage of stem tumours, plant hormones play key roles in the number of cells in the stem, which increase rapidly accompanied by cell splintering. pks expression was upregulated during the tumour inflation stage. to our knowledge, pks is phosphorylated in a phytochrome-dependent manner and negatively regulates phytochrome light signalling in arabidopsis. other reports have shown that pks regulates root phototropism and gravitropism and leaf flattening and positioning, and that it affects the state of phytochrome a in etiolated arabidopsis seedlings. chen et al. found that illumination and temperature were the main factors affecting the formation of stem tumours, but they did not determine whether pks plays key roles in illumination. is pks related to stem inflation in tumourous stem mustard? whether or not the overexpression or knockout/down of abr and pks alters the trend of inflation requires further research.  although the molecular functions of individual tumourous mustard genes and the associated signal transduction pathways remain largely unknown, the present transcriptome analysis provides valuable information regarding tumourous stem development, which may facilitate future investigations of the detailed regulatory mechanisms and pathways. additionally, we obtained about gb of raw sequence data and predicted assembled unigenes, and annotated unigenes using illumina sequencing technology. this is the first large-scale transcriptome and expression study of the tumourous stem mustard, brassica juncea var. tumida tsen et lee for which the data has been deposited in genbank. we believe that our data not only provide more molecular information for researchers, but will also help to accelerate gene expression and function research in brassica juncea as well as that on the evolutionary relationships of the cruciferae. authors' contributions cyf conceived the study. cyf, sq, zxy participated in experiment materials preparation. zxy participated in rna extraction. sq analyzed data and performed qrt-pcr materials. cyf, sq wrote the paper. zgf provided the plant varieties and participated in experiment materials preparation, fyh, lyh, hxh, sjj, jhz, hdw, pz, xlx and yjp participated in part of experiment. all authors read and approved the final manuscript. supplementary material additional file 1 table s primers used in qrt-pcr for validating differentially expressed genes. click here for file additional file 2 table s top hits obtained by blastx for the unigenes. blastx against the nr, swiss-prot, kegg and cog protein database was used with a cutoff e-value < 1e-. click here for file additional file 3 table s ko's assigned of unigenes. click here for file additional file 4 table s expression level of unigenes in each library. click here for file additional file 5 table s discription of thirteen screened out genes. expression ratios absolute values > in six comparison groups were used for rule. click here for file additional file 6 table s ko's assigned of screening different genes. click here for file additional file 7 table s the unigenes' id of screened out genes. click here for file acknowledgements we acknowledge the beijing genomics institute at shenzhen for its assistance in original data processing and related bioinformatics analysis. this work was supported by the natural science foundation of china and chongqing natural science foundation.
tumourous stem mustard is an economically and nutritionally important vegetable crop of the cruciferae family that also provides the raw material for fuling mustard. a total of genes with significantly different expression levels occurring simultaneously in the six comparison groups were screened out. finally, the altered expression levels of a number of randomly selected genes were confirmed by quantitative real-time pcr. to gain greater insight into the molecular mechanisms underlying stem swelling in this vegetable and to provide additional information for molecular research and breeding, we sequenced the transcriptome of tumourous stem mustard at various stem developmental stages and compared it with that of a mutant variety lacking swollen stems. our data provide comprehensive gene expression information at the transcriptional level and the first insight into the understanding of the molecular mechanisms and regulatory pathways of stem swelling and development in this plant, and will help define new mechanisms of stem development in non-model plant organisms. the genetics breeding, physiology, biochemistry and classification of mustards have been extensively studied, but little information is available on tumourous stem mustard at the molecular level. in our analysis, we assembled genetic information for tumourous stem mustard at various stem developmental stages. transcriptometumourous stem mustarddevelopmentrna sequencingbrassica juncea var. in addition, we constructed five dge libraries, which covered the strains yong'an and dayejie at various development stages. using illumina short-read technology with a tag-based digital gene expression system, we performed de novo transcriptome assembly and gene expression analysis. we also compared the gene expression profiles of three swollen stem samples with those of two non-swollen stem samples.
 tumourous stem mustard is an economically and nutritionally important vegetable crop of the cruciferae family that also provides the raw material for fuling mustard. the genetics breeding, physiology, biochemistry and classification of mustards have been extensively studied, but little information is available on tumourous stem mustard at the molecular level. to gain greater insight into the molecular mechanisms underlying stem swelling in this vegetable and to provide additional information for molecular research and breeding, we sequenced the transcriptome of tumourous stem mustard at various stem developmental stages and compared it with that of a mutant variety lacking swollen stems.  using illumina short-read technology with a tag-based digital gene expression system, we performed de novo transcriptome assembly and gene expression analysis. in our analysis, we assembled genetic information for tumourous stem mustard at various stem developmental stages. in addition, we constructed five dge libraries, which covered the strains yong'an and dayejie at various development stages. illumina sequencing identified unigenes, including clusters and singletons. the unigenes were subjected to a blast search and annotated using the go and ko databases. we also compared the gene expression profiles of three swollen stem samples with those of two non-swollen stem samples. a total of genes with significantly different expression levels occurring simultaneously in the six comparison groups were screened out. finally, the altered expression levels of a number of randomly selected genes were confirmed by quantitative real-time pcr.  our data provide comprehensive gene expression information at the transcriptional level and the first insight into the understanding of the molecular mechanisms and regulatory pathways of stem swelling and development in this plant, and will help define new mechanisms of stem development in non-model plant organisms. transcriptometumourous stem mustarddevelopmentrna sequencingbrassica juncea var. tumida 