Skip to main content

Transcriptome analysis reveals the potential lncRNA-mRNA modules involved in genetic male sterility and fertility of Chinese cabbage (brassica rapa L. ssp. pekinensis)

Abstract

Background

Long non-coding RNAs (lncRNAs) play a crucial role in regulating gene expression vital for the growth and development of plants. Despite this, the role of lncRNAs in Chinese cabbage (Brassica rapa L. ssp. pekinensis) pollen development and male fertility remains poorly understood.

Results

In this study, we characterized a recessive genic male sterile mutant (366–2 S), where the delayed degradation of tapetum and the failure of tetrad separation primarily led to the inability to form single microspores, resulting in male sterility. To analyze the role of lncRNAs in pollen development, we conducted a comparative lncRNA sequencing using anthers from the male sterile mutant line (366–2 S) and the wild-type male fertile line (366–2 F). We identified 385 differentially expressed lncRNAs between the 366–2 F and 366–2 S lines, with 172 of them potentially associated with target genes. To further understand the alterations in mRNA expression and explore potential lncRNA-target genes (mRNAs), we performed comparative mRNA transcriptome analysis in the anthers of 366–2 S and 366–2 F at two stages. We identified 1,176 differentially expressed mRNAs. Remarkably, GO analysis revealed significant enrichment in five GO terms, most notably involving mRNAs annotated as pectinesterase and polygalacturonase, which play roles in cell wall degradation. The considerable downregulation of these genes might contribute to the delayed degradation of tapetum in 366–2 S. Furthermore, we identified 15 lncRNA-mRNA modules through Venn diagram analysis. Among them, MSTRG.9997-BraA04g004630.3 C (β-1,3-glucanase) is associated with callose degradation and tetrad separation. Additionally, MSTRG.5212-BraA02g040020.3 C (pectinesterase) and MSTRG.13,532-BraA05g030320.3 C (pectinesterase) are associated with cell wall degradation of the tapetum, indicating that these three candidate lncRNA-mRNA modules potentially regulate pollen development.

Conclusion

This study lays the foundation for understanding the roles of lncRNAs in pollen development and for elucidating their molecular mechanisms in regulating male sterility in Chinese cabbage.

Peer Review reports

Introduction

Chinese cabbage is an important vegetable crop in the Brassica family of Cruciferae. It undergoes cross-pollination through bisexual flowers and exhibits significant heterosis. Utilizing a male sterile line serves as an optimal method for hybrid seed production [1, 2]. Understanding the genes and molecular mechanisms associated with male sterility is pivotal for advancing Chinese cabbage hybrid breeding and thus holds substantial economic implications.

Male sterility primarily involves anther development, a delicate and complex process encompassing anther wall formation and pollen grain generation. The anther wall comprises layers such as the epidermis, endoderm, middle layer, and tapetum, progressing from the outer layers to the inner layers. Tapetum cells transform into secretory cells post-formation, while haploid microspores become enveloped by callose. As the anther grows, the microspores released from the tetrads gradually mature into pollen grains within the anther chamber [3,4,5]. Throughout anther development, processes like tapetum formation and apoptosis, meiosis, callose production and degradation, and pollen wall formation intricately influence pollen development. These tightly regulated processes, if disrupted, lead to varying degrees of reduced male fertility [6].

Long non-coding RNAs (lncRNAs), extending beyond 200 nucleotides, wield significant influence across diverse life processes, modulating crucial biological pathways, including plant flowering, photomorphogenesis, reproduction, and responses to abiotic/biological stress [7,8,9]. Increasingly, studies highlight the involvement of numerous lncRNAs in male reproductive development in plants. In maize, Zm401 plays a role in anther development by functioning as short open reading frame mRNA (sORF mRNA) and/or non-coding RNA (ncRNA). Knockout of Zm401 significantly affects the expression of key pollen development genes (ZmMADS2, Zm3-3, and ZmC5), resulting in abnormal development of microspore and tapetum and ultimately male infertility [10,11,12]. In rice, the long-day-specific male-fertility-associated RNA (LDMAR) has been identified as an lncRNA capable of regulating photosensitive male sterility. Under extended daylight conditions, highly expressed LDMAR maintains normal pollen development [13]. Conversely, mutations in ldmar alter the secondary structure of LDMAR, causing increased methylation in the LDMAR promoter region. This reduction in LDMAR transcription levels initiates premature programmed cell death of anther tapetum cells under extended daylight, leading to photosensitive male sterility [14]. Furthermore, the photosensitive sterility gene PMS1 generates the transcription of lncRNA PMS1T in rice. This lncRNA, targeted and activated by miR2118, generates 21-nt phasiRNAs that preferentially accumulate in photosensitive male sterility under prolonged sunlight conditions, actively participating in rice reproductive development [15].

Similarly, several lncRNAs have been implicated in pollen development within Chinese cabbage. LncRNA BcMF11 is specifically expressed in anthers. Suppression of BcMF11 results in delayed tapetum degradation, ultimately affecting pollen development [16, 17]. A recent study of Brassica highlighted the expression of 12,051 lncRNAs during pollen development and fertilization [18]. Another investigation identified 4,347 lncRNAs and 2,045 lncNATs across five stages of pollen development, utilizing genome-wide identification methods [19]. The utilization of diverse male sterile materials in identifying lncRNAs might reveal corresponding anther development processes. Moreover, lncRNAs often act as regulators of the transcription factors or protein-coding genes that participate in the formation of male sterility. Comprehensive studies are undoubtedly important for interpreting lncRNA-mediated gene expression networks linked to anther development and male sterility.

In this study, a new Chinese cabbage genic male sterile line (366–2 S) was characterized. During anther development, the tapetum of 366–2 S remained intact for a long time and presented delayed degradation. Additionally, the tetrads of 366–2 S remained unseparated and clustered together, resulting in the absence of pollen release and subsequent male sterility. This study was carried out to identify differentially expressed lncRNAs and mRNAs involved in anther development through lncRNA and mRNA sequencing using fertile line 366–2 F and sterile line 366–2 S. Furthermore, to find out the potential lncRNA-mRNA modules associated with male sterility, target gene prediction analysis of lncRNAs and GO enrichment analysis of mRNAs were performed. This study will broaden our knowledge of the role of lncRNAs and the molecular mechanisms underpinning male sterility in Chinese cabbage.

