Comparative transcriptome analysis and ChIP-sequencing reveals stage-specific gene expression and regulation profiles associated with pollen wall formation in Brassica rapa

Background Genic male sterility (GMS) line is an important approach to utilize heterosis in Brassica rapa, one of the most widely cultivated vegetable crops in Northeast Asia. However, the molecular genetic mechanisms of GMS remain to be largely unknown. Results Detailed phenotypic observation of ‘Bcajh97-01A/B’, a B. rapa genic male sterile AB line in this study revealed that the aberrant meiotic cytokinesis and premature tapetal programmed cell death occurring in the sterile line ultimately resulted in microspore degeneration and pollen wall defect. Further gene expression profile of the sterile and fertile floral buds of ‘Bcajh97-01A/B’ at five typical developmental stages during pollen development supported the result of phenotypic observation and identified stage-specific genes associated with the main events associated with pollen wall development, including tapetum development or functioning, callose metabolism, pollen exine formation and cell wall modification. Additionally, by using ChIP-sequencing, the genomic and gene-level distribution of trimethylated histone H3 lysine 4 (H3K4) and H3K27 were mapped on the fertile floral buds, and a great deal of pollen development-associated genes that were covalently modified by H3K4me3 and H3K27me3 were identified. Conclusions Our study provids a deeper understanding into the gene expression and regulation network during pollen development and pollen wall formation in B. rapa, and enabled the identification of a set of candidate genes for further functional annotation. Electronic supplementary material The online version of this article (10.1186/s12864-019-5637-x) contains supplementary material, which is available to authorized users.

fundamental structure, comprised of an outer exine and an inner intine [3]. The exine, whose synthesis is regulated by both the sporophytic tapetum and the microspore [4], is constructed primarily of sporopollenin, a robust biopolymer comprised predominantly of polyhydroxylated aliphatic compounds and phenolics [5]. The intine, which is initiated during the early stages of male gametogenesis and is controlled gametophytically [6], is consisting of pectin, cellulose, hemicellulose, hydrolytic enzymes, hydrophobic proteins. As the major event after microspore release from the tetrad, the entire dynamic complex and well-coordinated process of pollen wall development requires the precise spatial and temporal cooperation of gametophytic and sporophytic tissues and metabolic events [7][8][9]. Recently, well-characterized genes in which mutations cause impaired exine and male sterile phenotype have enriched our understanding in pollen wall development [4,8,10,11]. One of the most striking discoveries is the conversed exine regulation pathway that formed by five transcription factors, DYT1, TDF1, AMS, MYB80, and MS1. This pathway regulates tapetum development and function and thereby influences the developing microspores by controlling callose dissolution, pollen exine formation and tapetal programmed cell death [12][13][14][15][16][17]. Many other genes, such as lipid transfer protein family members and genes related to lipid and phenolic metabolism are involved in pollen exine formation [18][19][20][21][22][23][24][25][26][27][28]. Previous studies on male sterility focusing on pectin degrading enzymes such as pectin methylesterases (PMEs), polygalacturonases (PGs) and pectate lyases like proteins (PLLs) have also emphasized the important roles of cell wall modification-related genes during the regulation of intine formation and male fertility [29][30][31][32]. However, the molecular mechanisms underlying pollen wall patterning remain largely elusive.
Genic male sterility (GMS) is a main type of crop male sterility. In our previous study, a Chinese cabbage (B. rapa ssp. chinensis cv. Aijiaohuang) genic male sterile A/B line system, named as 'Bcajh97-01A/B' , was established. The progenies of the A/B line were segregated into sterile and fertile types during reproduction at a 1:1 ratio. In the sterile plant, an aberrant meiotic cytokinesis at the early pollen developmental stage resulted in the degeneration of microspore content and consequently the formation of aborted mature pollen grains only coated with a defective exine wall [33,34]. Thus, this GMS system serves as a good material to study the pollen wall formation and pollen development.
To better understand the mechanism of pollen wall development and male sterility, in this work, using the sterile and fertile floral buds of 'Bcajh97-01A/B' as materials, a detailed gene expression profile at five typical developmental stages, namely, pollen mother cells, tetrad, uninucleate pollen, binucleate pollen, and mature pollen stage, were further examined by RNA sequencing (RNA-seq). Together with the detailed phenotype investigation on the difference of dynamic pollen callose change and tapetum degradation between the sterile and fertile anther, candidate genes that associate with these two developmental events were identified. What's more, by using ChIP-sequencing (ChIP-seq), we mapped the genomic and gene-level distribution of trimethylated histone H3 lysine 4 (H3K4) and H3K27, two histone modifications associated with gene activation and silencing, respectively, on the fertile floral buds, to explore the epigenetic control on gene expression during pollen development. These studies provided a deeper understanding into the gene expression network during pollen wall development process in B. rapa, and enabled the identification of a set of candidate genes for further functional annotation.

