Genome-wide identification of the restorer-of-fertility-like (RFL) gene family in Brassica napus and expression analysis in Shaan2A cytoplasmic male sterility

Background Cytoplasmic male sterility (CMS) is very important in hybrid breeding. The restorer-of-fertility (Rf) nuclear genes rescue the sterile phenotype. Most of the Rf genes encode pentatricopeptide repeat (PPR) proteins. Results We investigated the restorer-of-fertility-like (RFL) gene family in Brassica napus. A total of 53 BnRFL genes were identified. While most of the BnRFL genes were distributed on 10 of the 19 chromosomes, gene clusters were identified on chromosomes A9 and C8. The number of PPR motifs in the BnRFL proteins varied from 2 to 19, and the majority of BnRFL proteins harbored more than 10 PPR motifs. An interaction network analysis was performed to predict the interacting partners of RFL proteins. Tissue-specific expression and RNA-seq analyses between the restorer line KC01 and the sterile line Shaan2A indicated that BnRFL1, BnRFL5, BnRFL6, BnRFL8, BnRFL11, BnRFL13 and BnRFL42 located in gene clusters on chromosomes A9 and C8 were highly expressed in KC01. Conclusions In the present study, identification and gene expression analysis of RFL gene family in the CMS system were conducted, and seven BnRFL genes were identified as candidates for the restorer genes in Shaan2A CMS. Taken together, this method might provide new insight into the study of Rf genes in other CMS systems. Supplementary Information Supplementary information accompanies this paper at 10.1186/s12864-020-07163-z.

mitochondrial proteins [10]. To date, many other Rf genes have been identified in different CMS systems. Most of these Rf genes encode pentatricopeptide repeat (PPR) proteins. Examples of such Rf genes include Rf-PPR592 in petunia (Petunia hybrida) [11], Rfo [12] and orf687 [13] in radish (Raphanus sativus), Rf4 [14], Rf5 [15] and Rf6 [16] in rice (Oryza sativa) and Rfp [17] and Rfn [18] in rapeseed. The PPR proteins were first identified as tandem repeats of degenerate 35-amino-acid motifs (PPR motifs) in Arabidopsis thaliana [19] and were classified into PLS and P subfamilies, according to the PPR motif structure [20]. The PPR gene family is a large family comprising 441 members in Arabidopsis [21], 491 members in rice [22] and 626 members in poplar (Populus alba) [23]. Except for the PPR13 in sorghum (Sorghum bicolor), most of the RF-related PPR proteins belong to the P subfamily and lack the catalytic sites for RNA editing or binding [15]. Two partner proteins of the RF-related PPR proteins have been reported, including GRP162, which associates with RF5 [15], and hexokinase 6 (HXK6), which functions together with RF6 [16] to rescue CMS in rice.
To date, all reported Rf genes have been identified via genetic mapping, which is a time-consuming method and takes several years to narrow down the genomic region of interest. However, the concept of restorer-of-fertility-like (RFL) gene was put forward in 2011, and 212 RFL genes were identified based on BLAST searches using the Rf-PPR592 and Rf5 sequences against the genomes of 13 different dicot and monocot species, including Arabidopsis, soybean (Glycine max) and sorghum [24]. AtRFL2, together with RNase P, regulates the processing of mitochondrial orf291 RNA [25]. AtRFL4 is needed for processing the 5′-end of nad4 mRNA in mitochondria [26]. AtRFL9, also known as RNA PROCESSING FAC-TOR 4 (RPF4), participates in the generation of extra 5′ termini of ccmB transcripts in Arabidopsis [27]. These results enhanced our understanding of mitochondrial RNA processing in plants and provided novel insights into the function of RFL proteins.
In the present study, we performed BLAST searches using the sequences of Rf-PPR592 and AtRFLs against the genome of rapeseed and identified 53 BnRFL genes. Based on these 53 BnRFL genes, candidate Rf genes were analyzed in the Shaan2A CMS system by RNA-seq and tissue-specific expression analyses. Our data provide a strong foundation for the study of Rf genes in other CMS systems.