Results

Characterization of abnormal anther development in 366–2 S

During the pollination phase, the sterile line 366–2 S was distinguished from the progenies of the heterozygous fertile line 366-2. Anther appearance within flowers exhibited notable differences between the 366–2 F and 366–2 S lines (Fig. 1). Both 366–2 F and 366–2 S showed typical flower traits, blooming normally with four petals (Fig. 1, C, J), four sepals (Fig. 1, D, K), and one pistil adhering to six stamens (Fig. 1, E, L). However, distinct variations were observed in the anthers of these lines. In 366–2 F, the anthers appeared yellow and plump (Fig. 1, G), exhibiting normal dehiscence and the release of active pollens. Conversely, the anthers of 366–2 S appeared gray and withered (Fig. 1, N), failing to release pollens, indicating a lack of pollen development within 366–2 S.

Fig. 1
figure 1

Morphological characteristics of floral organs of Chinese cabbage at full flowering stage of 366–2 S and 366–2 F lines. (A–G) A floret at the anthesis stage with normal flower organs in 366–2 F. (H–N) A floret at the anthesis stage with shorter filaments and anthers without pollen in 366–2 S. Bar = 5 mm

To understand the abnormal pollen development in 366–2 S, we meticulously observed the anthers at five consecutive developmental stages (meiosis, tetrad, uninucleate, bicellular, and tricellular) utilizing HE staining. At the meiosis and tetrad stages (Fig. 2, A, B, F, G), no obvious differences were discernible between 366 and 2 F and 366–2 S. Both exhibited similar and normal tapetums, microspore mother cells (MMCs), and tetrads. However, notable deviations emerged at the uninucleate and bicellular stages. In 366–2 F, the tetrads dissociated normally, releasing microspores (Fig. 2, C, D). Conversely, in 366–2 S, the tetrads remained unseparated, clustering together (Fig. 2, H, I), resulting in the absence of microspore formation. At the tricellular stage, microspores were spherical, and the anther chambers were filled in 366–2 F (Fig. 2, E). In contrast, anther chambers in 366–2 S showed atrophy, with only a few contents remaining (Fig. 2, J), indicating signs of male sterility. At the bicellular and tricellular stages, the tapetum underwent gradual degradation, ceasing the provision of necessary molecules for microspore development. Nevertheless, the tapetum of 366–2 S appeared more robust compared to that in 366–2 F (Fig. 2). This delayed degradation of the tapetum might contribute to male sterility in 366–2 S.

Fig. 2
figure 2

Paraffin sections of Chinese cabbage anthers of 366–2 F and 366–2 S. (A–E) 366–2 F pollen microspore development process. (F–J) 366–2 S pollen microspore development process; the tetrad period began to develop abnormally. T: Tapetum; MMC: Microspore mother cell; Td: Tetrad; M: Microspore. Bar = 20 μm

Identification of lncRNAs and prediction of their target genes in anthers

To analyze the role of lncRNAs in pollen development of both 366–2 S and 366–2 F, a comparative lncRNA-seq was conducted using anthers harvested from unopened buds of varying lengths. A substantial quantity of clean reads were obtained, totaling 92,287,338 for 366–2 S and 92,214,922 for 366–2 F (Table 1). These reads exhibited high quality, surpassing Q20 and Q30 thresholds of over 99%, with mapped rates exceeding 83% (Table 1). By aligning with the reference genome and conducting CPC coding potential prediction, 1,781 lncRNAs were identified in 366–2 S and and 2,057 lncRNAs were identified in 366–2 F. Ultimately, a total of 2 075 lncRNAs were obtained (Fig. 3, A and Supplementary Table S2). The target genes were predicted through upstream and downstream gene sequence alignment and base pairing. Among these, 1,164 lncRNAs had both upstream and downstream target genes, and 103 lncRNAs exhibited base pairing target genes, while 808 lncRNAs had no corresponding target genes (Fig. 3, B and Supplementary Table S2), resulting in the prediction of a total of 2,376 target genes (Supplementary Table S2). Differential expression analysis of lncRNAs between 366 and 2 S and 366–2 F revealed significant differential expression in 385 lncRNAs within 366–2 S, with 299 and 86 showing upregulation and downregulation, respectively (Fig. 3, C and Supplementary Table S2). Correspondingly, the predicted target genes associated with these DE-lncRNAs were identified (Fig. 3, D and Supplementary Table S2). These DE-lncRNAs and their corresponding predicted target genes are implicated in the abnormal pollen development observed in 366–2 S.

Table 1 LncRNA sequencing data of sterile and fertile Chinese cabbage lines
Fig. 3
figure 3

LncRNA sequencing in anthers of 366–2 S and 366–2 F. (A) Identification of lncRNAs in 366–2 S and 366–2 F. (B) Target gene prediction of lnRNAs. (C) Identification of DE-lncRNAs between 366–2 S and 366–2 F. (D) Target gene prediction of DE-lncRNAs between 366–2 S and 366–2 F

Differential expression of mRNAs in anther development

