Next Article in Journal
Inhibitory Effect of Opuntia humifusa Fruit Water Extract on Solar Ultraviolet-Induced MMP-1 Expression
Next Article in Special Issue
Genome-Wide Association Studies of 39 Seed Yield-Related Traits in Sesame (Sesamum indicum L.)
Previous Article in Journal
Co-Expression Network Analysis Identifies miRNA–mRNA Networks Potentially Regulating Milk Traits and Blood Metabolites
Previous Article in Special Issue
Genome-Wide Characterization and Analysis of Metallothionein Family Genes That Function in Metal Stress Tolerance in Brassica napus L.
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
IJMS
Volume 19
Issue 9
10.3390/ijms19092501
ijms-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Normal and Abortive Buds Transcriptomic Profiling of Broccoli ogu Cytoplasmic Male Sterile Line and Its Maintainer

by
Jinshuai Shu
,
Lili Zhang
,
Yumei Liu
*,
Zhansheng Li
,
Zhiyuan Fang
,
Limei Yang
,
Mu Zhuang
,
Yangyong Zhang
and
Honghao Lv
Institute of Vegetables and Flowers, Chinese Academy of Agricultural Sciences, Key Laboratory of Biology and Genetic Improvement of Horticultural Crops, Ministry of Agriculture, 12 Zhongguancun Nandajie Street, Beijing 100081, China
*
Author to whom correspondence should be addressed.
Int. J. Mol. Sci. 2018, 19(9), 2501; https://0-doi-org.brum.beds.ac.uk/10.3390/ijms19092501
Submission received: 30 June 2018 / Revised: 12 August 2018 / Accepted: 14 August 2018 / Published: 24 August 2018
(This article belongs to the Special Issue Plant Genetics and Molecular Breeding)
Browse Figures
Versions Notes

Abstract

:
Bud abortion is the main factor affecting hybrid seeds’ yield during broccoli cross breeding when using ogura cytoplasmic male sterile (ogu CMS) lines. However, the genes associated with bud abortion are poorly understood. We applied RNA sequencing to analyze the transcriptomes of normal and abortive buds of broccoli maintainer and ogu CMS lines. Functional analysis showed that among the 54,753 annotated unigenes obtained, 74 and 21 differentially expressed genes in common were upregulated and downregulated in ogu CMS abortive buds compared with ogu CMS normal buds, maintainer normal, and abortive buds, respectively. Nineteen of the common differentially expressed genes were enriched by GO terms associated with glycosyl hydrolases, reactive oxygen species scavenging, inhibitor, and protein degradation. Ethylene-responsive transcription factor 115 and transcriptional factor basic helix-loop-helix 137 were significantly upregulated; transcription factors DUO1 and PosF21/RF2a/BZIP34 were downregulated in ogu CMS abortive buds compared with the other groups. Genes related to polygalacturonase metabolism, glycosyl hydrolases, oxidation reduction process, phenylalanine metabolism, and phenylpropanoid biosynthesis were significantly changed in ogu CMS abortive buds. Our results increase our understanding of bud abortion, provide a valuable resource for further functional characterization of ogu CMS during bud abortion, and will aid in future cross breeding of Brassica crops.

Graphical Abstract

1. Introduction

Bud abortion is a very common biological phenomenon in Brassica species. During bud abortion, the buds stop growing, convert from green into yellow progressively from the base to the top of the bud, and eventually wither off before flowering [1,2], which is detrimental for hybrid seeds yield during cross breeding. In recent years, bud abortion has attracted considerable attention in Brassica species cross breeding [3,4,5], as it reduces the hybrid seed yield production of many Brassica crops and also affects the efficiency of Brassica crops cross breeding. In the past decade, to solve the bud abortion problem, many studies have investigated the factors that cause bud abortion [1,2,3,4,5,6,7,8]. However, most of the studies focused on morphology and anatomy, and few studies have reported the factors contributing to Brassica bud abortion at the molecular level. Therefore, the fundamental mechanism of bud abortion is poorly understood.
High-throughput sequencing technologies have played important roles in revealing the molecular mechanisms of various organismal biological processes. RNA sequencing (RNA-Seq) technology is an important high throughput sequencing technology that produces functional genomic data for non-model plants that lack genomic sequence data. RNA-Seq technology has been widely applied to assist in determining differentially expressed genes (DEGs) involved in different biological processes in many species and may be a promising method to address the genes associated with bud abortion in Brassica species. Surprisingly, no studies on bud abortion using RNA-Seq technology have been documented so far.
Brassica species display obvious heterosis, and using cytoplasmic male sterile (CMS) lines to produce hybrid seeds is an important approach that uses this heterosis. Broccoli (Brassica oleracea var. italica) is an important vegetable crop and plays a vital role in the vegetable production industry worldwide. More importantly, broccoli is also used to produce health products and develop anti-cancer drugs [9,10,11]. Thus, broccoli is the most important economical Brassica vegetable crop. However, a large problem exists in broccoli cross breeding: the hybrid seed purity cannot reach 100% when using inbred lines or self-incompatible lines to produce hybrid seeds. To improve the purity of broccoli hybrid seeds, many breeders execute wide crosses between broccoli cultivars and other Brassica CMS materials to transfer CMS genes. Our group used multiple broccoli inbred lines to with the cabbage CMS material ogura CMSR3629 (ogu CMS), which was introduced by the Asgrow Seed Co. (USA), which has been used to breed a plurality of broccoli CMS lines with a 100% sterile rate. However, the obtained CMS lines often show serious bud abortion, which leads to delayed flowering time and sharply reduced amounts of flowers, some serious bud abortion CMS lines cannot flower [5]. These problems limit the use of the CMS lines considerably and have reduced their breeding efficiency dramatically. So far, several studies have focused on the causes of bud abortion, e.g., the RsVPE1 gene encodes a vacuolar processing enzyme that is involved in radish floral bud abortion under heat stress [1]. The expression of stress response, energy metabolism, amino acid synthesizing and processing, signal transduction, disease resistance and senescence, transcription and translation, and transmembrane transport-related genes were different between normal and abortive buds, as detected by cDNA-amplified length polymorphism technology [6,7,8]. However, as far as we know, no study has investigated the molecular events occurring in bud abortion based on high throughput sequencing during broccoli cross breeding, or any other crop, thus little information is available concerning the genes involved in bud abortion. Thus, it is necessary to determine the elements involved in broccoli bud abortion at the molecular level, which will increase our understanding of the molecular events involved in bud abortion of broccoli and other crops. The objectives of this study were to identify the functional genes involved in bud abortion and to determine the genes expression characteristics related to broccoli bud abortion, based on RNA-Seq technology.

2. Results

2.1. High-Throughput Transcriptome Sequencing and Unigene Assembly

To explore gene expression and gene networks that control bud abortion of ogu CMS lines in broccoli, we performed RNA-Seq analysis of normal and abortive buds from ogu CMS line CMS93219 and its maintainer line ML93219. The cytoplasmic male sterility original of CMS93219 was cabbage CMS material ogu CMSR3629 and the backcross paternal was broccoli ML93219. CMS93219 was stabilized via sixteen generations of backcrossing before this study. There were no visible differences between CMS93219 and ML93219 in terms of the plants and the heading morphologies of the harvested material, except that CMS93219 showed significant bud abortion and ML93219 showed only slight bud abortion before or at the early anthesis stage (Figure 1A–D). The average number of abortive buds per branch of CMS93219 was 35.73 throughout the whole growth phase, accounting for 83.88% of the total bud number, which is extremely significantly higher than the 19.60 per branch (33.33%) in ML93219 (t-test, p < 0.01) (Figure 1E,F). Neither the abortive nor the normal buds of CMS93219 contained pollen, while the aborted buds of ML93219 contained less pollen than normal buds. Therefore, bud abortion and pollen abortion character were closely connected between ML93219 and CMS93219, and pollen abortion may promote bud abortion. Gene expression changes are usually associated with morphologic changes; therefore, we chose the maximum buds stage (around 3 days before anthesis) for transcriptome analyses (Figure 1B-I,B-II,D-I,D-II).
Twelve cDNA libraries from normal buds and abortive buds of the ML93219 and CMS93219, containing three biological replicates for each sample, were subjected to Illumina sequencing. After filtering invalid reads and data cleaning, 663,701,918 clean reads containing 99,555,287,700 nucleotides (99.6 Gb), with a mean length of 150 bp, were obtained. The Q20 and GC percentage were 92.93–94.33% and 46.46–47.12%, respectively (Table 1).
After assembling the clean reads, 97,347 transcripts (91,137,323 nucleotides) and 66,050 unigenes (51,896,834 nucleotides) were obtained. The average length of the transcripts and unigenes were 936 and 786 bp, respectively (Table 2).

2.2. Gene Annotation and Functional Classification

