Genomic identification, characterization and differential expression analysis of SBP-box gene family in Brassica napus

Identification and annotation of SBP-box genes in the B. napus genome

Firstly, the HMM profiles of the SBP domains (PF03110) in the Pfam database (http://pfam.xfam.org/) were downloaded and used to search the genome database of B. napus (http://www.genoscope.cns.fr/brassicanapus/) using HHMER search program. All non-redundant sequences were submitted to Interpro (http://www.ebi.ac.uk/interpro) to confirm the presence of the SBP domain. Sequences without complete SBP domain were excluded from the result. We also performed HHMER search against Brassica rapa and Brassica oleracea genome databases to identify SBP proteins. Secondly, Arabidopsis SBP protein sequences were downloaded from TAIR (http://www.arabidopsis.org/) to use as query to perform the BLASTP against B. napus genome. SBP-box gene accession numbers in B. napus genome database were extracted. The nomenclature of putative SBP-box genes in B. napus was in accordance with the homologous gene IDs in Arabidopsis. For one SBP-box gene in Arabidopsis, the orthologous SBP-box genes in oilseed rape were drawn up alphabetically. As the sequence of AtSBP1 and AtSBP12 shows high similarity, only BnaSBP1 genes were named in oilseed rape. SBP-box genes in rice were downloaded from rice genome project (http://rice.plantbiology.msu.edu/).

Gene structure, chromosomal location, duplication and phylogenetic analysis of BnaSBP genes

All the BnaSBP genes were mapped to the B. napus genome chromosomes according to the approximate position information. The exon/intron structure of each BnaSBP genes was displayed in Gene Structure Display Server program (http://gsds.cbi.pku.edu.cn/index.php) by comparing the coding sequence and genomic sequence. MCScanX software (http://chibba.pgml.uga.edu/mcscan2/) was used to analyze the duplication pattern of BnaSBP genes in oilseed rape genome. The local blast + software was used to perform the BLASTP analysis of B. napus with the e-value under 1e-5. The position of SBP-box genes and the blast output were imported into MCScanX software to generate a circle plot under a default criterion. Multiple sequence alignment of SBP-box protein sequence from Oryza sativa, Arabidopsis thaliana and Brasscia napus was performed using ClustalX2.0 with the default parameters [41]. Phylogenetic trees were constructed in MEGA6.0 software using the neighbor-joining (NJ) method and maximum likelihood (ML) method with 1000 bootstrap replications.

Conserved motif identification and miR156 target site prediction

The conserved motifs were identified using the MEME online tool (http://meme-suite.org/) with parameter setup as following: maximum number of motifs, 20; number of repetitions, any; the range of motif width was from 6 to 80. All the identified motifs were searched in InterPro database (http://www.ebi.ac.uk/interpro/) and sequence logos were created using Weblogo online software (http://weblogo.threeplusone.com/). To predict the putative target sites of miR156, full length of BnaSBP genes including exon, intron and UTR sequences were analyzed using psRNATarget tool (http://plantgrn.noble.org/psRNATarget/?function). The conserved target sequences were modified by Genedoc software.

Plant materials and growth condition

Plant samples used for expression pattern analysis and RNA-seq were collected from B. napus var. Zhongshuang 11 at the Oil Crops Research Institute of the Chinese Academy of Agricultural Sciences (OCRI-CAAS). The RNA-seq data were generated from twelve different tissues (root, leaf, bud, silique, stamen, new petal, blooming petal, wilting petal, stem, sepal, ovule and pericarp). The high resolution RNA-seq data of BnaSBP genes were kindly provided by Professor Shengyi Liu from OCRI-CAAS (data not published). The detailed FPKM value (Fragments Per Kilobase of exon model per Million mapped reads) was list in the supplemental data (Additional file 3: Table S2). The FPKM value was log2-transformed and the euclidean distances of all genes were calculated. Clustering tree was constructed and displayed by hierarchical cluster method of “complete linkage clustering” through R package.

To analyze the expression pattern of miR156 and BnaSBP genes, twelve tissue samples were also collected from the same tissue site at the same developmental stage as the sample for RNA-seq. All samples were collected and frozen in liquid nitrogen quickly and stored at the −80 °C. B. napus lines Purler and 6098B, harboring large and small branch angle respectively, were used for expression analysis. Results from different years showed that the branch angle of 6098B was 30−32° larger than that of Purler at the mature stage [42]. Tissue samples at the branch sites were collected at the bolting and early flowering stages for RNA-seq analysis. RNA-seq data were analyzed as described for Zhongshuang 11. Other tissue samples from 6098B and Purler were taken as those from Zhongshuang 11 to perform RT-PCR to verify the RNA-seq result. All plant materials were grown at the field in OCRI-CAAS, Wuhan, China.

RNA extraction and quantitative real-time RT-PCR analysis

Total RNA from diverse tissues at different growth stage was extracted with Trizol Reagent (Invitrogen, America). Before reverse transcription, total RNA was treated with RNase-free DNase I (Promega, America) for 15 min to degrade genomic DNA. Stem-loop RT-PCR was used to examine miR156 expression level in different tissues following the procedure reported previously [43]. miRNA sequences in B. napus were downloaded from miRBase Sequence Database [44]. Primers used for stem-loop RT were designed according to Zhao et al. (2012) [45]. U6 specific primer was added simultaneously as reference for accurate normalization in each reaction. As the mature sequence of miR156 family varies in the 5′ region, five different forward primers were designed for realtime qPCR. qRT-PCR was run in CFX96 Real Time System (Bio-Rad, Hercules, California, USA) using SYBR Green (Tiangen, China) according to the instructions. Briefly, 12.5 μl SYBR mixture, 1 μl universal reverse primer and 1 μl specific primer were added for each reaction. The U6 reaction as a control was conducted using the specific primer. Three replicate reactions were performed for each sample using following program: 10 min at 95 °C, 40 cycles of 5 s at 95 °C, and 30 s at 60 °C. The specificity of the amplification for each primer pair was verified by melting curve analysis. For RT-PCR, two μg of RNA was used for first strand cDNA synthesis with a Transcript First Strand cDNA Synthesis Kit (Tiangen, China) according to manufacturer’s instructions. The reaction was conducted using following program: 5 min at 95 °C, 31–37 cycles of 30 s at 95 °C, 40 s at 54–60 °C and 1 min at 72 °C. Primers used in the qPCR and RT-PCR were listed in Additional file 1: Table S1. The U6 and actin genes were selected as internal reference genes as described previously [45].

Identification of SBP genes in B. napus

Chromosome localization and gene duplication analysis

To determine chromosome distribution and gene duplication of SBP genes in B. napus, all the SBP genes except four located on unanchored scaffolds, were mapped to approximate chromosome positions (Fig. 1). These fifty-four SBP genes were unevenly distributed on the Brassica chromosomes. Except for A8 and C8, all chromosomes harbor at least one of the SBP genes. On chromosome A1, A3 and C1, only one SBP gene was found. Four chromosomes contain the maximum number of SBP genes, i.e., A5, A7, C4 and C5 each has five SBP genes. Four clusters each with two SBP genes were identified by the criteria that the distance of adjacent SBP genes is less than 50 kb. Twenty-six and thirty-two SBP genes were found to located at the A genome and C genome respectively. It was interesting to find that the number of SBP genes located at the A genome of B. napus was equal to the number of SBP gene found in the B. rapa genome. However, only 19 SBP genes were identified in B. oleracea genome, which is the progenitor of the C genome in B. napus, indicating that the SBP gene expansion may have occurred in the polyploid C genome.

The tandem and segmental duplication of Brassica SBP genes were also analyzed. Among all the SBP genes, eight members (13.8 %) showed tandem repeats, which include four clusters of tandem repeat genes (Fig. 1). In addition, 49 (84.5 %) of the fifty-eight BnaSBP genes were found to be segmentally duplicated genes. These genes were located at seventeen different chromosomes (Fig. 2).

Structural organization and conserved domain identification

To understand the evolutionary relationship among SBP protein in B. napus, we constructed the unrooted tree based on the alignments of full-length SBP protein sequences using neighbor-joining (NJ) method in MEGA 6.0. The fifty-eight SBP proteins in B. napus were divided into eight distinct groups (from Ito VII). Group I consist of the maximum number (14) of BnaSBPs, while group Vcontains only three BnaSBPs. The entire tandem duplicated BnaSBPs were assigned to one group, in accordance with the results reported in other species, such as tomato, Populus trichocarpa [40, 46]. The genomic sequence of the BnaSBP genes ranged from 510 bp to about 5 kb. To obtain further gene structure information, we compared the coding sequence with the genomic sequence of all BnaSBP genes (Fig. 3a). Different introns (from 0 to 10) were observed among the BnaSBP genes. Except BnaSBP6d, all BnaSBP genes contain at least one intron. The genes possess maximum number of introns were in group IV and VII. The BnaSBP gene clusters that were divided into the same group exhibited similar structure. Several motifs were identified among SBP proteins in B. napus (Fig. 3b). One motif (S) containing the SBP-domain was detected in all BnaSBP proteins except BnaSBP8b which contains a similar SBP-domain that could not be detected due to missing of a few amino acids. The BnaSBP protein in the same group exhibited similar motif composition.

All the BnaSBP proteins were aligned by the ClustalX 2.0 and the conserved SBP domain was created by the Weblogo online tools. Fifty-eight BnaSBP proteins contained the complete SBP domain with two Zinc motifs and one nuclear localization signal (Fig. 4). The first zinc finger motif was C3H type in all the SBP proteins except BnaSBP5 group. All the SBP proteins contain the second CCHC type zinc motif. As SBP proteins possess the character of transcription factors, all the SBP proteins contain the conserved nuclear localization signal.

Phylogenetic analysis of SBP genes in B. napus, Arabidopsis and rice

The phylogenetic relationship among BnaSBP genes and other SBP genes with known functions from other species is useful for predicting their roles in oilseed rape development. Sixteen SBP genes from Arabidopsis and nineteen SBP genes from rice, which are model plants for dicot and monocot species respectively, were extracted from the public gene pool. Fifty-eight SBP genes from oilseed rape together with the Arabidopsis and rice genes were used for the construction of an unrooted phylogenetic tree (Fig. 5, Additional file 2: Figure S2). According to phylogenetic analysis, SBP genes from these three plant species can be classified into seven groups (SBP-a to SBP-h). The largest group (SBP-e) contains 21 members which account for 23 % of the total SBPs, whereas group SBP-a forms the smallest group containing only five members. As shown in Fig. 5, genes in group SBP-a were more diverged than those in other groups. BnaSBP genes showed a high similarity to their orthologs from Arabidopsis and were classified into the same group. Among the groups revealed by phylogenetic analysis, group SBP-f only contain SBPs from Arabidopsis and oilseed rape, indicating the diversification of SBP genes between monocot and dicot plants.

MiR156 family in B. napus and their target site to BnaSBP genes

Seven putative members of miR156 (BnaMiR156a-g) in oilseed rape were found after querying the miRBase database. Recently, thirty-two putative pre-mature structures of miR156 were predicted in B. napus by high throughput small RNA deep sequencing [47]. Previous results showed that miR156 complementarily bind to SBP genes either at the coding or 3′UTR region and reduced gene activity by translation suppression or cleavage [27, 29]. It was shown that 44 SBP proteins have miR156 binding site, with 30 and 14 at coding and 3′UTR regions, respectively (Fig. 6). According to previous results, 11 out of 17 SBP genes in Arabidopsis are targeted by miR156. The homologous genes in oilseed rape are also predicted to be target of miR156. These results suggest that relationship between miR156 and SBP genes is conserved across species. However, three BnaSBP genes targeted by miR156 differed from other genes. BnaSBP5c possesses the binding site within the coding region, while the other three BnaSBP5 genes are targeted by miR156 in 3′UTR. MiR156 was predicted to bind to 3′UTR sequence of BnaSBP6d and BnaSBP10a, while the relative homologous gene in Arabidopsis were bound by miR156 at the coding region. The distinct regulation pattern of the homologous genes between B. napus and Arabidopsis reveals the divergence of the SBP-box genes in oilseed rape.

Expression profile of BnaSBP

A wide range of SBP genes play important roles in plant development process. In the absence of SBP gene mutants, the expression pattern may provide a clue to elucidate the potential role of the different SBP genes in B. napus. The expression level of BnaSBP genes in twelve tissues were shown by heat map representation (Fig. 7, Additional file 3: Table S2). Transcript of BnaSBP6c was zero in all twelve tissue samples and only very low expression level of BnaSBP4c in leaf was detected. Based on the hierarchical clustering analysis, the BnaSBP genes could be divided into eight categories. The transcription of a large number of BnaSBP genes was enriched in bud, stamen and pericarp. By contrast, most of BnaSBP genes exhibit low expression level in ovule and petal. Eight BnaSBP genes, BnaSBP1a, 1b, 11e, 14a, 14b, 14c, 16a and 16b seemed to be expressed constitutively, from root to pericarp. It should be noted that all these genes, excluding BnaSBP11e, are not predicted to be targeted by the miR156. BnaSBP4c, 4d, 5c, 5d, 10d and 13d sustained low expression level in most tissues. The expression level of BnaSBP3a and 3d was not detected in most tissue samples, but reached clearly higher levels in pericarp. A relative higher expression level of BnaSBP2b and 11d could also be discerned in root tissue. Compared with the SBP genes not bound by miRNA, the BnaSBP genes have the target site represent more divergent expression pattern. We also performed RT-PCR to confirm the expression levels of some BnaSBPs in eight different tissues (Fig. 8). Thirty-nine BnaSBPs were selected to verify the result of RNA-seq data. Results showed that RT-PCR data was generally consistent with RNA-seq data for relative expression of BnaSBPs in most of the tissues. For example, expression level of BnaSBP1a, 1b and 11e could be detected in most tissues (Fig. 8). Though BnaSBPs were expressed at least in one of the tissues, distinction of expression patterns were observed across the gene groups. Some BnaSBPs belongs to a same group exhibited similar expression pattern, such as BnaSBP1a and 1b in group IV, BnaSBP15a and 15b in group III, indicating redundant roles of BnaSBPs in the same group. Therefore, the oilseed rape SBP transcription factors have diverse expression patterns and may be redundant in biological function with each individual in charge of certain physiological processes.

To investigate the putative genes involved in branch angle regulation, the expression profile of two B. napus material (6098B and Purler) with different branch angle was conducted (Additional file 4: Figure S1). Sample of branch site from two materials at bolting and early flowering stage was harvested to perform DEGs (Different Expression Genes). The transcription level of all SBP genes was extracted from expression profile (Additional file 5: Table S3). Heat maps representing expression levels in the lines at two developmental stages are shown in Fig. 8. Many BnaSBP genes showed different expression patterns between the two lines at the two development stages. BnaSBP5c, 8a and 7d showed high expression at bolting stage but no or little expression at early flowering stage in the two materials. Ten and thirteen BnaSBP genes were found differentially expressed between the two lines at the two development stages, respectively. Among them, six BnaSBP genes were differentially expressed at the two development stages (Fig. 9). Further studies may focus on the role of these genes on branch angle regulation. RT-PCR was performed to confirm the expression level of BnaSBPs in the same tissues used for RNA-seq. A large number of BnaSBPs in Purler expressed at higher level than those in 6098B (Fig. 10). This RT-PCR result was generally consistent with that from RNA-seq data.

Expression profile of miR156

Several BnaSBP genes carry the complementary sequences to miR156. MiR156 was thus expected to be an important determinant for the expression of these BnaSBP genes. The expression level of miR156 was mostly abundant in bud and silique of Zhongshuang 11 at different developmental stages (Fig. 11a). Relative low levels were found in leaf sample. Meanwhile, the expression level of miR156 in 6098B and Purler was also determined. It was showed that the abundance of miR156 decreased significantly at early flowering time compared to bolting time (Fig. 11b). Besides the stem sample of two materials, the transcription of miR156 was stronger in Purler than in 6098B of the other tissues.

SBP-box genes in Brasscia and their evolution

The SBP-box proteins are characterized by a conserved SBP domain with 76 amino acids and constitute one large family of transcription factors in plants. Plant specific SBP-box transcription factors were only detected in green plants suggesting that it might originate predating the divergence of green algae and the ancestor of land plants [5, 48]. Different numbers of SBP-box genes have been characterized in various land plants [39, 40, 49]. In present study, 58 SBP-box genes in B. napus genome were identified, which is about four times the number of Arabidopsis SBP-box genes. B. napus contains 13 more SBP-box genes than the sum of B. rapa (AA, 2n = 20) and B. oleracea (CC, 2n = 18), which are two immediate progenitor species of B. napus (AACC, 2n = 38) [50]. For one gene family, tandem and segmental duplication events are the main reasons for gene expansion. SBP-box genes are unevenly distributed on 17 of the 19 chromosomes of B. napus, and four clusters each with two BnaSBPs were identified (Fig. 1). Uneven and cluster distribution of SBP-box gene family genes was also found in rice and peach [6]. There are seven and 49 BnaSBP genes which were found to be tandem and segmental duplications respectively. Diversification of BnaSBP genes was observed from many aspects, including phylogenesis, genomic structure, as well as location of miR156 target site. This diversity of SBP-box gene structure is likely to be trigged by gene duplication followed by intron and exon loss.

Functional divergence of SBP-box genes

As the SBP-box genes possess the character of transcription factors, their expression pattern is expected to be correlated with their function on plant development. The expression profile of BnaSBP-box genes showed distinct expression patterns among different tissues. In Arabidopsis, some SPL genes are constitutively expressed, while the transcription level of others is under developmental control [5]. Expression analysis of SBP-box genes in other organisms also presented diverse spatiotemporal expression patterns [39, 40, 49, 51]. SBP transcription factors in B. napus showed diverse expression patterns across tissues, indicating their possible functions in various biological processes. The transcription of a large number of BnaSBP genes was enriched in bud, stamen and pericarp, suggesting most of the SBP-box genes in oilseed rape may be involved in the development of reproductive organs.

SBP-box genes in many species, especially in rice and Arabidopsis, have been demonstrated to play essential roles in diverse developmental processes. The microRNA regulated SBP-box genes SPL9 and SPL15, which are the most close orthologous genes in Arabidopsis, was proven to control shoot maturation [52]. Further support of possible roles for BnaSBP in development comes from the rice genes SPL14 in panicle development and ideal rice plant architecture regulation [22, 23]. We identified four BnaSBP9 genes in oilseed rape genome. Although the BnaSBP9 genes possess similar gene structure, diverse expression patterns were observed. It should be noted that the expression of BnaSBP9d in the compact material Purler is higher than in the loose material 6098B (Figs. 9 and 10). The expression of BnaSBP9d visibly decreased from bolting to early flowering. Further study should be performed to verify whether BnaSBP9d might play a role in regulating branch angle in oilseed rape.

Arabidopsis gene SPL8 affects pollen sac development and also controls gynoecium patterning [18]. Three BnaSBP genes, BnaSBP8a, 8b and 8c showed most similarity to AtSBP8,joining the same group through phylogenetic analysis. BnaSBP8a and BnaSBP8b were highly expressed in the stamen. Further study may focus on the potential role of BnaSBP8 in flower development.

Constitutive expression of AtSPL3 resulted in early flowering [53]. The SPL3 homologous genes in Antirrhinum majus and Silver birch also regulate flower development by binding to the MADS-box genes [16, 54]. Tomato LeSPL-CNR, which is most similar to AtSPL3 gene, is crucial for normal fruit development and ripening [55]. In Arabidopsis, miR156-SPL3 module controls FT expression to regulate ambient temperature-responsive flowering [56]. Among the five genes homologous to AtSPL3 identified in B. napus in our study, BnaSBP3c showed much higher expression level in bud, stamen, silique and pericarp, indicating a possible role in the reproduction phase. Arabidopsis gene AtSPL2, AtSPL10 and AtSPL11 were shown to play important roles in determining leaf shape and embryonic morphogenesis [20, 57]. All the BnaSBP2, 10 and 11 genes were classified into a same group of SBP-e. It would be interesting to explore the exact role of these group SBP-box genes by functional characterization.

Conservation of miR156 target site in SBP-box genes

A larger number of miRNAs targets are transcription factors, such as SBP, MYB, NAC, ARF, GRAS, and AP2 [27]. MiRNAs play important roles in regulating the transcription of target genes. Previous results showed that overexpression of miR164, miR159a, and miR319 affected members of the NAC, MYB, and TCP families of transcription factor genes, respectively [58–60]. In present study, target prediction showed that 44 of the 58 BnaSBP genes were regulated by miR156. The complementary sites of miR156 locate in the coding region of 30 BnaSBP genes, and in the 3′ UTR of the other 14 BnaSBP genes. In Arabidopsis, 10 (AtSBP2, 3, 4, 5, 6, 9, 10, 11, 13, 15) out of 17 SBP genes were predicted or verified to be targeted by miR156. The other six AtSBP genes including (AtSBP1, 7, 8, 12, 14, 16) are not targets of miR156. AtSPL7 has been demonstrated to bind directly to the Cu-response element (CuRE) containing a core sequence of GTAC and regulate Cu homeostasis [3]. The 44 BnaSBP genes predicted to be targeted by miR156 are the homologous genes in Arabidopsis, which also formed 10 gene clusters. Therefore, the miR156 target site in SBP-box genes is conserved across plant species.

Over-expression of miR156 in Arabidopsis significantly represses the SPL transcription and thus reduces apical dominance, leading to dwarfism and increases in total leaf number and plant biomass [28]. The transcripts of the target SBP genes were also suppressed in other miR156 over-expression plants [29, 56]. In present study, the transcript level of miR156 was abundant in bud and silique (Fig. 11). By contrast, most putative target SBP genes with predicted miR156 target sites showed lower expression level in these tissues (Figs. 9 and 10). Among the floral organs, most BnaSBP genes showed a low expression level in petal and ovule, though transcript was relatively high in pericarp, which is a main component of silique. These results suggested that the transcript of miR156 is negatively correlated with the expression of most BnaSBP genes. The level of miR156 was declined with a concomitant rise in SPL levels during the aging time in Arabidopsis [61]. SPL9 and SPL10 mediated the transition from high levels of miR156 to high levels of miR172 through direct activation of miR172 expression, thereby promoting the juvenile to adult phase transition [57, 62]. Our results showed that the lower expression level of miR156 in 6098B with bigger branch angle than in Purler with smaller branch angle (Fig. 11) is negatively correlated with the expression difference of many SBP-box genes, eg. BnaSBP2a, 2d, 3d, 3e, 5d, 8b, 9a, 9b, 10b, 11a, 11c, 13d and 15c (Figs. 9 and 10), indicating that the SBP/miR156 module is likely involved in regulating plant architecture in B. napus.