The anthers of small buds (about 1  1.5 mm) in 366–2 S showed signs of sterility. To understand the differences in mRNA expression in anther development of 366–2 S and 366–2 F, comparative transcriptome analysis was performed at the early-stage (when bud length was approximately 1  1.5 mm) and the late-stage (when bud length was approximately 3.5  5 mm) of anther development. The sequencing yielded an average of 42,154,243; 41,243,279; 41,333,469 and 41,328,565 clean reads with 85.76%, 86.67%, 88.03%, and 87.86% of reads successfully mapped to the reference genome for E-366–2 S, E-366–2 F, L-366–2 S, and L-366–2 F, respectively (Table 2). A comprehensive total of 40,482 mRNAs were identified and annotated (Supplementary Table S3). The expression levels of mRNAs exhibited a highly significant correlation among the three repetitions in each treatment (Supplementary Fig. S1), affirming the reliability of the sequencing data. Subsequently, 6,982 mRNAs (4,055 upregulated and 2,927 downregulated) were identified in the early group (E-366–2 S vs. E-366–2 F), along with 13,224 mRNAs (8,625 upregulated and 4,599 downregulated) in the late group (L-366–2 S vs. L-366–2 F) (Fig. 4, A). Furthermore, 1,176 consistent mRNAs were identified from the two groups (Fig. 4, B). The observation indicated that the pollen of 366–2 S was abnormal compared to that of 366–2 F at both the early and late stages. Thus, these consistent DE-mRNAs are anticipated to be associated with pollen development and male fertility.

Table 2 mRNA sequencing data of 12 samples
Fig. 4
figure 4

mRNA sequencing in anthers of 366–2 S and 366–2 F. (A) Identification of DE-mRNA between 366–2 S and 366–2 F. (B) Venn analysis of DE-mRNAs at the early and late stages

Go enrichment analysis of DE-mRNAs

Based on the BP, GO enrichment analysis was performed for fertility-related DE-mRNAs. Five GO terms, i.e., GO: 0045229 (external encapsulating structure organization), GO: 0071555 (cell wall organization), GO: 0042545 (cell wall modification), GO: 0005975 (carbohydrate metabolic process), and GO: 0071554 (cell wall organization or biogenesis), were significantly enriched (Fig. 5, A, Supplementary Table S4). Significantly enriched GO terms contained a total of 87 DE-mRNAs (Fig. 5, B). Among them, in GO: 0045229, GO: 0071555, GO: 0042545, and GO: 0071554, 20 mRNAs encode pectinesterases; in GO: 0005975, 17 mRNAs encode polygalacturonases, all of which belong to pectinase regulating the degradation of pectins in the cell wall. In addition, many other mRNAs in GO: 0005975 encode some glycosidases, such as galactosidase (BraA03g003510.3 C, BraA06g027250.3 C, BraA01g003370.3 C, BraA03g059130.3 C, BraA03g059130.3 C, and BraA10g019920.3 C) and glucosidase (BraA02g016980.3 C, BraA02g016980.3 C, BraA02g016980.3 C, BraA02g016980.3 C, BraA02g016980.3 C, BraA05g017770.3 C, BraA06g002000.3 C, BraA05g033960.3 C, BraA05g033960.3 C, BraA09g049970.3 C, and BraA09g049950.3 C). They might be involved in the hydrolysis of polysaccharides in the cell wall. Notably, most of these mRNAs showed a downregulation pattern in 366–2 S (Fig. 5, B), which potentially impaired the ability of cell wall degradation.

Fig. 5
figure 5

GO enrichment analysis of differentially expressed mRNAs. (A) Five GO terms associated with the anther cell wall were significantly enriched. (B) A total of 87 DEGs were contained in the significantly enriched GO.terms

Major lncRNA-mRNA modules associated with pollen fertility

The role of lncRNAs in regulating the expression of target genes is pivotal in various BPs. A Venn diagram (Fig. 6, A, Supplementary Table S5) revealed 15 DE-mRNAs that corresponded to the target genes of DE-lncRNAs. Among these, certain lncRNAs exhibited consistent expression trends with their target mRNAs, exemplified by MSTRG.9997-BraA04g004630.3 C; while others showed opposing trends, such as MSTRG.8514-BraA03g042940.3 C (Fig. 6, B). This indicates the potential of lncRNAs to either activate or inhibit their gene expressions. Notably, three target mRNAs, namely, BraA04g004630.3 C (β-1,3-glucanase), BraA02g040020.3 C (pectinesterase), and BraA05g030320.3 C (pectinesterase), belonged to the significantly enriched GO terms as previously identified (Fig. 6, B). They were downregulated along with their corresponding lncRNAs in 366–2 S. Therefore, the MSTRG.9997-BraA04g004630.3 C, MSTRG.5212-BraA02g040020.3 C, and MSTRG.13,532-BraA05g030320.3 C were considered the major lncRNA-mRNA modules associated with pollen fertility/sterility.

Fig. 6
figure 6

Identification and analysis of differential lncRNAs and target mRNAs. (A) A total of 15 lncRNA-mRNA modules were identified. (B) FPKM expression levels and corresponding relationships of lncRNAs and mRNAs in the identified module. * indicates that mRNAs belong to significantly enriched GO terms

Validation of lncRNA and mRNA expression by qPCR

To verify the accuracy of the lncRNA-seq and mRNA-seq results, the relative expression levels of nine lncRNA-mRNA modules were investigated by qPCR. The expression patterns of most selected lncRNAs and genes analyzed by qPCR were similar to the expression levels that were computed using sequencing data (Fig. 7). Although the magnitude of changes in the expression level of the selected lncRNAs and mRNAs varies, the expression trends were consistent with expectations in sequencing, indicating that the sequencing data were reliable.

Fig. 7
figure 7

Expression validation of mRNAs and lncRNAs by qPCR. The relative expression of E-366–2 F was used as the control in qPCR analysis

Discussion

Abnormalities of the tapetum resulted in the male sterility of 366–2 S

From the formation of sporogenous cells, pollen mother cells are produced by mitosis; then, they form tetrads by meiosis. After the tetrads are dissociated, microspores are released and gradually develop into pollen grains [20]. The abnormality of these processes will lead to pollen abortion or a total absence of pollens [21]. In this study, Chinese cabbage mutant 366–2 S indicated abnormal anthers and pollens, which are typical characteristics of male sterility. Compared with the wild type 366–2 F, tetrads of 366–2 S could not be dissociated or they could not release microspores, which were directly responsible for male sterility. The degradation of tetrads depends on the callase enzyme (β-1,3-glucanase), which typically hydrolyzes the callose wall encompassing four microspores in a tetrad to promote the release of microspores [22]. Many studies have shown that β-1,3-glucanase is secreted by the tapetum cells [23,24,25,26]. Some abnormalities in the tapetum cells may affect the secretion or the function of β-1,3-glucanase and thus lead to the occurrence of male sterility. Cytological observations of anther development also showed that the tapetum in 366–2 S remains intact for a long time and presents delayed degradation. Therefore, these results suggested that some abnormalities of the tapetum delay self-degradation and subsequently affect the action of β-1,3-glucanase on dissociating the tetrads and releasing the microspores, which brings about male sterility of 366–2 S.