Identification of BnRFL genes
A total of 53 BnRFLs were identified in this study, based on the homology with the RFL genes in Arabidopsis (AtRFL1-26) and petunia (Rf-PPR592) ( Table 1). First, sequences of all 26 AtRFLs were searched in the database one at a time. Of the 26 AtRFL genes, nine showed no homologs in B. napus, including AtRFL5, AtRFL6, AtRFL9, AtRFL10, AtRFL14, AtRFL15, AtRFL16, AtRFL25 and AtRFL26. Then, using Rf-PPR592 as a reference [24], 26 BnRFLs were identified (E-value < 1e − 100 ). Taken together, there should be a total of 53 BnRFLs genes in B. napus. We also identified two known restorer genes, BnRFL6 (Rfn) and BnRFL13 (Rfp) (previously identified in the nap and pol CMS systems) and four candidate restorer genes (BnRFL2, BnRFL10, BnRFL11 and BnRFL42; previously identified in B. napus by fine genetic mapping) [31].
The number of PPR motifs in the BnRFL proteins varied from 2 to 19, although most of the BnRFL proteins contained at least 10 PPR motifs and the average number of PPR motifs was 12 (Table 1). Approximately one-fifth of the BnRFLs showed relatively low pI (< 6), whereas nearly half of the BnRFLs showed relatively high pI (> 8). The molecular weight of these RFL proteins ranged from 11.7-92.4 kDa. Additionally, the GRAVY value of nearly two-fifth BnRFLs and most of the selected Rf genes was less than 0, indicating that these RFL proteins were hydrophilic. Most of the BnRFLs were predicted to localize to the mitochondria, which is consistent with the subcellular localization of RF proteins (Table 1).

Chromosomal location and structural analysis of BnRFL genes
First, we downloaded the chromosomal distribution of AtRFLs from The Arabidopsis Information Resource (TAIR) (Fig. 1a). All 26 AtRFLs were located in a cluster on chromosome 1. Of the 53 BnRFLs identified in this study, 46 BnRFLs were distributed unevenly on 10 of the 19 chromosomes, and 18 and 10 BnRFLs formed highly dense clusters on chromosomes A9 and C8, respectively. The remaining seven BnRFLs were located on the unmapped scaffold (Fig. 1b and c).
Next, we determined the exon-intron structure of the BnRFL genes and a few known restorer genes (Additional file 1). Most of the BnRFLs were intron-less, similar to the restorer genes, such as Rf4, Rf5 and Rf6, in rice CMS line. Of the 53 BnRFL genes, 10 contained a single intron, similar to the Rf genes, Rfk1, Rfob and orf687, in radish. Notably, the intron in BnRFL23 was more than 4 kb in length, unlike other BnRFLs.
Because PPR proteins generally contain tandem repeats of PPR motifs, we searched for the PPR motifs in the BnRFL proteins and a few known restorer proteins (Additional file 2). To investigate whether BnRFLs contained additional motifs, 53 BnRFLs and 9 known restorer proteins from other species were submitted to MEME. The results showed 20 motifs in the BnRFL proteins (Fig. 2). Interestingly, all of the identified RFL