Phenotypic characterization of anther and pollen development in the fertile and sterile lines
The previous morphological observation showed the only difference between the sterile and fertile lines was that the anthers of the sterile line lacked normal mature pollen grains [33]. Here, scanning electron microscopy (SEM) showed that globular remnants with rough surface layers filled in pollen sac of the sterile anther instead of ellipsoidal pollen grains with reticulate exine structure and distinct apertures filling in that of the fertile anther ( Fig. 1b and f ). Brightfield microscopy revealed that numerous yellowish oily droplets were deposited on the surface of the globular remnants, demonstrating abnormal accumulation of exine-held materials in the sterile anther ( Fig. 1c and g). The 4′, 6-diamidino-2-phenylindole (DAPI) staining showed no nucleus formed in the globular remnants ( Fig. 1d and h). Previous cytological observation have determined that the aberrant cytokinesis at the end of meiosis caused failure of tetrads formation leading to the male sterility of the sterile line [34]. Here, aniline blue staining showed that the dyads of the sterile and fertile lines have identical callose deposition at the periphery and in the cell plate. However, in contrast to the normally deposited callose in the cell plate of tetrads in the fertile line, the intersporal walls were absent in the tetrads of the sterile line, while, the karyokinesis was demonstrated to be normal by DAPI staining (Fig. 2).
An in-depth cytological study and comparison was carried out on tapetum development of the sterile line utilizing semi-thin section and transmission electron microscopy (TEM) observation. No obvious difference was observed prior to tetrad stage. After meiosis, the sterile anther exhibited premature tapetal programmed cell death (PCD). The tapetal cells had indistinct cytoplasm and were loosely arranged. While, the tapetal cells of the fertile anther were tightly arranged and appeared metabolically active with abundant organelles. It was not until middle uninucleate stage, the tapetal cells of the fertile anther began to undergo PCD. At the binucleate stage, the tapetal cell in the fertile anther contained differentiating tapetosomes and elaioplasts, and the tapetum began to be thin. But in the sterile anther, no typical tapetosome was found. Elaioplasts with vague outlines differed slightly from the distinct elaioplast globules of the fertile anther which contained numerous round globuli within their stroma. At the late binucleate stage, most space between tapetal cells was consumed by elaioplasts and tapetosomes in the fertile anther, while the tapetal cytoplasm had been completely degraded in the sterile anther. However, the tapetal cell wall retained in the anther locule of the sterile line until the mature pollen stage when the whole tapetum layer was absent from the locule in the fertile anther (Fig. 3).

Transcriptome assembly and annotation of Unigenes
To understand the pollen abortion-causing mechanism in 'Bcajh97-01A/B' and identify candidate genes contributing to anther and pollen development in B. rapa, via RNA-seq, we performed a detailed gene expression profiling on both of the sterile and fertile floral buds at five typical pollen developmental stages. After filtering out low quality data in both fertile and sterile libraries of each stage, over 14,500,000 reads (designated herein as "clean" reads) were remained. All clean reads were assembled by running Trinity, and 25,509,284 contigs (including 101,292 contigs > 200 bp) were generated. The contig length distribution was shown in Additional file 1: Figure S1. After clustering, 207,932 transcripts and 72,168 Unigenes were obtained (Additional file 2: Table  S1), and 15,353 Unigenes (21.27%) were greater than 1000 bp in length and no Unigenes were shorter than 200 bp (Additional file 3: Figure S2).
For annotation, 72,168 Unigenes were subjected to BLASTX searches against the sequences in the NCBI non-redundant protein sequences (NR), Swiss-Prot, Gene Ontology (GO), the Clusters of Orthologous Groups (COG) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases. As a result, a total of 45,530 Unigenes (63.09% of all Unigenes) provided a significant BLAST result (Additional file 2: Table S2). Among the 45,530 Unigenes, approximately 32.83% could be annotated in COG classification. 14,946 Unigenes were classified into 24 function classifications (Additional file 4: Figure S3), which means that the identified genes are involved in various biological processes. GO classifications were also obtained to investigate the function of the Unigenes. In total, 31,340 annotated Unigenes were further classified into 49 functional groups (Additional file 5: Figure S4).