Differentially expressed lncRNAs between 366 and 2 F and 366–2 S associated with male sterility

Several lncRNAs have been identified and cloned, and they have been proven to be involved in anther development in different plants, such as maize (Zm401) [10,11,12], rice (LDMAR [13], PMS1T [14, 15], and eTM160 [27]), Brassica campestris (BcMF11) [16, 17], and Brassica rapa (eTM160) [18]. Moreover, it was found that these lncRNAs play a dominant role in a specific stage, such as the microspore mother cell stage, tetrad stage, mononuclear stage, and mature pollen stage, corresponding to the male sterile characteristics of their respective sterile lines. Thus, the use of different male sterile materials is conducive to the identification of new lncRNAs with male sterility. Huang et al. (2018) used the Chinese cabbage genic male sterility (GMS) A/B line (‘Bcajh97-01 A/B’) for lncRNA sequencing at five pollen developmental stages and identified 2,320 and 1,693 lncRNAs specifically expressed in fertile or sterile floral buds, respectively. Through the target gene prediction of lncRNAs, many pollen development-related genes were matched, most of which had been identified as male sterility genes. Thus, these results preliminarily provided some lncRNA databases associated with male sterility in Chinese cabbage. In this study, the Chinese cabbage fertile/sterile lines (366–2 F/S) were used for analyzing the roles of lncRNAs in pollen development. A total of 385 DE-lncRNAs were identified in 366–2 S compared with 366–2 F. The male sterility of 366–2 S primarily resulted from the delayed degradation of the tapetum and the non-separation of tetrads. Conversely, in the sterile line Bcajh97-01 A, disruptions occurred during the meiosis of the pollen mother cells, leading to the inability of the anther to form tetrads [28]. The use of different male sterile materials for identifying lncRNAs may provide insights into distinct developmental processes of pollen. In addition, compared to previously validated lncRNAs, a homology comparison revealed that none of them were identified both in 366–2 F and 366–2 S, which might be due to the different sampling stages or methods or due to the low expression level. Therefore, the DE-lncRNAs identified in this study were potentially different and would enrich lncRNA datasets associated with male sterility in Chinese cabbage.

A large number of downregulated mRNAs related to cell wall degradation potentially delayed the degradation of the tapetum

In this study, a total of 1,176 DE-mRNAs were screened by comparative transcriptome analysis of anthers between 366 and 2 S and 366–2 F. GO enrichment analysis showed that carbohydrate metabolic and cell wall-related processes, such as GO: 0071555 and GO: 0005975, are significantly enriched. The mRNAs involved in these processes encode pectinesterases, polygalacturonases, galactosidase, and glucosidase, all of which are involved in cell wall degradation. Most of these mRNAs were downregulated in 366–2 S, potentially inhibiting cell wall degradation during certain developmental processes. Coincidentally, delayed degradation of the tapetum was observed in the cell sections of 366–2 S. The cytological features of tapetal degradation primarily include the shrinkage of the cell, the rupture of the liquid bubble, the cavity of the cell, and the cracking of the cell wall [29]. Tapetal cells encircle the innermost anther wall layer (named tapetum), which surrounds the male reproductive cells, and undergo programmed cell death-triggered degradation after the meiosis of pollen mother cells [30, 31]. Tapetal degradation is essential for healthy pollen development to provide many molecules including nutrients, proteins, lipids, and polysaccharides; as a result, premature or delayed tapetal degradation causes male sterility [32]. Therefore, we concluded that these downregulated mRNAs related to cell wall degradation potentially hinder tapetal degradation, which then causes male sterility in 366–2 S.

Candidate lncRNA-mRNA regulated the pollen development in Chinese cabbage

LncRNAs target the mRNAs as lncRNA-mRNA modules to play a role in multiple processes of growth and development. For example, lncRNA HID1, which is associated with photomorphogenesis, can assemble into large nuclear protein-RNA complex(es) and bind to Phytochrome Interacting Factor 3 (PIF3) to inhibit its transcription, thereby regulating the development of Arabidopsis hypocotyl [33]. In this study, three major lncRNA-mRNA modules were identified, including one module (MSTRG.9997-BraA04g004630.3 C) associated with callose degradation and two modules (MSTRG.5212-BraA02g040020.3 C and MSTRG.13,532-BraA05g030320.3 C) associated with cell wall degradation, which might contribute to the male sterility of 366–2 S. Indeed, BraA04g004630.3 C encodes a β-1,3-glucanase, which can specifically catalyze the hydrolysis of callose wall during pollen development. In rice, the silencing of Osg1, a β-1,3-glucanase gene, resulted in the disrupted callose degradation in the tetrads and delayed release of young microspores and caused male sterility [34]; the scenario is similar to the results of our study. In 366–2 S, MSTRG.9997 was downregulated and then targeted BraA04g004630.3 C to bring about a low level of β-1,3-glucanase, probably causing non-separation of tetrads and non-production of microspores. BraA02g040020.3 C and BraA05g030320.3 C encode pectinesterase, by which pectin can be demethylated and then cleaved by endo-polygalacturonases, resulting in cell wall loosening. Many studies found that pectinesterase plays an important role in pollen and pollen tube development and its functional deficiency may lead to male sterility [35,36,37]. Thus, it seems that these two lncRNA-mRNA modules might not be causal factors but downstream results along with the male sterility occurrence in 366–2 S. However, cytological observations showed that tapetum in 366–2 S presents delayed degradation, which must involve the process of cell wall degradation. Probably, certain pectinesterases (such as BraA02g040020.3 C and BraA05g030320.3 C) also played a role in the degradation of tapetum, and these two lncRNA-mRNA modules had some key downstream links. Therefore, these three candidate lncRNA-mRNA modules might potentially be involved in regulating the pollen development in Chinese cabbage.

