Skip to main content

Deficiencies in the formation and regulation of anther cuticle and tryphine contribute to male sterility in cotton PGMS line



Male sterility is a simple and efficient pollination control system that is widely exploited in hybrid breeding. In upland cotton, CCRI9106, a photosensitive genetic male sterile (PGMS) mutant isolated from CCRI040029, was reported of great advantages to cotton heterosis. However, little information concerning the male sterility of CCRI9106 is known. Here, comparative transcriptome analysis of CCRI9106 (the mutant, MT) and CCRI040029 (the wild type, WT) anthers in Anyang (long-day, male sterile condition to CCRI9106) was performed to reveal the potential male sterile mechanism of CCRI9106.


Light and electron microscopy revealed that the male sterility phenotype of MT was mainly attributed to irregularly exine, lacking tryphine and immature anther cuticle. Based on the cytological characteristics of MT anthers, anther RNA libraries (18 in total) of tetrad (TTP), late uninucleate (lUNP) and binucleate (BNP) stages in MT and WT were constructed for transcriptomic analysis, therefore revealing a total of 870,4 differentially expressed genes (DEGs). By performing gene expression pattern analysis and protein-protein interaction (PPI) networks construction, we found down-regulation of DEGs, which enriched by the lipid biosynthetic process and the synthesis pathways of several types of secondary metabolites such as terpenoids, flavonoids and steroids, may crucial to the male sterility phenotype of MT, and resulting in the defects of anther cuticle and tryphine, even the irregularly exine. Furthermore, several lipid-related genes together with ABA-related genes and MYB transcription factors were identified as hub genes via weighted gene co-expression network analysis (WGCNA). Additionally, the ABA content of MT anthers was reduced across all stages when compared with WT anthers. At last, genes related to the formation of anther cuticle and tryphine could activated in MT under short-day condition.


We propose that the down-regulation of genes related to the assembly of anther cuticle and tryphine may lead to the male sterile phenotype of MT, and MYB transcription factors together with ABA played key regulatory roles in these processes. The conversion of fertility in different photoperiods may closely relate to the functional expression of these genes. These findings contribute to elucidate the mechanism of male sterility in upland cotton.


Cotton (Gossypium spp.), an important economic crop, is widely used mainly to provide natural fibers for the textile industry [1]. Heterosis utilization is an important way to improve the productivity and quality of cotton fiber [2]. At present, the production of cotton hybrids is still dominated by artificial pollination, which has the disadvantages of large labor demand, high raw material consumption, and high cost. The emergence of male sterile lines provides a new approach to the utilization of heterosis. Male sterile lines mainly involving cytoplasmic male sterility (CMS), genic male sterility (GMS) and environment-sensitive genic male sterility (EGMS) [3]. The application of CMS hybrids is based on a stringent genetic relationship among CMS, maintainer, and restorer lines, and is limited by the inaccessible restorer lines [4]. In the two-line hybrid technology, EGMS lines are used as the male sterile lines and maintainer lines under contrary environmental conditions, and have the advantage of broad restorer lines. The EGMS lines include photoperiod-sensitive genic male sterility (PGMS) lines, thermo-sensitive genic male sterility (TGMS) lines, and the recently reported humidity-sensitive genic male sterility (HGMS) lines [5]. In upland cotton, CCRI9106, a PGMS mutant isolated from CCRI040029 was controlled by one pair of recessive genes, and showed male sterile under long-day conditions and fertile under short-day conditions [6, 7]. Therefore, it will be an ideal genetic material for studying cotton anther development and heterosis utilization.

Anther cuticle is a hydrophobic barrier covering the surface of the land plant anthers and is essential to protect microspores / pollen grains from environmental and biological stresses. Cuticular wax, the major component of the cuticle, predominantly contain alcohols, aldehydes, ketones, alkanes, and esters derived from very long chain fatty acids (VLCFAs) [8]. However, some evidences demonstrates that carotenoids, flavonoids, isoprenoids (terpenoids) and sterol esters are also precursors for the anther cuticle formation [9, 10]. Several TTP-specific genes have been showed to be essential to the early formation of anther cuticle and pollen exine, such as DPW, DPW2, OsGPAT3, DWP3, OsNP1, CYP704B2 and CYP703A3 in rice [11,12,13,14,15,16,17]. In plants, the CER family, also known as the Glossy (GL) family, have been reported to be involve in the accumulation of epicuticular waxes. Such genes include CER1, CER6/KCS6 and LACS1/CER8 in Arabidopsis, the GL1 gene in maize, Wax Crystal-Sparse Leaf 2 (WSL2), OsGL1–1 and OsGL1–3 in rice [18,19,20,21,22,23,24,25,26].

Pollen coat (tryphine), the outermost layer of the pollen grain, protects the released pollen grain from desiccation and damage, enabling the completion of pollination, is reported assembly by compounds such as carotenoids, flavonoids, fatty acids and isoprenoids [27, 28]. Given the similar lipid-soluble precursors, as previously reported, mutation of anther cuticle synthesis genes often leads to the defects of tryphine. In Arabidopsis, CER1, CER6/KCS6 and LACS1/CER8 have been reported to participate in the pathways of cuticular wax and cutin biosynthesis, and their mutants display reduced wax deposition on stems and fruits, also characterized with abortion pollens due to the deletion of the pollen coat [22,23,24,25,26]. WDA1, the homologous gene of CER1-like in rice, is involved in cuticle and wax formation of rice anther walls and is essential for pollen development. Mutant plants of wda1 appear to lack the pollen wall and the epicuticular wax crystals on the outer layer of the anther [29]. Recently, a HGMS line is isolated in rice and the corresponding gene (HMS1) is characterized as a putative β-ketoacyl-CoA synthase (KCS). The hsm1 mutants show abnormal anther cuticle without white crystals and significant reduction of bacula and tryphine in pollen wall [5]. A glossy family protein, OsGL1–4, is require for the synthesis of very long chain alkanes, and is essential to the pollen adhesion and hydration. Mutant plants generate low-humidity inviable pollen due to the deficient pollen coat with a significant reduction in C25 and C27 alkanes [30]. Further, these genes are preferentially expressed during the later stages of pollen development, mainly abundant at S10-S11(late uninucleate pollen stage). Taken together, anther cuticle and tryphine may share the similar lipid-soluble precursors synthesized during the late uninucleate pollen stage, and mutations of these genes may cause abnormal anther surface and defective pollen wall, thus resulted in male sterility.

The regulation mechanism of the cuticle synthesis is complex, involving multiple signal networks, and the regulation process is mainly achieved through a variety of transcription factors and plant hormones [9]. One of the transcription factors that regulate cuticle formation is WIN1 /SHN1, which directly regulates cutin synthesis-related genes, such as LACS2, GPAT4, CYP86A4 and CYP86A7,and is regulated by GA and may be completed by DELLA protein [31]. In addition, the AP2/ERF transcription factor DEWAX negatively regulates wax synthesis during photoperiod regulation, while CER1, LACS2 and ECR genes are its direct targets [32]. Previous studies have shown that MYB96 promotes the accumulation of cuticular wax under drought conditions or ABA treatment, which may be achieved by actively regulating the expression of wax synthesis genes such as KCS1, KCS2, KCS6, KCR1 and CER3 [33, 34]. MYB94, a paralogous gene of MYB96, which is located upstream of WSD1, KCS2/DAISY, CER2 and ECR, is a positive regulator related to wax biosynthesis caused by drought and ABA treatment [35]. MYB30 has been identified to improve disease resistance by altering the composition and surface properties of the wax in apple and Arabidopsis [36]. The transcription level of MYB30 in apple increase sharply after 12-h of ABA treatment, but do not show any changes under GA and NaCl treatments [37]. MYB106 and MYB16 were found to participate in the regulation of cuticle biosynthesis through partially cooperated with WIN1/SHN1 [38]. MYB16 was showed act as a major regulator of cuticle formation by regulating the cuticle biosynthesis genes CYP86A8/LCR and CER1 [39]. These studies indicate that TFs, particularly MYBs, play a curial role for the regulation of cuticle formation by activating the expression of synthesis genes. Moreover, given the inducible expression to TFs, we consider ABA may also mediate cuticle formation during organ development and stress response. Similar mutant phenotypes of ABA synthetic genes (ABA2 and ABA3) and LACS2 confirm this view [40].

In the present study, morphological observation of PGMS line CCRI9106 showed that CCRI9106 display relatively delay-growth plant, pale-green leaves and abortion of pollens when compared with its wild type CCRI040029. Moreover, the MT exhibited abnormal anthers with a deleted cuticle surface, aborted pollen grains covered with irregular exine, lacking tryphine and thinner intine. Comparative transcriptome profiling analyses of anthers at three different stages were performed to reveal the key genes and mechanism of MT male sterility. By analyzing transcript expression trends and constructing co-expression networks, we identified several candidate genes which involving in lipidic precursors synthesis of anther cuticle and tryphine. Further, we speculated that MYBs and ABA may be potential regulators that lead to the phenogenesis of MT. The qRT-PCR expression analysis of several lUNP-specific down-regulated genes in MT under two photoperiods revealed that these genes could be activated under short-photoperiod conditions. Our data also provide new insights into the molecular events underlying male sterility.


Phenotypic characteristics and cytological observations of MT and WT

Compared with WT, the PGMS line CCRI9106 (MT) was accompanied by a relatively delay-growth plant and pale green leaf phenotype (Fig. 1a). Further, the chlorophyll and carotenoid content were analyzed to determine the pigment difference between MT and WT leaves. The chlorophyll contents of MT leaves were significantly reduced except the chlorophyll b content in the first leaf. Additionally, the second and fourth leaves of MT showed significant reduction of carotenoid (Fig. 1b and Supplementary Figures S1). The decrease in pigment content indicates that the pigment synthesis pathway of MT may be defective.

Fig. 1

Morphological characteristics of the PGMS MT line CCRI9106 and its WT line CCRI040029. a Intact plants of MT (right) and WT (left). b Chl (a + b) and Car content of leaves in MT and WT (Student’s t-test, *P < 0.05, **P < 0.01). L1-L5, the first to the fifth leaf. c Flower and pollen phenotype of WT and MT. d Paraffin sections in virous development stages of WT and MT anthers. En, endothecium; ePG, empty pollen grain; Mc, mother cells; PG, pollen grain; Td, tetrad; T, tapetum. The red arrows point out the major differences of pollen grains between WT and MT. Bar = 100 μm

Particularly, the MT flower displayed shorter filaments and shriveled anthers (Fig. 1c), with the anthers did not dehisce, and lacked mature pollen grains at the late development stages of anther when stained with 1% I2-KI solution (Fig. 1c).