A total of 54,753 unigenes (82.90%) were annotated in at least one database and 6070 unigenes (9.19%) were annotated in seven databases. In the NCBI non-redundant protein sequences (Nr), NCBI non-redundant nucleotide sequences (Nt), KEGG Ortholog (KO), a manually annotated and reviewed protein sequence database (SwissProt), protein family (Pfam), Gene Ontology (GO), eukaryotic ortholog groups (KOG) databases, 44,294 (67.06%), 50,157 (75.93%), 12,403 (18.77%), 29,861 (45.2%), 24,660 (37.33%), 30,285 (45.85%) and 12,492 (18.91%) unigenes were annotated, respectively (Table 3).
We used the KOG functional annotation system for the assembled unigenes to obtain clues to the functions of the unigenes. As a result, 12,492 unigenes have defined, specific protein functions, accounting for 22.82% of the total annotated unigenes and involving 26 KOG functional classes. The five largest categories were “General function prediction only”, “Posttranslational modification, protein turnover, chaperones”, “Signal transduction mechanisms”, “Transcription”, “Intracellular trafficking, secretion, and vesicular transport”. Two categories, “Unnamed protein” and “Cell motility”, accounted for the lowest percentages (Figure 2).
To classify the functions of the unigenes, GO assignments were used based on the annotations from the Nr and Pfam databases. The results showed that 30,285 unigenes could be categorized into 57 functional groups, which were separated into three main categories: biological process (25 subcategories), cellular component (18 subcategories), and molecular function (14 subcategories) (Figure 3). The “cellular process,” “metabolic process,” and “single-organism process” were the major subcategories of the biological processes category; “cell” and “cell part,” “binding,” and “catalytic activity” were remarkable in the cellular component and molecular function categories, respectively. The classification result indicated that biological processes play a notable role during broccoli bud development, while the terms “biological phase,” “cell aggregation,” “synapse,” “synapse part,” “receptor regulator activity,” and “metallochaperone activity” were rare.
To predict the activated biochemical pathways in broccoli buds, the unigenes were annotated against the KEGG database (e-value = 1 × 10−10). A total of 264 KEGG pathways represented by 12,403 unigenes were obtained (Table S1). The main pathways were “Biosynthesis of amino acids [ko01230] (454 unigenes),” “Plant hormone signal transduction [ko04075]” (447 unigenes), “Carbon metabolism [ko01200]” (433 unigenes), and “Starch and sucrose metabolism [ko00500]” (366 unigenes). Moreover, the pathways such as “Endocytosis [ko04144],” “RNA degradation [ko03018],” “ubiquitin mediated proteolysis [ko04120],” “Fatty acid degradation [ko00071],” “Valine, leucine and isoleucine degradation [ko00280],” “Proteasome [ko03050],” “Lysine degradation [ko00310],” and “Apoptosis [ko04210],” which contained 235, 187, 184, 82, 79, 77, 49, and 31 unigenes, respectively, which indicated that protein degradation and cell death were associated with bud abortion.

2.3. Genes Express Differences and DEGs Clustering

The 97,347 assembly unigenes sequences obtained from de novo assemble of a merged set of 663,701,918 clean reads were used as the reference. We acquired the number of readcounts that aligned to each unigene and calculated the expected number of Fragments Per Kilobase of transcript sequence per Millions base pairs sequenced (FPKM) value to estimate the expression levels of the unigenes. The FPKM distribution and heatmap indicated the presence of many DEGs among the four samples, and the expression patterns of aborted buds of ML93219 (ML_AB) and aborted buds of ogu CMS93219 (CMS_AB), and normal buds of ML93219 (ML_NB) and normal buds of ogu CMS93219 (CMS_NB) were more similar, respectively (Figure 4 and Figure 5). Changes in the DEGs indicated that the abundances of transcripts were different at the same stage of normal and abortive buds development between ML93219 and CMS93219.

2.4. DEGs in Normal and Abortive Buds of ogu CMS

To explore the reference sample CMS_AB gene expression levels, the gene expression variations were determined using three comparisons, between CMS_AB and CMS_NB, between CMS_AB and ML_NB and between CMS_AB and ML_AB. Compared with CMS_NB, there were 6575 and 5482 genes up- and downregulated in CMS_AB, respectively. Compared with ML_NB, in CMS_AB, there were 6192 and 6321 genes that were up- and downregulated, respectively. Compared with ML_AB, there were 182 and 825 genes that were up- and downregulated in CMS_AB, respectively (Figure 6a–c).
To determine the key genes involved in ogu CMS line bud abortion process, we compared the DEGs of normal and abortive buds from ML93219 and CMS93219 to define the differentially expressed genes in common. In CMS_AB, 74 and 21 genes were upregulated and downregulated, respectively, compared with CMS_NB, ML_NB and ML_AB (Figure 7; Tables S2–S5), which are most likely associated with ogu CMS bud abortion in broccoli.

2.5. Validation of RNA-Seq Data by Quantitative Real-Time Reverse Transcription PCR (qRT-PCR)

We performed qRT-PCR to verify the DEGs identified by RNA-Seq, using the same samples that were used for the RNA-Seq analysis. Among the 21 randomly selected DEGs based on the expression level fold-change and differences in expression, 11 genes displayed higher expression quantification and 10 genes showed lower expression quantification in CMS_AB, compared with CMS_NB, ML_NB, and ML_AB. All 21 genes showed the same expression patterns in qRT-PCR as they did in the RNA-Seq databases experiment, indicating the high reliability of our RNA-Seq data (Figure 8).

2.6. Ogu CMS Bud Abortion-Related Genes in Broccoli

To further understand the function and biological process of the DEGs, GO term (corrected p < 0.05) and KEGG pathway (corrected p < 0.05) enrichment were performed to analyse the 95 differentially coexpressed genes. Nineteen were enriched and the most significantly enriched GO terms were “cell wall organization” and “external encapsulating structure organization” (corrected p = 3.26 × 10−5) in the biological process (BP) group, “extracellular region” (corrected p = 0.003937) in the cellular component (CC) group and “polygalacturonase activity” (corrected p = 0.015471) in the molecular function (MF) group (Table 4). The GO terms of “cell wall organization or biogenesis” and “extracellular region” contained 10 and 13 genes, respectively, and were the biggest categories in the biological process and cellular component groups (Table 4), respectively. Only “polygalacturonase activity” GO terms in molecular function group was enriched and contained 4 genes (Table 4). Furthermore, there were 8 same genes were enriched in “sucrose metabolic process,” “starch metabolic process,” “disaccharide metabolic process,” “cellular glucan metabolic process,” “glucan metabolic process” and “oligosaccharide metabolic process” (Table 4); 5 same genes were enriched between “external encapsulating structure” and “cell wall” (Table 4). “Phenylalanine metabolism” and “Phenylpropanoid biosynthesis” (K00430) contained the same three genes (Gene ID: c20440_g1, c20661_g1, c9433_g1) that were significantly enriched in KEGG pathways; meanwhile, the three genes were enriched by GO enrichment (Table 4 and Table 5).
Based on the enrichment results and functional annotation of the 19 enriched differential expressed genes in common by GO term, we found that “Glycosyl hydrolases” and “Reactive oxygen species (ROS) scavenger” related genes accounted for a high proportion (Figure 9). In addition, “Inhibitor,” “Plant defense,” and “Cell division and expansion” related genes were also significantly upregulated in CMS_AB and “Transporter” related gene Kinesin-4 downregulated in CMS_AB (Figure 9). Furthermore, five of the differentially expressed genes in common (GI: c11939_g1, c14539_g1, c25960_g1, c52977_g1, c9433_g1) enriched by GO terms were verified by qRT-PCR, and had significantly different expression between the abortive buds and normal buds of CMS93219 (Figure 8). These results suggested that most of the differentially expressed genes in common found in our study are required for ogu CMS bud abortion in broccoli.

2.7. Transcription Factors Are Involved in Broccoli ogu CMS Bud Abortion Control

Among the 95 common differentially expressed genes, four transcription factors were identified, including predicted transcription factor basic helix-loop-helix (bHLH) 137 (GI: c13818_g1), ethylene-responsive transcription factor (ERF) 115 (GI: c21802_g2), transcription factor DUO1 (GI: c18372_g1), and PosF21/RF2a/BZIP34 (GI: c14038_g1), which were distributed in four different gene families. Compared with CMS_NB, ML_NB, and ML_AB, transcription factors bHLH137 ERF115 were upregulated in CMS_AB, suggesting that these transcription factors may function as positive regulators in ogu CMS bud abortion in broccoli. Conversely, transcription factors DUO1 and PosF21/RF2a/BZIP34 were downregulated in CMS_AB, suggesting that they are negative regulators. Moreover, the expression quantifications of the four transcription factors were confirmed by qRT-PCR (Figure 8), suggesting that the four transcription factors play key roles in ogu CMS bud abortion in broccoli.

3. Discussion