Conclusions

This study revealed that the abnormal pollen development in the Chinese cabbage sterile line 366–2 S might start at the tetrad stage, including the delayed degradation of the tapetum and the non-separation of tetrads. Through lncRNA sequencing, 385 DE-lncRNAs were identified between 366 and 2 F and 366–2 S, which are probably related to male sterility of 366–2 S. To understand exactly the alterations in mRNA expression and identify the potential target genes of lncRNAs, comparative transcriptome analyses were performed at two pollen developmental stages. A total of 1,176 DE-mRNAs were screened. GO enrichment analysis showed that carbohydrate metabolic and cell wall-related processes are significantly enriched, involving a large number of downregulated mRNAs related to cell wall degradation, which potentially play a role in delaying the tapetal degradation. Furthermore, a total of 15 lncRNA-mRNA modules were identified. Among them, three major lncRNA-mRNA modules were considered to be closely related to the male sterility of 366–2 S. Among them, MSTRG.9997-BraA04g004630.3 C (β-1,3-glucanase) was probably associated with callose degradation and separation of tetrads; MSTRG.5212-BraA02g040020.3 C (pectinesterase) and MSTRG.13,532-BraA05g030320.3 C (pectinesterase) were probably associated with cell wall degradation of tapetums, which is typically involved in regulating pollen development. Overall, this study provides some important clues and evidence that lncRNAs participate in pollen development, and it also lays the basic genomic foundation for further investigations.

Methods

Plant materials and cultivation

The genic male sterile mutant line 366–2 S originated as a natural mutant derived from the Chinese cabbage breeding material 366-2 population, sourced from the National Vegetable Germplasm Resource Intermediate Bank, Institute of Vegetable and Flower, Chinese Academy of Agricultural Sciences. The homozygous fertile line was identified as the wild type (366–2 F). Upon crossing 366–2 S with 366–2 F, the resulting F1 population presented a fertile type, while the F2 population exhibited a segregation ratio of fertile to sterile types at 3:1 due to a recessive single gene mutation present in 366–2 S. Additionally, in contrast with the tapetum of the fertile 366–2 F, the tapetum observed in the sterile 366–2 S consistently persisted without undergoing degradation, presenting a notable deviation from typical male sterile mutants. The cultivation of all plant materials took place in the standardized experimental field situated at the Henan Modern Agriculture Research and Development Base in Yuanyang, Henan Province, China.

Cytological observation

For identifying the period of pollen abortion, the flower buds of 366–2 S and 366–2 F were fixed by using FAA solution (50% ethanol, 5% glacial acetic acid, and 10% formalin) and then dehydrated in a graded ethanol series. The prepared anthers were embedded in paraffin and stained with hematoxylin and eosin and observed under an optical microscope (Nikon ECLIPSE 80i; Nikon, Japan). The buds were subsequently divided into five groups considering the length of the buds (BUD1 < 1 mm, BUD2: 1  1.5 mm, BUD3: 1.5  2.5 mm, BUD4: 2.5  3.5 mm, and BUD5 > 3.5 mm).

Sample preparation and transcriptome sequencing

A single biological replicate was prepared for lncRNA sequencing, aiming for accuracy and sequencing efficacy. Anthers from the 366–2 S line were combined from multiple buds of different lengths (ranging between 1 and 5 mm) before flowering, ensuring a diverse representation. Similar types of anthers from the 366–2 F line were utilized as the control group. The extraction of total RNA was conducted using the TRIzol reagent method (Invitrogen, Carlsbad, CA, USA). The concentration and quality of RNA were evaluated by using a NanoDrop 2000 system (Thermo Fisher Scientific, Wilmington, DE) and an Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA). Ribosomal RNA (rRNA) and poly-A enriched RNA (most coding sequences) were removed from the total RNA. The remaining RNA was randomly fragmented into short fragments for constructing libraries of sequencing. First-strand cDNA was synthesized by using random hexamers and M-MuLV reverse transcriptase; second-strand cDNA was subsequently synthesized by using DNA polymerase I and RNase H; after using AMPure XP beads to purify double-stranded cDNA, they were subjected to end repair, addition of poly-A tail, ligation of sequencing linker, and fragment size selection. Finally, two cDNA libraries were subjected to PCR experiments and sequenced in Illumina HiSeq 2500 platform (150 bp paired-end sequencing; 15 Gb data size).

Anthers of parallel buds at the early stage (E; about 1  1.5 mm) and parallel buds at the late stage (L; about 3.5  5 mm) were selected for mRNA sequencing. In total, 12 samples including four treatments (E-366–2 S, E-366–2 F, L-366–2 S, and L-366–2 F) with three replicates were used. The mRNA (coding sequences) was enriched with a poly-A tail from total RNA by using magnetic beads with Oligo dT. Subsequently, the sequencing libraries of all samples were constructed. Finally, 12 cDNA libraries were sequenced in the Illumina HiSeq 2500 platform (150 bp paired-end sequencing; 6 Gb data size).

LncRNA identification and target gene prediction

FastQC software [38] was used to confirm the quality (Q20, Q30) of lncRNA-seq data. Clean reads were obtained by removing reads with an adapter and containing poly-N (an ambiguous sequence content exceeding 10%) and low-quality reads (Q20 < 80%). After aligning the clean reads to the B. rapa reference genome (V3.0) [39], using HISAT2 software [40], all mapped reads were assembled into transcripts by using the StringTie software [41]. Considering the incomplete transcript assembly due to gaps in coverage, artificial reads were constructed and used for assembly based on reference transcripts. The transcripts with lengths less than 200 nt or possessing protein-coding ability were filtered out. Furthermore, the protein-coding potential of the remaining transcripts was predicted by using the CPC software [42]. The transcripts without protein-coding ability were identified as lncRNAs.