To understand the developmental abnormalities of anther and pollen in MT, light microscopic examination was applied. During the early developmental stages, there were no remarkable morphological differences between the WT and MT. Microspore mother cells (MMCs) underwent meiosis and successfully formed the tetrads of haploid microspores (Fig. 1d). Subsequently, at the early uninucleate pollen stage, young microspores were freely released from the tetrads (Fig. 1d). With the development of microspores, the WT anther chambers gradually expanded. In contrast, at the same stage, MT had relatively narrow anther chambers (Fig. 1d). At binucleate pollen stage, the WT displayed round pollen grains with abundant inclusions, whereas most of MT pollen grains failed to accumulate storage inclusions, and the microspores were irregular in shape (Fig. 1d). At the mature pollen stage, the anther endothecium layer of WT began to thicken, preparing for the dehiscence of the anther wall. Consistent with the previous observations, the MT anther wall started to collapse, resulting in unvital pollens. These results showed that the abortive phenotype of MT may be caused by abnormal anther and immature pollens.

Defects of anther structure and pollen wall development cause male sterility in MT

Transmission electron microscopy (TEM) was performed to further observe the ultrastructural changes in MT at various pollen development stages. As revealed via bright-field microscopy, the pollen defects of MT were first observed at uninucleate stage. MT displayed chaotic anther cell layers structure and a thinner pollen nexine when compared with WT (Fig. 2C-D, c-d). At binucleate pollen stage, the WT pollen grains had evenly stained cytoplasm which contained numerous small storage bodies (Fig. 2E). At higher magnification, the WT nexine was considerable thickened with the tectum and bacula regularly assembled on nexine surface (Fig. 2F). In contrast, the MT pollen grains showed severe cytoplasmic degradation (Fig. 2e) and abnormal exine assembly with irregular accumulation of tectum and bacula (Fig. 2f). At mature pollen stage, these differences were more apparent. Whereas WT pollen grains had accumulated starch granules as well as a mature exine filled with abundant tryphine (Fig. 2G, H), MT pollen had a completely empty antrum and absent tryphine (Fig. 2g, h). Since it is well known that the synthesis of intine is done under the control of the microspore, it is no wonder that MT aborted pollen had an abnormal intine (Fig. 2h). At later stage, the endothecium of WT anther particularly expanded for dehiscence (Fig. 2I), yet, MT anther showed a thin endothecium which failed to expand (Fig. 2i), resulting in the indehiscence phenotype. Combined with the finds of light microscopy of transverse sections, TSM results confirmed that the male sterility phenotype of MT may be attributed to the irregular exine, immature intine, lacking tryphine and abnormal anther wall development.

Fig. 2

TEM analysis of developing anthers in WT (A-I) and MT (a-i) plants. Cross-sections of WT and MT anthers at the tetrad stage (A, B and a, b), uninucleate pollen stage (C, D and c, d), binucleate pollen stage (E-H and e-h) and mature pollen stage (I and i) are shown. Anther wall is showed in B, C, I (WT) and b, c, i (MT). Pollen is showed in E, G (WT) and e, g (MT). Pollen wall is showed in D, F, H (WT) and d, f, h (MT). Ba, bacula; E, epidermis; En, endothecium; Ex, exine; In, intine; ML, middle layer; Ne, nexine; T, tapetum; Te, tectum; Tr, tryphine. Bar = 10 μm in (A-C, and a-c), 20 μm in (E, G, I, and e, g, i), 2 μm in (D, F, H, and d, f, h). The ellipses and arrows point out the major differences between WT and MT in the same development stages

RNA-seq for anthers of three developmental stages based on cytological observation

Cytological analysis demonstrated that the abortion of MT pollen occurs during the uninucleate stage (Figs. 1,2). Based on this observation, RNA-seq was performed using anther samples in the tetrad (TTP), late uninucleate (lUNP), and binucleate (BNP) to understand the gene expression profile alterations, investigate the underlying mechanisms for the male sterility phenotype. A total of 18 libraries were constructed for the anthers of WT and MT. High throughput sequencing generated more than 50 million clean reads for each sample, which was sufficient for the quantitative analysis of gene expression (Supplementary Table S1). Subsequently, by mapping the clean reads to the upland cotton TM-1 reference genome [41], we obtained exceeded 94% mapped reads in all samples, with up to 85% were unique match reads. Additionally, Q30 ≥ 90.86% and GC contents approximately 43% were found in each sample (Supplementary Table S1). The results above showed that the sequencing data is of high quality.

Fragments PerKilo base of transcript per Million mapped reads (FPKM) was used to calculate the gene expression levels. A total of 52,143 unigenes were obtained with the threshold that FPKM> 1 in at least one of the samples, of which including 44,993 original genes and 7150 new genes. Among the original genes, 8704 differentially expressed genes were filtered with the |log2 fold change| ≥ 1 and FDR ≤ 0.01. 159 (52 up-regulated, 107 down-regulated), 2239 (1334 up-regulated, 905 down-regulated) and 7614 (2971 up-regulated, 4643 down-regulated) DEGs were found in the MT line at three stages, respectively (Supplementary Figure S2A). The significantly increased number of DEGs from TTP to BNP was consistent with the more serious pollen abortion phenotype observed in the cytological analysis of MT anthers. The number of down-regulated genes were larger than up-regulated genes, indicating that genes related to the anther development were affected. Among these genes, 25 DEGs restricted to TTP and lUNP, 1183 to UNP and BNP, 24 to TTP and BNP, and only 38 were differently expression in all three stages, pointing that most DEGs may exhibited specific expression at different stages of anther development (Supplementary Figure S2B).

Down-regulated DEGs probably cause the male sterility phenotype of MT based on the analysis of gene expression patterns

To further provide insights into the functional transitions along anther development, we analyzed the 8704 DEGs by K-means clustering algorithm, and 11 co-expression clusters were finally obtained. Of these, cluster 3(2561 DEGs, accounting for 29.42% of 8704 genes), in which DEGs specific expressed at BNP of WT, were mostly enriched. The BNP-MT specific expression cluster (cluster 9) was the second enriched cluster with 2218 DEGs accounting for 25.48% of all DEGs (Fig. 3a and Supplementary Table S2). These results indicated that the late pollen maturation process requires precise regulation of large numbers of genes, and the down-regulation of these DEGs may cause the arrest of pollen grains development.

Fig. 3

K-means clustering and functional analysis of clusters. a Eleven clusters representing 8704 genes are shown with distinct expression patterns in different stage. Z-score of all DEGs were used to generate the heat map. TTP-W, lUNP-W, BNP-W together with TTP-M, lUNP-M, BNP-M represent different stages of WT and MT, respectively. b The number of genes from different TF families of down- or up-regulated clusters. c Significantly over-represented gene ontology (GO) terms of down- or up-regulated TFs. The size of the circles shows the significance of that node. Nodes that were specifically represented of down-regulated TFs are indicated as blue circles and as red circles for up-regulated TFs. Grey nodes were over-represented non-specifically. d GO functional categories enriched in cluster 1–6. Only top five terms are displayed

Cluster 1–6 displayed a relatively high transcript level in WT, while DEGs in cluster 7–9 were up-regulated in MT. Cluster 10 was special, for the DEGs showed an earlier expression in MT than WT. DEGs in cluster 11 were particularly enriched at TTP, but exhibited differentially expression at other stages with pretty low levels (Fig. 3a and Supplementary Table S2).

To obtain an overview of the potential functions for DEGs, GO enrichment analysis was performed on all down-regulated (cluster 1–6) and up-regulated (cluster 7–9) DEGs, respectively. The analysis of down-regulated DEGs showed representation of genes related to “cell tip growth”, “pollen tube development”, “cell wall organism process”, “actin cytoskeleton organization” and “carbohydrate metabolic process” (Supplementary Figure S3A). These processes were known to be involved in anther development and pollen maturation. A total of 150 TF-encoding genes belonging to 28 families were found from cluster 1–6 (Supplementary Table S3). The members of MYB, NAC, bHLH, C2H2, LBD and MADS families were highly represented (Fig. 3b), and were mainly related to “microsporogenesis”, “floral meristem determinacy”, “pollen maturation” and “secondary cell wall biogenesis” (Fig. 3c). In short, enrichment analyses of down-regulated DEGs and TFs indicated that these functional categories are closely related to the male sterile phenotype. Up-regulated DEGs were marked by processes such as “response to stimulus”, “response to hormone”, “response to oxygen-containing compound” and “response to chemical” (Supplementary Figure S3B), in which 394 TFs belonging to 37 families were found (Table S3), and ERF, WRKY, C2H2, MYB, bHLH and NAC families were overrepresented (Fig. 3b). GO analysis showed that these TFs mainly response to hormones (ethylene, salicylic acid, jasmonic acid) and stimulus (heat, cold, light or water) (Fig. 3c). Taken together, the appearance of these exceptionally GO terms and higher proportions of TFs in cluster 7–9 suggested that the upregulated gene expression may be due to the male sterility phenotype, in other words, the mutant might cause cellular stresses that lead to the general upregulated expression of these genes. Similar founding was reported in previous study [42].

Functional characteristic of DEGs in down-regulated clusters

Based on the above results, the down-regulated DEGs in cluster 1–6 were more likely to be related to male sterility and were therefore selected for further analysis.

The cluster 1 (95 genes, 0 TFs) corresponds to DEGs that were specifically expressed at TTP and then appeared to be low levels or not expressed at the later stages (Fig. 3a and Supplementary Table S2). Previous researches in Arabidopsis and rice supported the idea that pollen exine development, the meiosis of PMC and the degradation of callose walls were processed at the stage around TTP [43]. To explore whether genes in this cluster are beneficial to the above biological process, we performed GO analysis and found the representative significant enrichment terms are “lipid transport”, “external encapsulating structure organization”, “ubiquinone metabolic process”, “carbohydrate metabolic process” and “pollen exine formation” (Fig. 3d). These enriched terms suggested that the genes in cluster 1 may be critical for early anther development, which is consistent with previous researches.

Cluster 2 showed peak expression at lUNP, with 565 genes, including 36 TFs (Fig. 3a and Supplementary Table S2). GO enrichment showed that DEGs in this cluster mainly associated with “oxidation-reduction”, “secondary metabolic”, “lipid biosynthesis”, “single-organism and terpenoid metabolic process” (Fig. 3d). Genes in this cluster may respond to the sterile pollen because the exine and tryphine were assembly with lipidic precursors, and the secondary metabolites and terpenoids are often essential components for tryphine. In addition, nine (25% of 36) MYB TFs (MYB3, MYB7, MYB17, MYB36, MYB48, MYB68, MYB73, MYB85 and MYB105) were found in this cluster (Supplementary Table S2), reflecting their important role in lUNP stage during anther development. Taken together, genes in this cluster may contribute to exine and tryphine formation, and MYB TFs likely take part in these biological progresses as important regulators.

DEGs in cluster 3 were expressed specifically at BNP in WT. 2561 genes, including 68 (2.66% of 2561) TFs were found in this cluster (Fig. 3a and Supplementary Table S2), and the appearance of such a large number of DEGs reflected the criticality of this period for anther development. GO terms such as “cell tip growth”, “pollen tube development”, “actin cytoskeleton organization”, “polysaccharide catabolic process”, “pectin catabolic process” and “plant-type secondary cell wall biogenesis” were significantly enriched in this cluster (Fig. 3d), indicating that these genes may facilitated the pollen maturation, anther dehiscence and pollen tube growth to some extent.