The bud abortion phenomenon is a very complex bioprocess and demands many molecular events during bud development. In this study, we used RNA-Seq technology to explore the genes involved in ogu CMS bud abortion process and to provide a comprehensive analysis of the genes involved in ogu CMS bud abortion control in broccoli. Compared with CMS_NB, ML_NB, ML_AB, there were 6575, 182, and 6192 genes that were upregulated and 5482, 825 and 6321 genes that were downregulated in CMS_AB, respectively (Figure 6a–c), among which 74 genes were significantly upregulated and 21 genes were significantly downregulated equally in CMS_AB with serious bud abortion (Figure 7; Tables S2–S5). qRT-PCR proved that our RNA-Seq data was highly reliable (Figure 8). Functional categories of the common differentially expressed genes by GO term enrichment analysis showed that gene associated with cell wall composition and metabolism, such as cell wall organization, cell wall organization or biogenesis; extracellular region, such as external encapsulating structure organization; sugar metabolism, such as sucrose metabolic process, starch metabolic process, disaccharide metabolic process, cellular glucan metabolic process, glucan metabolic process, oligosaccharide metabolic process related genes were strongly induced in the buds abortion of broccoli (Table 4). These results indicated that bud abortion may be closely related to cell wall organization, external encapsulating structures and sugar metabolism related genes play an important role in bud abortion in broccoli. Our results were consistent with the differential expression genes obtained by cDNA-AFLP technique in radish and cabbage between normal bud and dead bud [6,7,8].

3.1. Genes Related to Programmed Cell Death Are Involved in Bud Abortion

Programmed cell death (PCD) is an important physiological process in single cells and penetrates the whole plant life cycle, which can help plants to control and organize the destruction of non-functional or redundant damaged cells [13,14,15]. Although PCD is a natural result of ageing, it may be switched on by environmental stress or irregular development in plants [15]. Increases in caspase-like proteases [16] and ROS [17,18] activities of metacaspase gene family related genes [19] are the important features of PCD. In this study, several genes associated with PCD were significantly changed in ogu CMS abortive buds compared with normal buds and the maintainer abortive buds (CMS_NB, ML_NB and ML_AB). Among the 8833 differential expression genes between ML_AB vs. ML_NB and CMS_AB vs. CMS_NB, there was 11.75% participate in the redox process, 11.08% with redox enzyme activity and 55.38% with catalytic activity, indicating that redox process involved in bud abortion in broccoli. ROS scavengers related genes, such as: peroxidase 27-like and peroxidase 45, l-lactate dehydrogenase, laccase-5-like and fatty acyl-CoA reductase 7 (Table S5), which were all upregulated in abortive buds compared with normal buds, suggesting that higher levels of ROS were produced in abortive buds and ROS scavenger-related genes were closely related with broccoli bud abortion. Moreover, 13 caspase-like and metacaspase activity genes involved in cell apoptosis were also discovered, such as metacaspase-1, metacaspase-3, metacaspase-5, metacaspase-6, metacaspase-7, metacaspase-9, pentapeptide repeats protein, caspase recruitment domain, BTB/POZ domain-containing protein POB1. The expression of these genes in abortion buds was significantly higher than that of normal buds suggesting that bud abortion was related with PCD.

3.2. Glycosyl Hydrolases, Inhibitors and Plant Defence Related Genes Are Implicated in Bud Abortion

In this study, several glycosyl hydrolase-related genes were determined as significantly upregulated in ogu CMS abortive buds, such as endoglucanase 20, endoglucanase 15-like, polygalacturonase ADPG2-like and polygalacturonase-like. However, polygalacturonase plays an important role during the life cycle of cell separation, being involved in cell wall modification, abscission and dehiscence in Arabidopsis thaliana [20,21], and endoglucanase is involved in cell wall biogenesis or degradation, cellulose degradation and polysaccharide degradation [22], which suggested that glycosyl hydrolase-related genes were required for ogu CMS bud abortion in broccoli. Pectinesterase/pectinesterase inhibitors modify cell walls via demethylesterification of cell wall pectin, which negatively regulates catalytic activity. As a voltage-gated inward-rectifying Ca2+ channel (VDCC) across the vacuole membrane, the calcium channel is an essential components of the slow vacuolar (SV) channel and is the major ROS-responsive Ca2+ channel, which is the possible target of Al-dependent inhibition and is involved in the regulation of stomatal movement [23,24]. In our study, pectinesterase/pectinesterase inhibitor 6 and 54 probable and calcium channel inhibitor-related genes were highly expressed in ogu CMS abortive buds, indicating that these genes may be important in controlling ogu CMS bud abortion. Moreover, E3 ubiquitin-protein ligase and peptidoglycan-binding LysM domain-containing protein both play important roles in the plant defence response [23,25,26,27,28,29,30] and nicotianamine synthase is involved in the cellular response to ethylene stimulus [31]: these genes were distinguished significantly upregulated in ogu CMS abortive buds, suggesting their likely involvement in ogu CMS bud abortion.

3.3. Transcription Factors Associated with Bud Development

Transcription factors are critical to regulate gene expression during plant development and in response to biotic and abiotic stresses [32,33,34,35,36,37,38]. In the present study, predicted transcription factors bHLH137, ERF115, DUO1 and PosF21/RF2a/BZIP34 were obviously changed in abortive broccoli buds compared with normal buds. Interestingly, the genes encoding stress-responsive transcription factors ERF115 and gibberellin-responsive transcription factors bHLH137 showed similar expression patterns: both were significantly upregulated in abortive buds (ML_AB and CMS_AB) compared with normal buds (ML_NB and CMS_NB), and in ogu CMS abortive buds (CMS_AB), they were significantly upregulated compared with maintainer abortive buds (ML_AB) (Figure 8). Conversely, the genes encoding transcription factors DUO1 and PosF21/RF2a/BZIP34, which are positive regulators of transcription showed similar expression patterns: both were significantly downregulated in abortive buds, and in ogu CMS abortive buds they were significantly downregulated compared with maintainer abortive buds (Figure 8). ERF 115 as a transcriptional activator of the phytosulfokine PSK5 peptide hormone family that binds to the GCC-box pathogenesis-related promoter element and limit quiescent center cell division activity when surrounding stem cells are damaged and is also a proteolytic target of the APC/C-FZR1 complex [39]. bHLH transcription factors belong to a family of transcriptional regulators and have a range of different roles in plant cells and tissue development, as well as in plant metabolism [40]. Transcription factor DUO1 could be involved in pollen sperm cell differentiation [41]. PosF21/RF2a/BZIP34 is a transcription factor with an activatory role [42,43], which might be involved in the sporophytic control of cell wall patterning and gametophytic control of pollen development, and play a role in the control of metabolic pathways regulating cellular transport and lipid metabolism [44,45]. Therefore, transcription factors play a key role in the complex regulatory networks of ogu CMS bud abortion.

3.4. Molecular Mechanisms Associated with CMS

Previous studies have shown that the genes involved in reactive oxygen species (ROS) homeostasis or antioxidative system balance may be important factors contributing to pollen abortion in cotton [46] and wheat [47]. Carbohydrate and energy metabolisms, oxidation-reduction system and phenylpropanoid metabolism pathways related genes may be important factor for CMS in soybean [48], rapeseed [49], cabbage [50], onion [51] and wheat [52]. Male sterility might be related to energy metabolism turbulence, excessive ethylene synthesis, and suffocation of starch synthesis in pepper [53]. In addition, pentatricopeptide repeat proteins, heat shock proteins, stress proteins, MYB, bHLH and heat shock transcription factors and anther development related genes may be the candidates for pollen abortion in Brassica crops [49,50,54]. In this study, we found that genes related to polygalacturonase metabolism, glycosyl hydrolases, oxidation reduction process, phenylalanine metabolism, phenylpropanoid biosynthesis were significantly changed in ogu CMS abortive buds compared with the other groups. Ethylene-responsive transcription factor 115 and transcriptional factor basic helix-loop-helix 137 were both significantly upregulated in ogu CMS abortive buds. Therefore, our results were basically consistent with the results of previous studies [46,47,48,49,50,51,52,53,54], and the genes discovered related to energy metabolism, oxidation reduction process and phenylpropanoid biosynthesis, ethylene-responsive transcription factor 115 and transcriptional factor basic helix-loop-helix 137 may be important factors contributing to Broccoli ogu CMS pollen abortion and bud abortion. Further experiments are needed to elucidate the molecular mechanisms of these genes that lead to broccoli CMS and bud abortion.

4. Materials and Methods

4.1. Plant Materials

Broccoli (Brassica oleracea var. italica) maintainer ML93219 (showing slight bud abortion) and ogu CMS93219 (serious levels of bud abortion) were bred by the Institute of Vegetables and Flowers, Chinese Academy of Agricultural Sciences. The backcross paternal line of ogu CMS93219 was ML93219 and the number of backcross generations was sixteen. In the spring of 2015, the plants were grown in an experimental greenhouse at the Institute of Vegetable and Flowers, Chinese Academy of Agricultural Sciences, Changping (Beijing, China). We handled the main bouquet using the approach proposed by Shu et al. [55]. When the plants began to flower, four kinds of bud samples, ML_NB, ML_AB, CMS_NB and CMS_AB, were collected and labelled with three biological replicates. To ensure the integrity of the sample RNA, isolated buds were immediately frozen in liquid nitrogen and stored at −80 °C before RNA extraction.

4.2. RNA Extraction and Quality Testing

Total RNA was extracted using an EASYspin Plus Plant RNA-38 Kit, according to the manufacturer’s instructions (Juhuatech Co., Ltd., Beijing, China). The integrity and purity of the RNA samples were determined by 1% agarose gels electrophoresis, and the RNA concentration was measured by Qubit® RNA Assay Kit in Qubit® 2.0 Flurometer (Life Technologies Corporation, Carlsbad, CA, USA) and the integrity of the RNA was assessed by an RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies Inc., Santa Clara, CA, USA).

