Morphological and cytological comparison of anthers and pollens in CHA treated materials and control floral buds
It was revealed that the petals of control flowers were significantly wider than those in CHA treated materials, and the anthers and filaments in flowers with no CHA treatment were obviously longer than those in CHA treated materials (Fig. 1a and b). The anthers in sterile buds gradually shrivel from the middle stage and didn’t split and release pollen at last, while the pistil with CHA treatment was normal.
To further understand the ultrastructure difference of anthers between CHA treated materials and control, SEM and TEM observation for pollens were conducted. For SEM analysis, no distinct differences were observed in pollens before stage 7 between CHA treated materials and control. The pollen grains that treated with CHA were appeared loss of inclusion, and presented shriveled compared to control at stage 8 (Fig. 1e and i). In addition, the tryphine of control was fulfilled with some substances, while it does not exist on CHA (Fig. 1u and v).
For TEM analysis, no obvious difference was observed before stage 8 between the control (Fig. 2a-c) and CHA treated materials (Fig. 2e-g), the abnormal development was observed in most of CHA pollen grains from stage 9, including plasmolysis phenomenon, deletion of organelles and cytoskeleton, and totally empty pollens at last (Fig. 2m-p). What’s more, the tryphine was not deposited completely in the pollen wall compared to normal pollens (Fig. 2q and r). In addition, no obvious difference was observed in the tapetum before stage 7, but the tapetosomes and elaioplasts disassembled earlier in tapetal cells of CHA treated materials than those in control until stage 8 (Fig. 2t and x). The thickness of CHA tapetal cells were 10.3 ± 0.21 μm (Fig. 2u and v), while only 1.02 ± 0.12 μm in control (Fig. 2y and z), which implied that the tapetal cells were undergoing normal PCD process in control, while this process was delayed with the tapetal cells stacked together in CHA treated materials.
Comparative proteomic analysis between the anthers that treated with CHA and control
In order to detect differentially expressed proteins (DEPs) at different stages, two-dimensional electrophoresis (2-DE) with three biological triplicates were conducted for total protein of anthers from CHA treated materials and control. As a result, 1 024, 723, 1 207 and 1 120 protein spots at SA (stage 6–7), MA (stage 8–9), LA (stage 10–11) and LA2 (stage 12–13) of control, and 847, 409, 940 and 771 protein spots at SA, MA, LA and LA2 stage of CHA treated materials were detected, respectively (Fig. 3). Compared to control, 87, 25, 74 and 60 increased expression spots were identified that changed in abundance over two-fold and exceeded the least significant difference (p < 0.05) at SA, MA, LA and LA2 stages respectively, and the decreased spots were 130, 220, 329 and 366, respectively.
One hundred nineteen DEPs that showed difference expression more than two stages were chosen for MALDI-TOF-MS-MS analysis and 101 DEPs were successfully identified (Additional file 1). These DEPs were grouped into 16 categories, including protein metabolism (29.7%), energy metabolism (14.9%), amino acid metabolism (8.9%), defense (6.9%), carbohydrate metabolism (5.9%), and so on. The enlarged area of 2-DE gel of some representative DEPs were showed in Fig. 3i, including protein tapetum 1 (spot 204; ATA1, AT3G42960), which was involved in amino acid metabolism, protein disulfide isomerase precursor-like (spot 8806; AT5G60640), cysteine protease (spot 8234; AT5G60360) and carboxypeptidase (spot 5604; AT3G45010) which were related to protein metabolism, chalcone synthase 3 protein (spot 2417; AT1G02050) which was involved in metabolism of terpenoids and polyketides, auxin responsive-like protein (spot 3712; AT5G13370) which was related to signal transduction, and one cytoskeleton protein, Actin-12 (spot 6407; AT3G46520). All of these proteins were down-regulated at CHA treated anther.
To further facilitate the biological interpretation of the identified DEPs, the hierarchical clustering analysis of the quantitative DEPs were conducted and four major clusters were recognized (Fig. 4a), which revealed that the DEPs between fertile and sterile anther at different development stages showed great changes in their expression patterns. Further research indicated that Cluster A accounted for 26 proteins and most of them down-regulated at MA stage, with energy metabolism and protein metabolism dominant (Fig. 4b), such as transketolase-like protein (spot 5702; AT3G60750) and chloroplast HSP70–1 (spot 8811; AT5G49910). Cluster B included 18 proteins with the down-regulation at MA and LA stage, with protein metabolism and amino acid metabolism dominant (Fig. 4b), for example, serine carboxypeptidase-like 48 (spot 5604; AT3G45010), s-adenosylmethionine synthetase (spot 627; AT2G36880), and so on. Cluster C, the largest cluster, almost decreased expressed at LA2 stage, represented 32 proteins with protein metabolism and energy metabolism proteins as the dominant classification (Fig. 4b), for instance, HSP70 (spot 6903; AT1G79920) and soluble inorganic pyrophosphatase (spot 5106; AT1G01050). Cluster D, accounted for 25 proteins that almost showed significant down-regulation at YB stage, in which the protein metabolism and amino acid metabolism were the majority classification (Fig. 4b), such as protein disulfide-isomerase A6 (spot 2314; AT2G47470) and glutathione transferase (spot 1101; AT2G30860).
Comparative transcriptomic analysis between the anthers of CHA treated materials and control
Comparative transcriptomics was conducted to identify genes associated with male sterility induced by SX-1, three biological replicates of RNA-seq of different development stages was profiled and good repetition among three biological replicates in both CHA treated materials and control were obtained, which produced 139.7, 107.3, 151.1, 114.7 million raw reads for the control and 194.2, 147.6, 251.6 and 149.4 million raw reads for the CHA treated materials at YB stage (stage 1–5), SA, MA and LA stage, respectively. These raw data were filtered to get clean reads respectively, and then average 85.59% clean reads perfectly matched to the reference genome [28] (Additional file 2).
The differential gene expression analysis was conducted, and a total of 998, 2 194, 2 428 and 10 627 up-regulated differentially expressed genes (DEGs) and 1 177, 2 488, 827 and 9 745 down-regulated DEGs (P < 0.05) were identified at YB, SA, MA and LA stage, respectively (Fig. 5a and b). The number of up- and down-regulated DEGs shared at all stages was only 60 and 9 respectively, while large number of DEGs in these adjoining stages were different. The hierarchical cluster analysis of all DEGs was performed to show the global gene expression pattern at different stages (Fig. 5c). As a result, CHA_YB and CHA_SA, YB and SA, MA and LA stage were located in the same cluster respectively, while CHA_MA and CHA_LA were far away from MA and LA respectively. Moreover, the differentially expression levels of several selected DEGs for further discussion were showed in Additional file 3. For example, two copies of CSR1, the target gene of SX-1, were down-regulated at YB stage.
As TFs usually play an important role in regulating at early stage [29], 68 up- and 144 down-regulated TFs were identified at YB stage (Additional file 4), and the top three TF families (ranked by DEG numbers) were “ERF” (41 DEGs), “MYB” (34 DEGs) and “NAC” (23 DEGs). Previous study had proposed a putative regulatory network model for Arabidopsis anther development, and DYT1 (a putative bHLH TF) played the key role in the model, and NAC025, MYB80 and WRKY33 were the downstream genes of DYT1 [30, 31]. DYT1 didn’t show differentially expression at this stage (Additional file 3), while two copies of NAC025, four copies of MYB80, and six copies of WRKY33 were significantly down-regulated at YB stage (Additional file 4). In Arabidopsis, MYB80 (formerly named as MYB103) was restrictedly expressed at the tapetum of developing anthers and trichomes, and its down-regulation leaded to abnormal tapetum and pollen development [32, 33].
The up- and down-regulated DEGs at different stages between control and CHA treated materials were assigned for KEGG pathway analysis (P < 0.05). The results revealed that “flavonoid biosynthesis” and “protein processing in ER” were enriched at YB stage, “Phenylpropanoid biosynthesis” and “Galactose metabolism” were enriched at SA stage (Additional file 5), as several genes in these pathways were significantly differentially expressed, such as Hsp70, CCOAMT, TT4 and UGT72E3. In Arabidopsis, UGT72E1- E3 encode glycosyltransferases that glucosylated phenylpropanoids in vitro [34], and UGT72E3 was predicted as the target of bra-miR9556b-5p in the present study. For MA stage, the down-regulated DEGs were enriched in “Protein export”, “Pyruvate metabolism” and “Fatty acid elongation” (Additional file 5). Among these genes, IPTs and BRI1 were identified. IPTs mediated the rate-limiting step of cytokinin biosynthesis [35], and were the predicted target of sample_miRNA_520 in the present study. At the last stage, many pathways related to amino acid metabolism were identified, including valine, leucine and isoleucine degradation, lysine degradation, alanine, aspartate and glutamate metabolism, etc. (Additional file 5). In the meanwhile, “Fatty acid degradation” and “Steroid biosynthesis” were enriched. Some important genes in these pathways were up- or down-regulated, for example, LACSs, AIM and FK. In Arabidopsis embryogenesis, FK was a sterol C-14 reductase, which was required for organized cell division and expansion [36], and was the predicted target of sample_miRNA_201 in the present study. More importantly, the down-regulated DEGs at these stages shared one same pathway, protein processing in ER, which might imply the male sterility induced by SX-1 largely because of the disordered protein processing.
The down-regulated DEGs between CHA treated materials and control were further enriched by Cytoscape EnrichmentMap. As shown in Fig. 5d, YB, SA, MA and LA stage generated 64, 35, 52 and 251 nodes, respectively. These nodes were classified into different categories. Interestingly, the common terms among the four stages was “Response to stimulus”, which might largely due to the continued effect by SX-1 from early stage. Very importantly, some down-regulated DEGs at different stages were enriched in different categories, which might because of the continuous effect of SX-1 at different anther development stages. For instance, “Cellular component organization” was only enriched at YB and SA stage and “Catabolic process” was only enriched at SA stage. What’s more, “Localization”, “Cellular process”, “Cellular metabolic process” and so on were also clustered.
Differentially expressed (DE) miRNAs between the CHA treated materials and control
A total of 125.9, 79.3, 97.9, 114.4 million raw reads and 86.3, 72.1, 73.7, 51.9 million raw reads were obtained at YB, SA, MA and LA stage in control and CHA treated materials, respectively. The length distribution of the small RNAs was analyzed and the 21–24 nt small RNAs were with the majority (the 24 nt small RNA was the most dominant) (Additional file 6). A total of 115 known miRNAs and 557 novel miRNAs were identified (Additional file 7). The hierarchical cluster analysis of all identified miRNAs was performed to show the different expression pattern at different stages (Fig. 6a). As a result, CHA_SA, CHA_MA and CHA_LA, YB and CHA_YB, SA and MA were located in the same cluster respectively, while LA stage was alone. Among the 672 miRNAs, 9, 14, 22 and 6 DE-miRNAs at YB, SA, MA and LA stage were identified, respectively (Additional file 8). Interestingly, these stages didn’t share even one common DE-miRNA, which suggested that miRNAs regulations induced by SX-1 at different anther development stages were obviously not alike.
Among the total 51 DE-miRNAs, 80, 101, 156 and 57 targets were identified for 7, 12, 20 and 5 DE-miRNAs at YB, SA, MA and LA stage, respectively (Additional file 9). Some miRNAs might be the important regulators, for example, hvu-miR156b, bra-miR5654b and one novel miRNA Samples_miRNA_520 targeted more than 20 genes, bra-miR9556-5p and mtr-miR169i target 12 and 5 genes in the miRNA-mRNA network in B. napus (Fig. 6b). Interestingly, most targets of these DE-miRNAs were TFs. For instance, hvu-miR156b targeted SPL, mtr-miR169i targeted WRKY and ata-miR164c-5p targeted NAC, which suggested the timely miRNAs/TFs regulation modules might be essential for cellular reprogramming early in the response to SX-1. In addition, bra-miR5654b and bra-miR158-5p targeting PPR superfamily proteins and a novel miRNA Samples_miRNA_520 targeting isopentenyl transferases (IPTs) were identified.
For the total 394 DE-miRNA/mRNA modules, 16 modules were inversely related, for example, Samples_miRNA_500, bra-miR5654b, hvu-miR156b and bra-miR5654b at YB, SA, MA and LA stage respectively, which targeted unknown protein, PPR, SPL and PPR respectively, whereas 20 modules were positively related (Additional file 10). Another 358 modules had no significant differential expression.
The correlation analysis of transcript-to-protein between DEPs and DEGs
As B. napus is an allotetraploid plant, its genomes contain a large number of homologous sequences and repeats [37]. The 101 DEPs identified above, were corresponding encoded by 423 genes in total. Subsequently, correlation analysis between DEPs and DEGs was performed among different development stages. According to the expression in CHA treated materials compared to control, DEPs and DEGs can be divided into 4 groups: Group I, 23 (12 up-regulated and 11 down-regulated), 14 (10 up-regulated and 4 down-regulated), 53 (27 up-regulated and 26 down-regulated) DEGs showed the same expression trend with DEPs at SA, MA and LA stage respectively. For example, Hsp70, BIP2 and SHD, which are associated with protein processing in ER, SAM and GLN1.3, proteins related to amino acids biosynthesis, and KASI, a fatty acids biosynthesis protein. Group II, 17 at SA stage, 24 at MA stage and 46 at LA stage, these DEGs showed opposite expression trend with DEPs, such as TLP-3, a thaumatin-like protein 3 and CDC48, cell division cycle 48, which was functioned in ATPase activity and involved in response to cadmium ion; Group III, 369, 366 and 240 DEGs at SA, MA and LA stage respectively were not differentially expressed, for example, GRF10 and FBR12, which were involved in protein phosphorylated amino acid binding; Group IV, 26 DEPs at SA stage, 34 DEPs at MA stage and 36 DEPs at LA stage didn’t show different expression, for instance, CAB1, a chlorophyll A/B binding protein and functioned in chlorophyll binding, and FIB, functioned in structural molecule activity and involved in photoinhibition. All the correlated DEGs and DEPs were displayed in Additional file 11. These results revealed that both independent and parallel correlations existed between the mRNA and protein expression profiles among different stages, which suggested there were a highly complex post-transcriptional regulatory network in the male sterility induced by SX-1.
Validation of selected DEPs, DEGs, DE-miRNA by qRT-PCR and free amino acids, lipids and flavonoids between CHA treated and control materials
To validate the reliability of the RNA-seq and miRNA-seq, 15 DEGs and 8 DE-miRNAs potentially related to CIMS were selected for qRT-PCR (Additional file 12). The expression patterns of 12 DEGs and 6 DE-miRNAs detected by qRT-PCR were basically consistent with the RNA-seq and miRNA-seq, these results indicated that the RNA-seq and miRNA-seq in the present study were reliable. What’s more, among 6 DEPs encoded by 8 DEGs selected above, only the relative expression levels of BnaC04g45570D, BnaA04g21720D and BnaA05g22420D were consisted with the proteomic results, and the low correlation largely agreed with transcript-to-protein analysis mentioned above.
The comparative proteomic and transcriptomic analysis showed the male sterility induced by SX-1 might involve in amino acid (AA), fatty acid (FA) and flavonoid biosynthesis. To confirm these findings, the contents of total free AA, FA and flavonoids were analyzed at YB, SA, MA and LA stage of both CHA treated and control materials (Additional file 13). As a result, total free AA at YB, MA and MA stages were significantly reduced after SX-1 treatment. Moreover, the concentration of BCAA at YB and MA stages were also largely reduced after SX-1 treatment. What’s more, total flavonoids (including quercetin, kaempferol and isorhamnetin) were obviously reduced at YB and MA stage, while was increased at SA. The total FAs was significantly reduced at MA stage, further analysis revealed that palmitic acid (C16:0) was slightly decreased at YB and largely decreased at MA and LA, and linolenic acid (C18:3) was also reduced at YB, MA and LA, except for a slightly increase at SA. These data indicated that, amino acid, fatty acid and flavonoid biosynthesis might be destroyed to a great extent by SX-1, especially at MA stage.
Interaction analysis for candidate genes and the putative metabolic pathway
The results of proteome, transcriptome and miRNA-seq implied that male sterility induced by SX-1 should be controlled by a complex mechanism. To further elucidate this mechanism, we investigated the known and predicted interactions among the candidate genes corresponding to DEPs, DEGs, target genes of DE-miRNAs and genes reported to be involved in the development of anther or pollen wall based on A. thaliana orthologues (Additional file 14) by STRING 10.0, and then visualized by Cytoscape 3.6.1, and an interaction network associated with anther and pollen wall development were constructed (Fig. 7). It’s interesting that, some biologic processes enriched by KEGG analysis, including protein processing in ER, plant hormone signal transduction, protein export, galactose metabolism, biosynthesis of flavonoid, phenylpropanoid, amino acids, fatty acids, steroid, participated in the anther and pollen wall development through some pivotal proteins, such as EMS1, ACOS5, ERL2, ATA1 and ABCG31.
According to the interaction regulation network and the KEGG analysis (Additional file 5), an enormous metabolic pathway related to the male sterility induced by SX-1 at different stage was constructed (Additional file 15). Firstly, protein processing in ER and flavonoid biosynthesis were destroyed. Subsequently, phenylpropanoid biosynthesis, the upstream pathway of flavonoid biosynthesis was down-regulated at SA stage, might be for the tight interaction between 4CL5, CYP84A1, At1g80820, At4g05160, CAD5 and At4g19010 in phenylpropanoid biosynthesis and all the flavonoid biosynthesis proteins showed in Fig. 7. Meanwhile, galactose metabolism was also found down-regulated. DIN10, reported to responding to the level of sugar in the cell, was exhibited slight link with AT1G24320, a six-hairpin glycosidases superfamily protein in protein processing in ER. What’s more, one up-regulated bra-miR9556-5p were identified located at phenylpropanoid biosynthesis pathway by repressing the expression of UGT72E3. Afterwards, at MA stage, fatty acid elongation was affected, and the LACS3, LACS6, LACS8 showed strong interaction with AT4G19010, AT4G05160 and 4CL5 in phenylpropanoid biosynthesis, and CYP98A3, C4H and TT7 in flavonoid biosynthesis. As well as the protein export and plant hormone signal transduction displayed strong interaction with other pathway by some important proteins, such as BIP2, At1g29310, At2g34250, ATERDJ2B, BRI1 and BRU6. Interestingly, one novel miRNA, samples_miRNA_520 targeting IPTs supposedly, which were of importance in zeatin biosynthesis, was significant up-regulated at MA stage. Finally, amino acids, fatty acids and steroid biosynthesis at the LA stage were down-regulated, which were the upstream of pyruvate metabolism, fatty acid elongation and brassinosteroid biosynthesis respectively, while these pathways were all down-regulated at MA stage. Meanwhile, the fatty acid degradation was up-regulated. Additionally, one novel miRNA, samples_miRNA_201 was supposed to target FK exhibiting up-regulated at steroid biosynthesis.
A potential regulation model for male sterility induced by SX-1
According to the results mentioned above, we proposed a possible regulation model for male sterility induced by SX-1 in B. napus (Fig. 8). After spraying on leaves, SX-1 was mainly absorbed by mesophyll, and then polar-transported to the anthers largely through the phloem. Firstly, CSR1, the target gene of SX-1 was down-regulated. At the same time, the accumulation of SX-1 in the anthers might affect the expression of several important TFs, such as SPL, MYB80 and NAC025, partially through the regulation of TFs/miRNAs modules, for example, miR156/SPL at YB stage, when the protein processing in ER and flavonoid biosynthesis were destroyed. Secondly, probably because of the disordered flavonoid biosynthesis ahead, its upstream pathway, phenylpropanoid biosynthesis which associated with bra-miR9556-5p were down-regulated at SA stage, as well as galactose metabolism. Thirdly, besides the galactose metabolism, due to the early cellular reprogramming in the response to SX-1, plant hormone signal transduction, pyruvate metabolism and fatty acid elongation were all down-regulated. Moreover, samples_miRNA_520, predictively targeting IPTs was identified up-regulated. In the meanwhile, bra-miR5654b/PPR modules were disordered. Then the delayed PCD process in tapetum during mitosis was happened, and the tapetum cells stacked together subsequently at LA stage. At the same time, some upstream pathways of the down-regulated ones at MA stage, such as steroid, amino acids and fatty acids biosynthesis were also imbalanced. Meanwhile, the function of some transport proteins, such as LTPs, ABCGs were impeded. Finally, the tapetum couldn’t timely provide young microspores nutrients, such as flavonoids, amino acids and lipids, resulting in the defective in tryphine of pollen wall and completely crinkled pollen grains with little cytoplasm and abnormal high content of vacuole, and male sterility at last.