DEGs in cluster 4 (200 genes, 8 TFs), cluster 5 (160 genes, 4 TFs) and cluster 6 (844 genes, 33 TFs) were expressed at more than one of the three stages (Fig. 3a), and all work properly until the down-regulated expression at BNP stage. GO enrichment analysis suggested that these genes may play unique roles for the development of anthers (Fig. 3d). These results once again proved that MT anthers at the genetic level showed obvious defects during the BNP stage, and the abnormal expression of these genes during this period may be important cause of male sterility.

PPI network construction to identified candidate genes in cluster 1–3

DEG expression pattern clustering and GO enrichment analysis implied that MT male sterility should be controlled by a complex mechanism, and DEGs in clusters 1–3 may be closely related to it. To elucidate this mechanism and further identify key genes and pathways that contribute to MT male sterility, we investigated the known and predicted interactions among genes in cluster 1–3, respectively. Genes in the networks were displayed by their Arabidopsis homologous genes. Functional analysis of those interacting genes was performed using ClueGO software and visualized by Cytoscape 3.3.0.

In network 1, ATA1, AT3G23770, and MEE48 were shown to interact with each other (Supplementary Table S4). ATA1 expressed specifically in tapetal cells and was peaked at the early anther development stage together with AT3G23770 [44]. Even though their exact role has not been fully elucidated, the interactions with the proven male sterile gene MEE48 in the network also suggest their probably functions to the anther development.

Eighty-four genes were found in network 2 and these genes are related to lipid biosynthetic process and the synthesis pathways of several types of secondary metabolites such as terpenoids, flavonoids and steroids (Fig. 4a and Supplementary Table S4). Notably, the productions of these process are exactly precursors for the formation of tryphine, indicating that the formation processes of MT tryphine are indeed abnormal, which is consistent with the cytological defects, thus confirming that the infertility of MT is closely related to genes in this cluster. Additionally, several male sterility genes are present in this network, such as CAS1, ABCD1, KCS6 and LACS1 (Fig. 4a). In terms of the abortive phenotype caused by the mutants, we found these genes not only affect the formation of pollen exine and tryphine, but are also required for the formation of anther cuticle [22,23,24,25, 45, 46]. Besides, anther cuticle shares similar precursor components with tryphine [9, 10]. Therefore, we considered whether the genes in cluster 2 also response to the formation of anther cuticle, and their down-regulation may affect the anther surface of MT. In order to testify this hypothesis, scanning electron microscope (SEM) analysis was carried out. At pollen maturation stage, the outer surface of the wild-type anther was covered by well-formed cuticle, while the MT anther surface was quite smooth and cuticle seems absent, which is consistent with our inference (Fig. 4b). In addition, CHLI, CYP97A3, LUT2, SBPASE, UGT78D2, ABCD1, ECHID and FAB1 were identified as key genes in this network. Of these, CYP97A3 and ABCD1 were reported involved in pollen development [46, 47], indicating that they may be key genes responsible for the male sterility phenotype of MT. Interestingly, several pigment synthesis genes like CHLI1, HEME1, PDS3, CYP97A3, CYP97B3 and LUT2 were found in this network, and involved in multiple biological processes (Fig. 4a). The obvious period-specific characteristics in anthers suggesting that they may also play an important role in the development of anthers (Table S4). Further, we used qRT-PCR to detect the expression of several genes in this network, and found the results were similar to those of RNA-seq data (Fig. 4c and Supplementary Table S7). These results indicated that genes in this cluster performed various and significant functions on the formation of anther cuticle and tryphine, and their abnormal expression might be the direct cause of the male sterility phenotype of MT.

Fig. 4

PPI network of cluster 2 DEGs. a GO modules enriched of PPI DEGs in cluster 2 visualized by the ClueGO plug-in in Cytoscape. Genes with blue font are MS genes reported in Arabidopsis. Trigonal node represents key genes in the network. b Appearance of the anther surface in WT and MT under scanning electron microscopy. The red arrow points out the major differences between WT and MT anther. c qRT-PCR expression analysis of DEGs in the network. TTSET is performed only in lUNP stage. **P < 0.01, ***P < 0.001

Two hundred five genes were found in network 3, and their involvement are consistent with the GO terms mentioned in cluster 3, confirming their potential roles to the pollen maturation (Supplementary Figure S4 and Supplementary Table S4). Based on the reported MS genes in Arabidopsis [48], several genes which associated with these biological processes were found, such as FIM5, PLIM2a, PLIM2b, PLIM2c, WLIM2, RGP1, UGP2, PGA4, NST1 and CESA4. Numbers of 7 key nodes with high degree were determined and were included in 4 processes (Supplementary Figure S4). Of these, CAP1 may have a central role for the interconnection of these biological processes (Figure S4). CAP1 encodes an actin monomer binding protein that accelerates the exchange of ADP for ATP, is key intermediate between actin-depolymerizing factor (ADF) mediated disassembly and the profilin-based nucleation and elongation machinery, and was reported as essential regulator of pollen tube growth [49]. No TFs were found in this network (Figure S4), indicating that DEGs in this cluster were structural genes located downstream. Down-regulation of these genes may be the results of microspore abortion, in other words, they may be needed for pollen grains maturation and pollen tube growth rather than microspore development.

Hub genes identification through WGCNA

Based on the analysis above, we speculated that the down-regulation of genes specifically expressed at lUNP stage are critical to the phenotype of MT, and several key genes associated with virous biological process were identified. But the regulatory mechanism is poor understood. Therefore, weighted correlation networks were constructed for the DEGs based on the pairwise correlations between genes in their common expression trends across all samples to identify key genes and the potential mechanism that are highly associated with the male sterile phenotype of MT, and 10 distinct modules were classified. The module eigengenes for 10 modules were correlated with different samples (Fig. 5a).

Fig. 5

The co-expression modules of the valid DEGs analyzed via WGCNA. a The correlation coefficients between modules and samples. Ten module eigengenes are generated and named as the left lane presented. b The eigengene expression levels in the yellow module

Of particular interest to us is the lUNP anther specific expression module (MEyellow) (Fig. 5b), which contains 526 genes, of these including 31 TFs (Supplementary Table S5). For these genes, 461 (87.6% of 526) were interacting genes in the PPI network of cluster 2 (Fig. 4a and Supplementary Table S5), the closely related community formed by these genes indicate their important role to lUNP pollen development. To further detected hub genes of the yellow module, a co-expression network was obtained and then visualized by Cytoscape 3.3.0, and the genes with higher connectivity were showed in Fig. 6a. Furthermore, several genes were identified as hub genes. Strikingly, most of them are directly related to lipids. Ghir_D07G012910 (NPC2), a non-specific phospholipase C2 protein with the central location in the network, is involved in gametophyte development. In Arabidopsis, double mutants of npc2–1 npc6–2 exhibited reduced viability of ovules and pollens [50]. A small rubber particle protein3 (Ghir_D11G032630) which associates with lipid droplet surfaces, was reported to be a positive regulator of tissue growth and development and was induced by ABA [51]. Other highly connected hub genes include 2 glycosylphosphatidylinositol-anchored lipid protein transfer 1 (Ghir_D04G019930/Ghir_A04G015210) which bind to lipids and function as a component of the cuticular lipid export machinery that performs extensive export of intracellular lipids from epidermal cells to the surface to build the cuticular wax layer and silique walls [52], a LTP family protein Ghir_A08G018060 with unclear function, and a isoprenoid biosynthesis enzyme DXR (Ghir_D09G017510) which caused a seedling lethal, albino phenotype when were knockout or strongly silenced [53]. Interestingly, Ghir_D13G017030 (CYP97A3) and Ghir_D10G003560 (CHLI1) were also identified as hub genes, which is consistent with the results of the PPI network 2 (Figs. 4a, 6a), thus highlighting them as possible key regulators of the developing anther.

Fig. 6

Information of yellow module network and DEGs expression characteristics. a The co-expression network of the DEGs in yellow module. Central genes are referred to as hub genes. The color gradation from turquoise to red and the increase of node size represents betweenness and degree from low to high, respectively. b and c qRT-PCR expression analysis of several hub genes or ABA-related genes in this network. TTSET is performed only in lUNP stage. *P < 0.05, **P < 0.01, ***P < 0.001. d ABA content of anthers in virous development stages of MT and WT. ***P < 0.001

In this network, 10 MYB, 4 NAC, 4 LBD, 3 C2H2, 3 bHLH, 2 MIKC_MADS, 2 ERF, 1 GRF, 1 Dof and 1 B3 were found (Supplementary Table S5). Therefore, MYB TFs which with the highest proportion (32.26%) were considered to have an important role in UNP pollen development. In Fig. 6a, MYB3 (Ghir_A08G003880), MYB7 (Ghir_A07G002060) and MYB48 (Ghir_A13G024180) which included in the PPI network of cluster 2 were also found in this module, indicating these potential roles to the lUNP anther development. In addition, Ghir_A12G017450 (MYB16) was showed a higher connection with hub genes (Fig. 6a). Taken together, MYB TFs may be crucial upstream key genes regulating the expression of other genes in this module, and further became the facilitators to the achievement of anther cuticle and pollen tryphine during lUNP stage. Additionally, these hub genes were taken to perform quantitative RT-PCR. The significantly down-regulation of these genes indicated their potential roles to lUNP anther development (Fig. 6b and Supplementary Table S7).

ABA signal pathway was altered in MT anthers

Previous studies demonstrated that ABA plays an important role in the regulation of plant cuticle formation, and ABA biosynthesis and signaling transduction were impaired in mutants with disrupted cuticle biosynthesis [54].

In this study, we found several hub genes were shown to be regulated by ABA, such as Ghir_A01G008890 (AT1G16310), Ghir_D01G009320 (AT1G16310), Ghir_D07G021480 (AT3G25280), Ghir_D08G014410 (MAKR6), Ghir_D11G032630 (AT3G05500), together with Ghir_A07G002060 (MYB7), Ghir_D10G003560 (CHLI1) and Ghir_A12G017450 (MYB16) (Fig. 6a). Further, we analyzed the functions of all DEGs in this module based on the homologous annotations in Arabidopsis and found a total of 28 genes in this module have been shown to be involved in the biosynthesis, transport or response of ABA (Supplementary Table S6), of these, Ghir_A07G002060 (MYB7), Ghir_D10G003560 (CHLI1), Ghir_D12G021400 (CYP707A3), Ghir_D13G004890 (FAAH), Ghir_A01G017870 (KCS6) and Ghir_D01G019360 (KCS6) were also found in the PPI network of cluster 2 (Fig. 4a), indicating that ABA signal may be a key regulator to lUNP anther development. Meanwhile, the expression trends of several genes assessed by qPCR (Fig. 6c and Supplementary Table S7) were consistent with the results of RNA-seq, supporting our deduction. Furthermore, the ABA content in WT and MT anthers was measured during different development stages, and the results suggested that the ABA content of MT anthers decreased across all stages (Fig. 6d). These results suggested that ABA signaling process was affected during MT pollen development, and the reduced endogenous ABA levels might lead to the inhibition of genes involved in the formation of anther cuticle, pollen exine and tryphine.