Transcripts differentially expressed in the fertile and sterile floral buds
With the restrictive conditions of False Discovery Rate (FDR) < 0.01 and log2 ratio > = 1.0, Unigenes that were differentially expressed in the fertile and sterile floral buds at five stages were identified. In total, 8288 genes (11.48% of Unigenes) were differentially expressed in at least one stage of the sterile floral buds compared with the fertile ones, and these genes were designated herein as DEGs ( Fig. 4a and Additional file 2: Table S3). Among these DEGs, down-regulated genes accounted for the majority in each floral buds except A1 floral buds, i.e. the sterile floral buds at stage I ( Fig. 4b and c). And as growth progresses, the number of down-regulated genes increased dramatically (Fig. 4c). Of the down-regulated genes, 4503 (77.70%) genes showed the lowest expression in stage IV or stage V in the sterile floral buds, the pollen maturation stage, reflecting the lack of mature pollen grains in the sterile line. More remarkable, 73.04% of down-regulated DEGs (764 genes) in stage II showed their differential expression only at this stage, suggesting that many genes are specific for meiosis or tetrad formation (Fig. 4c). In contrast, the number of the up-regulated genes changed gently, and 2030 (67.40%) genes were up regulated in the sterile floral buds at stage III or stage IV (Fig. 4b). Among the 8288 DEGs, 519 genes showed opposite trend at different stages (Fig. 4d). Hereinto, 217 genes were down-regulated at stage II, while up-regulated at stage III in the sterile floral buds compared with fertile ones. Further analysis showed that, as the fertile floral buds growth, these genes were down-regulated from stage II to stage III (Fig. 4e), however, the expression reached peak at stage III in the sterile floral buds (Fig. 4f ). This result demonstrated that the expression of these genes were delayed during the developmental process of the sterile floral buds.
GO annotation information of the DEGs was selected for further functional analysis. To reveal significantly enriched GO terms in DEGs comparing to the transcriptome generated by this study, GO enrichment analysis of functional significance was performed using AgriGO, and the GO term with FDR ≤ 0.05 was defined as significantly DEGs enriched GO term. This analysis allowed us to identify the major biological processes (BP), molecular functions (MF) and cellular components (CC) with which DEGs were involved at each stage. For enriched BP, more GO terms were significantly enriched by the development of anther in down-regulated DEGs. At each stage, 13, 11, 35, 45 and 86 terms were revealed, respectively ( Fig. 5b and Additional file 2: Table S4), which means that more biological processes were disturbed through the development. Several aspects are worthy of note. First, pollen development, pollen wall assembly, and pollen exine formation were GO terms significantly enriched at stage II through to stage V. Second, some GO terms related to cell wall organization or biogenesis were enriched at all stage III, IV and V. Third, a large number of GO terms finally pointing to pollen tube growth, such as pollination, cell differentiation, cell morphogenesis involved in differentiation, and cell growth, were also enriched at stage III through to stage V (Fig.  5c, Fig. 6 and Additional file 2: Table S4). In addition, pollen sperm cell differentiation were also enriched at stage IV and stage V (Additional file 2: Table S4). While at stage I, the enriched GO terms were different from the later four stages (Additional file 2: Table S4). For up-regulated DEGs, 69, 106, 55, 74 and 29 GO terms were significantly enriched at the five stages, respectively  Table S4). Most of these GO terms were associated with responding to stimulus, especially in the earlier four stages. There were some enriched GO terms which were related to protein folding, morphogenesis, and cell death. It is worth mentioning that at stage II, the significantly enrichment GO terms involved programmed cell death, host programmed cell death induced by symbiont and toxin metabolic process (Additional file 2: Table S4), indicating that these DEGs may relate to the phenotype of premature tapetum degradation in the sterile line. For enriched MF, it is worth pointing out that one of hydrolase activity, polygalacturonase activity were remarkably enriched at stage III, IV, and V (Additional file 2: Table S4), with the gene number increasing through development. Another enriched MF and CC GO terms were also showed in Additional file 2: Table S4. Using the KEGG database, the complex biological behaviors of genes were further studied. A total of 56, 83, 100, 99 and 103 biological pathways were identified by KEGG pathway analysis for the DEGs at five stages, and pathway enrichment analysis showed that there were 11, 11, 12, 11 and 18 pathways were significantly enriched, respectively (Additional file 2: Table S5). There were 10, 8.16, 8.36, 9.02 and 5.28% DEGs belonging to the enriched pathway "plant hormone signal transduction" at the five stages, respectively. These DEGs were involved in complicated regulatory network of plant hormone signal transduction, including auxin, cytokinine, abscisic acid, ethylene, brassinosteroid, jasmonic acid and salicylic acid. While, gibberellin signal transduction was not included. "Starch and sucrose metabolism" was significantly enriched at the later four stages. According to the KEGG maps, the up-or down-regulated DEGs belonging to this pathway mainly took part in the regulation of saccharometabolism, such as the biosynthesis of pectin, the production of monosaccharide and the degradation of pectin or cellulose, which suggests that early aberrant meiotic cytokinesis may finally affects the pollen wall modification.
Callose metabolism-related genes showed dramatically altered expression in the sterile line Callose, a polymer of β-1, 3 glucans, serves as a temporary wall to separate newly formed microspores in the tetrad and acts as the mold for primexine. Callose defects can affect pollen wall formation and pollen viability. The expression changes of 15 callose synthase genes were analyzed in the sterile line. None of these 15 genes displayed changed expression at stage II, suggesting normal synthesis of callose required for tetrad formation.
The accurate timing of callase activation is critical for normal callose degradation. Therefore, we analyzed the expression of all B. rapa genes encoding putative β-1, 3-glucanases. Among the 52 genes that were detected in the sequencing, 13 genes showed different expression between the fertile and sterile lines at least one floral developmental stage. Two major groups were distinguished according to the dominant expression pattern of these genes in the fertile line. One group includes seven early-expressed genes, and the other one contains six genes that expressed during the late pollen maturation stages. All of the six late-expressed genes showed reduced expression levels in the sterile floral buds. Interestingly, seven early-expressed genes showed consistent changes in gene expression between the fertile and sterile lines with a decrease first at stage II followed by a remarkable up-regulation at stage III in the sterile floral buds. For example, the expression levels of Bra032758 and Bra037057, two orthologs of Arabidopsis A6, were nearly reduced by half first but then up-regulated over 30-fold, and 60-fold respectively. Bra001918 showed a sharp decrease (< 5% remained) at stage II followed by an over 70-fold increase at stage III. While, further analysis of these genes revealed that, in the fertile floral buds, they were almost expressed specifically at one developmental stage, which is stage II. As regards for regulation on callose metabolism, Bra004288, the ortholog of Arabidopsis CDM1 was expressed at stage II specifically in the fertile floral buds and decreased dramatically in the sterile floral buds ( Table 1).
The expression of genes presumed to be involved in the formation of pollen exine was changed in the sterile line GO analysis revealed that GO term "pollen exine formation" was enriched at all the late three developmental stages and "sporopollenin biosynthetic process" was specially enriched at stage III, suggesting seriously impaired pollen exine formation. The Brassica and Arabidopsis genera share about 85% exon sequence similarity [35], and since genes regulating anther and pollen development have been well established in Arabidopsis by genetic and molecular biological studies, here, we analyzed the expression of B. rapa orthologs of Arabidopsis genes that are known to be involved in the formation of pollen exine. It was known that the biosynthesis and transport of the lipidic and phenolic precursors was important for the formation of pollen exine [8]. Therefore, we examined the expression of B. rapa orthologs of Arabidopsis genes that were associated with lipid metabolism and phenolic metabolism during pollen wall development.
Our results indicated that all the B. rapa orthologs of those tapetum-specific genes including MALE STERIL-ITY 2 (MS2), CYP703A2, CYP704B1, ACYL-COA SYN-THASE 5 (ACOS5), POLYKETIDE SYNTHASE A (PKSA), POLYKETIDE SYNTHASE B (PKSB) and TET-RAKETIDE α-PYRONE REDUCTASE 2 (TKPR2) were dramatically affected in the sterile floral buds. It was worth mentioning that they showed a decrease in Furthermore, our analyses revealed that six orthologous genes related to transport of precursors required for exine development showed different expression levels between the sterile and fertile lines. Among them, Bra039378 and Bra005048, the orthologs of Arabidopsis ATP-BINDING CASSETTE G26 (ABCG26) and ABCG1, performed with a similar change in expression as with those lipid and phenolic metabolism-related genes mentioned above. And other four orthologs of Arabidopsis ABCG transporters were all down-regulated in the sterile line. However, four orthologs of Arabidopsis lipid transfer protein genes exhibited increased expression levels ranging from 3.71-fold to 537.30-fold in the sterile floral buds at the late pollen developmental stages. The above data suggested that transport of precursor across anther tissue during pollen development was partially affected in the sterile floral buds.
We also examined the expression of several orthologs of Arabidopsis transcription factors (TFs) which had been reported to form regulatory networks to control pollen wall development. The expression of BrDYT1 (Bra013519) was up-regulated at stage II in the sterile floral buds while some genes essential for post-meiotic tapetal function including BrAMS (Bra002004 and Bra013041), BrMYB80 (Bra002847 and Bra035604) and BrMS1 (Bra002401) decreased first in expression at stage II and then increased at stage III. In addition, Bra004689, the ortholog of Arabidopsis bZIP34 exhibited reduced expression level in the late-stage sterile floral buds ( Table 2).