4.3. RNA-Seq Library Construction and Illumina Sequencing

Twelve strand-specific RNA-Seq libraries were constructed with cDNA fragments of 250–300 bp in length. An Illumina TruSeq PE Cluster Kit v3-cBot-HS on a cBot Cluster Generation System was then used to cluster the samples, according to the manufacturer’s instructions. After cluster generation, the libraries were sequenced on an Illumina Hiseq™ 4000 system and reads were generated.

4.4. RNA-Seq Data Quality Control and Transcriptome de Novo Assembly

We obtained clean reads by removing the adaptor reads, unknown sequences “N” (reads containing more than 10% unknown nucleotides), low quality reads (reads containing more than 50% bases with Q-value ≤ 5) from the raw data. The Q20, Q30 and GC-content were then calculated based on the clean reads. The high quality clean reads were used for downstream analyses. Transcriptome de novo assembly was executed using Trinity [56] with min_kmer_cov set to 2 by default and other parameters set at their defaults. After assembly, the longest transcripts of each gene were selected as the unigenes.

4.5. Unigene Function Annotation

We annotated the unigenes based on seven databases, NCBI blast (2.2.28+) was used to search against the Nr (E-value = 1 × 10−5), Nt (E-value = 1 × 10−5), Swiss-Prot (E-value = 1 × 10−5) and KOG databases (E-value = 1 × 10−3). The unigenes were divided into 26 groups and their participation in different metabolic pathways based on KOG annotation was assessed. Pfam annotated was determined using the HMMER 3.0 package [57], hmmscan (e-value = 0.01). GO annotations for the unigenes were determined by Blast2GO v2.5 [12] with the self-write script (e-value = 1 × 10−6) based on the annotation result of Nr and Pfam, which has three ontologies: molecular function, cellular component and biological process [58]. KEGG [59] related annotations were identified by the KAAS and KEGG Automatic Annotation Server [60] (E-value = 1 × 10−10) to determine the metabolic pathway of unigenes.

4.6. Analysis of DEGs

Alignment results of bowtie were counted by RSEM [61]. FPKM [62] values were used to calculate the gene expression levels of the four groups of normal and abortive buds from the maintainer and ogu CMS lines. FPKM has become the most commonly used method to estimate the level of gene expression and takes into account the effects of sequencing depth and gene length on the calculation of gene expression [62]. There were three biological replicates; therefore, the calculated gene expression could be used directly to compare the gene expression between samples. Referring to the statistical method of Storey and Tibshirani [63], |log2Fold change| > 1 and p-adjusted < 0.05 were set as the threshold for significantly differential expression. p-Values were adjusted to control the false discovery rate, referring to Benjamini and Hochberg’s approach [64]. Then, based on the Wallenius non-central hyper-geometric distribution [65], GO and KEGG functional enrichment analysis of the DEGs was executed by the GOseq [66] and KOBAS software [60], respectively.

4.7. qRT-PCR Validation