Analysis of lUNP-specific down-regulated genes in MT under two photoperiods

Compared with the male sterile phenotype in Anyang, the PGMS line CCRI9106 is male fertile under short-day conditions when planted in Sanya, China [55]. To examine whether the different phenotypes were reflected at the transcriptional level, the mRNA expression levels of 34 lUNP-specific down-regulated genes, which held higher connectivity in the PPI network and WGCNA module, were analyzed by qRT-PCR. Anther samples of the same stages (termed as TTP-S, lUNP-S and BNP-S) of MT in Sanya were collected for comparison analysis. As shown in Fig. 7, most of these genes were up-regulated in Sanya condition, indicating their activation by short photoperiod. HEME1, CYP97B3_A07, KCR1 and AAO3 displayed higher expression levels in Sanya than in Anyang at all stages and shed significant activation in TTP-S stage. While CHLI1, FLS1, MYB7, FAB1 and TT4 exhibited a constantly increasing expression to a peak at the BNP-S stage. DXR, Ghir_D11G032630, UGT78D2_D04, UGT78D2_A05 and LUT2 showed down-regulation in Sanya condition, especially under the last two stages. The remaining genes, which accounted for a large proportion, were specifically up-regulated at lUNP-S stage, even though some of them showed slightly down-regulation during BNP-S stage. Therefore, the differences in the expression levels of these genes under the two conditions indicated that the conversion of fertility was closely related to the functional expression of these genes, and the lUNP stage may be crucial for PGMS line CCRI9106 to respond to photoperiod regulation during pollen development.

Fig. 7

qRT-PCR expression analysis of 34 lUNP-specific down-regulated genes in MT under two photoperiods. TTP-A, lUNP-A and BNP-A represent the expression levels of MT under three development stages in Anyang. TTP-S, lUNP-S and BNP-S represent the expression levels of MT under three development stages in Sanya


Understanding the molecular mechanisms controlling anther development and male sterility is of great importance for heterosis utilization in many crop species, including cotton. In recent years, male sterility lines have been widely studied by transcriptomic analysis in different species, such as wheat [56, 57], maize [58], eggplant [59] and cotton [60,61,62,63]. However, most studies tend to carry out extensive enrichment of DEGs without systematic classification analysis based on different development stages. Therefore, the identification of key genes would be hindered because of the huge number of DEGs. In this study, DEGs expression pattern analysis was performed between MT and WT across three anther development stages to reveal the characteristics and possible functions of each cluster, and their contributions to the MT phenotype was evaluated in conjunction with the morphologic observation. Further, pivotal genes were identified via PPI network and WGCNA. These results provide comprehensive information on genes involved in the determination of MT phenotypes.

Defects at lUNP stage lead to the failure of microsporogenesis

Normal microsporogenesis is essential to the maturation of pollen grain and pollination, and it is closely related to anther cell structure development and lipidic pollen wall formation. Almost all male sterility lines display abnormal anther structure and/or defective exine. Functional anther required timely degeneration or thickening of cells in various layers, as well as the normal formation of the cuticle [12, 13]. According to the cytological observation, the anther cell layers seem irregular at UNP stage, and SEM results confirm the defective cuticle covering the surface of MT anther (Figs. 1d, 2C and 4b). In rice, several genes were reported essential to the formation of anther cuticle, such as DPW, CYP704B2 and CYP703A3 [14, 15, 17]. However, their homologous genes, MS2, CYP704B1 and CYP703A2 in Arabidopsis have no effects on the anther cuticle [64,65,66]. Functional differentiation of these genes between monocotyledons and dicotyledons may explains that phenomenon. Therefore, in dicotyledons, like cotton, anther cuticle formation may regulate by other genes. DEGs expressed at the lUNP stage were functional for synthesize or transport fatty acids, carotenoids, flavonoids, isoprenoids (terpenoids) and sterol esters (Figs. 3d and 4a). Interestingly, these components are consistent with the precursors of anther cuticle, which indicate that lUNP-specific genes may be related to the anther cuticle formation in cotton. In terms of microspores, the pollen wall seems abnormal with immature exine assembly and absented tryphine (Figs. 1d and 2). As the previous studies, exine-related genes are mainly expressed at early anther development stage, represented by AMS/TDR, MS1/PTC1, MS2/DPW, CYP703A2/CYP703A3, and CYP704B1/CYP704B2 in Arabidopsis or rice [64,65,66,67,68,69,70], and knock out of these genes often lead to serious defects on pollen exine rather than just irregular assembly like CCRI9106 in this study. This result suggested that exine-related genes expressed at early stage may be normal, which is consistent with the results of RNA-seq. Therefore, because of the similar lipidic precursors of pollen exine and tryphine with anther cuticle, we considered the disruptive microsporogenesis of MT mainly due to the down-regulation of lUNP genes.

Pigment-related genes are associated with the phenotype of MT

In cluster 2, several pigment related genes were found in this network, such as CHLI, HEME1, PDS3, CYP97A3, CYP97B3, DXR and LUT2, some of them even served as hub genes (Fig. 4a). Studies in Arabidopsis indicated that these genes may be involved in the synthesis of carotenoids or chlorophylls [47, 53, 71,72,73,74], and their subcellular localization to chloroplast and the specific expression in lUNP anthers indicated that they may essential to plastid morphogenesis and play an important role in the development of anthers.

Similar results have been found in previous studies, like OsMGD2 and OsDGD2 in rice, which are all localized to chloroplast [75, 76]. Their gene knockout plants have reduced fertility. Additionally, down-regulation of OsMGD2 also caused a decrease of chlorophyll content and plant height [76]. PUR4, another chloroplast localized synthase which involved in the purine biosynthetic pathway, is essential to the male gametophyte development [77]. However, silencing of its homologous genes in cotton leads to variation in leaf color [78]. In addition, there are several male sterility genes in this network, and their loss-of-function plants have color variations in addition to reduced fertility. Mutant plants of CAS1 carrying albino inflorescence shoots because of the low amounts of carotenoids and chlorophylls, and showed male sterility [45]. Plants of pxa1(abcd1) display blue-greenish leaves when exposed to extended night conditions [79]. The lacs1 plants have bright green stems and siliques, reduced fertility under conditions of low humidity, somewhat reduced plant height [24]. Mutant plants of kcs6 have very bright green stems and siliques with no abundant epicuticular wax, also show humidity-sensitive male sterility [22, 23]. The phenotype of these male sterility lines suggested that lUNP-specific genes may involve in anther development as well as plant color variations. In the present study, the male sterility line CCRI9106 displays pale-green leaves during the whole growth period (Fig. 1a, b and Supplementary Figure S1). Combined with the above analyses, we suppose that pigment-related genes may associated with anther development and the leaf phenotype of CCRI9106 may due to the down-regulated of pigment-related genes mentioned above.

Formation of anther cuticle and tryphine may rely on the regulatory effect of MYBs on lipidic genes

MYB TFs comprise one of the largest TF families and were reported involved in controlling anther development and cuticle formation. MYB103 is expressed in the tapetum and microspores during the stage 6 to 10, and play a key role to tapetal and pollen development. MYB103 functions downstream of AMS to activate CYP703A2 in sporopollenin biosynthesis [80]. Double mutants of MYB33 with MYB65 are male sterile with small anthers and defective pollens [81]. Other MYBs, like MYB106, MYB30, MYB41, MYB16 and MYB96 are regulated cuticle formation in virous tissues by activating the biosynthesis genes of epidermal wax in Arabidopsis [34, 37, 39, 82].

In this study, several secondary metabolism-related MYBs were found to be down-regulated at lUNP stage of MT anthers, and may be essential to the synthesis of lipidic precursors for cutin and tryphine. Based on the DEGs expression pattern analysis, nine (25% of 36) MYB TFs were found in cluster 2 (Supplementary Table S3), such as MYB3, MYB7, MYB17, MYB36, MYB48, MYB 68, MYB 73, MYB 85 and MYB105, indicating their potential role to lUNP anther development. In the WGCNA network, several MYBs interacted with lipidic genes and may serve as facilitators to the downstream genes (Fig. 6a). In Arabidopsis, MYB7 and MYB48 were showed to regulate flavanol biosynthesis, and MYB3 was involved in phenylpropanoid biosynthesis [83,84,85]. MYB16 (Ghir_A12G017450) was showed closely connected with hub genes, and its orthologous gene in Arabidopsis was previously reported to be an important role of epidermal cell morphogenesis, acting in the upstream of cuticle biosynthesis genes (such as CYP86A8/LCR and CER1) to regulate cuticle biosynthesis and wax accumulation in reproductive organs and trichomes [39]. Mutant plants of KCS6 and LACS1 exhibited the similar phenotype with cer1 [22,23,24,25], indicating that MYB16 may also be their regulator, even the regulator of more genes in this module. As mentioned above, we propose that MYBs act as global regulators of cuticular substances that regulate anther surface coating and male sterility.

ABA signal deficiency is responsible for the male sterility of MT

Several previous reports demonstrated that ABA plays a key role in the regulation of cuticular wax formation, and numerous wax biosynthesis-related genes were reported to respond to ABA treatment [54]. For example, the expression of ECERIFERUM 1(CER1) gene, which associated with production of stem epicuticular wax and pollen fertility, is induced by ABA treatments [86]. In addition, considerably greater KCS6 transcript was accumulated in ABA-treated plants, and the ABA-responsive elements in the promoters suggested that ABA directly regulates KCS6 expression [87]. Moreover, MYB96, MYB94, MYB30, MYB16 and MYB106, which directly binding to the consensus motifs of some cuticle synthesis genes, are regulated by ABA signaling and are part of the ABA-mediated signaling pathway [34, 37, 38, 82]. These studies suggested that ABA induced the formation of cuticle by activating the expression of cuticular wax synthesis genes.

In this study, many lUNP specific-expressed genes were predicted to be associated with ABA signaling pathway, of these including lipid precursor synthesis genes such as KCS6, LACS1, FAAH, FIB and CHLI1, and MYB TFs like MYB7 and MYB16 (Supplementary Table S6). These results indicated that ABA may exhibit close relationship with lUNP anther development and anther cuticle formation in cotton. Moreover, the decrease in endogenous ABA content in MT anthers further confirmed that the deficiency of ABA signaling may be a key factor to the male sterility phenotype of MT.

Deficiencies at BNP stage may be the results of microspore abortion

Based on the analysis of differential gene expression patterns, we found that many DEGs showed peak expression during the BNP stage (Fig. 3a and Supplementary Table S2), which indicates that genes in this period are required for the normal production of pollen grains. GO enrichment analysis demonstrated that larger numbers of genes in this stage are involved in pollen maturation, pollen tube growth and anther dehiscence (Fig. 3d and Supplementary Figure S4). During pollen maturation, pollen intine gradually forms and thickens, while the pollen cytoplasm is filled with inclusions such as starch and lipids. The aberrant expression of genes involved in these processes may cause the reduction of fertility rather than complete sterility.