Given the small sample size of sequencing, the trans-regulation of lncRNAs was not effectively analyzed. Thus, the target gene prediction of lncRNAs was carried out by co-located and base pairing analysis. An lncRNA located within 10 kb upstream or downstream of a coding gene was considered as co-location genes, and these genes were predicted as the lncRNA’s co-located target genes. The co-located target genes were identified by custom Perl scripts. In addition, the RNAplex package [43] was used to predict the complementary binding between lncRNA and mRNA (coding gene) to obtain the base pairing target genes.

Differential expression analysis of lncRNAs

Cufflinks software [43] was used to calculate the fragments per kilobase per million (FPKM) value of exons as the expression level of each lncRNA. Differential expression analysis of lncRNAs was performed between 366 and 2 S and 366–2 F (as control) by using the R package (DESeq) [44]. The differentially expressed lncRNAs (DE-lncRNAs) were screened out with a |fold change| ≥ 2 and p < 0.01. Further, their corresponding predicted target genes were acquired.

Identification and differential analysis of mRNA expression profiles

For mRNA-seq data, clean reads from all treatments (E-366–2 S, E-366–2 F, L-366–2 S, and L-366–2 F) were obtained and aligned to the B. rapa reference genome (v3.0). Based on the known genome annotation information, the mapped reads were assembled into annotated mRNAs and unannotated mRNAs (potentially representing novel genes). The BLAST software [45] was used to obtain the annotation information for these new genes through sequence alignment with NR (NCBI non-redundant protein sequences) and Gene Ontology (GO) databases. FPKM values were computed for all assembled mRNAs within each sample. Correlation analysis across three replicates for each treatment was performed to assess repeatability. The average of these three replicates was used as the expression profile for each treatment. Differential expression analysis of mRNA profiles was performed within two distinct groups: the early-stage group (E-366–2 S vs. E-366–2 F) and the late-stage group (L-366–2 S vs. L-366–2 F). Differentially expressed mRNAs (DE-mRNAs) were identified for each group using the criteria of |fold change| ≥ 2 and p < 0.01. Moreover, overlapped DE-mRNAs between the two groups were determined using a Venn diagram approach to explore commonalities and distinctions in gene expression patterns.

Go enrichment analysis

Overlapping DE-mRNAs underwent GO enrichment analysis using Blast2GO software [46]. A hypergeometric test served to pinpoint significantly enriched GO terms (biological process, BP). A significance level of 0.01 was selected for the enrichment analysis (p < 0.01). In addition, the analysis employed a rich factor to quantify the degree of enrichment. The rich factor gauges the ratio between the proportion of GO terms annotated within the overlapped DE-mRNAs and those present within the annotated mRNAs.

Screening of lncRNA-mRNA modules associated with male sterility

Some target genes of DE-lncRNAs associated with male sterility were predicted as above. Moreover, the overlapped DE-mRNAs associated with male sterility were identified as above. Through Venn analysis, the consistent mRNAs were selected between these predicted target genes and these overlapped DE-mRNAs. The consistent DE-mRNAs, particularly those associated with significantly enriched GO terms, were potentially paired with corresponding DE-lncRNAs, forming potential lncRNA-mRNA modules that could work together in regulating anther fertility.

Quantitative real-time PCR (qPCR) analysis