Phylogenetic and Syntenic analysis
To identify the homologs of BnRFLs in different monocot and dicot species, multiple sequence alignments were performed and sequence similarity was determined. The rice RF5 protein was used for BLAST searches against the rice and maize genomes, and Rf-PPR592 was used for BLAST searches against the radish genome. An additional 16 OsRFLs (including 4 reported restorer genes), 9 ZmRFLs, and 22 RsRFLs (including 4 known restorer genes) were identified (E-value <1e − 100 ). Phylogenetic analysis revealed that RFLs mainly formed two separate clusters, and RFLs in monocot and dicot species were clustered together, respectively (Fig. 3a). Additionally, two reported restorer genes (BnRFL6 and BnRFL13) and four candidate restorer genes (BnRFL2, BnRFL10, BnRFL11 and BnRFL42) clustered together. Six additional BnRFLs (BnRFL3, BnRFL4, BnRFL5, BnRFL8, BnRFL15 and BnRFL41) clustered together with the reported and candidate restorer genes. These 12 BnRFLs have been deeply investigated in the following study. Next, we examined the synteny of BnRFLs with their homologs in Arabidopsis, B. rapa and B. oleracea (Fig. 3b). The results showed syntenic relationships between AtRFL genes on chromosome 1 and RFL genes on chromosomes BraA8, BraA9, BolC8, BnaA8, BnaA9, BnaC3, BnaC8 and BnaC9. The AtRFL gene on chromosome 3 showed synteny with RFL genes on BraA1 and BnaA1. The AtRFL genes on chromosome 4 showed synteny with RFL genes on BraA1, BolC1, BnaA1 and BnaC1, and the AtRFL genes on chromosome 5 showed synteny with RFL genes on BraA4, BnaA4, BnaC4 and BnaC9.
The Ks and Ka values indicate the evolutionary pressure on species. A Ka/Ks ratio < 1indicates functional constraint, whereas Ka/Ks ratio > 1 indicates positive selection [32]. To explore the selection pressure on BnRFLs, we calculated the Ka/Ks ratios (Additional file 3).

Interaction analysis of AtRFL proteins
Most of the RFL proteins belong to the P subfamily and need to interact with other proteins to perform RNA processing [15]. To predict the interacting partners of RFL proteins (no B. napus information in STRING database), an interaction network for AtRFLs were constructed based on STRING 10.5 and Cytoscape 3.6.1. Except AtRFL10, which did not have interaction information in the STRING database, 25 AtRFLs and their predicted partners are shown in Fig. 4 and Additional file 4. Interestingly, AtRFL11, AtRFL12 and AtRFL13 and three HXKs, including HXK1, HXK2 and HXK3, shared a common interacting protein, namely replication factor C2 (RFC2), a multi-subunit complex critical for high-speed ATP-dependent DNA synthesis [33]. No homologs of AtRFL25 were identified in B.