In Arabidopsis, FLA3 is found to be specifically expressed in pollen grains and tubes, and the mutant plants have approximately 50% abnormal pollen grains which lacked viability [88]. Two pectate lyase-like gene, OsPLL3 and OsPLL4 exhibit strong and preferential expression in late development anthers, and the knockdown plants display defective pollen maturation process with a partial male sterile phenotype [89]. Mutants of other genes like UPEX1 in Arabidopsis, BcMF2, BcMF8, BcMF9 and BcMF18 in Brassica have the similar reduced male sterile phenotype [90,91,92,93,94]. Meanwhile, mutant lines of these genes often have normal microsporogenesis, but resulted in aberrant pollen grains with abnormal cellulose distribution, lacking intine and cytoplasm. These abortion characteristics indicate that genes expressed in the late development stages are mainly involved in pollen grains maturation, but have few effects on microspores development. In fact, the pollen intine is reported to be formed by microspores, unlike the exine, which is formed mainly by materials released from tapetum cells [90]. Contents such as starch are also synthesized after being transported to the microspores in the form of precursors such as monosaccharides [95]. This suggests that the survival of microspores is a prerequisite for these processes.

In the present study, the anther cuticle and pollen wall were found to be defective (Figs. 2 and 4b), and male gametophytes mainly aborted at uninucleate stage, leading to the failure of microsporogenesis. Thus, the defects of MT pollen grain during BNP stage are likely to be the results of microspore abortion, not the cause.

Short-photoperiod activated the expression of fertility-related genes

Though the importance for the theory and practice to research and insight the fertility conversion of P/TGMS, there are few achievements in revealing its underlying mechanism. Several P/TGMS lines in plants are widely used in the breeding of hybrid plants and their functional genes are cloned, but by what means photoperiod/temperature as a key factor reverses the fertility remains unclear [3, 95,96,97]. An Arabidopsis TGMS line, rvms, was reported with the characteristic that the sterility-fertility conversion was overcomed by slowing of development, highlighting a potential mechanism applicable to TGMS lines [98].

As a PGMS line, CCRI9106 is proven to be male sterile under long-photoperiod conditions and fertile under short-photoperiod conditions [55], but the underling mechanism of fertility conversion is largely unknown. According to our analysis in the present study, it is speculated that genes related to the formation of the anther cuticle and pollen tryphine were essential for the male sterile phenotype of MT. To study if these genes were response to the photoperiod conditions, the expression features of 34 down-regulated genes were explored under 2 environments (Anyang and Sanya, China). As a result, many genes of MT exhibited activated expression levels at the short-photoperiod region (Sanya) during 3 different development stages, especially during the lUNP stage (Fig. 7), which was concordant with the fertility restoration shown by MT in Sanya. Notably, pigment-related genes that involved in the biosynthesis of xanthophylls (or chlorophylls) (PDS3, CHLI1, HEME1, LUT2, CYP97A3 and CYP97B3_A07) and flavonoids (or anthocyanins) (FLS1, TT4 and DFR) showed significant higher expression in Sanya than in Anyang (Fig. 7). It have been reported that the expression of some of these genes (PDS3, CHLI1, HEME1, FLS1 and TT4) was regulated by light [99,100,101,102]. Therefore, we speculated that the expression activation of these genes was caused by the conversion of the photoperiod. Taken together, the down-regulation of lUNP-specific genes in MT in Anyang and their recovery in Sanya suggested that they might be involved in the fertility transition of MT, and the lUNP stage may be a key period for the regulation of photoperiod.

A potential regulation model for the male sterility phenotype of MT

According to the results mentioned above, we produced a probably regulation model for the male sterility phenotype of MT. MYB transcription factors regulate the formation of anther cuticle and tryphine by promoting the expression of lipidic precursor synthesis genes at lUNP stage, and may as the core of the module. Further, ABA serve as a key regulator of these processes by inducing the expression of MYBs and/or lipid-related genes. Under long-photoperiod conditions, the functional expression and precise regulation of these genes in MT are affected, and the microsporogenesis will be abnormal, resulting in immature pollen grains covering irregular exine together with thinner intine and few inclusions. Finally, vital pollen grains are completely absent and male sterility (Fig. 8).

Fig. 8

A possible gene regulation network in the male sterile line CCRI9106. Abnormal expression of lipidic genes disrupts the microsporogenesis and causes pollen abortion


To elucidate the underling mechanism of male sterility of MT, a comparative histological and transcriptome analysis were conducted between the male sterile line CCRI9106 and its wild type line CCRI040029. Cytological observation showed that the defects at lUNP stage may be responsible for the MT phenotype. At the transcriptional level, a total of 8704 DEGs were obtained between MT and WT. Gene expression pattern, PPI network and WGCNA analysis revealed that the down-regulated of genes involved in lipidic precursors synthesis for cutin and tryphine during lUNP stage were significantly associated with MT phenotype. Hub genes such as Ghir_D07G012910 (NPC2), Ghir_D04G019930 (LTPG1), Ghir_A04G015210 (LTPG1), Ghir_D11G032630 (AT3G05500), Ghir_A08G018060 (LTP), Ghir_D13G017030 (CYP97A3) and Ghir_D10G003560 (CHLI) were predicted regulated by MYBs and induced by ABA signaling pathway. In addition, expression levels of lUNP-specific down-regulated genes in MT could be activated under short-photoperiod conditions. Their changes suggested that they might be involved in the fertility transition of MT. Therefore, our study elucidated the possible genes and pathways related to MT phenotype, and provided important insights into the underling mechanism underlying male sterility. Lastly, our work offered a theoretical basis and foundation for further researches on anther development and heterosis utilization in cotton.

Materials and methods

Plant materials

As in our previous reports, upland cotton PGMS line CCRI9106 was generated from CCRI040029 trough space mutation in 2010 [6, 7]. In this study, CCRI9106 and its wild type CCRI040029 were used as the materials, and cultivated in conventional field of Institute of Cotton Research of CAAS (Anyang, Henan, P. R. China) from April to October.

Morphological and cytological observations

The morphological characteristics of grown plants were analyzed and the leaves were placed into mixture of ethanol and acetone (ratio 1:1) to detect absolute pigment content. The quantities of chlorophyll a (Chla), chlorophyll b (Chlb), Car and total Chl were calculated as the previous method described [78]. Pollen grains from anthesis buds of MT and WT were sampled, dissolved into mixed acids (chromic acid/nitric acid/hydrochloric acid, 15/10/5, v/v/v) and then stained by 2% iodine/potassium iodide solution (I2-KI). Flower buds of different lengths were collected in formalin-aceto-alcohol (FAA) for section observation to identify the pollen developmental stages and asses the cytological characteristics of MT. They were then photographed using an Olympus DP72 light microscope. Based on this identification, the anthers of the corresponding periods are selected for SEM and TEM. The detailed processes of paraffin sectioning, SEM and TEM were performed following the method described by Liu et al. (2015) [6]. Additionally, anthers from MT and WT during tetrad stage (bud, ~ 4.5 mm), late uninucleate stage (bud, ~ 8 mm), binucleate stage (bud, ~ 11 mm) were collected, frozen in liquid nitrogen immediately, and then stored at − 80 °C for RNA-seq and real-time PCR analysis.

RNA extraction and illumina sequencing

Anthers of MT and WT at three different stage were collected with three biological replicates for RNA-seq. Total RNA was extracted from anthers using DP441 Kit (Tiangen, China) according to the manufacturer’s instructions. RNA integrity was assessed on a 1% agarose gel, and a Nanodrop 2000 spectrophotometer (Thermo, USA) was used for the measure of RNA concentration and contamination. Then, a total of 18 libraries of 3 stages in WT and MT were sequenced using Illumina HiSeqTM 2000 to generate paired-end reads. The reads files of this study are accessible from the NCBI Sequence Read Archive (SRA) database under Accession Number PRJNA598402.

Low quality reads were first removed, and the resulting clean reads were mapped to the upland cotton TM-1 genome using HISAT2 [41]. StringTie was used to calculate the fragments per kilobase of transcript per million mapped reads (FPKM) for estimating gene expression levels. Differential expression analysis was performed using the DESeq R package. Genes with FDR < 0.01 and normalized Fold Change≥2 were considered differentially expressed.

Bioinformatics analysis

K-means clustering was performed to cluster DEGs with similar expression patterns by Multiexperiment Viewer v4.9.0. Further, Gene Ontology (GO) enrichment analysis of clusters was performed using OmicShare tools ( and ClueGO plugin in Cytoscape 3.3.0. By screen against the STRING database, protein interactions were identified in the 3 clusters. ClueGO plugin was then used in Cytoscape 3.3.0 to visualize the combination of genes and GO terms. Highly co-expressed gene modules were deduced from DEGs using weighted gene co-expression network analysis (WGCNA) [103], and typical genes in the yellow module were represented using Cytoscape 3.3.0.


To determine the expression of DEGs selected in this study, qRT-PCR analyses were carried out as described previously [104]. Gene-specific primers were designed depend on the reference gene sequences with Oligo 7( The cotton actin 7 was used to as an endogenous reference gene to normalize the amount of gene-specific PCR products, and the gene relative expression levels were calculated with the ΔΔCt algorithm method. All qRT-PCR analyses were performed with three biological replicates.

ABA determination

Endogenous ABA levels were estimated using 100 mg (fresh weight) of MT and WT anthers in TTP, lUNP, BNP and MP stages. Three biological replicates from each plant group were analyzed. ABA levels were determined as described by Min et al. (2013) [105].

Availability of data and materials

Sequence data of 18 RNA-Seq have been uploaded to NCBI Sequence Read Archive (SRA) database under the accession PRJNA598402 (









Cytoplasmic male sterility


Differentially expressed genes


Environment-sensitive genic male sterility




False discovery rate


Fragments PerKilo base of transcript per Million mapped reads


Genic male sterility


Gene Ontology


Humidity-sensitive genic male sterility


Iodine/potassium iodide solution


β-ketoacyl-CoA synthase


Late Uninucleate


Microspore mother cells


The Mutant


Photosensitive genetic male sterile


Protein-protein interaction


Quantitative Real-time polymerase chain reaction


RNA sequencing technology


Scanning electron microscope


Sequence Read Archive


Transmission electron microscopy


Transcription factor


Thermo-sensitive genic male sterility




Very long chain fatty acids


Weighted gene co-expression network analysis


The Wild type


  1. 1.

    Ma J, Geng Y, Pei W, Wu M, Li X, Liu G, Li D, Ma Q, Zang X, Yu S, et al. Genetic variation of dynamic fiber elongation and developmental quantitative trait locus mapping of fiber length in upland cotton (Gossypium hirsutum L.). BMC Genomics. 2018;19(1):882.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  2. 2.

    Zhu W, Liu K, Wang X-D. Heterosis in yield, fiber quality, and photosynthesis of okra leaf oriented hybrid cotton (Gossypium hirsutum L.). Euphytica. 2008;164(1):283–91.

    CAS  Article  Google Scholar 

  3. 3.

    Zhou H, Liu Q, Li J, Jiang D, Zhou L, Wu P, Lu S, Li F, Zhu L, Liu Z, 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.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  4. 4.

    Horn R, Gupta KJ, Colombo N. Mitochondrion role in molecular basis of cytoplasmic male sterility. Mitochondrion. 2014;19 Pt B:198–205.

    PubMed  Article  CAS  PubMed Central  Google Scholar 

  5. 5.

    Chen H, Zhang Z, Ni E, Lin J, Peng G, Huang J, Zhu L, Deng L, Yang F, Luo Q, et al. HMS1 interacts with HMS1I to regulate very-long-chain fatty acid biosynthesis and the humidity-sensitive genic male sterility in rice (Oryza sativa). New Phytol. 2019;255(5):2077–93.

  6. 6.

    Liu J, Pang C, Wei H, Song M, Meng Y, Ma J, Fan S, Yu S. iTRAQ-facilitated proteomic profiling of anthers from a photosensitive male sterile mutant and wild-type cotton (Gossypium hirsutum L.). J Proteome. 2015;126:68–81.

    CAS  Article  Google Scholar 

  7. 7.

    Liu J, Pang C, Wei H, Song M, Meng Y, Fan S, Yu S. Proteomic analysis of anthers from wild-type and photosensitive genetic male sterile mutant cotton (Gossypium hirsutum L.). BMC Plant Biol. 2014;14:390.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  8. 8.

    Li H, Zhang D. Biosynthesis of anther cuticle and pollen exine in rice. Plant Signal Behav. 2010;5(9):1121–3.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  9. 9.

    Nawrath C, Schreiber L, Franke RB, Geldner N, Reina-Pinto JJ, Kunst L. Apoplastic diffusion barriers in Arabidopsis. Arabidopsis Book. 2013;11:e0167.

    PubMed  PubMed Central  Article  Google Scholar 

  10. 10.

    Kunst L, Samuels AL. Biosynthesis and secretion of plant cuticular wax. Prog Lipid Res. 2003;42(1):51–80.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  11. 11.

    Mondol PC, Xu D, Duan L, Shi J, Wang C, Chen X, Chen M, Hu J, Liang W, Zhang D. Defective Pollen Wall 3 (DPW3), a novel alpha integrin-like protein, is required for pollen wall formation in rice. New Phytol. 2020;225(2):807–22.

  12. 12.

    Men X, Shi J, Liang W, Zhang Q, Lian G, Quan S, Zhu L, Luo Z, Chen M, Zhang D. Glycerol-3-phosphate Acyltransferase 3 (OsGPAT3) is required for anther development and male fertility in rice. J Exp Bot. 2017;68(3):513–26.

    CAS  PubMed  PubMed Central  Google Scholar 

  13. 13.

    Liu Z, Lin S, Shi J, Yu J, Zhu L, Yang X, Zhang D, Liang W. Rice no pollen 1 (NP1) is required for anther cuticle formation and pollen exine patterning. Plant J. 2017;91(2):263–77.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  14. 14.

    Li H, Pinot F, Sauveplane V, Werck-Reichhart D, Diehl P, Schreiber L, Franke R, Zhang P, Chen L, Gao Y, et al. Cytochrome P450 family member CYP704B2 catalyzes the {omega}-hydroxylation of fatty acids and is required for anther cutin biosynthesis and pollen exine formation in rice. Plant Cell. 2010;22(1):173–90.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  15. 15.

    Shi J, Tan H, Yu XH, Liu Y, Liang W, Ranathunge K, Franke RB, Schreiber L, Wang Y, Kai G, et al. Defective pollen wall is required for anther and microspore development in rice and encodes a fatty acyl carrier protein reductase. Plant Cell. 2011;23(6):2225–46.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  16. 16.

    Xu D, Shi J, Rautengarten C, Yang L, Qian X, Uzair M, Zhu L, Luo Q, An G, Wassmann F, et al. Defective Pollen Wall 2 (DPW2) encodes an acyl Transferase required for Rice pollen development. Plant Physiol. 2017;173(1):240–55.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  17. 17.

    Yang XJ, Wu D, Shi JX, He Y, Pinot F, Grausem B, Yin CS, Zhu L, Chen MJ, Luo ZJ, et al. Rice CYP703A3, a cytochrome P450 hydroxylase, is essential for development of anther cuticle and pollen exine. J Integr Plant Biol. 2014;56(10):979–94.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  18. 18.

    Zhou X, Li L, Xiang J, Gao G, Xu F, Liu A, Zhang X, Peng Y, Chen X, Wan X. OsGL1-3 is involved in cuticular wax biosynthesis and tolerance to water deficit in rice. PLoS One. 2015;10(1):e116676.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  19. 19.

    Qin BX, Tang D, Huang J, Li M, Wu XR, Lu LL, Wang KJ, Yu HX, Chen JM, Gu MH, et al. Rice OsGL1-1 is involved in leaf cuticular wax and cuticle membrane. Mol Plant. 2011;4(6):985–95.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  20. 20.

    Mao B, Cheng Z, Lei C, Xu F, Gao S, Ren Y, Wang J, Zhang X, Wang J, Wu F, et al. Wax crystal-sparse leaf2, a rice homologue of WAX2/GL1, is involved in synthesis of leaf cuticular wax. Planta. 2012;235(1):39–52.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  21. 21.

    Sturaro M, Hartings H, Schmelzer E, Velasco R, Salamini F, Motto M. Cloning and characterization of GLOSSY1, a maize gene involved in cuticle membrane and wax production. Plant Physiol. 2005;138(1):478–89.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  22. 22.

    Fiebig A, Mayfield JA, Miley NL, Chau S, Fischer RL, Preuss D. Alterations in CER6, a gene identical to CUT1, differentially affect long-chain lipid content on the surface of pollen and stems. Plant Cell. 2000;12(10):2001–8.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  23. 23.

    Millar AA, Clemens S, Zachgo S, Giblin EM, Taylor DC, Kunst L. CUT1, an Arabidopsis gene required for cuticular wax biosynthesis and pollen fertility, encodes a very-long-chain fatty acid condensing enzyme. Plant Cell. 1999;11(5):825–38.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  24. 24.

    Lu S, Song T, Kosma DK, Parsons EP, Rowland O, Jenks MA. Arabidopsis CER8 encodes LONG-CHAIN ACYL-COA SYNTHETASE 1 (LACS1) that has overlapping functions with LACS2 in plant wax and cutin synthesis. Plant J. 2009;59(4):553–64.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  25. 25.

    Jessen D, Olbrich A, Knufer J, Kruger A, Hoppert M, Polle A, Fulda M. Combined activity of LACS1 and LACS4 is required for proper pollen coat formation in Arabidopsis. Plant J. 2011;68(4):715–26.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  26. 26.

    Aarts MG, Keijzer CJ, Stiekema WJ, Pereira A. Molecular characterization of the CER1 gene of arabidopsis involved in epicuticular wax biosynthesis and pollen fertility. Plant Cell. 1995;7(12):2115–27.

    CAS  PubMed  PubMed Central  Google Scholar 

  27. 27.

    Hernandez-Pinzon I, Ross JH, Barnes KA, Damant AP, Murphy DJ. Composition and role of tapetal lipid bodies in the biogenesis of the pollen coat of Brassica napus. Planta. 1999;208(4):588–98.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  28. 28.

    Xue Z, Xu X, Zhou Y, Wang X, Zhang Y, Liu D, Zhao B, Duan L, Qi X. Deficiency of a triterpene pathway results in humidity-sensitive genic male sterility in rice. Nat Commun. 2018;9(1):604.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  29. 29.

    Jung KH, Han MJ, Lee DY, Lee YS, Schreiber L, Franke R, Faust A, Yephremov A, Saedler H, Kim YW, et al. Wax-deficient anther1 is involved in cuticle and wax production in rice anther walls and is required for pollen development. Plant Cell. 2006;18(11):3015–32.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  30. 30.

    Yu B, Liu L, Wang T. Deficiency of very long chain alkanes biosynthesis causes humidity-sensitive male sterility via affecting pollen adhesion and hydration in rice. Plant Cell Environ. 2019;42(12):3340–54.

  31. 31.

    Shi JX, Malitsky S, De Oliveira S, Branigan C, Franke RB, Schreiber L, Aharoni A. SHINE transcription factors act redundantly to pattern the archetypal surface of Arabidopsis flower organs. PLoS Genet. 2011;7(5):e1001388.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  32. 32.

    Go YS, Kim H, Kim HJ, Suh MC. Arabidopsis cuticular wax biosynthesis is negatively regulated by the DEWAX gene encoding an AP2/ERF-type transcription factor. Plant Cell. 2014;26(4):1666–80.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  33. 33.

    Lee SB, Suh MC. Recent advances in cuticular wax biosynthesis and its regulation in Arabidopsis. Mol Plant. 2013;6(2):246–9.

    CAS  PubMed  Article  Google Scholar 

  34. 34.

    Seo PJ, Lee SB, Suh MC, Park MJ, Go YS, Park CM. The MYB96 transcription factor regulates cuticular wax biosynthesis under drought conditions in Arabidopsis. Plant Cell. 2011;23(3):1138–52.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  35. 35.

    Lee SB, Suh MC. Advances in the understanding of cuticular waxes in Arabidopsis thaliana and crop species. Plant Cell Rep. 2015;34(4):557–72.

    CAS  PubMed  Article  Google Scholar 

  36. 36.

    Raffaele S, Vailleau F, Leger A, Joubes J, Miersch O, Huard C, Blee E, Mongrand S, Domergue F, Roby D. A MYB transcription factor regulates very-long-chain fatty acid biosynthesis for activation of the hypersensitive cell death response in Arabidopsis. Plant Cell. 2008;20(3):752–67.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  37. 37.

    Zhang YL, Zhang CL, Wang GL, Wang YX, Qi CH, Zhao Q, You CX, Li YY, Hao YJ. The R2R3 MYB transcription factor MdMYB30 modulates plant resistance against pathogens by regulating cuticular wax biosynthesis. BMC Plant Biol. 2019;19(1):362.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  38. 38.

    Oshima Y, Shikata M, Koyama T, Ohtsubo N, Mitsuda N, Ohme-Takagi M. MIXTA-like transcription factors and WAX INDUCER1/SHINE1 coordinately regulate cuticle development in Arabidopsis and Torenia fournieri. Plant Cell. 2013;25(5):1609–24.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  39. 39.

    Oshima Y, Mitsuda N. The MIXTA-like transcription factor MYB16 is a major regulator of cuticle formation in vegetative organs. Plant Signal Behav. 2013;8(11):e26826.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  40. 40.

    L'Haridon F, Besson-Bard A, Binda M, Serrano M, Abou-Mansour E, Balet F, Schoonbeek HJ, Hess S, Mir R, Leon J, et al. A permeable cuticle is associated with the release of reactive oxygen species and induction of innate immunity. PLoS Pathog. 2011;7(7):e1002148.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  41. 41.

    Wang M, Tu L, Yuan D, Zhu SC, Li J, Liu F, Pei L, Wang P, Zhao G, et al. Reference genome sequences of two cultivated allotetraploid cottons, Gossypium hirsutum and Gossypium barbadense. Nat Genet. 2019;51(2):224–9.

    CAS  PubMed  Article  Google Scholar 

  42. 42.

    Feng B, Lu D, Ma X, Peng Y, Sun Y, Ning G, Ma H. Regulation of the Arabidopsis anther transcriptome by DYT1 for pollen development. Plant J. 2012;72(4):612–24.

    CAS  PubMed  Article  Google Scholar 

  43. 43.

    McCormick S. Control of male gametophyte development. Plant Cell. 2004;16(Suppl):S142–53.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  44. 44.

    Lebel-Hardenack S, Ye D, Koutnikova H, Saedler H, Grant SR. Conserved expression of a TASSELSEED2 homolog in the tapetum of the dioecious Silene latifolia and Arabidopsis thaliana. Plant J. 1997;12(3):515–26.

    CAS  PubMed  Article  Google Scholar 

  45. 45.

    Babiychuk E, Bouvier-Nave P, Compagnon V, Suzuki M, Muranaka T, Van Montagu M, Kushnir S, Schaller H. Albinism and cell viability in cycloartenol synthase deficient Arabidopsis. Plant Signal Behav. 2008;3(11):978–80.

    PubMed  PubMed Central  Article  Google Scholar 

  46. 46.

    Footitt S, Dietrich D, Fait A, Fernie AR, Holdsworth MJ, Baker A, Theodoulou FL. The COMATOSE ATP-binding cassette transporter is required for full fertility in Arabidopsis. Plant Physiol. 2007;144(3):1467–80.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  47. 47.

    Kim J, DellaPenna D. Defining the primary route for lutein synthesis in plants: the role of Arabidopsis carotenoid beta-ring hydroxylase CYP97A3. Proc Natl Acad Sci U S A. 2006;103(9):3474–9.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  48. 48.

    Cui X, Wang Q, Yin W, Xu H, Wilson ZA, Wei C, Pan S, Zhang D. PMRD: a curated database for genes and mutants involved in plant male reproduction. BMC Plant Biol. 2012;12:215.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  49. 49.

    Jiang Y, Chang M, Lan Y, Huang S. Mechanism of CAP1-mediated apical actin polymerization in pollen tubes. Proc Natl Acad Sci U S A. 2019;116(24):12084–93.

    CAS  PubMed  PubMed Central  Google Scholar 

  50. 50.

    Ngo AH, Lin YC, Liu YC, Gutbrod K, Peisker H, Dormann P, Nakamura Y. A pair of nonspecific phospholipases C, NPC2 and NPC6, are involved in gametophyte development and glycerolipid metabolism in Arabidopsis. New Phytol. 2018;219(1):163–75.

    CAS  PubMed  Article  Google Scholar 

  51. 51.

    Chi YH, Kim SY, Lee ES, Jung YJ, Park JH, Paeng SK, Oh HT, Melencion SM, Alinapon CV, Lee SY. AtSRP1, SMALL RUBBER PARTICLE PROTEIN HOMOLOG, functions in pollen growth and development in Arabidopsis. Biochem Biophys Res Commun. 2016;475(2):223–9.

    CAS  PubMed  Article  Google Scholar 

  52. 52.

    Kim H, Lee SB, Kim HJ, Min MK, Hwang I, Suh MC. Characterization of glycosylphosphatidylinositol-anchored lipid transfer protein 2 (LTPG2) and overlapping function between LTPG/LTPG1 and LTPG2 in cuticular wax export or accumulation in Arabidopsis thaliana. Plant Cell Physiol. 2012;53(8):1391–403.

    CAS  PubMed  Article  Google Scholar 

  53. 53.

    Xing S, Miao J, Li S, Qin G, Tang S, Li H, Gu H, Qu LJ. Disruption of the 1-deoxy-D-xylulose-5-phosphate reductoisomerase (DXR) gene results in albino, dwarf and defects in trichome initiation and stomata closure in Arabidopsis. Cell Res. 2010;20(6):688–700.

    CAS  PubMed  Article  Google Scholar 

  54. 54.

    Martin LBB, Romero P, Fich EA, Domozych DS, Rose JKC. Cuticle biosynthesis in tomato leaves is developmentally regulated by abscisic acid. Plant Physiol. 2017;174(3):1384–98.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  55. 55.

    Ma J, Wei H, Liu J, Song M, Pang C, Wang L, Zhang W, Fan S, Yu S. Selection and characterization of a novel photoperiod-sensitive male sterile line in upland cotton. J Integr Plant Biol. 2013;55(7):608–18.

    CAS  PubMed  Article  Google Scholar 

  56. 56.

    Li J, Zhang J, Li H, Niu H, Xu Q, Jiao Z, An J, Jiang Y, Li Q, Niu J. The major factors causing the microspore abortion of genic male sterile mutant NWMS1 in wheat (Triticum aestivum L.). Int J Mol Sci. 2019;20(24):6252.

  57. 57.

    Wang Y, Yu H, Tian C, Sajjad M, Gao C, Tong Y, Wang X, Jiao Y. Transcriptome association identifies regulators of wheat spike architecture. Plant Physiol. 2017;175(2):746–57.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  58. 58.

    Li C, Zhao Z, Liu Y, Liang B, Guan S, Lan H, Wang J, Lu Y, Cao M. Comparative transcriptome analysis of isonuclear-alloplasmic lines unmask key transcription factor genes and metabolic pathways involved in sterility of maize CMS-C. PeerJ. 2017;5:e3408.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  59. 59.

    Yang Y, Bao S, Zhou X, Liu J, Zhuang Y. The key genes and pathways related to male sterility of eggplant revealed by comparative transcriptome analysis. BMC Plant Biol. 2018;18(1):209.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  60. 60.

    Li Y, Qin T, Wei C, Sun J, Dong T, Zhou R, Chen Q, Wang Q. Using transcriptome analysis to screen for key genes and pathways related to cytoplasmic male sterility in cotton (Gossypium hirsutum L.). Int J Mol Sci. 2019;20(20):5120.

  61. 61.

    Zheng J, Kong X, Li B, Khan A, Li Z, Liu Y, Kang H, Ullah Dawar F, Zhou R. Comparative transcriptome analysis between a novel allohexaploid cotton progeny CMS line LD6A and its maintainer line LD6B. Int J Mol Sci. 2019;20(24):6127.

  62. 62.

    Kong X, Liu D, Liao X, Zheng J, Diao Y, Liu Y, Zhou R. Comparative analysis of the cytology and transcriptomes of the cytoplasmic male sterility line H276A and its maintainer line H276B of cotton (Gossypium barbadense L.). Int J Mol Sci. 2017;18(11):2240.

  63. 63.

    Hamid R, Tomar RS, Marashi H, Shafaroudi SM, Golakiya BA, Mohsenpour M. Transcriptome profiling and cataloging differential gene expression in floral buds of fertile and sterile lines of cotton (Gossypium hirsutum L.). Gene. 2018;660:80–91.

    CAS  PubMed  Article  Google Scholar 

  64. 64.

    Morant M, Jorgensen K, Schaller H, Pinot F, Moller BL, Werck-Reichhart D, Bak S. CYP703 is an ancient cytochrome P450 in land plants catalyzing in-chain hydroxylation of lauric acid to provide building blocks for sporopollenin synthesis in pollen. Plant Cell. 2007;19(5):1473–87.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  65. 65.

    Dobritsa AA, Shrestha J, Morant M, Pinot F, Matsuno M, Swanson R, Moller BL, Preuss D. CYP704B1 is a long-chain fatty acid omega-hydroxylase essential for sporopollenin synthesis in pollen of Arabidopsis. Plant Physiol. 2009;151(2):574–89.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  66. 66.

    Chen W, Yu XH, Zhang K, Shi J, De Oliveira S, Schreiber L, Shanklin J, Zhang D. Male Sterile2 encodes a plastid-localized fatty acyl carrier protein reductase required for pollen exine development in Arabidopsis. Plant Physiol. 2011;157(2):842–53.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  67. 67.

    Wilson ZA, Morroll SM, Dawson J, Swarup R, Tighe PJ. The Arabidopsis MALE STERILITY1 (MS1) gene is a transcriptional regulator of male gametogenesis, with homology to the PHD-finger family of transcription factors. Plant J. 2001;28(1):27–39.

    CAS  PubMed  Article  Google Scholar 

  68. 68.

    Sorensen AM, Krober S, Unte US, Huijser P, Dekker K, Saedler H. The Arabidopsis ABORTED MICROSPORES (AMS) gene encodes a MYC class transcription factor. Plant J. 2003;33(2):413–23.

    CAS  PubMed  Article  Google Scholar 

  69. 69.

    Li N, Zhang DS, Liu HS, Yin CS, Li XX, Liang WQ, Yuan Z, Xu B, Chu HW, Wang J, et al. The rice tapetum degeneration retardation gene is required for tapetum degradation and anther development. Plant Cell. 2006;18(11):2999–3014.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  70. 70.

    Li H, Yuan Z, Vizcay-Barrena G, Yang C, Liang W, Zong J, Wilson ZA, Zhang D. PERSISTENT TAPETAL CELL1 encodes a PHD-finger protein that is required for tapetal cell death and pollen development in rice. Plant Physiol. 2011;156(2):615–30.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  71. 71.

    Kobayashi K, Mochizuki N, Yoshimura N, Motohashi K, Hisabori T, Masuda T. Functional analysis of Arabidopsis thaliana isoforms of the mg-chelatase CHLI subunit. Photochem Photobiol Sci. 2008;7(10):1188–95.

    CAS  PubMed  Article  Google Scholar 

  72. 72.

    Pogson B, McDonald KA, Truong M, Britton G, DellaPenna D. Arabidopsis carotenoid mutants demonstrate that lutein is not essential for photosynthesis in higher plants. Plant Cell. 1996;8(9):1627–39.

    CAS  PubMed  PubMed Central  Google Scholar 

  73. 73.

    Foudree A, Aluru M, Rodermel S. PDS activity acts as a rheostat of retrograde signaling during early chloroplast biogenesis. Plant Signal Behav. 2010;5(12):1629–32.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  74. 74.

    Kim JE, Cheng KM, Craft NE, Hamberger B, Douglas CJ. Over-expression of Arabidopsis thaliana carotenoid hydroxylases individually and in combination with a beta-carotene ketolase provides insight into in vivo functions. Phytochemistry. 2010;71(2–3):168–78.

    PubMed  Article  CAS  Google Scholar 

  75. 75.

    Basnet R, Hussain N, Shu Q. OsDGD2beta is the sole digalactosyldiacylglycerol synthase gene highly expressed in anther, and its mutation confers male sterility in rice. Rice (N Y). 2019;12(1):66.

    Article  Google Scholar 

  76. 76.

    Basnet R, Zhang J, Hussain N, Shu Q. Characterization and mutational analysis of a monogalactosyldiacylglycerol synthase gene OsMGD2 in Rice. Front Plant Sci. 2019;10:992.

    PubMed  PubMed Central  Article  Google Scholar 

  77. 77.

    Berthome R, Thomasset M, Maene M, Bourgeois N, Froger N, Budar F. pur4 mutations are lethal to the male, but not the female, gametophyte and affect sporophyte development in Arabidopsis. Plant Physiol. 2008;147(2):650–60.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  78. 78.

    Mao G, Wei H, Hu W, Ma Q, Zhang M, Wang H, Yu S. Fine mapping and molecular characterization of the virescent gene vsp in upland cotton (Gossypium hirsutum). Theor Appl Genet. 2019;132(7):2069–86.

  79. 79.

    Kunz HH, Scharnewski M, Feussner K, Feussner I, Flugge UI, Fulda M, Gierth M. The ABC transporter PXA1 and peroxisomal beta-oxidation are vital for metabolism in mature leaves of Arabidopsis during extended darkness. Plant Cell. 2009;21(9):2733–49.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  80. 80.

    Xiong SX, Lu JY, Lou Y, Teng XD, Gu JN, Zhang C, Shi QS, Yang ZN, Zhu J. The transcription factors MS188 and AMS form a complex to activate the expression of CYP703A2 for sporopollenin biosynthesis in Arabidopsis thaliana. Plant J. 2016;88(6):936–46.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  81. 81.

    Millar AA, Gubler F. The Arabidopsis GAMYB-like genes, MYB33 and MYB65, are microRNA-regulated genes that redundantly facilitate anther development. Plant Cell. 2005;17(3):705–21.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  82. 82.

    Lee SB, Kim HU, Suh MC. MYB94 and MYB96 additively activate cuticular wax biosynthesis in Arabidopsis. Plant Cell Physiol. 2016;57(11):2300–11.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  83. 83.

    Zhou M, Zhang K, Sun Z, Yan M, Chen C, Zhang X, Tang Y, Wu Y. LNK1 and LNK2 corepressors interact with the MYB3 transcription factor in phenylpropanoid biosynthesis. Plant Physiol. 2017;174(3):1348–58.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  84. 84.

    Mehrtens F, Kranz H, Bednarek P, Weisshaar B. The Arabidopsis transcription factor MYB12 is a flavonol-specific regulator of phenylpropanoid biosynthesis. Plant Physiol. 2005;138(2):1083–96.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  85. 85.

    Zhou H, Lin-Wang K, Liao L, Gu C, Lu Z, Allan AC, Han Y. Peach MYB7 activates transcription of the proanthocyanidin pathway gene encoding leucoanthocyanidin reductase, but not anthocyanidin reductase. Front Plant Sci. 2015;6:908.

    PubMed  PubMed Central  Google Scholar 

  86. 86.

    Kosma DK, Bourdenx B, Bernard A, Parsons EP, Lu S, Joubes J, Jenks MA. The impact of water deficiency on leaf cuticle lipids of Arabidopsis. Plant Physiol. 2009;151(4):1918–29.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  87. 87.

    Hooker TS, Millar AA, Kunst L. Significance of the expression of the CER6 condensing enzyme for cuticular wax production in Arabidopsis. Plant Physiol. 2002;129(4):1568–80.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  88. 88.

    Li J, Yu M, Geng LL, Zhao J. The fasciclin-like arabinogalactan protein gene, FLA3, is involved in microspore development of Arabidopsis. Plant J. 2010;64(3):482–97.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  89. 89.

    Zheng Y, Yan J, Wang S, Xu M, Huang K, Chen G, Ding Y. Genome-wide identification of the pectate lyase-like (PLL) gene family and functional analysis of two PLL genes in rice. Mol Gen Genomics. 2018;293(6):1317–31.

    CAS  Article  Google Scholar 

  90. 90.

    Huang L, Cao J, Zhang A, Ye Y, Zhang Y, Liu T. The polygalacturonase gene BcMF2 from Brassica campestris is associated with intine development. J Exp Bot. 2009;60(1):301–13.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  91. 91.

    Huang L, Ye Y, Zhang Y, Zhang A, Liu T, Cao J. BcMF9, a novel polygalacturonase gene, is required for both Brassica campestris intine and exine formation. Ann Bot. 2009;104(7):1339–51.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  92. 92.

    Lin S, Dong H, Zhang F, Qiu L, Wang F, Cao J, Huang L. BcMF8, a putative arabinogalactan protein-encoding gene, contributes to pollen wall development, aperture formation and pollen tube growth in Brassica campestris. Ann Bot. 2014;113(5):777–88.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  93. 93.

    Lin S, Yue X, Miao Y, Yu Y, Dong H, Huang L, Cao J. The distinct functions of two classical arabinogalactan proteins BcMF8 and BcMF18 during pollen wall development in Brassica campestris. Plant J. 2018;94(1):60–76.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  94. 94.

    Suzuki T, Narciso JO, Zeng W, van de Meene A, Yasutomi M, Takemura S, Lampugnani ER, Doblin MS, Bacic A, Ishiguro S. KNS4/UPEX1: a type II arabinogalactan beta-(1,3)-galactosyltransferase required for pollen exine development. Plant Physiol. 2017;173(1):183–205.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  95. 95.

    Zhang H, Liang W, Yang X, Luo X, Jiang N, Ma H, Zhang D. Carbon starved anther encodes a MYB domain protein that regulates sugar partitioning required for rice pollen development. Plant Cell. 2010;22(3):672–89.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  96. 96.

    Ding J, Lu Q, Ouyang Y, Mao H, Zhang P, Yao J, Xu C, Li X, Xiao J, Zhang Q. 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.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  97. 97.

    Zhou H, Zhou M, Yang Y, Li J, Zhu L, Jiang D, Dong J, Liu Q, Gu L, Zhou L, et al. RNase Z(S1) processes UbL40 mRNAs and controls thermosensitive genic male sterility in rice. Nat Commun. 2014;5:4884.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  98. 98.

    Zhu J, Lou Y, Shi QS, Zhang S, Zhou WT, Yang J, Zhang C, Yao XZ, Xu T, Liu JL, et al. Slowing development restores the fertility of thermo-sensitive male-sterile plant lines. Nat Plants. 2020;6(4):360–7.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  99. 99.

    Hartmann U, Sagasser M, Mehrtens F, Stracke R, Weisshaar B. Differential combinatorial interactions of cis-acting elements recognized by R2R3-MYB, BZIP, and BHLH factors control light-responsive and tissue-specific activation of phenylpropanoid biosynthesis genes. Plant Mol Biol. 2005;57(2):155–71.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  100. 100.

    Karppinen K, Zoratti L, Sarala M, Carvalho E, Hirsimaki J, Mentula H, Martens S, Haggman H, Jaakola L. Carotenoid metabolism during bilberry (Vaccinium myrtillus L.) fruit development under different light conditions is regulated by biosynthesis and degradation. BMC Plant Biol. 2016;16:95.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  101. 101.

    Stephenson PG, Terry MJ. Light signalling pathways regulating the mg-chelatase branchpoint of chlorophyll synthesis during de-etiolation in Arabidopsis thaliana. Photochem Photobiol Sci. 2008;7(10):1243–52.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  102. 102.

    Mock HP, Trainotti L, Kruse E, Grimm B. Isolation, sequencing and expression of cDNA sequences encoding uroporphyrinogen decarboxylase from tobacco and barley. Plant Mol Biol. 1995;28(2):245–56.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  103. 103.

    Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9:559.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  104. 104.

    Qin Y, Sun H, Hao P, Wang H, Wang C, Ma L, Wei H, Yu S. Transcriptome analysis reveals differences in the mechanisms of fiber initiation and elongation between long- and short-fiber cotton (Gossypium hirsutum L.) lines. BMC Genomics. 2019;20(1):633.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  105. 105.

    Min L, Zhu L, Tu L, Deng F, Yuan D, Zhang X. Cotton GhCKI disrupts normal male reproduction by delaying tapetum programmed cell death via inactivating starch synthase. Plant J. 2013;75(5):823–35.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

Download references


We are grateful to all the colleagues in our laboratory who provided helpful advice and technical assistance.


This research was financed by the National Natural Science Foundation of China (31501344). The funding body was not involved in the design of the study, data collection and analysis, or interpretation of the manuscript.

Author information




SY and HLW designed the research program. HTW, MZ, LZ and HLW analyzed the data. JL, PC, LM, QM, XF and YQ revised the language and collected the data. MZ performed the experiment and wrote the manuscript. All authors have read and approved the final manuscript.

Corresponding authors

Correspondence to Hengling Wei or Shuxun Yu.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

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

Supplementary Information

Additional file 1:

Figure S1. Relative content of Cal a and Cal b in MT and WT leaves. (Student’s t-test, *P < 0.05, **P < 0.01). L1-L5, the first to the fifth leaf.

Additional file 2:

Figure S2. DEG analysis of male sterile line CCRI9106 compared to wild type at three anther development stages. (A) Number of DEGs that are up or down-regulated in the three development stages. (B) Venn diagrams showing the number of DEGs expressed over three stages.

Additional file 3:

Figure S3. GO enrichment analysis of down- and up-regulated clusters. Top 20 GO enrichment terms of down- and up-regulated DEGs are show in (A) and (B), respectively.

Additional file 4:

Figure S4. PPI network of cluster 3 DEGs. (A) GO modules enriched of PPI DEGs in cluster 3 visualized by the ClueGO plug-in in Cytoscape. Genes with blue font are MS genes reported in Arabidopsis. Trigonal node represents key genes in the network.

Additional file 5:

Table S1. Data Output Quality List.

Additional file 6:

Table S2. The DEGs classified by K-means clustering.

Additional file 7:

Table S3. TFs of down/up-regulated DEGs.

Additional file 8:

Table S4. DEGs in the PPI network 1/2/3.

Additional file 9:

Table S5. DEGs in yellow module of WGCNA.

Additional file 10:

Table S6. ABA-related DEGs in yellow module of WGCNA.

Additional file 11:

Table S7. Information of the primers used in qRT-PCR.

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 The Creative Commons Public Domain Dedication waiver ( 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

Verify currency and authenticity via CrossMark

Cite this article

Zhang, M., Liu, J., Ma, Q. et al. Deficiencies in the formation and regulation of anther cuticle and tryphine contribute to male sterility in cotton PGMS line. BMC Genomics 21, 825 (2020).

Download citation


  • Cotton
  • Photosensitive genetic male sterile (PGMS)
  • Exine
  • Tryphine
  • Anther cuticle
  • MYB TFs
  • ABA