All 12 anther samples from the two stages were re-prepared for a quantitative real-time PCR (qPCR) analysis. Total RNAs from all samples were extracted and reverse transcribed into cDNA using the TransScript one-step gDNA removal and cDNA Synthesis Kit (Trans, Beijing, China). A subset of potential lncRNA-mRNA modules was selected for verification via qPCR. The qPCR experiments were performed using the Roche Light Cycler 480-II system (Roche Applied Sciences, Beijing, China) following specific thermal profiles, with the BrGAPDH gene serving as an endogenous control. Each plate included three technical replicates. Specific primers for the selected lncRNAs and mRNAs were designed online via the NCBI website (https://www.ncbi.nlm.nih.gov/) (Supplementary Table S1).

For both lncRNAs and mRNAs, qPCR was performed using TB Green Premix Ex Taq II (Takara, Beijing, China) within a final qPCR reaction mixture containing 5.0 µL of 2×TB Green Premix Ex Taq, 0.8 µL of forward primers, 0.8 µL of reverse primers, 1.2 µL of cDNA, and 2.2 µL of ddH2O. Expression levels were determined using the 2−ΔΔCT method [47], with the BrGAPDH gene as the reference gene and internal control.

Data availability

The original contributions presented in the study are publicly available. This data can be found here: China National GeneBank DataBase (CNGBdb) under accession number CNP0004830.

References

  1. Sun X, Li X, Lu Y, Wang S, Zhang X, Zhang K, et al. Construction of a high-density mutant population of Chinese cabbage facilitates the genetic dissection of agronomic traits. Mol Plant. 2022;15(5):913–24.

    Article  CAS  PubMed  Google Scholar 

  2. Kamiya Y, Shiraki S, Fujiwara K, Akter M, Akter A, Fujimoto R, et al. The role of epigenetic modifications in the transcriptional regulation of long noncoding RNAs in Brassica rapa vegetables. Acta Hortic. 2023;1362:65–70.

    Article  Google Scholar 

  3. Wilson ZA, Zhang DB. From Arabidopsis to rice: pathways in pollen development. J Exp Bot. 2009;60(5):1479–92.

    Article  CAS  PubMed  Google Scholar 

  4. Ge X, Chang F, Ma H. Signaling and transcriptional control of reproductive development in Arabidopsis. Curr Biol. 2010;20(22):R988–97.

    Article  CAS  PubMed  Google Scholar 

  5. Gómez JF, Talle B, Wilson ZA. Anther and pollen development: a conserved developmental pathway. J Integr Plant Biol. 2015;57(11):876–91.

    Article  PubMed  PubMed Central  Google Scholar 

  6. Kurusu T, Koyano T, Hanamata S, Kubo T, Noguchi Y, Yagi C, et al. OsATG7 is required for autophagy-dependent lipid metabolism in rice postmeiotic anther development. Autophagy. 2014;10(5):878–88.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  7. Mishra A, Bohra A. Non-coding RNAs and plant male sterility: current knowledge and future prospects. Plant Cell Rep. 2018;37(2):177–91.

    Article  CAS  PubMed  Google Scholar 

  8. Kopp F, Mendell JT. Functional classification and experimental dissection of long noncoding RNAs. Cell. 2018;172(3):393–407.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  9. Wierzbicki AT, Blevins T, Swiezewski S. Long noncoding RNAs in plants. Annu Rev Plant Biol. 2021;72:245–71.

    Article  CAS  PubMed  Google Scholar 

  10. Dai X, Yu J, Ma J, Ao G, Zhao Q. Overexpression of Zm401, an mRNA-like RNA, has distinct effects on pollen development in maize. Plant Growth Regul. 2007;52:229–39.

    Article  CAS  Google Scholar 

  11. Ma J, Yan B, Qu Y, Qin F, Yang Y, Hao X, et al. Zm401, a short-open reading-frame mRNA or noncoding RNA, is essential for tapetum and microspore development and can regulate the floret formation in maize. J Cell Biochem. 2008;105(1):136–46.

    Article  CAS  PubMed  Google Scholar 

  12. Wang D, Li C, Zhao Q, Zhao L, Wang M, Zhu D, et al. Zm401p10, encoded by an anther-specific gene with short open reading frames, is essential for tapetum degeneration and anther development in maize. Funct Plant Biol. 2009;36(1):73–85.

    Article  PubMed  Google Scholar 

  13. Ding J, Lu Q, Ouyang Y, Mao H, Zhang P, Yao J, et al. A long noncoding RNA regulates photoperiod-sensitive male sterility, an essential component of hybrid rice. Proc Natl Acad Sci U S A. 2012;109(7):2654–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. Zhou H, Liu Q, Li J, Jiang D, Zhou L, Wu P, et al. Photoperiod- and thermo-sensitive genic male sterility in rice are caused by a point mutation in a novel noncoding RNA that produces a small RNA. Cell Res. 2012;22(4):649–60.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  15. Fan Y, Yang J, Mathioni SM, Yu J, Shen J, Yang X, et al. PMS1T, producing phased small-interfering RNAs, regulates photoperiod-sensitive male sterility in rice. Proc Natl Acad Sci U S A. 2016;113(52):15144–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  16. Song JH, Cao JS, Yu XL, Xiang X. BcMF11, a putative pollen-specific non-coding RNA from Brassica campestris ssp. chinensis. J Plant Physiol. 2007;164(8):1097–100.

    Article  CAS  PubMed  Google Scholar 

  17. Song JH, Cao JS, Wang CG. BcMF11, a novel non-coding RNA gene from Brassica campestris, is required for pollen development and male fertility. Plant Cell Rep. 2013;32(1):21–30.

    Article  CAS  PubMed  Google Scholar 

  18. Huang L, Dong H, Zhou D, Li M, Liu Y, Zhang F, et al. Systematic identification of long non-coding RNAs during pollen development and fertilization in Brassica rapa. Plant J. 2018;96(1):203–22.

    Article  CAS  PubMed  Google Scholar 

  19. Lohani N, Golicz AA, Allu AD, Bhalla PL, Singh MB. Genome-wide analysis reveals the crucial role of lncRNAs in regulating the expression of genes controlling pollen development. Plant Cell Rep. 2023;42(2):337–54.

    CAS  PubMed  Google Scholar 

  20. Liu L, Wang T. Male gametophyte development in flowering plants: a story of quarantine and sacrifice. J Plant Physiol. 2021;258–259:153365.

    Article  PubMed  Google Scholar 

  21. Hu Z, Zhang M, Wen Q, Wei J, Yi H, Deng X, et al. Abnormal Microspore Development leads to Pollen Abortion in a Seedless Mutant of ‘Ougan’ Mandarin (Citrus Suavissima Hort. Ex Tanaka). J Amer Soc Hort Sci. 2007;132(6):777–82.

    Article  Google Scholar 

  22. Wang B, Andargie M, Fang R. The function and biosynthesis of callose in high plants. Heliyon. 2022;8(4):e09248.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  23. Wang Y, Li X, Fan B, Zhu C, Chen Z. Regulation and function of defense-related callose deposition in plants. Int J Mol Sci. 2021;22(5):2393.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Stieglitz H, Stern H. Regulation of beta-1,3-glucanase activity in developing anthers of Lilium. Dev Biol. 1973;34(1):169–73.

    Article  CAS  PubMed  Google Scholar 

  25. Stieglitz H. Role of beta-1,3-glucanase in postmeiotic microspore release. Dev Biol. 1977;57(1):87–97.

    Article  CAS  PubMed  Google Scholar 

  26. Scott R, Hodge R, Paul W, Draper J. The molecular biology of anther differentiation. Plant Sci. 1991;80:167–91.

    Article  CAS  Google Scholar 

  27. Ding J, Shen J, Mao H, Xie W, Li X, Zhang Q. RNA-directed DNA methylation is involved in regulating photoperiod-sensitive male sterility in rice. Mol Plant. 2012;5(6):1210–6.

    Article  CAS  PubMed  Google Scholar 

  28. Huang L, Cao J, Ye W, Liu T, Jiang L, Ye Y. Transcriptional differences between the male-sterile mutant bcms and wild-type Brassica campestris ssp. chinensis reveal genes related to pollen development. Plant Biol (Stuttg). 2008;10(3):342–55.

    Article  CAS  PubMed  Google Scholar 

  29. Xie F, Vahldick H, Lin Z, Nowack MK. Killing me softly - programmed cell death in plant reproduction from sporogenesis to fertilization. Curr Opin Plant Biol. 2022;69:102271.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Solís MT, Chakrabarti N, Corredor E, Cortés-Eslava J, Rodríguez-Serrano M, Biggiogera M, et al. Epigenetic changes accompany developmental programmed cell death in tapetum cells. Plant Cell Physiol. 2014;55(1):16–29.

    Article  PubMed  Google Scholar 

  31. Ono S, Liu H, Tsuda K, Fukai E, Tanaka K, Sasaki T, et al. EAT1 transcription factor, a non-cell-autonomous regulator of pollen production, activates meiotic small RNA biogenesis in rice anther tapetum. PLoS Genet. 2018;14(2):e1007238.

    Article  PubMed  PubMed Central  Google Scholar 

  32. Wei S, Ma L. Comprehensive insight into tapetum-mediated Pollen Development in Arabidopsis thaliana. Cells. 2023;12(2):247.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  33. Wang Y, Fan X, Lin F, He G, Terzaghi W, Zhu D, et al. Arabidopsis noncoding RNA mediates control of photomorphogenesis by red light. Proc Natl Acad Sci U S A. 2014;111(28):10359–64.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  34. Wang Y, Zhang D, Huang L, Zhang Z, Gao J, et al. Research Progress of Pectin Methylesterase and its inhibitors. Curr Protein Pept Sci. 2022;23(10):684–96.

    Article  PubMed  Google Scholar 

  35. Francis KE, Lam SY, Copenhaver GP. Separation of Arabidopsis pollen tetrads is regulated by QUARTET1, a pectin methylesterase gene. Plant Physiol. 2006;142(3):1004–13.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  36. Niu F, Bu Y, Yang X, Wu Y, He M, Zhang L, et al. Rfd1, a restorer to the Aegilops juvenalis cytoplasm, functions in fertility restoration of wheat cytoplasmic male sterility. J Exp Bot. 2023;74(5):1432–47.

    Article  CAS  PubMed  Google Scholar 

  37. Cock PJ, Fields CJ, Goto N, Heuer ML, Rice PM. The Sanger FASTQ file format for sequences with quality scores, and the Solexa/Illumina FASTQ variants. Nucleic Acids Res. 2010;38(6):1767–71.

    Article  CAS  PubMed  Google Scholar 

  38. Zhang L, Cai X, Wu J, Liu M, Grob S, Cheng F, et al. Improved Brassica rapa reference genome by single-molecule sequencing and chromosome conformation capture technologies. Hortic Res. 2018;5:50.

    Article  PubMed  PubMed Central  Google Scholar 

  39. Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. 2015;12(4):357–60.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Pertea M, Pertea GM, Antonescu CM, Chang TC, Mendell JT, Salzberg SL. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat Biotechnol. 2015;33(3):290–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  41. Kong L, Zhang Y, Ye ZQ, Liu XQ, Zhao SQ, Wei L et al. CPC: assess the protein-coding potential of transcripts using sequence features and support vector machine. Nucleic Acids Res. 2007:W345–9.

  42. Tafer H, Hofacker IL. RNAplex: a fast tool for RNA-RNA interaction search. Bioinformatics. 2008;24(22):2657–63.

    Article  CAS  PubMed  Google Scholar 

  43. Trapnell C, Roberts A, Goff L, Pertea G, Kim D, Kelley DR, et al. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nat Protoc. 2012;7(3):562–78.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  44. Love MI, Huber W, Anders S. Moderated estimation of Fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550.

    Article  PubMed  PubMed Central  Google Scholar 

  45. Camacho C, Coulouris G, Avagyan V, Ma N, Papadopoulos J, Bealer K, et al. BLAST+: architecture and applications. BMC Bioinformatics. 2009;10:421.

    Article  PubMed  PubMed Central  Google Scholar 

  46. Götz S, García-Gómez JM, Terol J, Williams TD, Nagaraj SH, Nueda MJ, Robles M, Talón M, Dopazo J, Conesa A. High-throughput functional annotation and data mining with the Blast2GO suite. Nucleic Acids Res. 2008:3420–35.

  47. Rao X, Huang X, Zhou Z, Lin X. An improvement of the 2ˆ(-delta delta CT) method for quantitative real-time polymerase chain reaction data analysis. Biostat Bioinforma Biomath. 2013;3(3):71–85.

    PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgements

We are deeply indebted to the entire research team for their support during this research work.

Funding

This research work was supported by China Agricultural Research System (CARS-23-G15), Innovation Team of Henan Academy of Agricultural Sciences (2023TD06), Joint Research on Agricultural Variety Improvement of Henan Province (2022010504), Key Research & Development Project of Henan Province (221111110100).

Author information

Authors and Affiliations

Authors

Contributions

YXY and XWZ conceived the study. WF and XCW designed the experiments. ZYW, SJY and HNS performed the experiments and analyzed the data. YYZ, WJZ and WWC participated in the qRT-PCR and phenotype identification. XCW and XQW drafted the manuscript. BMT and UKN revised and finalized the manuscript. All authors have read and approved the publication of the manuscript.

Corresponding authors

Correspondence to Fang Wei, Yuxiang Yuan or Xiaowei Zhang.

Ethics declarations

Ethics approval

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Electronic supplementary material

Below is the link to the electronic supplementary material.

Supplementary Material 1: Fig. S1 Expression correlation heat map of samples.

12870_2024_5003_MOESM2_ESM.xlsx

Supplementary Material 2: Table S1 The information of the primers used in q-PCR. Table S2 The analysis of lncRNA sequencing. Table S3 The mRNA sequencing analysis. Table S4 GO enrichment analysis of differentially expressed mRNAs based on biological process.

Supplementary Material 3

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Wei, X., Wang, X., Zhao, Y. et al. Transcriptome analysis reveals the potential lncRNA-mRNA modules involved in genetic male sterility and fertility of Chinese cabbage (brassica rapa L. ssp. pekinensis). BMC Plant Biol 24, 289 (2024). https://doi.org/10.1186/s12870-024-05003-w

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12870-024-05003-w

Keywords