Expression analysis of BnRFL genes
Based on the results of phylogenetic analysis, 12 BnRFL genes, which were mentioned in the phylogenetic analysis, were selected for tissue-specific expression analysis in the sterile line Shaan2A, the maintainer line Shaan2B and the restorer line KC01 by qRT-PCR. Although BnRFL10 and BnRFL11 were located on different chromosomes, the coding sequence (CDSs) of these genes were highly similar (identity = 1926/1947; 99%), and it was difficult to distinguish them by qRT-PCR. Therefore, we finally analyzed the expression of 11 BnRFL genes. In the restorer line, the expression of BnRFL6, BnRFL13 and BnRFL42 was lower in leaves than in the perianth (Fig. 5a, Additional file 5). The majority of the selected BnRFLs showed higher expression level in MA when compared with leaves, except for BnRFL2, BnRFL3 and BnRFL4. However, the expression of the BnRFL genes, except BnRFL41, was lower in the gynoecium and LA when compared with leaves. Compared with Shaan2A tissues, the expression of 11 BnRFL genes was higher in KC01 tissues, especially in MA (Fig.  5b, Additional file 5). What's more, the expression of these genes in Shaan2B LA was higher than those in Shaan2A (Fig. 5b, Additional file 5). However, the expression of most of these BnRFLs was lower in the gynoecium of the restorer line than in that of Shaan2A.
To compare the expression level of all 53 BnRFL genes in Shaan2A vs. KC01, three biological replicates of RNA-seq were performed using RNA isolated from young buds (YB, < 1 mm, representing pre-meiosis) and small anthers (SA, sampled from buds 1-2 mm in length (representing tetrad stage to microspore release stage). A total of 320,892,232 raw sequence reads were generated, with approximately 50 million raw reads representing each tissue sample (SRA number: PRJNA511929). Additionally, to conduct comparative transcriptome analysis of the three lines in the Shaan2A CMS system, raw transcriptome reads representing the same stages of Shaan2A and Shaan2B were downloaded (SRA number: PRJNA502996). RNA-seq data analysis of Shaan2A and KC01 revealed only nine BnRFLs exhibited differential expression (Fig. 5c). These results provide important clues for analyzing the candidate restorer genes in the Shaan2A CMS system.

Transcriptomic analysis between Shaan2A and KC01
To investigate the differences between Shaan2A and KC01, possibly caused by the male sterile genes and restorer genes, we also identified the differentially expressed genes (DEGs) between Shaan2B and KC01, as Shaan2A and Shaan2B share the same nuclear genetic background. Thus, common DEGs identified based on Shaan2A vs. KC01 comparison and Shaan2B vs. KC01 comparison would represent the DEGs identified between different genetic backgrounds, i.e., Shaan2A (or Shaan2B) and KC01. A total of 2980 and 8243 DEGs were identified in YB and SA stage, respectively, based on the comparison between Shaan2A and KC01 (|log2 Ratio| > 1; Additional file 6). Based on GO analysis, only one GO term in the molecular function category, 'sequence-specific DNA binding', was significant at YB stage (Fig. 6a). By contrast, at SA stage, 8243 DEGs identified between Shaan2A and KC01 were categorized in three main categories, including molecular function, cellular component and biological process, which were further classified into many functional sub-categories (Additional file 7). The top 30 sub-categories, including 'disaccharide metabolic process', 'regulation of RNA biosynthetic process', 'regulation of RNA metabolic process' and 'regulation of transcription, DNA-dependent', are shown in Fig. 6b.
To validate the RNA-seq data, the expression of 11 DEGs, potentially involved in anther development, was analyzed by qRT-PCR (Additional file 8). The results of qRT-PCR analysis for most of these DEGs at the two stages were consistent with those of RNA-seq analysis, indicating that the reliability of our RNA-seq data.

Discussion
Many CMS systems have been used in B. napus, including pol, nap, Ogu and Shaan2A [4][5][6][7]. To date, Rf genes have been identified via genetic mapping only in pol and nap CMS systems [17,18]. In the present study, 53 BnRFL genes were identified in the Shaan 2A CMS system. Most of these genes contained more than 10 PPR motifs, which is consistent with the previously reported restorer proteins such as Rf-PPR592 (14 PPR motifs), Rf4 (18 PPR motifs) and Rfo (17 PPR motifs) [11,17,24]. Moreover, most of the BnRFLs identified in this study were predicted to localize to mitochondria, similar to the known restorer genes [15][16][17][18]. Since the Rf genes function with the toxic chimeric genes in mitochondria to rescue male sterility [14][15][16], the mitochondrial localization of the proteins seems appropriate. More importantly, we also identified BnRFL6 (Rfn) and BnRFL13 (Rfp), previously confirmed as restorer genes in nap and pol CMS systems, respectively [17,18], and four candidate restorer genes (BnRFL2, BnRFL10, BnRFL11 and BnRFL42), previously identified in B. napus via genetic mapping [31]. Taken together, analysis of the RFL gene family for the identification of candidate restorer genes were viable, which would also provide a new way to analysis the restorer genes in other CMS systems as one supplementary method except for the traditional genetic mapping to locate the candidate genes.
Nearly 7500 years ago, B. napus originated from the hybridization of B. rapa and B. oleracea [34], and the Brassica plants experienced the extra whole genome triplication (WGT) event when compared with Arabidopsis [35]. The Arabidopsis genome contains 26 RFL gene family members, so considering the WGT event there should be over 78 RFL genes in B. oleracea or B. rapa genome, and finally generate even more RFL genes in B. napus. While only 53 BnRFLs were identified in the present study, which implied that nearly 50% RFL genes were lost after the WGT event.
Most of the BnRFLs were unevenly distributed on 10 of the 19 chromosomes of B. napus, while a few formed gene clusters on chromosomes A9 and C8, similar to the gene cluster in Arabidopsis (chromosome 1; Table 1, Fig.  1), rice and barley (Hordeum vulgare) [36,37]. Gene clustering has also been observed in many other gene families, such as the LEA gene family in B. napus [38] and Phyllostachys edulis [39] and laccase gene family in Citrus sinensis [40]. The LEA gene clusters on B. napus chromosomes A9 and C4 probably resulted from chromosomal rearrangement during the evolution of Brassica species [38]. The RFL genes on Arabidopsis chromosome 1 showed synteny with the RFL genes on BraA9, BolC8, BnaA9 and BnaC8. Additionally, BnRFLs maintained a syntenic relationship with RFL genes in B. rapa and B. oleracea, suggesting that a conserved role of BnRFLs located on chromosomes A9 and C8. Moreover, AtRFL2, AtRFL4 and AtRFL9 were located within the gene cluster on Arabidopsis chromosome 1 and participated in the processing of the mitochondrial RNA [25][26][27]. The phylogenetic analysis revealed that the BnRFLs have the closer phylogenetic relationship with AtRFLs and RsRFLs (Fig. 3), and the structural analysis showed that all of the BnRFLs and the known restorer genes in radish share a conserved motif (Fig. 2), and all BnRFL genes showed a Ka/Ks ratio < 1 (Additional file 3), which indicated that there was no positive selection on the BnRFL genes during the evolution. What's more, BnRFL6 (Rfn) and BnRFL13 (Rfp) were located within the gene cluster on chromosome A9. These data suggest that the RFL genes within gene clusters on chromosomes A9 and C8 represent the restorer genes in the CMS system, as these likely exhibit a conserved role in mitochondrial RNA processing.
Tandem repeats of a degenerate 35-amino-acid PPR motif are the most prominent feature of the PPR family, and all of the 53 BnRFL proteins showed this trait. Although 212 RFL genes in 13 different species [24] and 26 RFL genes in barley [37] have been identified previously, the conserved domain of the RFL proteins has not yet been analyzed. Therefore, we investigated motifs other than PPR in the RFL proteins, revealing 20 motifs among the 53 BnRFLs and a few known RF proteins. Interestingly, all RFLs contained motif 1, comprising 80 amino acids. We propose motif 1 as the conserved motif in the RFL protein family. This motif will serve as a reference for RFL family analysis in other species.
Because RF-related PPR proteins belong to the P subfamily and do not exhibit endonuclease activity, these proteins form functional complexes with other proteins [20]. To date, only two RFL-interacting partner proteins have been identified, including GRP162 and HXK6 in the rice CMS system [15,16]. In the present study, we constructed an interaction network for AtRFLs (Fig. 4). Interestingly, AtRFL11, AtRFL12 and AtRFL13 and HXKs (HXK1, HXK2 and HXK3) shared a common interacting partner, RFC2, which was critical for highspeed DNA synthesis [22], whereas AtRFL25 showed interaction with GRP7 and GR-RBP2. Moreover, AT1G48510, SURF1, COX11 and COX15 were predicted to interact with most of the AtRFLs. AT1G48510 is a surfeit locus 1 cytochrome c oxidase biogenesis protein. SURF1 is associated with cytochrome c oxidase assembly [41]. Both COX11 and COX15 are mitochondrial proteins and belong to the cytochrome c oxidase protein family [42]. COX11 likely plays a key role as a mitochondrial chaperone in the assembly of the COX complex and regulates pollen germination and plant growth [43]. Overall, the interaction network indicates possible partner proteins of RFL proteins in Arabidopsis. These data provide important clues for the identification of interaction factors of RF proteins in other species.
Previously, it has been shown that Rf4 is constitutively expressed in different rice organs at relatively low levels [14]. Although Rf6 expression is detectable in various rice tissues, it is expressed to a higher level in the panicle than in other tissues [16]. In the pol CMS system, Rfp shows relatively high expression in flower buds and weak expression in opening flowers, leaves, stems and roots [17]. In the phylogenetic tree, two previously reported restorer genes (BnRFL6 and BnRFL13), four candidate restorer genes (BnRFL2, BnRFL10, BnRFL11 and BnRFL42) and another six BnRFLs (BnRFL3, BnRFL4, BnRFL5, BnRFL8, BnRFL15 and BnRFL41) clustered together, suggesting these genes as the more probable candidates of restorer genes in the rapeseed CMS system. Analysis of expression patterns revealed that most of these genes, except for BnRFL2, BnRFL3 and BnRFL4, were expressed to relatively higher levels in MA than in leaves in KC01. Additionally, these BnRFLs showed higher expression in KC01 tissues, especially MA, than in Shaan2A tissues. Expression profiling of BnRFL genes in Shaan2A vs. KC01 showed that BnRFL1, BnRFL6, BnRFL10, BnRFL13, BnRFL15, BnRFL38 and BnRFL49 were down-regulated in Shaan2A. However, BnRFL15 only harbored five PPR motifs, which was much lower than the number of PPR motifs in the known RF proteins. While BnRFL38 and BnRFL49 are located on chromosomes A1 and A10, respectively, BnRFL1, BnRFL5 and BnRFL8, BnRFL6, BnRFL11, BnRFL13 and BnRFL42 are located in gene clusters on chromosomes C8 and A9. Interestingly, almost all of the known rice Rf genes are located in the RFL gene cluster on chromosome 10 [36]. These data suggest BnRFL1, BnRFL5, BnRFL6, BnRFL8, BnRFL11, BnRFL13 and BnRFL42 as the more likely candidates of restorer genes in the Shaan2A CMS system. In O. sativa, RF6 with a characteristic duplication of PPR motifs in the restorer line of Honglian CMS can restore sterility, while the duplicated motifs are absent in rf6 of sterile line [16]. In the next steps, we will clone these candidate restorer genes in the restorer line and sterile line of Shaan2A CMS respectively, and compare the difference of sequences between them, for we wonder if there is the similar motif difference between these candidate restorer genes. Then we will narrow down the list of candidate genes, and conduct the transgenic work in sterile line to investigate their function.
Furthermore, DEGs identified in small anthers of Shaan2A vs. KC01 were annotated as involved in the 'regulation of RNA biosynthetic process', 'regulation of RNA metabolic process' and 'regulation of transcription, DNA-dependent'. The RF-related PPR proteins interact with their partner proteins to bind or to edit RNA [15]. Here, the regulation of RNA biosynthetic, RNA metabolic process and transcription was different between the sterile line and restorer line, which might be caused by the sterile genes in Shaan2A and restorer genes in KC01. However, the detailed mechanism needs further investigation.

Conclusions
In CMS, the Rf nuclear genes rescue the sterile phenotype and most of the Rf genes encode pentatricopeptide repeat (PPR) proteins. In the present study, a total of 53 BnRFL genes were identified in B. napus. Most of the BnRFL genes were distributed on 10 of the 19 chromosomes, and gene clusters were identified on chromosomes A9 and C8. The interaction network analysis was performed to predict the interacting partners of RFL proteins. Tissue-specific expression and RNA-seq analyses between the restorer line KC01 and the sterile line Shaan2A indicated that BnRFL1, BnRFL5, BnRFL6, BnRFL8, BnRFL11, BnRFL13 and BnRFL42 located in gene clusters on chromosomes A9 and C8 were highly expressed in KC01, which suggest these seven BnRFL genes as strong candidates for the restorer genes in Shaan2A CMS. Our results would provide new insight into the study of Rf genes in other CMS systems.

Plant materials
The sterile line Shaan2A, maintainer line Shaan2B and restorer line KC01 of B. napus, gifted by Professor Dianrong Li at the Hybrid Rape Research Center of Shaanxi Province, were used in this study. Plants were cultivated on the experimental field of the Huazhong University of Science and Technology (Wuhan, Hubei province, China). After harvest, plant samples were immediately frozen in liquid nitrogen and stored at − 80°C until needed for total RNA isolation.

Structural analysis of RFL genes
The exon-intron structure of BnRFL genes and a few known Rf genes were based on the alignments of the CDS with the corresponding genomic sequences, and the diagram was conducted using the Gene structure display server (GSDS, http://gsds.cbi.pku.edu.cn/) [50]. The PPR motifs in all BnRFL proteins and a few known RF proteins were analyzed using the NCBI Conserved Domain Search tool (http://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi) [51]. Conserved motifs in RFL proteins were analyzed using Multiple Expectation Maximization for Motif Elicitation (MEME, http://alternate.meme-suite.org/) [52].

Phylogenetic and Syntenic analysis of RFL genes
Multiple sequence alignment of the predicted amino acid sequences of BnRFLs, AtRFLs, RsRFLs, OsRFLs and ZmRFLs was performed using ClustalX [53]. A phylogenetic tree of these RFL proteins was constructed with MEGA 7 using the Neighbor Joining (NJ) method [54]. Analysis of synteny among BnRFL, AtRFL, BoRFL and BrRFL genes was performed using the syntenic gene tool in the Brassica database (BRAD, http://brassicadb.org/ brad/) [55]. The non-synonymous to synonymous nucleotide substitution ratio (Ka/Ks) was calculated using TBtools [56].

Interaction analysis
The interaction analysis of AtRFLs was based on the STRING 10.5 database, which included the known and predicted protein-protein interactions. First, the interaction proteins of AtRFLs were searched. After deleting the repeat proteins, the interaction network was visualized using Cytoscape 3.6.1.

RNA extraction, RNA-seq and qRT-PCR analysis
Gene expression was analyzed in various tissues of the sterile line Shaan2A and restorer line KC01 including leaves, perianths, gynoecium, medium anthers (MA) and large anthers (LA). MA were harvested from buds 2-4.5 mm in length and represented the uninuclear microspore stage, whereas LA were harvested from buds 4.5 mm in length and represented the mature pollen formation stage.
Total RNA extraction, RNA-seq analysis and qRT-PCR were conducted according to our previous protocols [57], with minor modifications. Briefly, approximately 100 mg plant samples were used for total RNA extraction using TRIzol Reagent (Invitrogen, Carlsbad, CA, USA), according to the manufacturer's instructions. Then, cDNA sequencing libraries were constructed using TruSeq™ RNA Sample Preparation Kit (Illumina, San Diego, CA, USA). RNA-seq was performed on the Illumina NovaSeq 6000 platform. The raw data were filtered using the NGSQC toolkit (v2.2.3), and the clean reads were mapped to the reference genome of the rapeseed cultivar 'ZS11'. The differentially expressed genes (DEGs) were evaluated using DESeq2, with normalized fold-change ≥2 and p-value < 0.05. Gene Ontology (GO) annotation was using the Web Gene Ontology Annotation Plot (WEGO) software.
To perform qRT-PCR analysis, RNA was reverse transcribed using the TaKaRa PrimeScript™ RT Reagent Kit with gDNA Eraser, according to the manufacturer's instructions. Actin was used as the internal reference gene [58]. The qRT-PCR experiments and transcript quantification were performed as described previously [57]. Primers used in this study are listed in Additional file 9.