Abstract 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 H3K4me^3 and H3K27me^3 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. Keywords: Brassica rapa, Male sterility, Pollen, Cytokinesis, Tapetum, Pollen wall, Transcription factor, Transcriptome, ChIP-sequencing Background Anther and pollen development in plants is a complex process which includes a series of extraordinary events such as anther cell division and differentiation, male meiosis, microspores being released from the tetrads, pollen wall development as well as pollen maturation and anther dehiscence. Any abnormality during this process will lead to male sterility, a trait utilized for hybrid breeding and crop yield increases [[41]1]. As a special structure of mature pollen grain, the pollen wall is important for reproduction not only because it provides protection for male gametophytes, but more importantly, because it functions in male-female interaction, fertilization, and seed production [[42]2]. The complex multi-layered pollen wall displays a variety of surface morphologies but a generally similar fundamental structure, comprised of an outer exine and an inner intine [[43]3]. The exine, whose synthesis is regulated by both the sporophytic tapetum and the microspore [[44]4], is constructed primarily of sporopollenin, a robust biopolymer comprised predominantly of polyhydroxylated aliphatic compounds and phenolics [[45]5]. The intine, which is initiated during the early stages of male gametogenesis and is controlled gametophytically [[46]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 [[47]7–[48]9]. Recently, well-characterized genes in which mutations cause impaired exine and male sterile phenotype have enriched our understanding in pollen wall development [[49]4, [50]8, [51]10, [52]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 [[53]12–[54]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 [[55]18–[56]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 [[57]29–[58]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 [[59]33, [60]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. Results 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 [[61]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. [62]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. [63]1c and g). The 4′, 6-diamidino-2-phenylindole (DAPI) staining showed no nucleus formed in the globular remnants (Fig. [64]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 [[65]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. [66]2). Fig. 1. Fig. 1 [67]Open in a new tab Comparison of pollen grains of the sterile line ‘Bcajh97-01A’ (A line) and the fertile line ‘Bcajh97-01B’ (B line) of Brassica rapa. a, e Scanning electron micrographs of mature pollen grains. b, f Magnified images show details of the pollen morphology. c, g Pollen morphology under bright field microscopy. d, h Pollen grains stained with 4′, 6-diamidino-2-phenylindole (DAPI) solution Fig. 2. Fig. 2 [68]Open in a new tab Cytochemical staining for chromatin with DAPI and for callose with aniline blue during the meiosis process of the sterile line ‘Bcajh97-01A’ (A line) and the fertile line ‘Bcajh97-01B’ (B line) of Brassica rapa 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. [69]3). Fig. 3. Fig. 3 [70]Open in a new tab Analysis of tapetum development in the sterile line ‘Bcajh97-01A’ (A line) and the fertile line ‘Bcajh97-01B’ (B line) of Brassica rapa through semi-thin section and transmission electron microscopy observation. Stage 4, mother cell at pachytene/diplotene stage. Stage 5, mother cell at meiosis telophase I. Stage 6, tetrad stage. Stage 7, early uninucleate microspore stage. Stage 8, middle uninucleate stage. Stage 9, late uninucleate stage. Stage 10, early binucleate stage. Stage 11, middle binucleate stage. Stage 12, late binucleate stage. Scale bars in semi-thin sections of anthers: 50 μm. Scale bars in transmission electron microscopy observation: 4 μm 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 [71]1: Figure S1. After clustering, 207,932 transcripts and 72,168 Unigenes were obtained (Additional file [72]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 [73]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 [74]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 [75]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 [76]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. [77]4a and Additional file [78]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. [79]4b and c). And as growth progresses, the number of down-regulated genes increased dramatically (Fig. [80]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. [81]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. [82]4b). Among the 8288 DEGs, 519 genes showed opposite trend at different stages (Fig. [83]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. [84]4e), however, the expression reached peak at stage III in the sterile floral buds (Fig. [85]4f). This result demonstrated that the expression of these genes were delayed during the developmental process of the sterile floral buds. Fig. 4. [86]Fig. 4 [87]Open in a new tab Venn diagrams and line charts showing differentially expressed genes (DEGs) between the sterile line ‘Bcajh97-01A’ and the fertile line ‘Bcajh97-01B’. a A summary of the numbers of DEGs at five developmental stages (I-V). b A summary of the numbers of DEGs up-regulated (U) in the sterile floral buds. c A summary of the numbers of DEGs down-regulated (D) in the sterile floral buds. d A summary of the numbers of DEGs showing opposite trend (U_D) at different developmental stages. e The expression (RPKM value) of 217 DEGs in the fertile floral buds at five developmental stages (B1-B5). Log[10] (RPKM) equal to − 2 represents no expression. The red line represents the average expression. f The expression (RPKM value) of 217 DEGs in the sterile floral buds at five developmental stages (A1-A5). Log[10] (RPKM) equal to − 2 represents no expression. The red line represents the average expression 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. [88]5b and Additional file [89]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. [90]5c, Fig. [91]6 and Additional file [92]2: Table S4). In addition, pollen sperm cell differentiation were also enriched at stage IV and stage V (Additional file [93]2: Table S4). While at stage I, the enriched GO terms were different from the later four stages (Additional file [94]2: Table S4). For up-regulated DEGs, 69, 106, 55, 74 and 29 GO terms were significantly enriched at the five stages, respectively (Fig. [95]5a and Additional file 2: 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 [96]2: Table S4), indicating that these DEGs may relate to the phenotype of premature tapetum degradation in the sterile line. Fig. 5. [97]Fig. 5 [98]Open in a new tab The numbers of significantly enriched GO terms by DEGs between the sterile line ‘Bcajh97-01A’ and the fertile line ‘Bcajh97-01B’ and some common Biological Processes GO annotation in some developmental stages.a.The numbers of significantly enriched Biological Processes GO terms by DEGs up-regulated (U) in the sterile floral buds. b. The numbers of significantly enriched Biological Processes GO terms by DEGs down-regulated (D) in the sterile floral buds. c. Eight and 27 Biological Processes GO terms were all significantly enriched by DEGs down-regulated in the sterile floral buds at stage II to V and stage III to V Fig. 6. [99]Fig. 6 [100]Open in a new tab Some Biological Processes GO terms of DEGs with down-regulated levels in sterile buds at stage II to V. The color in each cell indicates FDR of the GO enrichment according to the scale shown, and blank cells indicate not significant 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 [101]2: Table S4), with the gene number increasing through development. Another enriched MF and CC GO terms were also showed in Additional file [102]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 [103]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 [104]1). Table 1. Callose metabolism-related genes showed dramaticaly altered expression Brassica rapa Arabidopsis Locus B1/A1 B2/A2 B3/A3 B4/A4 B5/A5 Locus Gene name Description Bra037213 – – 4.9 5.6 9.4 AT2G13680 GSL2 Callose synthase Bra004288 – 1.8 – – – AT1G68200 CDM1 Zinc finger C-×8-C-×5-C-×3-H type family protein Bra026644 – – – – 1.8 AT2G01630 – Glycosyl hydrolases family 17 protein Bra036718 – – – 2.2 – AT1G64760 – Glycosyl hydrolases family 17 protein Bra010330 – – – – 2.0 AT4G29360 – Glycosyl hydrolases family 17 protein Bra028988 – – – 2.4 – AT5G55180 – Glycosyl hydrolases family 17 protein Bra037795 – – 7.0 7.1 8.7 AT5G64790 – Glycosyl hydrolases family 17 protein Bra031901 – – – – 6.2 AT5G64790 – Glycosyl hydrolases family 17 protein Bra003475 – 1.1 – 2.9 – AT3G61810 – Glycosyl hydrolases family 17 protein Bra001918 −2.7 5.3 −6.2 – – AT3G23770 – Glycosyl hydrolases family 17 protein Bra019084 – 4.4 −4.5 – 2.3 AT4G26830 – Glycosyl hydrolases family 17 protein Bra037057 – – −6.0 – – AT4G14080 A6 Glycosyl hydrolases family 17 protein Bra032758 −2.1 1.2 −5.1 1.8 – AT4G14080 A6 Glycosyl hydrolases family 17 protein Bra020110 – – −7.1 – – AT5G20330 – Glycosyl hydrolases family 17 protein Bra014979 – – −4.3 – – AT3G23770 – Glycosyl hydrolases family 17 protein [105]Open in a new tab All values are expressed in terms of log2FC (fold change), so that positive values indicate depression of gene expression in male sterile line ‘Bcajh97-01A’. Shot dashes represent either no significant difference or no expression 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 [[106]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 [[107]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 STERILITY 2 (MS2), CYP703A2, CYP704B1, ACYL-COA SYNTHASE 5 (ACOS5), POLYKETIDE SYNTHASE A (PKSA), POLYKETIDE SYNTHASE B (PKSB) and TETRAKETIDE α-PYRONE REDUCTASE 2 (TKPR2) were dramatically affected in the sterile floral buds. It was worth mentioning that they showed a decrease in expression at stage II and then a significant up-regulation at stage III in the sterile floral buds, compared with the fertile ones. 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 [108]2). Table 2. Pollen exine formation-related genes showed changed expression Classification Brassica rapa Arabidopsis Locus B1/A1 B2/A2 B3/A3 B4/A4 B5/A5 Locus Gene name Gene expression pattern Lipid metabolism Bra038691 − 2.7 4.0 − 2.8 – – AT3G11980 MS2 tapetum shortly after release of microspore from tetrad Bra033272 −2.6 3.6 −6.1 – 5.6 AT1G01280 CYP703A2 tapetum and microspore specific Bra004386 −1.4 1.7 −4.0 – – AT1G69500 CYP704B1 anther-specific Phenolic metabolism Bra036646 −2.4 1.4 −8.5 – – AT1G62940 ACOS5 tapetum-specific Bra011566 −1.7 2.5 −8.1 – – AT4G34850 PKSA tapetum-specific Bra017147 −1.7 3.3 −6.5 – – AT1G02050 PKSB tapetum-specific Bra004316 −3.0 – −4.4 – – AT1G68540 TKPR2 tapetum-specific Transporters Bra039378 – 1.5 −3.2 – – AT3G13220 ABCG26 tapetum-specific Bra005048 – 1.4 −5.2 – – AT2G39350 ABCG1 tapetum-specific Bra026352 – – 2.3 – – AT4G27420 ABCG9 tapetum-specific Bra021598 – – 4.2 4.6 4.7 AT2G29940 ABCG31 tapetum-specific Bra000469 – 1.9 3.9 3.3 3.3 AT2G29940 ABCG31 tapetum-specific Bra014776 – 3.2 3.5 2.0 2.7 AT3G55090 ABCG16 tapetum-specific Bra038907 −1.8 – – – −8.9 At3G51590 LTP12 anther wall tapetum Bra028294 −2.8 – −6.8 – – AT5G52160 – flower stage 9 Bra009282 – – −2.1 – – AT5G07230 – flower stage 9,10,11 Bra012819 −3.9 – −8.2 – – AT3G52130 – flower stage 9 Transcription factors Bra013519 – −3.0 4.3 5.3 – AT4G21330 DYT1 highly expressed in the tapetum, meiocytes and microspores Bra002004 – – −1.4 1.9 AT2G16910 AMS tapetum Bra013041 – – −1.2 – – AT2G16910 AMS tapetum Bra002847 – 3.8 −6.1 – – AT5G56110 MYB80 tapetum, middle layers and developing microspores Bra035604 – 2.2 −8.9 – – AT5G56110 MYB80 tapetum, middle layers and developing microspores Bra002401 – 5.3 −3.7 – – AT5G22260 MS1 tapetum shortly after microspore release from tetrad Bra004689 – – 1.7 – – AT2G42380 bZIP34 anthers and pistils [109]Open in a new tab All values are expressed in terms of log2FC (fold change), so that positive values indicate depression of gene expression in male sterile line ‘Bcajh97-01A’. Shot dashes represent either no significant difference or no expression 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 ([110]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 [111]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 function-known tapetum-related genes during distinct flower developmental stages in the ms1 mutant. Table 3. Tapetum development and function-related genes Classification Brassica rapa Locus Arabidopsis Locus Description Arabidopsis microarray data WT/ems1^a WT/spl1^a WT/tdf1^b WT/ams^c WT/ms1^d Known genes Bra002847 AT5G56110 AtMYB80 2.2 2.5 – – −2.2 Bra035604 AT5G56110 Bra002401 AT5G22260 MALE STERILITY 1 – – – – −2.7 Bra039378 AT3G13220 ABCG26 52.2 56.7 2.2 5.9 −1.6 Bra002004 AT2G16910 AMS 31.8 28.8 3.3 – −2.1 Bra013041 AT2G16910 Bra033272 AT1G01280 CYP703 47.9 43.6 12.0 9.9 −2.9 Bra004386 AT1G69500 CYP704B1 117.5 127.8 8.6 11.6 −2.7 Bra038691 AT3G11980 MALE STERILITY 2 42.8 50.8 31.7 17.5 −2.2 Bra017147 AT1G02050 LESS ADHESIVE POLLEN 6 17.4 14.8 6.2 6.5 – Bra011566 AT4G34850 LESS ADHESIVE POLLEN 5 – – 22.4 8.7 – Bra036646 AT1G62940 ACYL-COA SYNTHETASE 5 53.5 93.4 18.4 6.9 −4.0 Taptum specific genes Bra038803 AT4G20420 tapetum specific protein TAP35/TAP44 33.9 53.2 11.7 – −7.2 Bra020920 AT4G20420 Bra028324 AT3G42960 ATA1/TAPETUM 1 61.4 106.6 2.9 8.2 −1.6 Carbonhydrate metabolism and transport Bra020859 AT4G22080 pectin lyase-like superfamily protein – 7.5 – – 3.4 Bra013442 AT4G20050 QRT3 34.3 37.3 5.3 6.3 −1.5 Bra023600 AT5G17200 pectin lyase-like superfamily protein – 1.8 – – 2.1 Bra001918 AT3G23770 Glycosyl hydrolases family 17 protein 28.5 25.0 4.9 – 4.8 Bra032758 AT4G14080 Glycosyl hydrolases family 17 protein 82.9 204.5 54.0 8.6 −3.9 Bra022636 AT5G53190 nodulin MtN3 family protein 37.2 54.3 – 1.2 1.1 Oxidation reduction Bra015740 AT1G76470 NAD(P)-binding Rossmann-fold superfamily protein 18.7 18.0 11.4 13.1 22.7 Bra021337 AT3G57620 glyoxal oxidase-related protein 7.7 8.2 – – 2.4 Bra029151 AT5G51950 glucose-methanol-choline oxidoreductase family protein 12.2 14.8 – 3.6 −1.6 Bra022571 AT5G51950 Bra028286 AT5G51950 The others Bra020777 AT3G06100 NOD26-LIKE INTRINSIC PROTEIN 7 8.5 12.2 7.5 5.4 −2.1 Bra037505 AT5G48210 prolamin-like protein (DUF1278) 5.2 6.7 7.3 – 9.6 Bra020712 AT5G48210 Bra015870 AT1G75050 pathogenesis-related thaumatin superfamily protein 3.8 6.1 – 4.6 1.3 Bra032134 AT2G23945 eukaryotic aspartyl protease family protein – – – – – Bra011161 AT4G30030 eukaryotic aspartyl protease family protein – – – – – Bra009795 AT5G24820 eukaryotic aspartyl protease family protein 49.4 62.3 – – – Bra030647 AT1G06280 LOB DOMAIN-CONTAINING PROTEIN 2 – 18.4 19.5 – 8.2 Bra018697 AT1G47980 desiccation-like protein 19.8 17.5 32.6 17.7 6.3 Bra018660 AT1G08065 ALPHA CARBONIC ANHYDRASE 5 26.9 24.5 4.9 - - Bra016531 AT1G18960 myb-like transcriptional regulator family protein - 7.5 9.6 - −4.2 [112]Open in a new tab 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. ^aGenes showing expression change in spl and ems1 mutant [[113]36] ^bGenes showing expression change in tdf1 mutant [[114]15] ^cGenes showing expression change in ams mutant [[115]16] ^dGenes showing expression change in ms1 mutant [[116]37] 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. [117]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. [118]7b). Then, Bra016531 promoter-driven GUS expression was examined in anther cells using thin sections taken from floral buds at various developmental stages (Fig. [119]7c-h). Blue GUS signal was first detected both in sporocyte and tapetal cells during meiosis (Fig. [120]7d), and then in tetrads and tapetal cells (Fig. [121]7e). Later, the signal was detected mainly in tapetal cells at the early vacuolate stage (Fig. [122]7f). These temporal and spatial patterns of Bra016531 expression supports our hypothesis that Bra016531 is associated with tapetum development and function during pollen development. Fig. 7. [123]Fig. 7 [124]Open in a new tab Expression pattern analysis of Bra016531. a Quantitative RT-PCR analysis of Bra016531 transcripts in different tissues of Brassica rapa: roots (R), stems (Ste), leaves (L), inflorescences (Inf) and siliques (Si). b Bra016531 promoter-GUS activity. (c-h) Thin sections of anthers from Bra016531pro::GUS-expressing plants. c.pollen mother cell stage. d.meiosis stage. e.tetrad stage. f.early vacuolated stage.g.tapetum degeneration stage. h.mature pollen stage. Scale bars, 20 μm 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. [125]8). Fig. 8. [126]Fig. 8 [127]Open in a new tab Hierarchical cluster display of the differentially expressed cell wall hydrolytic enzyme-coding genes in Brassica rapa. The color scale bar shown under the cluster indicates the maximum and minimum brightness values that represent the values of log[2] (RPKM) Transcription factors showed expressional changes in the sterile line According to the information of transcription factors in Brassica database ([128]http://brassicadb.org/brad/) and Plant Transcription Factor Database ([129]http://planttfdb.cbi.pku.edu.cn/) as well as annotations of Unigenes in transcriptome analysis, we identified 2567 Unigenes that encode putative transcription factors (TFs). These TFs were divided into 65 gene families. Different percentages of Unigenes in 53 gene families displayed changed expression patterns in the sterile floral buds compared with the fertile floral buds. Among those over-represented were the HSF family (38%), the LBD family (36%), the NAC family (33%), the MADS family (32%), the MYB-related family (28%), the AP2 family (23%), the C2H2 family (22%) and the BHLH family (21%). In contrast, the E2F (8%), Aflin (7%), GRF (6%) and FHA (5%) families were all under-represented. It’s worth mentioning that DEGs belonging to the NAC (48), BHLH (41), AP2 (40), MYB (32), MADS (32), C2H2 (27) and WRKY (26) families accounted for half of all 493 putative TFs that with changed expression patterns (Additional file [130]2: 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. [131]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 pollen stage compared with the other four stages (Fig. [132]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. [133]9d). Fig. 9. [134]Fig. 9 [135]Open in a new tab 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 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. [136]10). GO classifications were obtained to investigate the functions of these genes (Additional file [137]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 [138]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 [139]7: Figure S7). Fig. 10. Fig. 10 [140]Open in a new tab Venn diagram detailing DEGs enriched for H3K4me^3 or H3K27me^3 at mature pollen stage 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 [141]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. Table 4. Pollen development-related genes covalently modified by H3K4me^3 and H3K27me^3 Classification Brassica rapa Type Arabidopsis Locus Fold enrichment FDR (%) Locus Gene name Description References