qRT-PCR analyses with the three biological replicates samples used for RNA-Seq were performed to verify the DGEs results. Twenty-one common differentially expressed genes were randomly selected that accounted for about 22.1% of the 95 common differentially expressed genes. Specific primers were designed using the Primer-BLAST tool (available online: http://0-www-ncbi-nlm-nih-gov.brum.beds.ac.uk/tools/primer-blast/index.cgi? LINK _LOC=BlastHome) in NCBI and synthesized by Sangon Biotech Co., Ltd. (Shanghai, China). cDNAs were reverse transcribed from total RNA using a PrimeScript RT reagent Kit (Takara, Dalian, China). qRT-PCR was carried out according to the SYBR PrimeScript RT-PCR Kit manufacturer specifications (Takara) on an ABI Prism®7900 Real-Time PCR System (Applied Biosystems, Foster City, CA, USA). To normalize the gene expression data, we used the broccoli β-actin gene as an internal standard [67]. The 2−ΔΔCt method [68] was used to determine the relative expression of genes. The standard deviation was calculated based on the three biological replicates. The specific primers sequences are listed in Table S6.

5. Conclusions

In this study, we found that buds abortion was related with polygalacturonase metabolism, glycosyl hydrolases, oxidation reduction process, phenylalanine metabolism, and phenylpropanoid biosynthesis. Moreover, 19 common differentially expressed genes associated glycosyl hydrolases, reactive oxygen species scavenging, inhibitor, plant defense, cell division and expansion, transporter, and four transcriptional factors—ethylene-responsive transcription factor 115, transcriptional factor basic helix-loop-helix 137, transcription factors DUO1, and PosF21/RF2a/BZIP34—may be associated with ogu CMS abortive buds. In conclusion, our results not only increased our understanding of ogu CMS bud abortive mechanisms and provided a valuable resource for further functional characterization of ogu CMS bud abortion, but also laid the foundation for molecular breeding to overcoming bud abortion in broccoli, as well as other Brassica crops in the future.

Supplementary Materials

Supplementary materials can be found at https://0-www-mdpi-com.brum.beds.ac.uk/1422-0067/19/9/2501/s1. Table S1: Sequences of the specific primers; Table S2: KEGG classification of the unigenes; Table S3: Unigenes expression fold change of CMS_ABvsCMS_NB; Table S4: Unigenes expression fold change of CMS_ABvsML_AB; Table S5: Unigenes expression fold change of CMS_ABvsML_NB; Table S6: FPKM and annotation of the 95 common differentially expressed genes.

Author Contributions

Conceptualization, J.S. and Y.L.; Data curation, J.S.; Formal analysis, J.S., L.Z., Z.L. and H.L.; Funding acquisition, Y.L.; Investigation, J.S.; Methodology, J.S.; Resources, Y.L., Z.F., L.Y., M.Z. and Y.Z.; Supervision, Y.L.; Writing—original draft, J.S. and L.Z.; Writing—review & editing, J.S. and Y.L. All authors read and approved the final manuscript.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (Grant No. 31372067), the China Agriculture Research System (Grant No. CARS-25-A), Key Projects in the National Science and Technology Pillar Program of China (Grant No. 2013BAD01B04), the Key Laboratory of Biology and Genetic Improvement of Horticultural Crops, Ministry of Agriculture, China, and the Science and Technology Innovation Program of the Chinese Academy of Agricultural Sciences (Grant No. CAAS-ASTIP-IVFCAAS). We thank Liwen Bianji, Edanz Group China (available online: www.liwenbianji.cn/ac), for editing the English text of a draft of this manuscript.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Zhang, J.; Li, Q.F.; Huang, W.W.; Xu, X.Y.; Zhang, X.L.; Hui, M.X.; Zhang, M.K.; Zhang, L.G. A vacuolar processing enzyme RsVPE1 gene of radish is involved in floral bud abortion under heat stress. Int. J. Mol. Sci. 2013, 14, 13346–13359. [Google Scholar] [CrossRef] [PubMed]
  2. Zhang, J.; Sun, X.L.; Zhang, L.G.; Hui, M.X.; Zhang, M.K. Analysis of differential gene expression during floral bud abortion in radish (Raphanus sativus L.). Genet. Mol. Res. 2013, 12, 2507–2516. [Google Scholar] [CrossRef] [PubMed]
  3. Liu, J.; Zhang, L.G.; Wang, F.M.; Hui, M.X.; Zhang, M.K. Observation of histocytological feature on radish flower bud during aborting. Acta Agric. Bor. Occid. Sin. 2008, 5, 272–276. [Google Scholar]
  4. Wang, Q.B.; Zhang, Y.Y.; Liu, Y.M.; Yang, L.M.; Zhuang, M.; Sun, P.T. The floral and seed setting characteristics in two types of male sterile lines of cabbage (Brassica oleracea L. var. capitata). Acta Hortic. Sin. 2011, 38, 61–68. [Google Scholar]
  5. Shu, J.S.; Liu, Y.M.; Li, Z.S.; Zhang, L.L.; Fang, Z.Y.; Yang, L.M.; Zhuang, M.; Zhang, Y.Y.; Li, Z.S.; Sun, P.T. Study on the floral characteristics and structure in two types of male sterile lines of broccoli (Brassica oleracea var. italica). J. Plant. Genet. Resour. 2014, 15, 113–119. [Google Scholar]
  6. Jia, J.; Zhang, L.G. mRNA differential display and EST sequence analysis of aborted bud and normal bud in radish (Raphanus sativus). Acta Agric. Nucl. Sin. 2008, 22, 426–431. [Google Scholar]
  7. Zhang, L.G.; Jia, J.; Zhang, S.L.; Zhang, Y. cDNA-AFLP differential expression analysis of genes related with aborting bud in Chinese cabbage. J. Agric. Biotechnol. 2010, 18, 489–492. [Google Scholar]
  8. Wang, Q.; Liu, Z.; Feng, H. cDNA-AFLP differential expression analysis of genes about bud aborting in genetic male sterile line of Brassica campestris L. ssp. chinensis (L.) Makino var. rosularis Tsen et Lee. Mol. Plant. Breed. 2014, 12, 118–126. [Google Scholar]
  9. Kensler, T.W.; Chen, J.G.; Egner, P.A.; Fahey, J.W.; Jacobson, L.P.; Stephenson, K.K.; Ye, L.X.; Coady, J.L.; Wang, J.B.; Wu, Y.; et al. Effects of glucosinolate-rich broccoli sprouts on urinary levels of aflatoxin-DNA adducts and phenanthrene tetraols in a randomized clinical trial in He Zuo township, Qidong, People’s Republic of China. Cancer Epidemiol. Biomark. Prev. 2005, 14, 2605–2613. [Google Scholar] [CrossRef] [PubMed]
  10. Canene-Adams, K.; Lindshield, B.L.; Wang, S.; Jeffery, E.H.; Clinton, S.K.; Erdman, J.W., Jr. Combinations of tomato and broccoli enhance antitumor activity in dunning r3327-h prostate adenocarcinomas. Cancer Res. 2007, 67, 836–843. [Google Scholar] [CrossRef] [PubMed]
  11. Jeffery, E.H.; Araya, M. Physiological effects of broccoli consumption. Phytochem. Rev. 2009, 8, 283–298. [Google Scholar] [CrossRef]
  12. Götz, S.; García-Gómez, J.M.; Terol, J.; Williams, T.D.; Nagaraj, S.H.; Nueda, M.J.; Robles, M.; Talon, M.; Dopazo, J.; Conesa, A. High-throughput functional annotation and data mining with the Blast2GO suite. Nucleic Acids Res. 2008, 36, 3420–3435. [Google Scholar] [Green Version]
  13. Thomas, S.G.; Franklin-Tong, V.E. Self-incompatibility triggers programmed cell death in Papaver pollen. Nature 2004, 429, 305–309. [Google Scholar] [CrossRef] [PubMed]
  14. Reape, T.J.; Molony, E.M.; McCabe, P.F. Programmed cell death in plants: Distinguishing between different modes. J. Exp. Bot. 2008, 59, 435–444. [Google Scholar] [CrossRef] [PubMed]
  15. Rybaczek, D.; Musiałek, M.W.; Balcerczyk, A. Caffeine-induced premature chromosome condensation results in the apoptosis-like programmed cell death in root meristems of Vicia faba. PLoS ONE 2015, 10, e0142307. [Google Scholar] [CrossRef] [PubMed]
  16. Bosch, M.; Franklin-Tong, V.E. Temporal and spatial activation of caspase-like enzymes induced by self-incompatibility in Papaver pollen. Proc. Natl. Acad. Sci. USA 2007, 104, 18327–18332. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  17. Wilkins, K.A.; Bancroft, J.; Bosch, M.; Ings, J.; Smirnoff, N.; Franklin-Tong, V.E. Reactive oxygen species and nitric oxide mediate actin reorganization and programmed cell death in the self-incompatibility response of papaver. Plant Physiol. 2011, 156, 404–416. [Google Scholar] [CrossRef] [PubMed]
  18. Zhang, F.J.; Wang, Z.Q.; Dong, W.; Sun, C.Q.; Wang, H.B.; Song, A.P.; He, L.Z.; Fang, W.M.; Chen, F.D.; Teng, N.J. Transcriptomic and proteomic analysis reveals mechanisms of embryo abortion during chrysanthemum cross breeding. Sci. Rep. 2014, 4, 6536. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  19. Fagundes, D.; Bohn, B.; Cabreira, C.; Leipelt, F.; Dias, N.; Bodanese-Zanettini, M.H.; Cagliari, A. Caspases in plants: Metacaspase gene family in plant stress responses. Funct. Integr. Genom. 2015, 15, 639–649. [Google Scholar] [CrossRef] [PubMed]
  20. Gonzalez-Carranza, Z.H.; Elliott, K.A.; Roberts, J.A. Expression of polygalacturonases and evidence to support their role during cell separation processes in Arabidopsis thaliana. J. Exp. Bot. 2007, 58, 3719–3730. [Google Scholar] [CrossRef] [PubMed]
  21. Ogawa, M.; Kay, P.; Wilson, S.; Swain, S.M. ARABIDOPSIS DEHISCENCE ZONE POLYGALACTURONASE1 (ADPG1), ADPG2, and QUARTET2 are polygalacturonases required for cell separation during reproductive development in Arabidopsis. Plant Cell 2009, 21, 216–233. [Google Scholar] [CrossRef] [PubMed]
  22. Libertini, E.; Li, Y.; McQueen-Mason, S.J. Phylogenetic analysis of the plant endo-β-1, 4-glucanase gene family. J. Mol. Evol. 2004, 58, 506–515. [Google Scholar] [CrossRef] [PubMed]
  23. Kawano, T.; Kawano, T.; Kadono, T.; Fumoto, K.; Lapeyrie, F.; Kuse, M.; Isobe, M.; Furuichi, T.; Muto, S. Aluminum as a specific inhibitor of plant TPC1 Ca2+ channels. Biochem. Biophys. Res. Commun. 2004, 324, 40–45. [Google Scholar] [CrossRef] [PubMed]
  24. Peiter, E.; Maathuis, F.J.; Mills, L.N.; Knight, H.; Pelloux, J.; Hetherington, A.M.; Sanders, D. The vacuolar Ca2+-activated channel TPC1 regulates germination and stomatal movement. Nature 2005, 434, 404–408. [Google Scholar] [CrossRef] [PubMed]
  25. Zhao, J.W.; Wang, J.L.; An, L.L.; Doerge, R.W.; Chen, Z.J.; Grau, C.R.; Meng, J.L.; Osborn, T.C. Analysis of gene expression profiles in response to Sclerotinia sclerotiorum in Brassica napus. Planta 2007, 227, 13–24. [Google Scholar] [CrossRef] [PubMed]
  26. Libault, M.; Wan, J.; Czechowski, T.; Udvardi, M.; Stacey, G. Identification of 118 Arabidopsis transcription factor and 30 ubiquitin-ligase genes responding to chitin, a plant-defense elicitor. Mol. Plant Microbe Interact. 2007, 20, 900–911. [Google Scholar] [CrossRef] [PubMed]
  27. Seo, D.H.; Ryu, M.Y.; Jammes, F.; Hwang, J.H.; Turek, M.; Kang, B.G.; Kwak, J.M.; Kim, W.T. Roles of four Arabidopsis U-Box E3 ubiquitin ligases in negative regulation of abscisic acid-mediated drought stress responses. Plant Physiol. 2012, 160, 556–568. [Google Scholar] [CrossRef] [PubMed]
  28. Wan, J.R.; Zhang, X.C.; Neece, D.; Ramonell, K.M.; Clough, S.; Kim, S.Y.; Stacey, M.G.; Stacey, G. A LysM receptor-like kinase plays a critical role in chitin signaling and fungal resistance in Arabidopsis. Plant Cell 2008, 20, 471–481. [Google Scholar] [CrossRef] [PubMed]
  29. Brotmana, Y.; Landau, U.; Pnini, S.; Lisec, J.; Balazadeh, S.; Mueller-Roeber, B.; Zilberstein, A.; Willmitzer, L.; Chet, I.; Viterbo, A. The LysM receptor-like kinase LysM RLK1 is required to activate defense and abiotic-stress responses induced by overexpression of fungal chitinases in Arabidopsis plants. Mol. Plant. 2012, 5, 1113–1124. [Google Scholar] [CrossRef] [PubMed]
  30. Tanaka, K.; Nguyen, C.T.; Liang, Y.; Cao, Y.; Stacey, G. Role of LysM receptors in chitin-triggered plant innate immunity. Plant. Signal. Behav. 2013, 8, e22598. [Google Scholar] [CrossRef] [PubMed]
  31. Garcia, M.J.; Lucena, C.; Romera, F.J.; Alcantara, E.; Perez-Vicente, R. Ethylene and nitric oxide involvement in the up-regulation of key genes related to iron acquisition and homeostasis in Arabidopsis. J. Exp. Bot. 2010, 61, 3885–3899. [Google Scholar] [CrossRef] [PubMed]
  32. Xu, Z.S.; Chen, M.; Li, L.C.; Ma, Y.Z. Functions and application of the AP2/ERF transcription factor family in crop improvement. J. Integr. Plant Biol. 2011, 53, 570–585. [Google Scholar] [CrossRef] [PubMed]
  33. Chen, M.K.; Shpak, E.D. ERECTA family genes regulate development of cotyledons during embryogenesis. FEBS Lett. 2014, 588, 3912–3917. [Google Scholar] [CrossRef] [PubMed]
  34. Nakano, T.; Fujisawa, M.; Shima, Y.; Ito, Y. The AP2/ERF transcription factor SlERF52 functions in flower pedicel abscission in tomato. J. Exp. Bot. 2014, 65, 3111–3119. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  35. Ito, Y.; Nakanot, T. Development and regulation of pedicel abscission in tomato. Front. Plant Sci. 2015, 6, 442. [Google Scholar] [CrossRef] [PubMed]
  36. Jisha, V.; Dampanaboina, L.; Vadassery, J.; Mithöfer, A.; Kappara, S.; Ramanan, R. Overexpression of an AP2/ERF type transcription factor OsEREBP1 confers biotic and abiotic stress tolerance in rice. PLoS ONE 2015, 10, e0127831. [Google Scholar] [CrossRef] [PubMed]
  37. Mishra, S.; Phukan, U.J.; Tripathi, V.; Singh, D.K.; Luqman, S.; Shukla, R.K. PsAP2 an AP2/ERF family transcription factor from Papaver somniferum enhances abiotic and biotic stress tolerance in transgenic tobacco. Plant Mol. Biol. 2015, 89, 173–186. [Google Scholar] [CrossRef] [PubMed]
  38. Babitha, K.C.; Vemanna, R.S.; Nataraja, K.N.; Udayakumar, M. Overexpression of EcbHLH57 transcription factor from Eleusine coracana L. in tobacco confers tolerance to salt, oxidative and drought stress. PLoS ONE 2015, 10, e0137098. [Google Scholar] [CrossRef] [PubMed]
  39. Heyman, J.; Cools, T.; Vandenbussche, F.; Heyndrickx, K.S.; Van Leene, J.; Vercauteren, I.; Vanderauwera, S.; Vandepoele, K.; De Jaeger, G.; Van Der, S.D.; et al. ERF115 controls root quiescent center cell division and stem cell replenishment. Science 2013, 342, 860–863. [Google Scholar] [CrossRef] [PubMed]
  40. Heim, M.A.; Jakoby, M.; Werber, M.; Martin, C.; Weisshaar, B.; Bailey, P.C. The basic helix–loop–helix transcription factor family in plants: A genome-wide study of protein structure and functional diversity. Mol. Biol. Evol. 2003, 20, 735–747. [Google Scholar] [CrossRef] [PubMed]
  41. Borg, M.; Brownfield, L.; Khatab, H.; Sidorova, A.; Lingaya, M.; Twel, D. The R2R3 MYB transcription factor DUO1 activates a male germline-specific regulon essential for sperm cell differentiation in Arabidopsis. Plant Cell 2011, 23, 534–549. [Google Scholar] [CrossRef] [PubMed]
  42. Aeschbacher, R.A.; Schrott, M.; Potrykus, I.; Saul, M.W. Isolation and molecular characterization of PosF21, an Arabidopsis thaliana gene which shows characteristics of a b-Zip class transcription factor. Plant J. 1991, 1, 303–316. [Google Scholar] [CrossRef] [PubMed]
  43. Yin, Y.H.; Zhu, Q.; Dai, S.H.; Lamb, C.; Beachy, R.N. RF2a, a bZIP transcriptional activator of the phloem-specific rice tungro bacilliform virus promoter, functions in vascular development. EMBO J. 1997, 16, 5247–5259. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  44. Shen, H.; Cao, K.; Wang, X. A conserved proline residue in the leucine zipper region of AtbZIP34 and AtbZIP61 in Arabidopsis thaliana interferes with the formation of homodimer. Biochem. Biophys. Res. Commun. 2007, 362, 425–430. [Google Scholar] [CrossRef] [PubMed]
  45. Gibalová, A.; Renák, D.; Matczuk, K.; Dupl’áková, N.; Cháb, D.; Twell, D.; Hongys, D. AtbZIP34 is required for Arabidopsis pollen wall patterning and the control of several metabolic pathways in developing pollen. Plant Mol. Biol. 2009, 70, 581–601. [Google Scholar]
  46. Yang, P.; Han, J.; Huang, J. Transcriptome sequencing and de novo analysis of cytoplasmic male sterility and maintenance in JA-CMS cotton. PLoS ONE 2014, 9, e112320. [Google Scholar] [CrossRef] [PubMed]
  47. Liu, Z.; Shi, X.; Li, S.; Hu, G.; Zhang, L.; Song, X. Tapetal-Delayed programmed cell death (PCD) and oxidative stress-induced male sterility of Aegilops uniaristata cytoplasm in wheat. Int. J. Mol. Sci. 2018, 19, 1708. [Google Scholar] [CrossRef] [PubMed]
  48. Li, J.; Han, S.; Ding, X.; He, T.; Dai, J.; Yang, S.; Gai, J. Comparative transcriptome analysis between the cytoplasmic male sterile line NJCMS1A and its maintainer NJCMS1B in soybean (Glycine max (L.) Merr.). PLoS ONE 2015, 10, e0126771. [Google Scholar] [CrossRef] [PubMed]
  49. Du, K.; Liu, Q.; Wu, X.; Jiang, J.; Wu, J.; Fang, Y.; Li, A.; Wang, Y. Morphological structure and transcriptome comparison of the cytoplasmic male sterility line in Brassica napus (SaNa-1A) derived from somatic hybridization and its maintainer line SaNa-1B. Front. Plant Sci. 2016, 7, 1313. [Google Scholar] [CrossRef] [PubMed]
  50. Wang, S.; Wang, C.; Zhang, X.X.; Chen, X.; Liu, J.J.; Jia, X.F.; Jia, S.Q. Transcriptome de novo assembly and analysis of differentially expressed genes related to cytoplasmic male sterility in cabbage. Plant Physiol. Biochem. 2016, 105, 224–232. [Google Scholar] [CrossRef] [PubMed]
  51. Liu, Q.; Lan, Y.; Wen, C.; Zhao, H.; Wang, J.; Wang, Y. Transcriptome sequencing analyses between the cytoplasmic male sterile line and its maintainer line in Welsh onion (Allium fistulosum L.). Int. J. Mol. Sci. 2016, 17, 1058. [Google Scholar] [CrossRef] [PubMed]
  52. Zhang, G.; Ye, J.; Jia, Y.; Zhang, L.; Song, X. ITRAQ-based proteomics analyses of sterile/fertile anthers from a thermo-sensitive cytoplasmic male-sterile wheat with Aegilops kotschyi cytoplasm. Int. J. Mol. Sci. 2018, 19, 1344. [Google Scholar] [CrossRef] [PubMed]
  53. Wu, Z.; Cheng, J.; Qin, C.; Hu, Z.; Yin, C.; Hu, K. Differential proteomic analysis of anthers between cytoplasmic male sterile and maintainer lines in Capsicum annuum L. Int. J. Mol. Sci. 2013, 14, 22982–22996. [Google Scholar] [CrossRef] [PubMed]
  54. Mei, S.; Liu, T.; Wang, Z. Comparative transcriptome profile of the cytoplasmic male sterile and fertile floral buds of radish (Raphanus sativus L.). Int. J. Mol. Sci. 2016, 17, 42. [Google Scholar] [CrossRef] [PubMed]
  55. Shu, J.; Liu, Y.; Li, Z.; Zhang, L.; Fang, Z.; Yang, L.; Zhuang, M.; Zhang, Y.; Lv, H. Organelle simple sequence repeat markers help to distinguish carpelloid stamen and normal cytoplasmic male sterile sources in broccoli. PLoS ONE 2015, 10, e0138750. [Google Scholar] [CrossRef] [PubMed]
  56. Grabherr, M.G.; Haas, B.J.; Yassour, M.; Levin, J.Z.; Thompson, D.A.; Amit, I.; Adiconis, X.; Fan, L.; Raychowdhury, R.; Zeng, Q.; et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat. Biotechnol. 2011, 29, 644–652. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  57. Zhang, Z.; Wood, W.I. A profile hidden Markov model for signal peptides generated by HMMER. Bioinformatics 2003, 19, 307–308. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  58. Ashburner, M.; Ball, C.A.; Blake, J.A.; Botstein, D.; Butler, H.; Cherry, J.M.; Davis, A.P.; Dolinski, K.; Dwight, S.S.; Eppig, J.T.; et al. Gene Ontology: Tool for the unification of biology. Nat. Genet. 2000, 25, 25–29. [Google Scholar] [CrossRef] [PubMed]
  59. Kanehisa, M.; Araki, M.; Goto, S.; Hattori, M.; Hirakawa, M.; Itoh, M.; Katayama, T.; Kawashima, S.; Okuda, S.; Tokimatsu, T.; et al. KEGG for linking genomes to life and the environment. Nucleic Acids Res. 2008, 36, D480–D484. [Google Scholar] [CrossRef] [PubMed]
  60. Mao, X.; Cai, T.; Olyarchuk, J.G.; Wei, L. Automated genome annotation and pathway identification using the KEGG Orthology (KO) as a controlled vocabulary. Bioinformatics 2005, 21, 3787–3793. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  61. Li, B.; Dewey, C. RSEM: Accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinform. 2011, 12, 323. [Google Scholar] [CrossRef] [PubMed]
  62. Trapnell, C.; Williams, B.A.; Pertea, G.; Mortazavi, A.; Kwan, G.; van Baren, M.J.; Salzberg, S.L.; Wold, B.J.; Pachter, L. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat. Biotechnol. 2010, 28, 511–515. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  63. Storey, J.D.; Tibshirani, R. Statistical significance for genomewide studies. Proc. Nat. Acad. Sci. USA 2003, 100, 9440–9445. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  64. Benjamini, Y.; Hochberg, Y. Controlling the false discovery rate: A practical and powerful approach to multiple testing. J. R. Stat. Soc. Ser. B Methodol. 1995, 57, 289–300. [Google Scholar]
  65. Wallenius, K.T. Biased Sampling: The Non-Central Hypegeometric Probability Distribution. Ph.D. Thesis, Stanford University, Stanford, CA, USA, 1963. [Google Scholar]
  66. Young, M.D.; Wakefield, M.J.; Smyth, G.K.; Oshlack, A. Gene ontology analysis for RNA-Seq: Accounting for selection bias. Genome Biol. 2010, 11, R14. [Google Scholar] [CrossRef] [PubMed]
  67. Gómez-Lobato, M.E.; Hasperuéb, J.H.; Civelloa, P.M.; Chaves, A.R.; Martínez, G.A. Effect of 1-MCP on the expression of chlorophyll degrading genes during senescence of broccoli (Brassica oleracea L.). Sci. Hortic 2012, 144, 208–211. [Google Scholar]
  68. Livak, K.J.; Schmittgen, T.D. Analysis of relative gene expression data using real-time quantitative PCR and the 2−ΔΔCT Method. Methods 2001, 25, 402–408. [Google Scholar] [CrossRef] [PubMed]
Ijms 19 02501 g001 550
Figure 1. Morphological characterization of two broccoli lines with different degrees of bud abortion. (A,B) ML93219, (C,D) CMS93219. B-I and B-II represent the normal and abortive bud of ML93219, respectively. D-I and D-II represent the normal and abortive bud of CMS93219, respectively. The bars in (E,F) represent the standard deviation (n = 15). Asterisks indicate that the average number of abortive buds and the abortive buds rate per branch are very significantly different between ML93219 and CMS93219 (unpaired t test, p < 0.01).
Figure 1. Morphological characterization of two broccoli lines with different degrees of bud abortion. (A,B) ML93219, (C,D) CMS93219. B-I and B-II represent the normal and abortive bud of ML93219, respectively. D-I and D-II represent the normal and abortive bud of CMS93219, respectively. The bars in (E,F) represent the standard deviation (n = 15). Asterisks indicate that the average number of abortive buds and the abortive buds rate per branch are very significantly different between ML93219 and CMS93219 (unpaired t test, p < 0.01).
Ijms 19 02501 g001
Ijms 19 02501 g002 550
Figure 2. Clusters of eukaryotic orthologous groups (KOG) functional classification of the broccoli bud transcriptome.
Figure 2. Clusters of eukaryotic orthologous groups (KOG) functional classification of the broccoli bud transcriptome.
Ijms 19 02501 g002
Ijms 19 02501 g003 550
Figure 3. Gene Ontology (GO) classifications of the broccoli bud transcriptome.
Figure 3. Gene Ontology (GO) classifications of the broccoli bud transcriptome.
Ijms 19 02501 g003
Ijms 19 02501 g004 550
Figure 4. Boxplot of Fragments per kb per million fragments (FPKM) distribution for the four samples. Five statistics are represented by different regions of the Boxplot; from the top down they are the maximum, upper quartile, median, lower quartile, and minimum, respectively.
Figure 4. Boxplot of Fragments per kb per million fragments (FPKM) distribution for the four samples. Five statistics are represented by different regions of the Boxplot; from the top down they are the maximum, upper quartile, median, lower quartile, and minimum, respectively.
Ijms 19 02501 g004
Ijms 19 02501 g005 550
Figure 5. Cluster analysis of differentially expressed genes among the four samples. Heatmap of differentially expressed genes among the four samples. Red indicates high expression, and blue indicates low expression. Color from red to blue represents descending log10 (FPKM + 1).
Figure 5. Cluster analysis of differentially expressed genes among the four samples. Heatmap of differentially expressed genes among the four samples. Red indicates high expression, and blue indicates low expression. Color from red to blue represents descending log10 (FPKM + 1).
Ijms 19 02501 g005
Ijms 19 02501 g006 550
Figure 6. Differentially expressed genes (DEGs) between normal and abortive buds from ML93219 and CMS93219. (a) DEGs between CMS AB and CMS NB. (b) DEGs between CMS AB and ML AB. (c) DEGs between CMS AB and ML NB. The x-axis indicates the log2 (fold change) between the two samples. The y-axis indicates −log10 (padj) (p-adjusted). The scatter points in the figure represent individual genes, the green dots indicate genes with no significant differences, the red dots indicate up-regulated genes with significant differences, and the blue dots indicate down-regulated genes with significant differences. The screening condition for DEGs is padj < 0.05.
Figure 6. Differentially expressed genes (DEGs) between normal and abortive buds from ML93219 and CMS93219. (a) DEGs between CMS AB and CMS NB. (b) DEGs between CMS AB and ML AB. (c) DEGs between CMS AB and ML NB. The x-axis indicates the log2 (fold change) between the two samples. The y-axis indicates −log10 (padj) (p-adjusted). The scatter points in the figure represent individual genes, the green dots indicate genes with no significant differences, the red dots indicate up-regulated genes with significant differences, and the blue dots indicate down-regulated genes with significant differences. The screening condition for DEGs is padj < 0.05.
Ijms 19 02501 g006
Ijms 19 02501 g007 550
Figure 7. DEGs between normal and abortive buds from broccoli maintainer and ogu cytoplasmic male sterile (CMS) lines.
Figure 7. DEGs between normal and abortive buds from broccoli maintainer and ogu cytoplasmic male sterile (CMS) lines.
Ijms 19 02501 g007
Ijms 19 02501 g008 550
Figure 8. Verification of the DEGs by qRT-PCR. Eleven DEGs with higher expression and ten DEGs with lower expression in CMS_AB were selected for qRT-PCR validation. The relative expression level of each gene was expressed as the FPKM among four samples in the RNA-Seq data (black line) and qRT-PCR data (blue bar). To normalize the expression data, the broccoli β-actin gene was used as the internal control. The bars represent the standard deviation.
Figure 8. Verification of the DEGs by qRT-PCR. Eleven DEGs with higher expression and ten DEGs with lower expression in CMS_AB were selected for qRT-PCR validation. The relative expression level of each gene was expressed as the FPKM among four samples in the RNA-Seq data (black line) and qRT-PCR data (blue bar). To normalize the expression data, the broccoli β-actin gene was used as the internal control. The bars represent the standard deviation.
Ijms 19 02501 g008
Ijms 19 02501 g009 550
Figure 9. Heatmap of the 19 enriched common differentially expressed genes in normal and abortive broccoli buds. The bar represents the expression levels for each gene (log10 (FPKM + 1)) in the ML_NB, ML_AB, CMS_NB, and CMS_AB groups, as indicated by red or green rectangles. Red means upregulation of genes and green means downregulation.
Figure 9. Heatmap of the 19 enriched common differentially expressed genes in normal and abortive broccoli buds. The bar represents the expression levels for each gene (log10 (FPKM + 1)) in the ML_NB, ML_AB, CMS_NB, and CMS_AB groups, as indicated by red or green rectangles. Red means upregulation of genes and green means downregulation.
Ijms 19 02501 g009
Table
Table 1. Summary of broccoli bud transcriptome sequencing data.
Table 1. Summary of broccoli bud transcriptome sequencing data.
Sample NameRaw ReadsClean Reads (Clean/Raw)Clean Bases (Gb)Q20 (%)GC Content (%)Mapped Reads (Mapped/Clean)
ML_NB157,767,24655,194,896 (95.55%)8.2894.3046.9341,252,694 (74.74%)
ML_NB244,577,57043,849,450 (98.37%)6.5892.9346.4632,614,176 (74.38%)
ML_NB361,251,85259,738,162(97.53%)8.9694.0346.9644,511,290 (74.51%)
ML_AB161,588,74060,493,582 (98.22%)9.0894.3246.7345,965,016 (75.98%)
ML_AB254,823,83653,530,342 (97.64%)8.0293.9046.7440,115,268 (74.94%)
ML_AB358,528,85256,475,844 (96.49%)8.4894.3346.8142,816,816 (75.81%)
CMS_NB154,455,25853,170,054 (97.64%)7.9894.2347.0239,776,946 (74.81%)
CMS_NB253,265,51452,284,140 (98.16%)7.8493.9747.1238,965,364 (74.53%)
CMS_NB351,816,26850,884,086 (98.20%)7.6493.7946.9137,800,172 (74.29%)
CMS_AB164,340,17663,146,694 (98.15%)9.4894.1647.0447,846,088 (75.77%)
CMS_AB256,566,28055,530,456 (98.17%)8.3294.2746.7842,582,890 (76.68%)
CMS_AB360,701,77059,404,212 (97.86%)8.9294.0946.7044,549,934 (74.99%)
Note: ML_NB and ML_AB represent the normal and abortive bud samples of ML93219, respectively. CMS_NB and CMS_AB represent the normal and abortive bud samples of CMS93219, respectively.
Table
Table 2. Summary of de novo transcriptome length distribution.
Table 2. Summary of de novo transcriptome length distribution.
Length200–500 bp500–1 kbp1 k–2 kbp>2 kbpTotalMin (bp)Mean (bp)Median (bp)Max (bp)N50Total Nucleotides
No. of transcripts44,03021,34821,46710,50297,34720193657816,361151091,137,323
No. of Unigenes37,34412,54510,691547066,05020178642416,361136351,896,834
Table
Table 3. Annotation of unigene sequences in broccoli buds.
Table 3. Annotation of unigene sequences in broccoli buds.
Sequence DatabaseNumber of Annotated UnigenesPercentage of Annotated Unigene Sequences (%)Software and Parameters
Annotated in Nr44,29467.06NCBI blast 2.2.28+, e-value = 1 × 10−5
Annotated in Nt50,15775.93NCBI blast 2.2.28+, e-value = 1 × 10−5
Annotated in KO12,40318.77KAAS, KEGG Automatic Annotation Server, e-value = 1 × 10−10
Annotated in SwissProt29,86145.2NCBI blast 2.2.28+, e-value = 1 × 10−5
Annotated in PFAM24,66037.33HMMER 3.0 package, hmmscan, e-value = 0.01
Annotated in GO30,28545.85Blast2GO v2.5 [12] and self-write the script, e-value = 1 × 10−6
Annotated in KOG12,49218.91NCBI blast 2.2.28+, e-value = 1 × 10−3
Annotated in all Databases60709.19-
Annotated in at least one Database54,75382.89-
Total Unigenes66,050100-
Table
Table 4. Results of DEGs enriched by Gene Ontology (GO) term.
Table 4. Results of DEGs enriched by Gene Ontology (GO) term.
GO AccessionDescriptionTerm Typep ValueCorrected p ValueDEG ItemGene Names
GO:0071555cell wall organizationBiological process9.34 × 10−103.25 × 10−69c23382_g2, c9513_g2, c9513_g1, c24866_g1, c25141_g1, c11939_g1, c54907_g1, c49924_g1, c26161_g1
GO:0045229external encapsulating structure organizationBiological process1.34 × 10−93.25 × 10−69c54907_g1, c26161_g1, c49924_g1, c23382_g2, c24866_g1, c25141_g1, c11939_g1, c9513_g1, c9513_g2
GO:0071554cell wall organization or biogenesisBiological process1.66 × 10−93.25 × 10−610c49924_g1, c26161_g1, c52977_g1, c54907_g1, c9513_g1, c9513_g2, c24866_g1, c25141_g1, c11939_g1, c23382_g2
GO:0005985sucrose metabolic processBiological process6.13 × 10−60.0054048c54907_g1, c54637_g1, c23382_g2, c25141_g1, c11939_g1, c24866_g1, c9513_g2, c9513_g1
GO:0005982starch metabolic processBiological process6.46 × 10−60.0054048c23382_g2, c25141_g1, c11939_g1, c24866_g1, c9513_g2, c9513_g1, c54907_g1, c54637_g1
GO:0005984disaccharide metabolic processBiological process8.99 × 10−60.0058498c54907_g1, c54637_g1, c23382_g2, c9513_g2, c9513_g1, c11939_g1, c25141_g1, c24866_g1
GO:0006073cellular glucan metabolic processBiological process1.2 × 10−50.0063958c24866_g1, c11939_g1, c25141_g1, c9513_g2, c9513_g1, c23382_g2, c54637_g1, c54907_g1
GO:0044042glucan metabolic processBiological process1.2 × 10−50.0063958c23382_g2, c9513_g1, c9513_g2, c25141_g1, c11939_g1, c24866_g1, c54907_g1, c54637_g1
GO:0009311oligosaccharide metabolic processBiological process1.66 × 10−50.0080848c23382_g2, c25141_g1, c11939_g1, c24866_g1, c9513_g2, c9513_g1, c54907_g1, c54637_g1
GO:0044264cellular polysaccharide metabolic processBiological process3.33 × 10−50.0130069c54637_g1, c24364_g1, c54907_g1, c24866_g1, c11939_g1, c25141_g1, c9513_g1, c9513_g2, c23382_g2
GO:0005976polysaccharide metabolic processBiological process5.27 × 10−50.0171349c25141_g1, c11939_g1, c24866_g1, c9513_g2, c9513_g1, c23382_g2, c54637_g1, c24364_g1, c54907_g1
GO:0044723single-organism carbohydrate metabolic processBiological process7.34 × 10−50.02263411c11939_g1, c25141_g1, c24866_g1, c9513_g1, c9513_g2, c562_g1, c23382_g2, c17871_g1, c54637_g1, c54907_g1, c24364_g1
GO:0044262cellular carbohydrate metabolic processBiological process0.000150.0438879c23382_g2, c9513_g1, c9513_g2, c25141_g1, c11939_g1, c24866_g1, c24364_g1, c54907_g1, c54637_g1
GO:0030312external encapsulating structureCellular component3.87 × 10−60.0054046c23382_g2, c54907_g1, c49924_g1, c11939_g1, c26161_g1, c24866_g1
GO:0005576extracellular regionCellular component7.52 × 10−60.00550711c56244_g1, c28395_g1, c55261_g1, c14539_g1, c49924_g1, c20661_g1, c24866_g1, c25141_g1, c20440_g1, c9513_g2, c9513_g1
GO:0005618cell wallCellular component2.31 × 10−50.0103975c11939_g1, c26161_g1, c24866_g1, c49924_g1, c54907_g1
GO:0071944cell peripheryCellular component4.68 × 10−50.0162578c49924_g1, c50518_g1, c26161_g1, c54907_g1, c11939_g1, c24866_g1, c23382_g2, c22601_g2
GO:0004650polygalacturonase activityMolecular function6.12 × 10−60.0054044c9513_g1, c9513_g2, c25141_g1, c24866_g1
GO:0004553hydrolase activity, hydrolyzing O-glycosyl compoundsMolecular function2.96 × 10−50.0123668c54637_g1, c25141_g1, c24866_g1, c9513_g1, c9513_g2, c23382_g2, c562_g1, c17871_g1
GO:0016798hydrolase activity, acting on glycosyl bondsMolecular function4.72 × 10−50.0162578c17871_g1, c562_g1, c23382_g2, c9513_g1, c9513_g2, c25141_g1, c24866_g1, c54637_g1
Table
Table 5. Results of DEGs enriched by Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway.
Table 5. Results of DEGs enriched by Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway.
#TermIDInput Numberp ValueCorrected p-ValueInput
Phenylalanine metabolismko0036030.0003370.006398357c9433_g1, c20661_g1, c20440_g1
Phenylpropanoid biosynthesisko0094030.0013980.013279387c9433_g1, c20661_g1, c20440_g1
Propanoate metabolismko0064010.0425470.179678515c17871_g1
Cutin, suberine and wax biosynthesisko0007310.0482130.179678515c24364_g1
Fatty acid elongationko0006210.0510330.179678515c21267_g1
Endocrine and other factor-regulated calcium reabsorptionko0496110.0603790.179678515c22601_g2
alpha-Linolenic acid metabolismko0059210.0833590.179678515c27465_g3
Photosynthesisko0019510.087890.179678515c28693_g1
Synaptic vesicle cycleko0472110.0905980.179678515c22601_g2
Parkinson’s diseaseko0501210.0950950.179678515c25708_g1
Glycerolipid metabolismko0056110.1040240.179678515c57011_g1
Peroxisomeko0414610.1320370.187274571c24364_g1
Pyruvate metabolismko0062010.1457270.187274571c17871_g1
Cysteine and methionine metabolismko0027010.1491170.187274571c17871_g1
Huntington’s diseaseko0501610.1617160.187274571c22601_g2
Glycerophospholipid metabolismko0056410.1642140.187274571c57011_g1
Oxidative phosphorylationko0019010.1749580.187274571c25708_g1
Glycolysis/Gluconeogenesisko0001010.1774180.187274571c17871_g1
Endocytosisko0414410.2087720.208772251c22601_g2

Share and Cite

MDPI and ACS Style

Shu, J.; Zhang, L.; Liu, Y.; Li, Z.; Fang, Z.; Yang, L.; Zhuang, M.; Zhang, Y.; Lv, H. Normal and Abortive Buds Transcriptomic Profiling of Broccoli ogu Cytoplasmic Male Sterile Line and Its Maintainer. Int. J. Mol. Sci. 2018, 19, 2501. https://0-doi-org.brum.beds.ac.uk/10.3390/ijms19092501

AMA Style

Shu J, Zhang L, Liu Y, Li Z, Fang Z, Yang L, Zhuang M, Zhang Y, Lv H. Normal and Abortive Buds Transcriptomic Profiling of Broccoli ogu Cytoplasmic Male Sterile Line and Its Maintainer. International Journal of Molecular Sciences. 2018; 19(9):2501. https://0-doi-org.brum.beds.ac.uk/10.3390/ijms19092501

Chicago/Turabian Style

Shu, Jinshuai, Lili Zhang, Yumei Liu, Zhansheng Li, Zhiyuan Fang, Limei Yang, Mu Zhuang, Yangyong Zhang, and Honghao Lv. 2018. "Normal and Abortive Buds Transcriptomic Profiling of Broccoli ogu Cytoplasmic Male Sterile Line and Its Maintainer" International Journal of Molecular Sciences 19, no. 9: 2501. https://0-doi-org.brum.beds.ac.uk/10.3390/ijms19092501

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop