Genome-wide identification and characterization of MADS-box family genes related to organ development and stress resistance in Brassica rapa

MADS-box transcription factors (TFs) are important in floral organ specification as well as several other aspects of plant growth and development. Studies on stress resistance-related functions of MADS-box genes are very limited and no such functional studies in Brassica rapa have been reported. To gain insight into this gene family and to elucidate their roles in organ development and stress resistance, we performed genome-wide identification, characterization and expression analysis of MADS-box genes in B. rapa. Whole-genome survey of B. rapa revealed 167 MADS-box genes, which were categorized into type I (Mα, Mβ and Mγ) and type II (MIKCc and MIKC*) based on phylogeny, protein motif structure and exon-intron organization. Expression analysis of 89 MIKCc and 11 MIKC* genes was then carried out. In addition to those with floral and vegetative tissue expression, we identified MADS-box genes with constitutive expression patterns at different stages of flower development. More importantly, from a low temperature-treated whole-genome microarray data set, 19 BrMADS genes were found to show variable transcript abundance in two contrasting inbred lines of B. rapa. Among these, 13 BrMADS genes were further validated and their differential expression was monitored in response to cold stress in the same two lines via qPCR expression analysis. Additionally, the set of 19 BrMADS genes was analyzed under drought and salt stress, and 8 and 6 genes were found to be induced by drought and salt, respectively. The extensive annotation and transcriptome profiling reported in this study will be useful for understanding the involvement of MADS-box genes in stress resistance in addition to their growth and developmental functions, which ultimately provides the basis for functional characterization and exploitation of the candidate genes for genetic engineering of B. rapa.

Background MADS-box genes play important roles in many aspects of plant development [1]. They are the major components in the well-known ' ABC' model that describes their roles in floral organ development [2]. MADS-box genes were identified initially as floral homeotic genes and are some of the most extensively studied transcription factors (TFs) involved in developmental control [3][4][5]. MADS-box proteins are characterized by the presence in the N-terminal region of a conserved MADSbox DNA-binding domain of approximately 58-60 amino acids that binds to so-called CArG boxes (CC[A/ T] 6 GG) [6].
Plant MADS-box genes have been subdivided into two main groups viz. M-type, also designated as type I, and MIKC, also known as type II [7]. The M-type MADSbox genes are grouped into Mα, Mβ and Mγ based on phylogenetic relationships within their MADS-box regions [4]. The MIKC genes are characterized by the presence of keratin-like (K) domain and are classified as either MIKC c or MIKC*-type [8]. The MIKC c genes are further partitioned into 14 clades based on phylogeny [9].
MIKC-type proteins generally contain four common domains. In addition to the MADS (M) domain, MIKC proteins contain intervening (I), K and C-terminal (C) domains [10,11]. The I domain is relatively less conserved, and contributes to the DNA binding specificity and dimerization of these proteins [12]. The K domain is characterized by a coiled-coil structure that mainly functions in the dimerization of MADS-box proteins. The K domain, which is present in MIKC MADS-box proteins but absent from M-type proteins, is more highly conserved than the I domain [4,13], and the MIKC* group has longer I domains and less conserved K domains than the MIKC c group [8]. The C domain, which is the least conserved, plays important roles in transcriptional activation and the formation of multimeric MADS-box protein complexes [14].
The biological functions of MIKC c genes in flower organogenesis can be grouped into five classes, A, B, C, D and E, which are required in different combinations to specify the identity of sepals (A + E), petals (A + B + E), stamens (B + C + E), carpels (C + E) and ovules (D + E) [20,24,25]. Expression of MIKC c genes has also been detected outside reproductive organs, e.g., of genes belonging to the AGL12 and AGL17 subfamilies [1,26]. This expression suggested a role for those genes in vegetative development, which was later demonstrated for some of them in root development. Nevertheless, AGL12 and AGL17 have been proposed to play roles as flowering promoters [27]. By contrast, M-type (type I) MADS-box genes in Arabidopsis appear to function exclusively during female gametophyte and seed development [28].
The genus Brassica includes a number of important crops that provide oil, vegetables, condiments, dietary fiber, and vitamin C [29]. Among Brassica species, Brassica rapa comprises several subspecies, including Chinese cabbage (B. rapa ssp. pekinensis), non-heading Chinese cabbage (B. rapa ssp. chinensis) and turnip (B. rapa ssp. rapifera). Chinese cabbage is one of the most important vegetables in Asia. In addition, B. rapa is used as the model species representing the Brassica ' A' genome and, therefore, was selected for genome sequencing [30,31]. This species has already proven a useful model for studying polyploidy, in part because it has a relatively small genome [approximately 529 megabase pairs (Mbp)] compared to other Brassica species. Comparative genomic analysis confirmed that B. rapa underwent genome triplication since its divergence from Arabidopsis [32]. MADSbox family genes have been thoroughly studied in its close relative Arabidopsis, but have not been characterized in the relatively large and complex genome of B. rapa. Over the course of evolution, the number of genes in this family steadily increased as the reproductive system became more complex; concomitant with this expansion of the lineage, MADS-box genes have been found to perform more diversified functions [33]. In addition to growth and development-related functions, some stress-responsive MADS-box genes have also been reported in wheat and rice [34,35]. As an important vegetable crop worldwide, Brassica species are subject to a variety of abiotic stresses. Identification of stress-resistance-related MADSbox genes in Brassica could be highly useful.
The recent sequencing of the Brassica rapa ssp. pekinensis genome [36] offers the possibility of genome-wide analysis of MADS-box genes. In this study, we analyzed the genomic localization, protein motif structure, phylogenetic relationships, and gene structure of all candidate MADS-box genes in B. rapa. We carried out extensive expression profiling for specific MIKC c subfamilies in vegetative and reproductive organs, as well as during flower developmental stages. Additionally, we investigated a considerable number of MADS-box genes, selected from whole-genome, low temperature-treated microarray data in the cold-tolerant and -susceptible inbred lines of B. rapa, Chiifu and Kenshin, respectively.