Screening of potential tapetum development and function-related genes
As the innermost sporophytic cell encasing the developing pollen, tapetum plays an essential role controlling pollen exine formation and pollen development. Previous expression analysis of DEGs has revealed that a part of tapetum development and function-related genes displayed a similar developmental expression change. As genes associated with the same metabolic pathway are perceived to be more highly coexpressed than genes from different pathways, here, we analyzed those genes which were down-regulated at stage II and then up-regulated at stage III in the sterile line to screen potential participants in tapetum development and function. The information of their homologous genes in Arabidopsis provided deeper understanding of functions these genes may perform during anther development.
Besides those pollen exine formation-involved genes we have mentioned above, twenty-five other genes encoding different kinds of proteins drew our attention. Their homologous genes in Arabidopsis were found to be specifically expressed or highly expressed at flower stage 9 according to the data of Arabidopsis eFP browser (http://bar.utoronto.ca/efp/cgi-bin/efpWeb.cgi), coinciding with those of tapetum development and function-related genes that have been reported to be involved in pollen exine formation. What's more, nearly all these genes showed variable expressional changes in several Arabidopsis floral mutants (Table 3). According to the microarray data, those genes with expressional changes in ems1, spl, tdf1 and ams were all down-regulated without exception. Interestingly, nearly two-thirds of those genes with expressional changes in ms1 exhibited up-regulated expression coinciding with expression analysis at stage III in our study. Further analysis revealed that Arabidopsis homologs of eight genes, Bra038803, Bra020920, Bra028324, Bra032758, Bra029151, Bra022571, Bra028286 and Bra016531, displayed absolutely similar expressional changes to those functionknown tapetum-related genes during distinct flower developmental stages in the ms1 mutant.
To verify the function of these newly identified genes that may be involved in tapetum development, we selected one of them, Bra016531, which encodes a single MYB domain-containing protein, to be further investigated. Quantitative RT-PCR analysis of Bra016531 transcripts in different tissues revealed that Bra016531 was specific to inflorescence (Fig. 7a). Its detailed expression pattern was then examined in Arabidopsis plants expressing a Bra016531 promoter-driven β-glucuronidase gene (Bra016531pro::GUS). In more than 10 independent transgenic lines, Bra016531 promoter-driven GUS activity was apparent only in the anthers of young flowers (Fig. 7b). Then, Bra016531 promoter-driven GUS expression was examined in anther cells using thin sections taken from floral buds at various developmental stages ( Fig. 7c-h). Blue GUS signal was first detected both in sporocyte and tapetal cells during meiosis (Fig.  7d), and then in tetrads and tapetal cells (Fig. 7e). Later, the signal was detected mainly in tapetal cells at the early vacuolate stage (Fig. 7f ). These temporal and spatial patterns of Bra016531 expression supports our hypothesis that Bra016531 is associated with tapetum development and function during pollen development.
A high proportion of cell wall modification-related genes were down regulated in the sterile line We found a high proportion of cell wall modification-related genes were down regulated in the sterile floral buds, especially at stage V with more than 12% of the total genes identified by GO analysis. Among these genes, some cell wall hydrolytic enzyme-encoding gene families that involved in polysaccharide metabolism captured our attention. Among PMEs family, 15 genes were down-regulated in the sterile floral buds. Most of these PMEs had very low or undetectable expression during the early stages but went up dramatically at the mature pollen stage. For example, Bra000438 (homolog of AtVGD1), Bra003491 (homolog of AtVGDH2) and Bra028699 (homolog of AtPPME1) were highly and specifically expressed in the fertile floral buds at stage V. The similar general tendency was apparent in another cell wall modification-associated gene family which is closely related to PMEs, that is pectin methylesterase inhibitor protein (PMEI) gene family. Twenty-two PMEIs were expressed specifically in the fertile floral buds including counterparts of AtPMEI1 (Bra014099 and Bra032239) and AtPMEI2 (Bra021235). As for PGs family, most PGs had very low or undetectable expression during the early stages then went up dramatically at stage V. Among these 14 PGs, some genes showed extremely high expression at the last stage, such as Bra029683 and Bra033347. However, in the sterile floral buds, these PGs showed almost undetectable expression. In addition, five genes encoding PLLs, namely Bra012756, Bra017412, Bra017412, Bra008721 and Bra016700, expressed at an extremely high level at the mature pollen stage. And except for Bra016700, the other four genes presented a specific expression in the fertile floral buds (Fig. 8).
Transcription factors showed expressional changes in the sterile line According to the information of transcription factors in Brassica database (http://brassicadb.org/brad/) and Plant   Table S6). Most of these differently expressed TFs were constitutively expressed in both of the sterile and fertile lines, but 70 genes were expressed specifically in the fertile floral buds, whereas 18 genes were specific to the sterile floral buds. Among those 70 genes, 24 genes were exclusively expressed at stage II and they belonged to the NAC, AP2, MADS, C2H2, C3H, PHD, ABI3, ARF, TAZ and Aflin families (Fig. 9a). It is noteworthy that NAC TFs accounted for one-thirds of these 24 genes revealed the important roles that NAC TFs play in regulating early anther development. We also found 20 genes including three orthologs of Arabidopsis DAZ1, DAZ2 and DAZ3, exhibited remarkably high expression at mature   The values in the column of "WT/ms1" are expressed in terms of log2FC (fold change) and the values in other columns are expressed in terms of the ration of wild type to mutant. Shot dashes represent either no significant difference or no expression. a Genes showing expression change in spl and ems1 mutant [36] b Genes showing expression change in tdf1 mutant [15] c Genes showing expression change in ams mutant [16] d Genes showing expression change in ms1 mutant [37] pollen stage compared with the other four stages (Fig. 9c). Eighteen sterile floral buds-specific genes which belonged to the NAC, AP2, MYB, MADS, bZIP and zf-HD families were mainly expressed at stage IV and stage V, suggesting that they were up-regulated at later pollen development processes in the sterile line (Fig. 9d).
Pollen development-related genes were covalently modified by H3K4me 3 and H3K27me 3 To get more information of gene expression regulation in pollen development, we characterized the epigenetic control during pollen development by performing ChIP-seq on the fertile floral buds at mature pollen stage. The specific antibodies against trimethylated H3K4, a typical histone modification pattern characterizing active chromatin, and trimethylated H3K27 which is proposed to inhibit transcription, were used. A total of 13,008 and 8091 genes were found enriched for H3K4me 3 and H3K27me 3 , respectively. Combined with RNA-seq transcriptome analysis, we found that 433 genes enriched for H3K4me 3 and 750 genes enriched for H3K27me 3 were down-regulated in the sterile floral buds, while 151 genes enriched for H3K4me 3 and 143 genes enriched for H3K27me 3 were up-regulated in the sterile floral buds at mature pollen stage. And there were 47 down-regulated and 11 up-regulated genes covalently modified by both H3K4me 3 and H3K27me 3 (Fig. 10). GO classifications were obtained to investigate the functions of these genes (Additional file 2: Table S7). It was worth mentioning that enriched GO terms in the group of down-regulated genes marked with H3K27me 3 were associated with pollen tube development, cell tip growth, cell wall modification, pollination and reproductive developmental process suggesting an important role of H3K27me 3 during pollen tube growth (Additional file 6: Figure S6). The molecular functions of these targets concentrated mainly on ATP binding, metal ion binding, for example, calcium ion binding and hydrolase activity, especially pectinesterase activity and polygalacturonase activity (Additional file 7: Figure S7). Besides, 41 and 25 pollen development-related genes covalently modified by H3K4me 3 and H3K27me 3 were selected respectively, according to the published reports of Arabidopsis mutants affecting anther or pollen development as well as pollen germination or pollen tube growth ( Table 4). The majority of these genes were involved in pollen maturation, pollen germination and pollen tube growth, over-represented by the genes encoding calmodulin-binding proteins. Those late-expressed genes enriched for H3K4 me 3 participated in biological processes including glucose catabolic process, protein phosphorylation, nucleotide-sugar transport, oxylipin biosynthetic process and so on. While for late-expressed genes enriched for H3K27me 3 , they were mainly involved in auxin signaling transduction, calcium ion transport, membrane trafficking and regulation of transcription. Several early-expressed genes during anther and pollen development were also found to be marked with H3K4me 3 or H3K27me 3 , over-represented by genes required for pollen exine formation. It is noticed that H3K27me 3 targeted three key tapetum development-related genes including two transcription factor-encoding genes, BrMS1, BrAMS and Bra000615, the ortholog of Arabidopsis CEP1. Among these 66 genes, Bra014776 was the only example of a bivalent gene displayed enrichment for both H3K4me 3 and H3K27me 3 .

Discussion
In our previous analysis, we have concluded that the defect in male meiotic cytokinesis leads to male sterility in 'Bcajh97-01A' [34]. As it is well known, a characteristic feature of male meiotic cytokinesis is the abundant presence of callose. In this study, cytochemical staining for callose during the meiosis process further supported our result. Among 12 genes encoding putative callose synthase identified in Arabidopsis, GSL1, GSL5 and GSL2 have been shown to be involved in callose synthesis during microsporogenesis. Mutants of these genes display defective callose synthesis and exine patterning [59,91]. On the other hand, callase mixture secreted by the tapetum is responsible for callose degradation, represented by the A6 gene encoding β-1, 3-glucanase. Knockout or knockdown of β-1, 3-glucanase genes in Arabidopsis, B. napus and rice frequently caused defective callose degradation and ultimately led to male sterility [92,93]. In addition, AtMYB80 and CDM1 have been reported to play an important role in regulation of callose metabolism during microsporogenesis [13,94]. Mutation in CDM1 disturbs the normal callose metabolism by down-regulating of AtGSL5 and AtGSL2 and activating in advance of A6. At last, the cdm1 mutant is completely male sterile resulting from delayed meiotic cytokinesis and microspore degeneration. In our analysis, although ortholog of Arabidopsis CDM1 was dramatically down-regulated and the phenotype of pollen abortion was similar to that of cdm1, the regulation of callose metabolic pathways may be different as no one callose synthase gene was affected in the sterile line at early pollen development. In addition, all the expression of early-expressed β-1, 3-glucanase genes seemed to be down-regulated first and then remarkably up-regulated. However, the expression change of β-1, 3-glucanase genes was different in cdm1. Although Arabidopsis has 55 genes encoding putative β-1, 3-glucanases, only A6 encoding protein was thought to be part of the callose enzyme complex [95]. However, a recent research has speculated that three other β-1, 3-glucanase genes (At3g24330, At3g55780, and At3g61810) might participate in callose dissolution during microsporogenesis and pollen development [94]. In this study, we identified two more groups of β-1, 3-glucanase genes showing changed expression in the sterile floral buds. All seven early-expressed genes displayed similar expressional Fig. 9 Expression of transcription factors (TFs) showing fertile and sterile floral buds-specific features throughout anther development. The intensities of the colors increase with increasing expression levels, as indicated at the bottom. a the fertile floral buds-specific TFs exclusively expressed at tetrad stage. b the fertile floral buds-specific TFs highly expressed at ninucleate stage and binucleate stage. c the fertile floral buds-specific TFs mainly expressed at mature pollen stage. d the sterile floral buds-specific TFs changes in the sterile floral buds. And their Arabidopsis orthologs (At4g14080, At3g61810 and At3g23770) were closely related with each other according to phylogenetic analysis and belong to the expression group K which includes β-1, 3-glucanase genes highly specific to anther [95]. All these suggested they might play a similar role in degradation of callose wall surrounding the tetrads.
Another particular phenotype of 'Bcajh97-01A' , compared with cdm1, was the disorder development of tapetum. In cdm1, the tapetum layer appeared normal and genes known to be important for tapetum development were not dramatically affected, while the tapetum layer of 'Bcajh97-01A' went on premature PCD followed by abnormal development of tapetosomes and elaioplasts. And a series of regulatory genes for tapetum development were altered, including BrDYT1, BrAMS, BrMYB80 and BrMS1. In Brassicaceae species, the main pollen exine materials are stored and transferred by tapetosomes, the lipid accumulating ER-derived spherical organelles in tapetal cells and elaioplasts, the specialized plastid derived from proplastids [96][97][98]. In Arabidopsis mutants, ams and ms1, as well as a B. napus male sterile mutant, abnormal tapetal development caused decrease or even absence of tapetosomes and elaioplasts leading to disordered pollen wall formation [17,99,100]. Studies have shown that either premature or retarded tapetal PCD could induce defects in typical tapetosomes and elaioplasts formation [45,101]. Here, combined with our phenotype observation, we presumed that abnormal development of tapetosomes and elaioplasts caused by premature tapetal PCD led to defective pollen exine formation in the sterile 'Bcajh97-01A' line.
Unlike exine development which is mainly controlled by tapetum cell, intine formation has long been assumed to be mainly contributed by microspores. Although the exact mechanism is not clear, pectin metabolism-associated genes have been reported to be essential for pollen wall development [8]. As with other transcriptomic analyses related with male sterility, a high proportion of cell wall modification-related genes involved in polysaccharide metabolism, including PMEs, PMEIs, PGs and PLLs, were significantly down-regulated in the sterile line in our study [102][103][104]. In Arabidopsis, 14 of 66 putative PMEs are specifically expressed in pollen grains or pollen tubes [105,106]. Two of these PMEs, AtVDG1 and AtPPME1 have been reported to be essential for pollen tube elongation and stability [29,107]. And, AtPME48, the second most expressed PME in dry and imbibed pollen grains, was found to be closely associated with pollen grain germination by influencing remodeling the homogalacturonan of the intine cell wall [30]. Here in our study, a total of 13 PMEs that showed particularly high expression at mature pollen stage, including orthologs of AtVDG1 and AtPPME1, were down-regulated in the sterile floral buds. These genes are presumed to be involved in pollen maturation, pollen germination and pollen tube growth affecting development of pollen intine. Recently, one of these PMEs has been proved to function in intine construction by our lab [108]. As for those 22 PMEIs highly expressed in pollen, they are presumed to be key regulators of cell wall stability during pollen tube growth by interacting with PMEs, although more details about the molecular mechanism of interaction needs to be clarified [109][110][111][112]. In B. rapa, several PG genes have been demonstrated to regulate different aspects of pollen development in our previous studies, with BcMF2 affecting pollen intine, BcMF6 affecting microsporogenesis and pollen maturation, and BcMF9 affecting both intine and exine formation [113][114][115]. In fact, two down-regulated PGs found in our study, that were Bra032138 and Bra037005, have been confirmed to be pollen-expressed and the latter one was demonstrated to participate in pollen wall construction by modulating intine information [31,116]. The overlapping but distinct expression patterns for most of these cell wall hydrolytic enzyme-coding genes implied complicated roles they may play during the plant reproductive processes.  Transcriptional regulation plays pivotal roles in the control of gene expression in plants. In Arabidopsis, 612 of 1350 predicted TFs were expressed in developing male gametophytes and a number of TFs have been identified as key regulators during plant reproductive development as Khurana et al. summarized [117]. Here in our study, 493 predicted TFs showed changed expression level in the sterile line compared with the fertile line representing candidates for transcriptional regulators of male gametophyte development. The analysis of stage-specific expression profiling of these TFs during floral bud development helps to make a prediction of their gene functions. One of these TFs, Bra016531, a novel MYB gene was indicated to be related with tapetum development and function in our study.
Comparative transcriptome analysis of the sterile and fertile floral buds of 'Bcajh97-01A/B' at five typical Genes were selected on the basis of previous reports of Arabidopsis mutants affecting anther or pollen development as well as pollen germination or pollen tube growth. The Brassica ID in bold means the gene showed differential expressional levels between the fertile and sterile floral buds at mature pollen development in our RNA-seq analysis pollen developmental stages identified a great deal of differentially expressed genes. On the basis of advances in Arabidopsis, and combined with detailed phynotype observation, we make a summary of genes involved in various processes of pollen wall development in B. rapa (Fig. 11). Besides proteins functioning as transcriptional regulators, epigenetic modifications of DNA and associated proteins have been identified as triggers for gene activity. Histone lysine methylation has emerged as a critical player in the regulation of gene expression by controlling the chromatin state in many organisms [118,119]. The deposition of H3K4me 3 have been reported to be closely related with plant reproductive development [120,121]. But how H3K4me 3 and H3K27me 3 deposition participate in pollen development still remains largely unknown. In this study, a large number of H3K4me 3 or H3K27me 3 enriched DEGs and pollen development-related genes were identified providing materials for further analysis of epigenetic control on gene expression during pollen development. We have also noticed that both H3K4me 3 and H3K27me 3 targeted several key genes involved in tapetum development and pollen exine formation suggesting vital roles of H3K4me 3 and H3K27me 3 deposition in regulation of pollen wall construction.

Conclusions
In this study, transcriptome profiling on the floral buds of a GMS line with aberrant meiotic cytokinesis allowed the generating of stage-specific gene expression profiles and identification of candidate genes underlying pollen wall formation and genic male sterility. With the application of ChIP-sequencing on floral buds, a new sight into epigenetic control on gene expression during pollen development and pollen wall formation was provided. The data presented here provides a powerful platform for future functional and molecular regulation research of pollen wall formation and pollen development in B. rapa.

Morphological and cytological observation
Pollen grains were observed by a scanning electron microscope as previously described [123]. Aniline blue staining was performed as described [13] and DAPI staining was made as previously described [31]. Both the micrographs of callose fluorescence and chromatin fluorescence were captured with a Leica DMLB fluorescence microscope under UV light. For transmission electron microscopy (TEM) analysis, the procedures were as described in reference [123] and the samples for TEM observation were also used for semi-thin section observation.

RNA sequencing, transcriptome assembly and annotation of Unigenes
Total RNA was isolated from both the sterile and fertile floral buds at five pollen developmental stages according to the instructions of the TRIzol kit (Invitrogen, USA) and sent for RNA-seq. After enrichment using NEBNext Poly (A) mRNA Magnetic Isolation Module (NEB, E7490), mRNA was used to construct cDNA libraries using NEBNext mRNA Library Prep Master Mix Set for Illumina (NEB, E6110) and NEBNext Multiplex Oligos for Illumina (NEB, E7500). The high-quality libraries were then sequenced by Biomarker Technologies (http:// www.biomarker.com.cn) on the Illumina HiSeqTM2500. After filtering out low quality data, over 14,500,000 clean reads were assembled by running Trinity, resulting 25,509,284 contigs. Then all contigs were clustering to transcripts and Unigenes according to the pair-end information and similarities between contigs.

GO and pathway analysis
GO annotations of the Unigenes and DEGs were determined using agriGO [124]. Then, the results were submitted to WEGO to obtain the GO classification graph [125]. GO analysis applied a hypergeometric test to identify significantly enriched GO terms in DEGs in comparison to the transcriptome generated by this study. We chose the Bonferroni to do the multi-test correction, set 5 as minimum number of mapping entries, and we used a corrected-p value ≤0.05 as the threshold value. GO term (P ≤ 0.05) was defined as significantly enriched GO term. For pathway analysis, the DEGs were submitted to KEGG Automatics Annotation Server (KAAS) and classified with the single directional best hit (SBH) method [126]. Then, the results were submitted to KEGG Mapper to obtain the KEGG map.

qRT-PCR analysis
Total RNA was extracted from five different tissues including roots, stems, leaves, inflorescences and siliques using TRIzol Reagent (nvitrogen, USA) and reverse transcribed into the first strand of cDNA using PrimerScript RT reagent Kit (Takara, Dalian, China). For qRT-PCR analysis of Bra016531, BcUBC10 was used as the reference gene and qRT-PCR reaction was performed using the SYBR® Premix Ex Taq™ Kit (TaKaRa, Dalian, China) in a CFX96 Real-Time System (Bio-Rad, California, USA) with the gene-specific pairs (Additional file 2: Table S8). For each sample, three biological replicates were conducted with three technical replicates, and the results were calculated using the 2 -ΔΔCt method [127].

Promoter analysis
A 996-bp genomic DNA fragment containing nucleotide sequence from 976-bp upstream of the Bra016531 start codon to the first 20 bp of the first exon was amplified (Additional file 8: Figure S5) and subcloned upstream of the GUS reporter gene in the pBI101 vector (Dalian,