Results
Identification and sequence analysis of MADS-box genes in B. rapa A set of 167 candidate MADS-box genes from the B. rapa genome was recovered using key word 'MADS-box' to search Swissprot annotations at the Brassica database (BRAD) (http://brassicadb.org/brad/) [37]. This number of candidates B. rapa (167) is higher than the number of MADS-box genes in Arabidopsis, rice, soybean, maize and sorghum (Additional file 1: Table S1) [4,35,38,39]. A domain search using EMBL (http://smart.embl.de/smart/ set_mode.cgi?GENOMIC=1) with the corresponding B. rapa candidate protein sequences confirmed 162 of them to contain a 'MADS' domain, whereas the other 5 did not. The five candidates (BrMADS85, 87, 89, 119 and 127) that lacked a 'MADS' domain shared considerable sequence similarity with MADS-box proteins of other crop species that also lack 'MADS' domains and are considered to be MADS-box proteins (4 published and 1 unpublished MADS-box genes; Additional file 1: Table S2). We classified all 167 putative B. rapa MADS-box proteins into five classes (i. e., MIKC c and MIKC* of type II and Mα, Mβ and Mγ of type I) in accord with the previously reported classification of the MADS-box family members in flowering plants [4]. We designated the 167 annotated MADSbox genes of B. rapa as BrMADS followed by Arabic numbers 1-167, consecutively following the five classes (MIKC c , MIKC*, Mα, Mβ and Mγ). Subsequent sequence analysis of the 167 genes showed open reading frame (ORFs) ranging from 180 to 2379 bp and predicted protein lengths from 59 to 792 amino acid (data not shown). Sequence analysis also revealed that B. rapa MIKC (type II) MADS-box genes usually contained multiple introns, with a maximum of 15 introns; the exceptions were BrMADS84, BrMADS86 and BrMADS88, which did not have any introns. Almost all of the M-type (type I) genes lacked introns or had only a single intron; however, M-type MADS-box genes BrMADS109 and BrMADS119 had 3 and 2 introns respectively (Table 1 and Additional file 2: Figure S2). These features are consistent with those of MADS-box genes in other flowering plants such as Arabidopsis, rice, grapevine, and soybean [4,13,35,38].

Phylogenetic analysis of MADS-box genes in B. rapa
Independent phylogenetic trees for M-type and MIKCtype MADS-box TFs were constructed using the B. rapa MADS-box proteins along with those from Arabidopsis and rice. There were 67 M-type members (i.e., Mα, Mβ and Mγ) from B. rapa, with the other 100 proteins belonging to MIKC-type (MIKC c and MIKC * ; Figure 1). Notably, the MIKC c family included 89 members of this latter group, more than in Arabidopsis, rice, and soybean (Additional file 1: Table S1). Among the 89 MIKC c genes, BrMADS84, 86, 87, 88 and 89 could not be assigned in the tree using the bootstrap method with 1000 replicates, possibly due to high sequence divergence in the conserved regions and sequence length. To test their relationships and relevance with other MADS-box genes, we generated an alternative phylogenetic tree without using bootstrap replications and found these five genes in the different clades of MIKC c (Additional file 2: Figure S1b).
In accordance with the known classes of Arabidopsis MADS-box genes, we found 13 MIKC c clades in B. rapa. Although most of the B. rapa MADS-box genes were consistent with Arabidopsis in terms of sequence similarity and grouping, we found some genes viz. BrMADS41, 47, 167, that were placed as close sisters of rice MADS-box genes in the tree. Interestingly, OsMADS59, instead of being included in the AGL15-like clade, paired with BrMADS47 in the TM3 clade. There was some disparity in the distribution of rice Mβ genes between the two phylogenetic trees prepared with the different methods ( Figure 1a and Additional file 2: Figure S1a). Among the 13 MIKC c clades, the TM3 clade contained the most B. rapa sequences (18). The FLC clade included three previously identified FLC genes of B. rapa viz. BrFLC1, BrFLC2, BrFLC3 [40] which showed 99.51, 100 and 100% similarity to BrMADS13, 12 and 14 respectively at the amino acid level. MIKC*/Mδ included 11 members, which is almost double that in Arabidopsis (6), rice (5) and soybean (5).
In case of type I MADS-box proteins, the Mα and Mγ groups had more members in B. rapa (29 and 22 respectively), than in Arabidopsis, rice and soybean. By contrast, the 16 Mβ genes found in B. rapa was less than that in Arabidopsis, but more than in rice and soybean (Additional file 1: Table S1) [4,35,38].
Analysis of conserved motifs in MADS-box proteins of B. rapa Ten conserved motifs among related proteins were identified from the 167 candidate MADS-box genes of B. rapa using the MEME (Multiple Em for Motif Elicitation) motif search tool (Figure 2 and Additional file 2: Figure S3       were found to contain the K-domain motifs (2, 5, and 7) less frequently than did MIKC c proteins ( Figure 2). Comparatively less conserved motifs 3 and 4 representative of the I domain were found in both M-type and MIKC MADS-box proteins. Mβ and Mγ type proteins contained I domains at lower frequencies as compared to members of the other groups. A considerable number of non-MIKC proteins, especially from the Mα group, showed partial K domain motifs. Finally, motifs 8, 9 and 10 representing the C-terminal domains were also weakly conserved among B. rapa MADS-box genes. Motif 9 was restricted to 14 MIKC c and 1 MIKC* proteins. All Mγ proteins except BrMADS161 and 162 consistently showed both the Cterminal-representing motifs 8 and 10. Motif 8 and 10 were limited to only M-type MADS-box proteins. The Mα group showed motif 8, but motif 10 was exclusively present in the Mγ proteins. The Mβ group showed an interesting pattern, wherein 7 genes contained only a single motif, specifically one representative of the 'MADS' domain. Only 4 Mβ genes out of 16 had more than two full or partial motifs (Additional file 2: Figure S3).

Syntenic relationships between MADS-box genes of B. rapa and Arabidopsis
Polyploidy [arising from whole-genome duplication (WGD)] has played a vital role in the evolution and genetic diversity of angiosperm genomes [41]. WGD events are generally followed by changes in gene expression and widespread gene loss [42]. The Brassica genus is closely related to the model species A. thaliana and both are members of the Brassicaceae family. Comparative genetic and physical mapping as well as genome sequencing studies have authenticated the syntenic relationships between the Arabidopsis genome and the triplicate genome of B. rapa, with subgenomes having evolved by genome fractionation [43,44]. Comparative analysis was conducted to identify homologous MADS-box transcription factors between B. rapa and Arabidopsis. Based on our phylogenetic results and BLASTX reconfirmation, we determined which Arabidopsis MADS-box genes were orthologous to the 167 MADS-box B. rapa homologs. Among the homologous gene sets, we found that most Arabidopsis MADS-box genes were represented by one to three copies of B. rapa MADS-box genes (Additional file 1: Table S3).

Chromosomal location of MADS-box genes and their genomic duplication in B. rapa
We mapped the physical locations of the MADS-box genes on the 10 chromosomes of B. rapa (except two genes mapped to scaffolds Scaffold000343 and Scaffold000385; Figure 3). The highest numbers of MADS-box genes were found on chromosomes 9 (26 genes; 15.8%) and 2 (24 genes; 14.5%), while chromosomes 8 and 10 contained the fewest (10 each). Among the five types of MADS-box genes, MIKC* and Mγ genes were clustered along chromosomes 1, 6, 7, 8, 9 and chromosomes 1, 2, 5, 6, 7, 9, 10, respectively. A high of 18 MIKC c genes was found on chromosome 3, but other than that there was no bias was observed in the distribution of MIKC c , Mα or Mβ genes ( Figure 3). Duplication analysis revealed that 67 out of 167 MADS-box genes (40.12%) were present in two or more copies. This gene duplication occurred as a result of tandem and segment duplications. A total of 63 MADS-box genes were found to have counterparts on duplicated segments. We observed, higher frequencies of segmental duplications generated many homologs of MADS-box genes along all chromosomes of B. rapa (black dotted lines in Figure 3). Conversely, lower frequencies of tandem duplications were evident among M-type B. rapa MADSbox genes. Only 4 tandemly duplicated genes (from Mβ and Mγ) were found on chromosomes 1 and 4. Evolutionary analysis of B. rapa also validated our findings, wherein only 14% of the B. rapa genes were tandem duplicates, compared with 27% of Arabidopsis genes in a 100-kbp window interval [45]. No large gene clusters or hot spots for B. rapa MADS-box genes were identified, possibly due to the very few tandem duplications.
Transcript analysis of B. rapa MADS-box genes during organ development MADS-box genes have been found to be involved primarily in floral organ specification; although some recent studies revealed their involvement in other processes as well. Specifically, MIKC c proteins among all the MADSbox groups have been found to have diverse functions related to plant growth and development [1,25,35,46]. We therefore examined the expression of all 89 B. rapa MIKC c genes in root, stem, leaf and flower buds. We also investigated these genes in the sepal, petal, stamen and pistil of B. rapa flower which had expressions only in the flower buds. And, we discussed the expression of all MIKC c genes here in accord with thirteen clades identified in our study. Additionally, we included all MIKC* genes in the four floral tissue expression study as they have been reported to be involved in the development of reproductive organs [47]. Finally, we conducted an expression study in six flower bud developmental stages (young to mature bud stage) for selected MIKC c genes (those expressed only in flower buds) and all MIKC* genes to justify their roles during the flower bud development (Figure 4).

AGL15-like genes
It has been reported that AGL15 in Arabidopsis strongly delays abscission and senescence in reproductive tissues [9]. The B. rapa genome has nine AGL15-like genes (BrMADS2, 3, 4, 5, 6, 84, 85, 86, 87) and their expression in different tissues was consistent with that of their closest Arabidopsis homologs. All of the genes had predominant expression in flower buds while a few of them were expressed at low levels in different vegetative tissues (Figure 4a).

FLC-like genes
FLC acts as an inhibitor of flowering and is a convergence point for environmental and endogenous pathways that regulate flowering time in Arabidopsis [9]. We found ten FLC homologs [BRMADS1, 7,8,9,10,11, and 15 in addition to the previously identified BrFLC1 (BrMADS12), BrFLC2 (BrMADS13), and BrFLC3 (BrMADS14)] in B. rapa with very similar expression patterns in most organs. BrMADS1 is a distant member of this subfamily and showed strong expression in the four tissues tested. Our root expression results for BrFLC1 and BrFLC2 contrast with those previously reported [40]. This might be due to varietal differences of B. rapa between the two studies. BrMADS9 is the only member of this subfamily that was not expressed in any of the organ tissues (Figure 4b).

AGL17-like genes
The AGL17-like genes show unusually diverse expression patterns, with members being expressed in roots (majority of genes), in pollen (DEFH125 in Antirrhinum), in both (ZmMADS2 in maize), or in leaf guard cells and trichomes (AGL16) [9]. We identified six AGL17-like genes (BrMADS16, 17,18,19,20,21) and found expression primarily in roots of B. rapa like their Arabidopsis counterparts. Additionally, they were expressed in flower buds like in other eudicots [9]. We also observed low expression in stem and leaf tissues (Figure 4c).

STMADS11-like genes
Genes of this clade perform contrasting roles in flower development. SVP (SHORT VEGETATIVE PHASE) functions as a floral repressor, whereas AGL24 belongs to the same subfamily but promotes flowering in Arabidopsis [48,49]. We identified four genes (BrMADS22, 23,24,25) in this subfamily and detected their widespread expression in the four organs of B. rapa (Figure 4d). This is in contrast to the expression of SVP in Arabidopsis, which is restricted to leaves and shoots [9].

GGM13-like genes
The GGM13-like genes are expected to represent a sister group of the B genes and hence are termed B sister (B s ) genes [9]. ABS/TT16 is the only Arabidopsis GGM13-like gene and has been shown to function in the specification of endothelial cells as well as in the control of flavonoid biosynthesis in the seed coat [23]. We identified three GGM13-like genes (BrMADS26, 27, 28), with expression exclusively in the flower buds like their Arabidopsis counterparts. All three were expressed in the female reproductive organ of B. rapa flowers, whereas BrMADS27 was also expressed in the male reproductive organ. Interestingly, transcript accumulation of all GGM13-like genes gradually decreased from early to mature bud stage of flower development (Figure 4e).

GLO and DEF-like genes
These genes are B class floral homeotic genes in eudicots and are involved in specifying petals and stamens during flower development [50]. We found three GLOlike genes (BrMADS29, 30, 31) and two DEF-like genes (BrMADS32, 33) that were expressed exclusively in the flower buds. Transcripts for these genes were abundant in the petals and stamens of B. rapa flowers. We also found low expression in sepals and pistils (Figure 4f & 4g).

AGL12-like genes
Three AGL12-like genes (BrMADS34, 35, 88) with preferential expression in roots were detected in B. rapa. BrMADS34 and 88 were also expressed in the flower buds, similar to their Arabidopsis counterpart AGL12 with the exception that AGL12 has also been detected in shoots (Figure 4h).

TM3-like genes
These genes are expressed preferentially in vegetative parts of other plant species [51,52]. SOC1 is an important member of this family expressed abundantly in the apical meristem and acting as a flowering time regulator [53]. We identified eighteen TM3-like genes (BrMADS36, 37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52 and 89) with variable expression patterns in vegetative and reproductive parts of B. rapa. BrMADS36, 37 and 38 are close homologs of SOC1 and were primarily expressed in stem, leaf and flower buds. Moreover, we found BrMADS39, 40 and 42 to be expressed primarily in roots, but unlike their Arabidopsis counterparts (AGL14 and AGL19), we detected their expression in other parts of the plant as well (Figure 4i).

AG-like genes
Genes of this clade are mainly involved in specifying stamen and carpel identity, and in providing floral determinacy [9]. We identified eight AGAMOUS-like (AG) genes (BrMADS53, 54, 55, 56, 57, 58, 59, 60) that were expressed exclusively in flower buds of B. rapa. Our results are consistent with those for the Arabidopsis AG subfamily, members of which specify stamen and carpel identity [54]. Some of these B. rapa genes were pistil specific (BrMADS53 and 54) and some were expressed in both male and female reproductive organs (BrMADS55, 56, 57, 58, 59 and 60) (Figure 4j).

SQUA-like genes
SQUA-like genes are typically expressed in inflorescence or floral meristems, and most of them function as meristem identity genes [9]. In addition, they are involved in specifying sepals and petals and thus are class ' A' floral organ identity genes [55]. We identified ten SQUA-like genes (BrMADS61, 62, 63, 64, 65, 66, 67, 68, 69, 70) that had variable transcript patterns, but were expressed mainly in flower buds like their Arabidopsis counterparts. Some BrMADS SQUA-like genes showed strong expression in the stem and leaf as well. Our results in this case are also consistent with the Gu et al. findings, where they detected the SQUA-like gene 'FRUITFULL' in stems and leaves of Arabidopsis [21]. BrMADS67 was the only member of this subfamily expressed in all tested organ tissues of B. rapa (Figure 4k).

AGL6-like genes
The functions of AGL6-like genes are not clear. We isolated three AGL6-like genes (BrMADS71, 72, 73) from B. rapa with expression in the flower buds, like their Arabidopsis counterparts AGL6 and AGL13. BrMADS72 and 73, unlike their close homolog AGL6, also showed expression in vegetative tissues (Figure 4l).

BrMIKC* genes
There were eleven MIKC* genes (BrMADS90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100) that were placed apart from the other MIKC genes in the phylogeny. Most of these genes were found to be expressed exclusively in the stamens, except in the case of BrMADS98 and 99, that were detected in the four floral organ tissues. Moreover, these genes showed differential expression in six flower bud developmental stages (young to mature bud stage). BrMADS96, 98, 99 and 100 were preferentially expressed in the young bud stage while their expression gradually decreased until to the mature bud stage. The rest of the genes exhibited widespread expression mainly in the early stages of bud development. However, two MIKC* genes (BrMADS90 and 91) appeared to be nonfunctional, as they were not expressed in any stage of bud development or in any floral organ tissues (Figure 4n).

Microarray expression against cold and freezing stress
Four weeks old seedlings of two inbred lines of B. rapa, Chiifu and Kenshin, were treated with cold and freezing stresses (4°C, 0°C, −2°C and −4°C) during 2 hours and the expression of the 167 MADS-box genes were subsequently analyzed using microarrays. Chiifu originated in temperate regions, whereas Kenshin originated in subtropical and tropical regions and therefore, these two lines are expected to respond differently against cold and freezing stresses. Only 19 MADS-box genes from different groups showed differential cold-or freezing-responsive expression between the two lines ( Figure 5), while the remaining 148 genes showed very low or no expression (Additional file 2: Figure S4). Among the 19 differentially expressed genes, 14 MIKC c genes showed varying levels of expression, with BrMADS7, 10, 24 and 39 displaying similar expression patterns in response to cold and freezing. BrMADS11, 12, 14, 20, 23, 36, 38 and 40 were expressed at different levels than the aforementioned four MIKC c genes in both lines of B. rapa. BrMADS43 and 44, two MIKC c genes, were expressed at low levels in Chiifu throughout the stress period, while in Kenshin they showed constitutive expression. By contrast, three genes from the Mα group (BrMADS103, 109 and 127) showed differential expression within and between the two lines, with Chiifu exhibiting higher expression than Kenshin. Notably, two Mγ genes (BrMADS146 and BrMADS155) showed higher responsiveness in Kenshin than in Chiifu upon exposure to cold and freezing temperatures ( Figure 5).

qPCR expression of MADS-box genes against abiotic stress
One of our main objectives was to identify MADS-box genes that might show stress responsiveness in addition to having different growth functions. At first, a qPCR experiment was conducted to validate the cold and freezing responsiveness of the 19 BrMADS genes which were selected from the microarray analysis. We observed their expression patterns and found them consistent with the microarray results in most of the cases. Only two genes (BrMADS43 and 44) were found to show their expressions differently from those in the microarray experiment ( Figure 5) [35]. The 19 differentially expressed MADS-box genes from the whole-genome low temperature-treated data set were selected for qPCR experiments (Additional remained static except at 24 h when it was induced more than 2 fold. By contrast, these six MADS-box genes in Kenshin were down-regulated soon after drought treatment and remained that way throughout the stress period. Though BrMADS11, 24 and 127 showed up-regulation at an early stage, they eventually became inactive for the rest of the period (Figure 6c). The error bars represent the standard error of the means of three independent replicates. Values denoted by the same letter did not differ significantly at P < 0.05 according to Duncan's multiple range tests.

Discussion
Duplication among MIKC genes seems to have played major role in the expansion of MADS-box genes in B. rapa In this study, we have reported 167 MADS-box genes of B. rapa, which is higher in number than the MADS-box genes in Arabidopsis (107) [4]. The whole genome of B. rapa underwent triplication events since its divergence from Arabidopsis [32]. Thus, evolutionary relationship between B. rapa and Arabidopsis is also supportive to our findings. On the other hand, we observed the expansion of MIKC and M-type genes in these two linages. We found some disparity on the duplication events between the MIKC and M-type genes of B. rapa and Arabidopsis. For example, duplication events took place with higher frequency among MIKC-type B. rapa MADS-box genes compared to M-type genes. And, in case of Arabidopsis this scenario was reverse, where more number of M-type genes than MIKC genes was found in the duplicated segments. More specifically, 57 MIKC genes were found in duplicated segments of B. rapa (black dotted lines in Figure 3). This might be related to the fact that there are more pseudogenes of M-type than of MIKC-type MADSbox genes in the Arabidopsis genome and they experienced faster birth and death rates than MIKC type [57]. Although the B. rapa genome is triplicated relative to that of Arabidopsis, the number of M-type genes in B. rapa is almost the same as in Arabidopsis (Additional file 1: Table S1). We speculate this might be due to the presence of many non-functional M-type genes (i.e., psuedogenes) that remained inactive and were not duplicated or were deleted from the B. rapa genome. MIKC-type genes have functioned in growth and development of plants since their evolution and after multiple duplication events in B. rapa, MIKC-type genes appear to have functionally differentiated in a relatively short time and been maintained as functional genes in the genome to perform more complex functions flower and organ development.

Involvement of MADS-box genes in organ development of B. rapa Role in reproductive organ development
Investigations regarding the genetic and molecular basis of floral development in the model eudicots Arabidopsis and Antirrhinum have revealed the involvement of a number of MADS-box genes in specifying floral organ identity [58]. The high degree of sequence identity and remarkably conserved genome structure between Arabidopsis and Brassica genomes enables comparison of crop genomics among the Brassica complex [45]. In this study, we investigated the Arabidopsis MADS-box homologs in B. rapa that play specific roles in flower development.
Consideration of the ABCDE model of flower development in B. rapa revealed extensive similarities with that of Arabidopsis and other higher plants.
All SQUA-like genes in B. rapa were typically expressed in the flower buds like their Arabidopsis counterparts. AP1 is involved in specifying sepals and petals as class A floral organ identity gene [53]. Our results also suggest that BrMADS61, 62, and 63 as putative orthologs of AP1 might play similar role, and they have sepal-and petalspecific expression in B. rapa flowers (Figure 4k).
Regarding the B class genes in B. rapa, we found five close homologs of Arabidopsis PISTILLATA (PI) and APETALA3 (AP3) that showed distinct expression in male reproductive organs but not female reproductive organs. Besides being involved in the male and female reproductive parts, these genes were also recruited for petal identity in Arabidopsis [59]. We also found petal expression for them in B. rapa flowers.
Genes involved in C and D functions are from the monophyletic AG subfamily. All AG family genes in B. rapa had higher expression in female organs than in male. C and D class genes like STK/AGL11, SHATTERPROOF1 (SHP1), and SHP2, are together required for ovule identity [52]. Close homologs of SEP (SEPALLATA) genes from the AGL2-like subfamily in B. rapa showed widespread expression mainly in the aboveground parts; this is suggestive of their involvement in organ development. Pelaz et al. studied triple mutants of Arabidopsis SEP family genes (SEP1, SEP2 and SEP3) and found that their redundant functions are required for petal, stamen and carpel development and to prevent indeterminate growth of the flower meristem [20]. Genes of this family have been identified in fruits during the ripening stage of grapevine [13]. Similarly, two tomato SEP genes, TM29 and LeMADSRIN, appear to play roles in tomato fruit development [60]. The AGL12 subfamily has three members in B. rapa, two in poplar and one each in Arabidopsis and grapevine. Genes from this subfamily have found to play roles in the regulation of cell cycle in root meristems and as promoters of flowering transition through up-regulation of SOC1, FLOWERING LOCUS T (FT) and LEAFY (LFY) [27].
We found both reproductive and vegetative expression of AGL15 subfamily genes in B. rapa, as in Arabidopsis, whereas they were restricted to the flower buds, flowers and fruits in grapevine [13]. AGL15 and AGL18 are proposed to function as repressors of floral transition, acting upstream of FT and probably in combination with other floral repressors like SVP or FLC [61]. Our results regarding AGL17-like genes correspond with their expression in Arabidopsis, where they are expressed primarily in roots, which indicate that they might function in B. rapa root development. The flower bud expression of the AGL17like genes in B. rapa is also consistent with the assumption of a flowering promoter role for AGL17, which could participate in the photoperiodic induction of AP1 and LFY independent of FT [62].
Predominant expression of B. rapa MIKC* genes in the young bud stage demonstrates their importance in male reproductive organ development. Our results contrast with those for AtMIKC*, for which Verelst el al. reported predominant expression during late stages (mature pollen grain stage) of pollen development [47].
Predominant expression of three TT16 homologs (GGM13-like genes) in the early stage of female reproductive growth demonstrates their importance in the development of this organ (Figure 4e). These findings are similar to that of a previous investigation in Arabidopsis, where GGM13-like gene expression was observed in female reproductive organs, especially in ovules, which is also consistent with the situation in gymnosperms and other angiosperms [63]. Moreover, TT16 from Arabidopsis is the only GGM13-like gene for which a mutant phenotype is known. Analysis of this mutant revealed that TT16 is involved in the specification of endothelial cells and control of flavonoid biosynthesis in seed coat [23].

Role of MADS-box genes in vegetative tissue development
Transcription of a number of MADS-box genes outside flowers and fruits as well as an increasing number of mutant and transgenic flowering plants suggest that members of this gene family play regulatory roles during vegetative development also, such as in embryo, root and leaf development [1,10]. The existence of MADSbox genes in gymnosperms, ferns, and mosses, which do not form flowers or fruits, further demonstrates the role of these genes in plants is not restricted to flower or fruit development [12,64].
All homologs from the AGL17-like clade in the B. rapa genome were predominantly expressed in roots and some of them were detected in stem and leaf tissues as well. Reports from different studies indicate that AGL17-like genes show unusually diverse expression patterns in roots, pollen, leaf guard cells and trichomes. It is likely that the ancestral AGL17-like gene had an expression domain restricted to vegetative tissues [1].
In Arabidopsis, AGL18 and AGL15 showed high expression in roots, flowers, siliques, and significant expression was also observed in stem and leaves. Moreover, AGL18 was detected up to the heart stage of embryo development but not in the developing embryos at any stage [1]. Accordingly, we can also predict that BrMADS2, 3, 4 and 85 in B. rapa, as putative orthologs of AGL18, might play roles in vegetative tissue development.
TM3-like genes in Arabidopsis (AGL14 and AGL19) have been reported to function in the roots (in the columella, lateral root cap, and epidermal cells of the meristematic region and in the central cylinder of the mature roots) [1,13]. SOC1, a floral pathway integrator, expressed most abundantly in aboveground parts, is repressed by another MADS-box gene, the floral transition repressor FLC, which is involved in vernalization [65,66].
The ubiquitous expression of some B. rapa FLC genes corresponds to that of their Arabidopsis homologs. Kim et al. reported that the expression of three BrFLC genes (BrFLC1, BrFLC2, BrFLC3) was associated with flowering time and concluded that BrFLC genes act similarly to AtFLC and ultimately help in controlling of flowering time in B. rapa and other crops as well to produce higher vegetative yields [40].
The ubiquitous expression of B. rapa STMADS11-like genes suggests that these might be good candidates to play regulatory roles. Reports on STMADS11 genes from different crops demonstrated that they play important roles in developing vegetative tissues. For example, JOINTLESS, a tomato (Solanum lycopersicum) MADS-box gene is required for the development of a functional abscission zone in tomato flowers [67]. Transcripts of the potato MADS-box genes STMADS11 and STMADS16 are present in all vegetative tissues of potato, including roots and new tubers, but are not detected in floral organs [68].
BrMADS SQUA-like genes expressed in the vegetative tissues might have some regulatory roles related to vegetative tissue development. Potato MADS-box 1 (POTM1) a potato SQUA-like gene, exhibited widespread expression in actively growing tissues such as meristems, roots, new leaves and new tubers [69].

Stress responsive MADS-box genes in B. rapa
MADS-box genes have already been identified to play roles under low temperature stress in tomato [70], while seven MADS-box genes have been demonstrated to take part in stress (cold, salt and drought) responses in rice [35]. Our qPCR analysis revealed differential expression of thirteen MADS-box genes (BrMADS11,12,14,20,23,24,36,38,39,40,44,103, and 127) in response to cold stress (Figure 6a). We observed, expression patterns some of these potential genes (BrMADS23, 24, 36, 38, 44 and 103) were not consistent with the microarray results. However, we identified some candidate stress-resistance and stress-susceptibility genes based on up-and downregulation of the genes between two inbred lines, Chiifu and Kenshin, of B. rapa. We found that Chiifu, as a coldresistant line, showed more up-regulation of MADS-box genes than did Kenshin in response to cold stress via qPCR analysis. The exceptions were BrMADS24 and 36, which exhibited much higher up-regulation in Kenshin than in Chiifu and these two genes might be related to cold susceptibility in Kenshin. The highly expressed MADS-box genes in Chiifu might be involved in cold resistance, while their inactivity or very low activity in Kenshin might play a role in the cold susceptibility of that line. We also identified six (BrMAD12, 14, 20, 39,103, and 127) and eight (BrMADS11, 12,14,24,38,44,103, and 127) MADS-box genes as differentially expressed in response to salt and drought, respectively (Figure 6b & 6c). Similar phenomena as in cold stress were also observed in case of resistance against salt and drought stresses between the two lines of B. rapa. Finally, we found BrMADS12, 14, 103 and 127 to be co-responsive against all three stresses, suggesting that these genes might have multiple stress resistance related functions in B. rapa. Among the stress-induced genes, eleven were from the important MIKC c group, which is well known for regulatory roles in growth and development of different higher plants. FLC is repressed by cold and others FLC-like genes are also responsive to temperature in different ways [71]. We also identified three cold responsive B. rapa FLC-like genes (BrMADS11, 12 and 14) from this clade. In rice, all seven stressresponsive genes were also from MIKC c [35]. Likewise, in wheat, a large number of genes involved in flower development are associated with abiotic stress responses [34]. Moreover, we found two Mα genes (BrMADS103 and 127) to show stress responsiveness in B. rapa, which has not been reported in any plant yet. Our findings here serve as an important resource guiding specific investigations on the stress resistance of B. rapa related to MADSbox genes.

Conclusion
This is a comprehensive and systemic analysis of MADSbox TFs in B. rapa wherein we demonstrated their expression patterns in different growth organs and examined their responses to various abiotic stresses as well. Our data set presented here, which includes likely B and C function genes that display male organ-specific expression, should be an important resource for study of male sterility in B. rapa. Furthermore, the stress-responsive genes described in this study might be exploited for molecular breeding of B. rapa. The results presented here also facilitate selection of appropriate candidate genes for further functional characterization.

Identification of MADS-box genes
A search of SWISSPROT annotations at the Brassica database (BRAD) was conducted using keyword 'MADSbox' (http://brassicadb.org/brad/) [37]. Protein and CDS of the resulting candidate B. rapa MADS-box genes were obtained from the Brassica database (http://brassicadb.org/brad/) [37]. To confirm the presence of a MADSbox domain, the web tool from EMBL (http://smart.embl.de/smart/set_mode.cgi?GENOMIC=1) and homology searches using the Basic Local Alignment Search Tool (BLAST; http://www.ncbi.nlm.nih.gov/BLAST/) were performed on the set of candidate MADS-box genes in B. rapa. The primary structure of the genes was analyzed using protParam (http://expasy.org/tools/protparam.html). The number of introns and exons was determined by manually aligning the CDS sequences with the genomic sequences using ClustalW [72] and with the 'Gene Structure Display Server' (GSDS) web tool [73].
Phylogenetic analysis of MADS-box proteins B. rapa MADS-box proteins were aligned using ClustalX with those of rice and Arabidopsis. [74]. The phylogenetic trees were generated with MEGA6.06 using the Neighbor -Joining (NJ) algorithm [75]. Bootstrap analysis with 1,000 replicates was used to evaluate the significance of the nodes. Pairwise gap deletion mode was used to ensure that the divergent domains could contribute to the topology of the NJ tree. For generating alternative phylogenetic trees all the protein sequences were aligned in ClustalW using default parameters [72] and the phylogenetic trees were constructed using MEGA6.06 [75].

Analysis of conserved motifs in MADS-box proteins
The MADS-box protein sequences were analyzed using the MEME software (Multiple Em for Motif Elicitation, V4.9.0) [76]. A MEME search was executed with the following parameters: (1) optimum motif width ≥6 and ≤200; (2) maximum number of motifs to identify =10.

Chromosomal locations and gene duplication of MADS-box genes
All MADS-box genes of B. rapa were BLAST searched (http://www.ncbi.nlm.nih.gov/BLAST/) against each other to identify duplicate genes, with the criteria that both the similarity and query coverage percentage of the candidate genes were > 80% [77]. Positional information for all candidate MADS-box genes along the 10 chromosomes of B. rapa were obtained from the Brassica database (http:// brassicadb.org/brad/) [37]. The map of all genes along the 10 chromosomes and duplication lines among genes were drawn manually.

Analysis of syntenic relationships
To identify Arabidopsis orthologues of MADS-box genes in B. rapa, each candidate MADS-box gene nucleotide sequence was employed in a BLASTX search of the NCBI database (http://blast.ncbi.nlm.nih.gov/Blast.cgi) using A. thaliana as reference organism and the best hit A. thaliana homologue was considered to be the orthologue of the B. rapa MADS-box gene.
Collection and preparation of plant material B. rapa 'SUN-3061' plants were grown in the Department of Horticulture, Sunchon National University, Korea. For the organ study, fresh roots, stems, leaves and flower buds were harvested, frozen immediately in liquid nitrogen, and stored at −80°C for RNA isolation. For the three abiotic stress treatments, two inbred lines of B. rapa ssp. pekinensis 'Chiifu' and 'Kenshin' were used. Chiifu originated in temperate regions, whereas Kenshin originated in subtropical and tropical regions [78]. Plants were cultivated under aseptic conditions in semisolid media for 10 d, after which plants were transferred into liquid media to minimize stress during the treatment time. Three stress treatments, cold, drought and salt, were administered over 8 time periods (0 h, 30 min, 1 h, 4 h, 8 h, 12 h, 24 h and 48 h). Plant samples were transferred to the incubator at 4°C to induce cold stress. Drought/desiccation stress was simulated by drying the plants on Whatmann 3 mm filter sheets. To induce salt stress, plant samples were transferred to rectangular petri dishes (72 × 72 × 100 mm) with medium containing 200 mM NaCl for the designed time courses [35]. In each stress experiment, leaves of treated samples were collected and processed to study the expression of different MADS-box genes.

Microarray expression analysis
Br135K microarray (Brapa_V3_microarray, 3'-Tiling microarray) is a high-density DNA array prepared with Maskless Array Synthesizer (MAS) technology by NimbleGen (http://www.nimblegen.com/). Probes are designed from 41,173 genes of B. rapa accession Chiifu-401-42, a Chinese cabbage [36]. For the microarray experiment four-weekold B. rapa inbred lines, Chiifu and Kenshin, were treated with cold or freezing stress (4°C, 0°C, −2°C and −4°C). Stress treatments were applied for 2 h and immediately after stress, total and polysomal RNA was extracted from the leaf tissues using the RNeasy Mini kit (Qiagen, USA). RNA protect reagent (Qiagen) and DNA was removed by on-column DNase digestion with the RNase-Free DNase set (Qiagen). Labeling was performed by NimbleGen Systems Inc. (Madison, WI USA), following their standard operating protocol (www.nimblegen.com). The raw data (pair file) was subjected to RMA (Robust Multi-Array Analysis) [79], quantile normalization [80], and background correction as implemented in the NimbleScan software package, version 2.4.27 (Roche NimbleGen, Inc.). To assess the reproducibility of the microarray analysis, we repeated the experiment three times with independently prepared total RNA. The complete microarray data have been deposited in Omics database of NABIC (http:// nabic.rda.go.kr) as enrolled number, NC-0024-000001 − NC-0024-000012.

RT-PCR expression analysis
RT-PCR was conducted using an AMV one step RT-PCR kit (Takara, Japan). Specific primers for all genes were used in RT-PCR, and Actin primers for B. rapa (FJ969844) were used as a control (Additional file 3: Table S4). PCR was conducted using 50 ng cDNA from the plant and flower organs as templates in master mixes composed of 20 pmol each primer, 150 μM each dNTP, 1.2 U Taq polymerase, 1x Taq polymerase buffer and double-distilled H 2 O diluted to a total volume of 20 μL in 0.5-mL PCR tubes. The samples were subjected to the following conditions: pre-denaturing at 94°C for 5 min, followed by 30 cycles of denaturation at 94°C for 30 s, annealing at 55°C for 30 s and extension at 72°C for 45 s, with a final extension for 5 min at 72°C.

qPCR expression analysis
Real-time quantitative PCR was performed using 1 μL cDNA in a 20-μL reaction volume employing iTaqTM SYBR® Green Super-mix with ROX (California, USA). The specific primers used for real-time PCR are listed in Additional file 4: Table S5. The conditions for real-time PCR were as follows: 10 min at 95°C, followed by 40 cycles at 95°C for 20 s, 58°C for 20 s, and 72°C for 25 s. The fluorescence was measured following the last step of each cycle, and three replicates were used for each sample. Amplification detection and data analysis were conducted using LightCycler96 (Roche, Germany).

Additional files
Additional file 1: Table S1. Total number of MADS-box genes within each group of Arabidopsis, Rice, Soybean, Maize, Sorghum and B. rapa. Table S2. Homology analysis of 167 MADS-box genes in B. rapa. Table S3. Synteny table showing A. thaliana orthologous MADS-box gene pairs in B. rapa.
Additional file 2: Figure S1. (a) Phylogenetic analysis of 138 type I MADS-box proteins from B.rapa (67), Arabidopsis (43) and Rice (28).  Figure S3. Distribution of Conserved motifs in Brassica rapa MADS-box type I proteins identified using MEME search tool. Schematic representation of motifs identified in B.rapa MADS-box type I proteins using MEME motif search tool for each group (Mα, Mβ and Mγ) given separately. Different motifs are indicated by different colors, and the names of all members are shown on the left side of the figure. The order of the motifs corresponds to the position of the motifs in individual protein sequences. Figure S4. Microarray expression analysis of MADS-box genes in B. rapa under different temperature treatment. Here C and K indicates Chiifu and Kenshin, were treated under five (5) temperatures as control (C1&K1), 4°c (C2&K2), 0°c (C3&K3), -2°c (C4&K4), and-4°c (C5&K5). Color bar at the top representing differential expression like purple representing medium level expression where pink to white showing low to no expression.