Skip to main content


Insights into the cotton anther development through association analysis of transcriptomic and small RNA sequencing



Plant anther development is a systematic and complex process precisely controlled by genes. Regulation genes and their regulatory mechanisms for this process remain elusive. In contrast to numerous researches on anther development with respect to mRNAs or miRNAs in many crops, the association analysis combining both omics has not been reported on cotton anther.


In this study, the molecular mechanism of cotton anther development was investigated with the employment of association analysis of transcriptome and small RNA sequencing during the predefined four stages of cotton anther development, sporogenuous cell proliferation (SCP), meiotic phase (MP), microspore release period (MRP) and pollen maturity (PM). Analysis revealed that the differentially expressed genes are increasingly recruited along with the developmental progress. Expression of functional genes differed significantly among developmental stages. The genes related with cell cycle, progesterone-mediated oocyte maturation, and meiosis are predominantly expressed at the early stage of anther development (SCP and MP), and the expression of genes involved in energy metabolism, flavonoid biosynthesis, axon guidance and phospholipase D signaling pathways is mainly enriched at the late stage of anther development (MRP and PM). Analysis of expression patterns revealed that there was the largest number of differentially expressed genes in the MP and the expression profiles of differentially expressed genes were significantly increased, which implied the importance of MP in the entire anther development cycle. In addition, prediction and analysis of miRNA targeted genes suggested that miRNAs play important roles in anther development. The miRNAs ghr-miR393, Dt_chr12_6065 and At_chr9_3080 participated in cell cycle, carbohydrate metabolism and auxin anabolism through the target genes, respectively, to achieve the regulation of anther development.


Through the association analysis of mRNA and miRNA, our work gives a better understanding of the preferentially expressed genes and regulation in different developmental stages of cotton anther and the importance of meiotic phase, and also the involvement of miRNAs in precise regulation for this process, which would be valuable for clarifying the mechanism of plant anther development in response to internal and external environments.


Cotton (Gossypium hirsutum L.) is one of the most economically important crops due to its fiber used as the principal natural source for the textile industry worldwide. The development of flowers directly determines the formation of cotton yield. So far, studies on male gametophytes have focused on model plants, such as Arabidopsis, [1,2,3], but there are less research reports on cotton. Though the basic mechanisms of anther development could be cross-referenced, each species has its own peculiarity [4]. The occurrence of flowers to the maturity of pollen grains is a systematic and complex process involving the regulation of a large number of genes [5,6,7]. Reportedly, genes involving in anther development concentrated on regulation of starch and sucrose metabolism, carbohydrate metabolism, antioxidant production, flavonoid biosynthesis, cell cycle, meiosis and plant hormone pathways [4, 8, 9]. In early stages of pollen development, large amounts of sugars are directed to the anther to support its development [10]. At late stage of pollen development and maturation, starch is required as energy storage which contributes to the secondary cell wall thickening and the swelling of the pollen grains, and ultimately promotes the dehiscent of the anther [5, 11]. Therefore, the disturbances in sugar and starch synthesis, transportation and metabolism during this stage could severely impair pollen development, resulting in male sterility [12]. Compared with the wild type, male sterile mutant shows declination in starch and sucrose production [8, 13]. Further investigations revealed that genes, including Ghck1, SnRK1, hexokinase and 14–3-3 protein genes are involved in the metabolism of sucrose in starch biosynthesis. Differential expression of these genes can lead to male sterility [9, 13,14,15]. Mitochondria, chloroplasts and peroxisomes are of great importance for energy-related metabolism, and also the site of reactive oxygen species (ROS) production. Pollen and tapetum cells are known to harbor a great number of mitochondria, and exhibit high respiration during pollen development [16]. Flavonoids are free radical scavengers and also are components of pollen coat [17]. During anther development of upland cotton, many genes related to the biosynthesis of quercetin, kaenpferol and myricetin were found with higher expression than in other tissues [4]. Plant hormones are important metabolic regulator, and the occurrence of male sterility in plants is usually accompanied by changes in the content of endogenous hormones in various reproductive organs [18]. GA is an absolute requirement for flower initiation [19], which involved in the development of tapetal and programmed cell death [20]. The synthesis and accumulation of GA in anthers is mainly regulated by the GA20ox and GA3ox1 [21, 22]. The differential expression of GA20ox and GA3ox1 usually leads to the occurrence of male sterility [13, 22]. IAA is another essential plant hormone that regulates plant growth and development [23]. Low content of IAA in anthers is the main reason of male sterility, and the low expression of ILR1 and indole-3-acetic acid-amido consequently lead to low levels of IAA [13]. Similarly, the high level of IAA is also detrimental to anther development, which can increase the indehiscent anther ratio [9]. MicroRNAs (miRNAs) are a class of endogenous small regulatory non-coding RNA which negatively regulates expression of target genes at the post-transcriptional through degrading target mRNAs or repressing gene translation [24]. Resent researches showed that miRNAs were involved in multiple developmental processes, including seed germination [25], root development [26] and floral organ identity [27]. Studies showed that miR159 participated in the development of male gametophytes through target genes DUO1 [28], and overexpression of miR159 lead to male sterility [29]. The genes ARFs play important roles in floral development [30], overexpression of targeted miRNA160 in cotton can induce anther indehiscence [31].

So far, a few studies on the molecular mechanism of anther development were carried out through separated applications of miRNAs and transcriptomes [32, 33], yet the association analysis combining both omics has not been reported. In our study, the association analysis of mRNA and miRNA in cotton anther at four continuous development stages was carried out, and a better understanding of the genes expression and regulation for this process was achieved.


DEGs between different stages during anther development

Based on the statistical analysis of sequencing results, the 6.8 Gb raw bases data can be measured for each sample on average, and the ratio of clean bases is above 95%. The comparative analysis of cotton reference genome showed that the average mapped ratio of clean reads is 81.5% (Additional file 1: Table S1). The distance between samples is calculated by the clustering method, and the results showed that the two biological repetitions samples which cluster together preferentially at different stages (Additional file 2: Figure S1).

In order to characterize the gene expression in anther development, the differentially expressed genes (DEGs) were screened between the four continuous stages, and 2910, 7054 and 10,386 DEGs were identified when SCP compared with MP, MRP and PM, respectively; When MP compared with the MRP and PM, 5326 and 9664 DEGs were screened respectively, and 4973 DEGs were obtained from the comparison between MRP and PM (Fig. 1). The results indicate that the number of DEGs is gradually increased with the progress of anther development, suggesting that more specific genes are recruited at different developmental stages to ensure the normal development and maturation.

Fig. 1

Histograms of DEGs with different control. The first three groups of bar charts were the number of DEGs obtained using SCP as control; the fourth and fifth groups of bar charts were based on MP as the control, and they were compared with MRP and PM respectively; the last of bar chart is based on MRP as control and compared with PM

Functional analysis of DEGs

To further understand the function of DEGs, we performed Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis. The DEGs which were screened out with SCP as a control and compared with MP, MRP and PM, respectively, were mainly enriched in cell cycle, progesterone-mediated oocyte maturation and meiosis pathways (Fig. 2a and Additional file 3: Table S2). These pathways were necessary for early stage of anther development. The vast majority of DEGs which are obtained from MP compared with MRP and PM were involved in energy metabolic pathways, such as starch and sucrose metabolism, pentose and glucuronate interconversions and glycolysis/gluconeogenesis pathways (Fig. 2b and Additional file 3: Table S2); The DEGs between MRP and PM were mainly enriched in energy metabolic pathways, in addition to the flavonoid biosynthesis, axon guidance and phospholipase D signaling pathways (Fig. 2c and Additional file 3: Table S2).

Fig. 2

Heat maps of the KEGG pathway enrichment statistics. a Analysis of differentially expressed genes with KEGG enrichment using SCP as control; b Analysis of differentially expressed genes with KEGG enrichment using MP as control; c Analysis of differentially expressed genes with KEGG enrichment using MRP as control; P-values range from 0 to 1, and its less value means greater intensiveness. Top 10 pathway terms enriched by KEGG database are dispalied

About 29 metabolic and signal pathways, which related with the development of pollen and anther [4, 34], such as cell cycle, meiosis, glutathione metabolism and homologous recombination were selected for further enrichment analysis. Among the DEGs which came from individual comparison between SCP as a control and other three stages, and 281, 688 and 1054 DEGs were enriched respectively in the above-mentioned 29 pathways. A total of 106 genes were shared by three sets of DEGs (Fig. 3a and Additional file 4: Table S3) and These genes were clustered into the eleven pathways, such as cell cycle, meiosis, starch and sucrose metabolism and so on. The pathway which enriched the largest number of DEGs is cell cycle, followed by the meiosis with 30 and 18 DEGs, respectively. This suggested the importance of cell cycle and meiosis during the whole anther development process [1]. Among the 106 DEGs, the expression profiles of genes enriched in pathways like cell cycle, homologous recombination, meiosis, flavonoid biosynthesis and peroxisome exhibited downward trend with the development of anther (Fig. 3b, c and d; Additional file 4: Table S3; Additional file 5: Figure S2). While the genes participated in carbohydrate digestion and absorption and starch and sucrose metabolism showed the upward expression trend continuously (Fig. 3e; Additional file 4: Table S3; Additional file 5: Figure S2), which serves as a marker of pollen maturity and provides an energy reserve for pollen germination [9, 11]. The genes involved in carbon metabolism and plant hormone signal transduction had no consistent expression trend (Fig. 3f and Additional file 4: Table S3).

Fig. 3

Expression profiles of genes related to metabolism and signaling pathways for anther development. a Venn diagram of anther development-related genes in three differentially expressed genes obtained from SCP as control; b-d The number of differential genes enriched in the cell cycle, meiosis and homologous recombination pathways, respectively; e Starch and sucrose metabolism was enriched in seven differentially genes; f Plant hormone signal transduction pathway was enriched five differential genes, and different hormone synthesis related genes showed different expression trends. The horizontal axis represents stages, and the vertical axis shows the stages series of gene expression for the gene after Log normalized transformation

Expression pattern of DEGs varied during anther development

Three sets of DEGs are obtained with SCP as a control and compared with MP, MRP and PM, respectively, and eight main kinds of expression pattern are analyzed (Fig. 4a and b; Additional file 6: Table S4).Cluster 8 comprised the largest number of genes (78 genes), followed by cluster 4 (35) and cluster 1 (30). The expression profile of DEGs which contained in cluster 8 reached the highest in the MP (Fig. 4a and Additional file 6: Table S4), indicating that meiosis normally proceeded and required a large number of highly expressed genes [35]. Gene Ontology (GO) enrichment found that DEGs were mainly enriched in the oxidation-reduction process, transport, membrane, carbohydrate metabolic process and ATP binding terms (Fig. 4c).

Fig. 4

Gene expression patterns of the DEGs. a The heat map shows the K-means clustering of transcription levels of the DEGs. Genes were clustered into eight groups with distinct expression patterns in four stages. b DEGs Venn diagram. Eight groups expression patterns displayed on Venn diagram. c The top five most significantly enriched GO terms in each clusters; Log10 (Num) range from 0 to 1, and its larger value means greater

In cluster 4, the genes expression reached the highest in the MP and MRP (Fig. 4a and Additional file 6: Table S4), which was consistent with the fact that MP require a large number of specific genes to participate. The GO enrichment results were mostly the same as cluster 8, and with part of DEGs enriched in the hydrolase activity and protein binding terms (Fig. 4c). The DEGs enriched in cluster 1 had the highest expression at the SCP, and then decreased with the process of anther development (Fig. 4a and Additional file 6: Table S4), GO annotation indicated that these genes related with DNA binding, nucleosome assembly and DNA replication were necessary for the proliferation of pollen mother cells (Fig. 4c). The remaining five clusters were enriched with a small number of differentially expressed genes and there was no significant GO enrichment.

Different expression of miRNAs during anther development

By deeply small RNA sequencing anther samples at four stages, an average of 21.1 M clean reads was obtained for each sample. The obtained clean reads were aligned to the reference genome and the average alignment rate was found to be 93.86% (Additional file 7: Table S5). The clean reads were homologous alignment in the miRBasa, meanwhile, new miRNA predictions were performed for reads that failed to undergo homology comparisons. It was found that 97 miRNAs were detected from all samples, including 56 known miRNAs and 41 new miRNAs, respectively.

When SCP compared with the MP, MRP and PM, respectively, 2, 14, and 13 differentially expressed miRNAs were obtained (Fig. 5). This trend was similar to the number of DEGs, presumably because more different genes were involved in the regulation of anther development, thus more miRNAs were required to regulate the mRNA expression [32]. When MP was used as a control and compared with the MRP and PM, 9 and 12 differentially expressed miRNAs were screened out (Fig. 5), respectively, and two difference expression miRNAs were obtained from the comparison between MRP and PM (Fig. 5).

Fig. 5

Histograms of differentially expressed miRNAs with different control. The first three groups of bar charts were the number of differentially expressed miRNAs obtained using SCP as control; the fourth and fifth groups of bar charts were based on MP as the control, and they were compared with MRP and PM respectively; the last of bar chart is based on MRP as control and compared with PM

Target gene prediction of differentially expressed miRNAs

The target genes of different expression miRNAs were predicted in the DEGs of same phase of transcriptome sequencing and enrichment on anther development related pathways. In the comparison of SCP and MRP, there were a total of nine target genes which were predicted for seven miRNAs. In the comparation of SCP and PM, eight out of 13 differentially expressed miRNAs predicted 13 target genes in total (Additional file 8: Table S6; Additional file 9: Table S7). However, no target gene was predicted when comparing with SCP and MP.

KEGG enrichment analysis of the predicted target genes found that target genes were mainly enriched in pathways like cell cycle, glutathione metabolism and plant hormone signal transduction (Additional file 10: Figure S3; Additional file 11: Figure S4), suggesting the importance of corresponding miRNA in pollen development [30, 36]. GO clustering analysis found that target genes mainly clustered into molecular function (MF) and biological process (BP), and few genes in cell component (CC).

Among the predicted target genes, there were three sharing target genes, termed CotAD_22690 (TIR1), CotAD_67461 (MPS1) and CotAD_13956 (GSTU24), and the corresponding miRNAs were At_chr9_3038, ghr-miR393 and ghr-miR396a, respectively. The expressions analysis showed that all three miRNAs and GSTU24 were up-regulated in the whole process of anther development, while the TIR1 and MPS1 were always down-regulated (Fig. 6). The interaction of these three target genes was predicted with string website (Fig. 7). As TIR1 is auxin receptor and mediates Aux/IAA proteins proteasomal degradation, it had been verified interaction with several IAAs, but its interaction with GSTU24 remained putative [37]. MPS1 had no interaction with the other two proteins.

Fig. 6

Gene expression analysis. qRT-PCR analysis of mRNAs (TIR1, MPS1, GSTU24, IAA14, AUX1, PYL) and miRNAs (At_chr9_3038, ghr-miR393, ghr-miR396a) expression at four developmental stages. Error bars represent SD, which were calculated from three technical repeats. Significant differences between treatments are indicated by letters above each bar (P ≤ 0.05). Data are present as means ± SE (n = 5)

Fig. 7

Protein-protein interaction (PPI) network. Differential miRNAs obtained with SCP as control group, three common target genes TIR1, MPS1 and GSTU24 PPI network illustrated by STRING (high confidence score = 0.7). Among the predicted genes, only IAA14 was significant differently expressed at different developmental stages

The other genes were not found expressed in our transcriptome sequencing, but the IAA14 was up-regulated. IAA14 functions as repressors of early auxin response at low auxin concentrations [38] and TIR1 functions as promoting the degradation of IAA14 [31, 39], which could explain the expression of TIR1 decreased while the expression of IAA14 increased during anther development.


Metabolism and signal pathways in anther development

The comparative analyses among four stages showed that more DEGs existed between distant stages (Fig. 1), which indicate that more specific genes started to participate in the regulation with the process of anther development [40, 41].

KEGG enrichment analysis of the DEGs, which were obtained with the SCP as control (Fig. 2a), the results indicate that anthers made material reserves for the following meiosis and ensured the transition from SCP to MP. In the comparison of MP as control (Fig. 2b), the results suggest that pollen mother cells had reached or been close to meiosis, and beginning to retain material and energy for subsequent pollen maturation [8]. In the comparison between MRP and PM (Fig. 2c), the results implied that in later stage of pollen development, beside continuing energy substance reserving, microspores also provided the basis for the subsequent pollination and fertilization of reserve materials [5, 8, 42, 43]. Flavonoids in PM stage are not only the components of pollen coat but also free radical scavengers, therefore, the up-regulation of flavonoids biosynthesis could be contributive to prevent or reduce pollen’s suffering the adversity environment in the pollination process, especially the abiotic stresses [44]. Phosphorylation signals contribute to the identification of pollen grains and stigmas [42]. The axon guidance allows the pollen tube to grow in the right direction, guiding the pollen tube to the location of the pearl hole so that the fertilization process can be successfully completed [8, 43].

Different gene expression pattern during anther development

The DEGs expression analysis showed that the genes which enriched in pathways like cell cycle, homologous recombination and meiosis, were down-regulated (Fig. 3b, c, d). In the SCP, these genes were the highest expression and reached the lowest in the PM, and the expression trend gradually decreased along with the process of anther development, which indicated that the three pathways played a decisive role in early stages of anther development and the importance of regulation gradually weakened along with the anther development [1, 45]. However, the expression of genes which were enriched in starch and sucrose metabolism and carbohydrate digestion and absorption pathways increased gradually along with the process of anther development (Fig. 3e and Additional file 5: Figure S2F), which indicated that in later phases of anther development the energy metabolism process was mainly carried out in the microsporocyte [5, 46].

There was no consistent expression profiling in the DEGs which enriched in plant hormone signal transduction and carbon metabolism pathways (Fig. 3f and Additional file 5: Figure S2E). The expression of auxin synthesis gene AUX1 (CotAD_18928) slowly increased from SCP to MP, and rapidly increased in the MRP, and then quickly decreased in the PM (Fig. 6). The expression profiling indicated that AUX1 played a crucial role in the MRP [9]. However, the abscisic acid synthesis gene PYL (CotAD_22336) showed stably up-regulated expression during the whole anther development (Fig. 6). These different expression patterns indicated that the PYL might play an important but different role with AUX1 during the whole anther development [47, 48].

Meiosis as the most important stage in pollen development process

The expression pattern analysis found that cluster 8 had enriched the most DEGs, followed by cluster 4 and cluster 1 in turn (Fig. 4a and Additional file 6: Table S4). Cluster 8 and cluster 4 were in the meiosis stage, where the DEGs reached the highest, indicating that a large number of specific genes were involved in MP and the complexity of gene regulation was existedduring meiosis. Therefore, it could be considered that MP was the most important stage in the whole anther development period [35]. In cluster 1, DEGs reached to the highest in the SCP, which indicating that the proliferation of pollen mother cells also requires a large number of different genes to participate in. SCP as an early stage of MP, a lot of specific genes are speculated to begin to be expressed in pollen mother cells, which reserved material and energy for the following meiosis to ensure the smooth progress of meiosis and promote the pollen development process [1]. This conjecture can be confirmed from the SCP as control and compared with the MP, MRP and PM, respectively, the DEGs mostly clustered in the pathways related to maturation and meiosis of oocytes.

GO enrichment analysis revealed that DEGs in cluster 8 were mainly enriched in oxidation-reduction process, transport, membrane, carbohydrate metabolic process and ATP binding terms (Fig. 4c). Presumably, pollen mother cells were required to obtain a lot of substances to satisfy the meiosis process. Such course is an intensive ATP consuming and reactive oxygen species (ROS)-producing process causing pollen mother cells to mobilize the active oxygen scavenging mechanism to prevent damaged by oxidation. In cluster 4, a part of DEGs enriched in the hydrolase activity and protein binding terms (Fig. 4c), which regulate the degradation of spindles and cytoskeletons restructuring. The DEGs GO enrichment in cluster 1 revealed that the DEGs were mainly enriched in terms of DNA binding, nucleosome assembly and DNA replication, which was consistent with the proliferation of pollen mother cells in the SCP through mitosis (Fig. 4c).

miRNA involved in regulation of anther development

MiRNA regulation of the target gene was mainly achieved by complementary pairing with the protein-coding region of the target gene and cutting the target gene or suppressing the translation of the target [49].

MPS1 (CotAD_67461), which regulated cell cycle function during anther development, was regulated by ghr-miR393 at gene expression. The MPS1 showed a sustained up-regulated expression and the ghr-miR393 showed an opposite trend (Fig. 8). Sequence alignment analysis in miRBase found that there was one base difference between the sequence of the ath-miR393 and ghr-miR393. In Arabidopsis thaliana, mTIR1 could promote the expression of ath-miR393, resulting in flowering delay [50]. The same trend was found between miRNA Dt_chr12_6065 and the target gene PME40 (CotAD_16305), this gene was associated with starch and sucrose metabolism (Fig. 8). Auxin is necessary for anther development. At different developmental stages, anther has strict requirements on auxin concentration, high or low auxin concentrations may lead to abnormal development of anthers. Studies have shown that too low concentrations of auxin can lead to male sterility [51], while excessive concentrations of auxin can cause anther indehiscence, which affecting the pollination and fertilization process [31].

Fig. 8

Real-time quantitative analysis of miRNAs and target genes. qRT-PCR analysis of the expression of miRNAs (ghr-miR393, At_chr9_3038 and Dt_chr12_6065) and their corresponding target genes (CotAD_67461, CotAD_22690 and CotAD_16305). Significant differences between treatments are indicated by letters above each bar (P ≤ 0.05). Data are present as means ± SE (n = 5)

Auxin signaling primarily depended on four partially redundant F-box proteins of the TIR1/AFB2 auxin receptor (TAAR) clade to trigger the degradation of AUX/IAA transcriptional repressors [52]. TIR1 is the target gene of At_chr9_3038 (the sequence of At_chr9_3038 is same as the ath-miR393), the expression of At_chr9_3038 increased gradually with the anther development process and the opposite expression profile of TIR1 (Fig. 8) suggesting that in the early stages of anther development, a higher concentration of auxin was required, and the auxin concentration decreased gradually with the anther development [51].

IAA14 inhibits the expression of auxin responsive genes [31]. In our study, the expression profiling of IAA14 was similar to that of At_chr9_3038 (Fig. 6), showing a gradual increase to inhibit the function of auxin in the later stages of anther development. It is speculated that during the anther development, a variety of ways of coordinated regulation were adopted to regulate the concentration of auxin, in order to reach the appropriate concentration of auxin in different stages of anther development [31, 53, 54].

Because of the stable and conserved characteristics of miRNAs, three miRNAs which could be used for the judgment of anther development stage according to the expression relationship were screened and confirmed by real-time quantitative RT-PCR. With the development of cotton anther, the order of expression amount is ghr-miR396a > Dt_chr7_8575 > ghr-miR2948-5p in SCP, Dt_chr7_8575 > ghr-miR396a > ghr-miR2948-5p in MP, ghr-miR2948-5p > Dt_chr7_8575 > ghr-miR396a in MRP, and Dt_chr7_8575 > ghr-miR2948-5p > ghr-miR396a in PM (Additional file 12: Figure S5). The relative expression amount of the three miRNAs provides a new method for the stages division of cotton anthers, and also provides a molecular basis at different period, but this reason and mechanism need to be further investigated.


The results of this study reveal that different genes are involved in the regulation of anthers at different developmental stages and at the same time promote the development of anthers. The gene expression profile analysis suggests that the importance of meiosis in the entire anther development, this conclusion can explain the MP is most susceptible to environmental influences, resulting in the occurrence of male sterility. Through the prediction of miRNA target genes, it was found that miRNAs participate in the regulation of anther development. Especially in the precise regulation of plant hormone auxin during the process of cotton anther development, miRNAs play a very important role. The results will enhance our understanding of the genes expression and regulation during cotton anther development and would be valuable for clarifying the mechanism of plant anther development in response to internal and external environments.


Plant materials

The genetically standard line of upland cotton was planted in experimental field of Hunan Agricultural University [55]. According to our earlier studies [56], the sporogenuous cell proliferation (SCP), meiotic phase (MP), microspore release period (MRP) and pollen maturity (PM) four stages buds were harvested, respectively, in the morning at the peak flowering period. Anthers collected from buds at four stages and then stored at freezer used to the transcriptomic and small RNA sequencing. The bud, stamen and anther size at four stages showed in Additional file 13: Figures S6A, S6B and S6C, and the each stage of cellular state were confirmed under optical microscope (Additional file 13: Figure S6D).

RNA extraction and sequencing

Twenty buds were collected for each period to separate anther tissue. Approximately 20 μg of total RNA from each sample was extracted from isolated anthers using TRIzol reagent kit (Invitrogen, Carlsbad, CA, US) according to the manufacturer’s specification. The yield of RNA was determined using a NanoDrop 2100 spectrophotometer (Thermo Scientific, USA), and the integrity was evaluated using agarose gel electrophoresis stained with ethidium bromide. The mRNA were enriched by magnetic beads coated with Oligo (dT), and then randomly fragmented by ultrasound. Using the mRNA fragments, first and second strand cDNAs were synthesized with six-base primers. Double stranded cDNAs were repaired and connected to sequencing heads. The total RNA was subjected to gel electrophoresis; fragments with a certain size were recovered and 3’and 5′-end linkers were added; the platform of Illumina HiSeq 2500 was adopted to have small RNA sequencing (transcriptome and small RNA sequencing were both performed by OE biotech company (Shanghai). The sequencing data are deposited in NCBI Sequence Read Archive (SRA, with accession number SRP143909.

Data analysis

Two biological repetitions were taken at each stage, respectively named as TC1–1 and TC1–2(SCP), TC2–1 and TC2–2 (MP), TC3–1 and TC3–2(MRP), TC4–1 and TC4–2(PM). The reference genome and transcriptome link of upland cotton were The comparative analysis of cotton reference genome is made for clean bases [57]. Gene expression was calculated and subsequently normalized to RPKM [58]. The uniform screening conditions for differential mRNA and differential microRNA (miRNA) were p < =0.05, and fold change > 2 or fold change < 0.5. The Bowtie software is adopted to compare the sequence obtained by sequencing with the miRNA mature body sequence in miRBase, which is considered to be a known miRNA [59]. Combined with homologous miRNA sequence of cotton and the RNA secondary structure prediction software like RNA fold, the miRDeep2 software is adopted to predict new mature body. The screened differentially expressed genes were made GO and KEGG analyses [60]. The K-means clustering algorithm was adopted to analyze the expression pattern, and the on-line prediction website ( was used to perform the target gene prediction, and the prediction parameters: Maximum Entropy = 4.0; Target accessibility = 50. The String website (https: //string. was used to predict protein interactions.

Real-time quantitative RT-PCR

The residual RNA of sequencing was used for real-time quantitative RT-PCR. The primer pairs (Additional file 14: Table S8) used for real-time quantitative RT-PCR were designed using Roche LCPDS2 software and synthesized by Generay Biotech (Generay, PRC).

The mRNAs real-time quantitative RT-PCR method was according to previous studies [61].The microRNAs quantification was performed with a two-step reaction process: reverse transcription (RT) and PCR. Each RT reaction consisted of 0.5 μg RNA, 2 μl of miScript HiSpec Buffer, 1 μl of Nucleics Mix and 0.5 μl of miScript Reverse Transcriptase Mix (Qiagen, Germany), in a total volume of 10 μl. Reactions were performed in a GeneAmp® PCR System 9700 (Applied Biosystems, USA) for 60 min at 37 °C, followed by heat inactivation of RT for 5 min at 95 °C. The 10 μl RT reaction mix was then diluted × 10 in nuclease-free water and held at − 20 °C. Real-time PCR was performed using LightCycler® 480 II Real-time PCR Instrument (Roche, Swiss) with 10 μl PCR reaction mixture that included 1 μl of cDNA, 5 μl of 2 × QuantiFast® SYBR® Green PCR Master Mix (Qiagen, Germany), 0.2 μl of universal primer (Qiagen, Germany), 0.2 μl of microRNA-specific primer and 3.6 μl of nuclease-free water. Reactions were incubated in a 384-well optical plate (Roche, Swiss) at 95 °C for 5 min, followed by 40 cycles of 95 °C for 10 s, 60 °C for 30 s. Each sample was run in triplicate for analysis. At the end of the PCR cycles, melting curve analysis was performed to validate the specific generation of the expected PCR product.

The expression levels of mRNAs and miRNAs were normalized to Actin (NM_001327051.1) and 5.8S gene, respectively. The expressin amounts were calculated using the 2-ΔΔCt methd [62].



Differentially expressed gene


Gene ontology


Kyoto encyclopedia of genes and genomes


Meiotic phase


Microspore release period


Pollen maturity


Reactive oxygen species


Sporogenuous cell proliferation


  1. 1.

    Rutley N, Twell D. A decade of pollen transcriptomics. Plant Reprod. 2015;28(2):73–89.

  2. 2.

    Guo SY, Sun B, Looi LS, Xu YF, Gan ES, Huang JB, Ito T. Co-ordination of flower development through epigenetic regulation in two model species: Rice and Arabidopsis. Plant Cell Physiol. 2015;56(5):830–42.

  3. 3.

    Jeong HJ, Kang JH, Zhao MA, Kwon JK, Choi HS, Bae JH, Lee HA, Joung YH, Choi D, Kang BC. Tomato male sterile 1035 is essential for pollen development and meiosis in anthers. J Exp Bot. 2014;65(22):6693–709.

  4. 4.

    Ma JH, Wei HL, Song MZ, Pang CY, Liu J, Wang L, Zhang JF, Fan SL, Yu SX. Transcriptome profiling analysis reveals that flavonoid and ascorbate-glutathione cycle are important during anther development in upland cotton. PLoS One. 2012;7(11):e49244.

  5. 5.

    Wilson ZA, Song J, Taylor B, Yang C. The final split: the regulation of anther dehiscence. J Exp Bot. 2011;62(5):1633–49.

  6. 6.

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

  7. 7.

    Ma J, Skibbe DS, Fernandes J, Walbot V. Male reproductive development: gene expression profiling of maize anther and pollen ontogeny. Genome Biol. 2008;9(12):R181.

  8. 8.

    Wei MM, Song MZ, Fan SL, Yu SX. Transcriptomic analysis of differentially expressed genes during anther development in genetic male sterile and wild type cotton by digital gene-expression profiling. BMC Genomics. 2013;14:97.

  9. 9.

    Min L, Li YY, Hu Q, Zhu LF, Gao WH, Wu YL, Ding YH, Liu SM, Yang XY, Zhang XL. Sugar and auxin signaling pathways respond to high-temperature stress during anther development as revealed by transcript profiling analysis in cotton. Plant Physiol. 2014;164(3):1293–308.

  10. 10.

    Oliver SN, Van DJT, Alfred SC, Mamun EA, Zhao X, Saini HS, Fernandes SF, Bianchard CI, Sutton BG, Geigenberger P, Dennis ES, Dolferus R. Cold-induced repression of the rice anther-specific cell wall invertase gene OSINV4 is correlated with sucrose accumulation and pollen sterility. Plant Cell Environ. 2005;28:1534–51.

  11. 11.

    Datta R, Chamusco KC, Chourey PS. Starch biosynthesis during pollen maturation is associated with altered patterns of gene expression in maize. Plant Physiol. 2002;130(4):1645–56.

  12. 12.

    Mamun EA, Alfred S, Cantrill LC, Overall RL, Sutton BG. Effects of chilling on male gametophyte development in rice. Cell Biol Int. 2006;30(7):583–91.

  13. 13.

    Fang WP, Zhao FA, Sun Y, Xie DY, Sun L, Xu ZZ, Zhu W, Yang LR, Zhao YM, Lv SP, Tang ZJ, Nie LH, Li W, Hou JN, Duan ZZ, Yu YB, Yang XJ. Transcriptomic profiling reveals complex molecular regulation in cotton genic male sterile mutant Yu98-8A. PLoS One. 2015;10(9):e0133425.

  14. 14.

    Zhang Y, Shewry PR, Jones H, Barcelo P, Lazzeri PA, Halford NG. Expression of antisense SnRK1 protein kinase sequence causes abnormal pollen development and male sterility in transgenic barley. Plant J. 2001;28(4):431–41.

  15. 15.

    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.

  16. 16.

    Selinski J, Scheibe R. Pollen tube growth: where does the energy come from? Plant Signal Behav. 2014;9(12):e977200.

  17. 17.

    Hsieh K, Huang AHC. Tapetosomes in Brassica tapetum accumulate endoplasmic reticulum-derived flavonoids and alkanes for delivery to the pollen surface. Plant Cell. 2007;19(2):582–96.

  18. 18.

    Shukla A, Sawhney VK. Abscisic acid: one of the factors affecting male sterility in Brassi canapus. Physiol Plant. 1994;91:522–8.

  19. 19.

    King RW, Moritz T, Evans LT, Martin J, Andersen CH, Blundell C, Kardailsky L, Chandler PM. Regulation of flowering in the long-day grass Lolium temulentum by gibberellins and the FLOWERING LOCUS T gene. Plant Physiol. 2006;141(2):498–507.

  20. 20.

    Aya K, Ueguchi-Tanaka M, Kondo M, Hamada K, Yano K, Nishimura M, Matsuoka M. Gibberellin modulates anther development in rice via the transcriptional regulation of GAMYB. Plant Cell. 2009;21(5):1453–72.

  21. 21.

    Lee DJ, Zeevaart JAD. Regulation of gibberellin 20-oxidase1 expression in spinach by photoperiod. Planta. 2007;226(1):35–44.

  22. 22.

    Sakata T, Oda S, Tsunaga Y, Shomura H, Kawagishi-Kobayashi M, Aya K, Saeki K, Endo T, Nagano K, Kojima M, Sakakibara H, Watanabe M, Matsuoka M, Higashitani A. Reduction of gibberellin by low temperature disrupts pollen development in rice. Plant Physiol. 2014;164(4):2011–9.

  23. 23.

    Zhao Y. Auxin biosynthesis and its role in plant development. Annu Rev Plant Biol. 2010;61:49–64.

  24. 24.

    Jones-Rhoades MW, Bartel DP, Bartel B. MicroRNAs and their regulatory roles in plants. Annu Rev Plant Biol. 2006;57:19–53.

  25. 25.

    Reyes JL, Chua NH. ABA induction of miR159 controls transcript levels of two MYB factors during Arabidopsis seed germination. Plant J. 2007;49(4):592–606.

  26. 26.

    Gutierrez L, Bussell JD, Pacurar DI, Schwambach J, Pacurar M, Bellini C. Phenotypic plasticity of adventitious rooting in Arabidopsis is controlled by complex regulation of AUXIN RESPONSE FACTOR transcripts and microRNA abundance. Plant Cell. 2009;21(10):3119–32.

  27. 27.

    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.

  28. 28.

    Zheng BL, Chen XM, McCormick S. The anaphase-promoting complex is a dual integrator that regulates both microRNA-mediated transcriptional regulation of cyclin B1 and degradation of cyclin B1 during Arabidopsis male gametophyte development. Plant Cell. 2011;23(3):1033–46.

  29. 29.

    Achard P, Herr A, Baulcombe DC, Harberd NP. Modulation of floral development by a gibberellin-regulated microRNA. Development. 2004;131(14):3357–65.

  30. 30.

    Yang XJ, Zhao YM, Xie DY, Sun Y, Zhu XL, Esmaeili N, Yang ZR, Wang Y, Lv SP, Nie LH, Tang ZJ, Zhao FA, Li W, Mishra N, Sun L, Zhu W, Fang WP. Identification and function analysis of microRNAs involved in anther development in cotton genic male sterile line Yu98-8A. Int J Mol Sci. 2016;17(10):1677.

  31. 31.

    Ding Y, Ma Y, Liu J, Xu J, Hu Q, Li YY, Wu YL, Xie S, Zhu LF, Min L, Zhang XL. microRNAs involved in auxin signaling modulate male sterility under high-temperature stress in cotton (Gossypium hirsutum). Plant J. 2017;91(6):977–94.

  32. 32.

    Seifert F, Bossow S, Kumlehn J, Gnad H, Scholten S. Analysis of wheat microspore embryogenesis induction by transcriptome and small RNA sequencing using the highly responsive cultivar “Svilena”. BMC Plant Biol. 2016;16:97.

  33. 33.

    Wei X, Zhang X, Yao Q, Yuan Y, Li X, Wei F, Zhao Y, Zhang Q, Wang Z, Jiang W, Zhang X. The miRNAs and their regulatory networks responsible for pollen abortion in Ogura-CMS Chinese cabbage revealed by high-throughput sequencing of miRNAs, degradomes, and transcriptomes. Front Plant Sci. 2015;6:894.

  34. 34.

    Liu J, Pang CY, Wei HL, Song MZ, Meng YY, Ma JH, Fan SL, Yu SX. iTRAQ-facilitated proteomic profiling of anthers from a photosensitive male sterile mutant and wild-type cotton (Gossypium hirsutum L.). J. Proteomics. 2015;126:68–81.

  35. 35.

    Jin L, Neiman AM. Post-transcriptional regulation in budding yeast meiosis. Curr Genet. 2016;62(2):313–5.

  36. 36.

    Zhang D, Yang L. Specification of tapetum and microsporocyte cells within the anther. Curr Opin Plant Biol. 2014;17:49–55.

  37. 37.

    Su Y, Luo W, Chen X, Liu H, Hu Y, Lin W, Xiao L. Auxin extraction and purification based on recombinant Aux/IAA proteins. Biol Proced Online. 2017;19:1.

  38. 38.

    Kepinski S, Leyser O. Auxin-induced SCFTIR1-Aux/IAA interaction involves stable modification of the SCFTIR1 complex. Proc Natl Acad Sci U S A. 2004;101(33):12381–6.

  39. 39.

    Paque S, Weijers D. Q&A: Auxin: the plant molecule that influences almost anything. BMC Biol. 2016;14(1):67.

  40. 40.

    Twell D, Brownfield L. Analysis of fluorescent reporter activity in the male germline during pollen development by confocal microscopy. Methods Mol Biol. 2017;1669:67–75.

  41. 41.

    Twell D. Male gametogenesis and germline specification in flowering plants. Sex Plant Reprod. 2011;24(2):149–60.

  42. 42.

    Lord EM, Russell SD. The mechanisms of pollination and fertilization in plants. Annu Rev Cell Dev Biol. 2002;18:81–105.

  43. 43.

    Gilles LM, Khaled A, Laffaire JB, Chaignon S, Gendrot G, Laplaige J, Bergès H, Beydon G, Bayle V, Barret P, Comadran J, Martinant JP, Rogowsky PM, Widiez T. Loss of pollen-specific phospholipase NOT LIKE DAD triggers gynogenesis in maize. EMBO J. 2017;36(6):707–17.

  44. 44.

    Li Q, Yu HM, Meng XF, Lin YJ, Hou BK. Ectopic expression of glycosyltransferase UGT76E11 increases flavonoid accumulation and enhances abiotic stress tolerance in Arabidopsis. Plant Biol. 2017;20(1):10–9.

  45. 45.

    Zhang DB, Luo X, Zhu L. Cytological analysis and genetic control of rice anther development. J Genet Genomics. 2011;38(9):379–90.

  46. 46.

    Hirose T, Hashida Y, Aoki N, Okamura M, Yonekura M, Ohto C, Terao T, Ohsugi R. Analysis of gene-disruption mutants of a sucrose phosphate synthase gene in rice, OsSPS1, shows the importance of sucrose synthesis in pollen germination. Plant Sci. 2014;225:102–6.

  47. 47.

    Storme ND, Geelen D. The impact of environmental stress on male reproductive development in plants: biological processes and molecular mechanisms. Plant Cell Environ. 2014;37(1):1–18.

  48. 48.

    Wang CS, Hsu SW, Hsu YF. New insights into desiccation-associated gene regulation by Lilium longiflorum ASR during pollen maturation and in transgenic Arabidopsis. Int Rev Cell Mol Biol. 2013;301:37–94.

  49. 49.

    Khraiwesh B, Zhu JK, Zhu JH. Role of miRNAs and siRNAs in biotic and abiotic stress responses of plants. Biochim Biophys Acta. 2011;1819(2):137–48.

  50. 50.

    Chen ZH, Bao ML, Sun YZ, Yang YJ, Xu XH, Wang JH, Han N, Bian HW, Zhu MY. Regulation of auxin response by miR393-targeted transport inhibitor response protein 1 is involved in normal development in Arabidopsis. Plant Mol Biol. 2011;77(6):619–29.

  51. 51.

    Cheng YF, Dai XH, Zhao YD. Auxin biosynthesis by the YUCCA flavin monooxygenases controls the formation of floral organs and vascular tissues in Arabidopsis. Genes Dev. 2006;20(13):1790–9.

  52. 52.

    Windels D, Bielewicz D, Ebneter M, Jarmolowski A, Szweykowska-Kulinska Z, Vazquez F. miR393 is required for production of proper auxin signaling outputs. PLoS One. 2014;9(4):e95972.

  53. 53.

    Baranwal VK, Negi N, Khurana P. Auxin response factor genes repertoire in mulberry: identification, and structural, functional and evolutionary analyses. Genes (Basel). 2017;8(9):202.

  54. 54.

    Yan Z, Liu X, Ljung K, Li S, Zhao W, Yang F, Wang M, Tao Y. Type B response regulators act as central integrators in transcriptional control of the auxin biosynthesis enzyme TAA1. Plant Physiol. 2017;175(3):1438–54.

  55. 55.

    Zhang XW, Ye ZW, Wang TK, Xiong HR, Yuan XL, Zhang ZG, Yuan YL, Liu Z. Characterization of the global transcriptome for cotton (Gossypium hirsutum L.) anther and development of SSR marker. Gene. 2014;551(2):206–13.

  56. 56.

    Deng JM, Xiong GS, Yuan YL, Jia F, Liu Z. Differences in SOD, POD, CAT activities and MDA content and their responses to high temperature stress at peak flowering stage in cotton lines with different tolerance to high temperature. Cotton Sci. 2010;22(3):242–7.

  57. 57.

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

  58. 58.

    Morrissy AS, Morin RD, Delaney A, Zeng T, McDonald H, Jones S, Zhao Y, Hirst M, Marra MA. Next-generation tag sequencing for cancer gene expression profiling. Genome Res. 2009;19:1825–35.

  59. 59.

    Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nat Methods. 2012;9(4):357–9.

  60. 60.

    Kanehisa M, Araki M, Goto S. KEGG for linking genomes to life and the environment. Nucleic Acids Res. 2008;36:480–4.

  61. 61.

    Huang Y, Mo YJ, Chen PY, Yuan XL, Meng FN, Zhu SW, Liu Z. Identification of SET domain-containing proteins in Gossypium raimondii and their response to high temperature stress. Sci Rep. 2016;6:32729.

  62. 62.

    Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2-CT method. Methods. 2001;25(4):402–8.

Download references


We thank Dr. Wenkai Ren for critical reading and editing the manuscript.


This work was supported by the National Natural Science Foundation of China (Nos. 31771847 and U1403284), and the Graduate student research innovation project of Hunan province (CX2017B360). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Availability of data and materials

The sequencing data are deposited in NCBI Sequence Read Archive (SRA, with accession number SRP143909.

Author information

JC and ZL conceived, organized and planned the research. JC and PS drafted the manuscript. JC, PC and QL analyzed the data of transcriptomic and microRNA sequencing. ZL and XY modified the manuscript. All authors read and approved the final manuscript.

Correspondence to Zhi Liu.

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.

Publisher’s Note

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

Additional files

Additional file 1:

Table S1. Data statistics for transcriptome sequencing. (XLSX 11 kb)

Additional file 2:

Figure S1. Sample-to-Sample clustering. (PDF 9 kb)

Additional file 3:

Table S2. Differentially expressed genes KEGG enriched TOP10. (XLSX 21 kb)

Additional file 4:

Table S3. Metabolism and signaling pathway gene expression levels. (XLSX 13 kb)

Additional file 5:

Figure S2. The expression profiles of gene related with metabolism and signal pathways during anther development. (PDF 221 kb)

Additional file 6:

Table S4. Differential genes expression. (XLSX 31 kb)

Additional file 7:

Table S5. Data statistics for microRNA sequencing. (XLSX 11 kb)

Additional file 8:

Table S6. Target gene prediction for differentially expressed miRNAs, which between SCP and MRP, and KEGG enrichment of target genes. (XLSX 11 kb)

Additional file 9:

Table S7. Prediction of target genes for differentially expressed miRNAs, which between SCP and PM and KEGG enrichment of target genes. (XLSX 11 kb)

Additional file 10:

Figure S3. Target gene enrichment and cluster analysis. A: KEGG enrichment of miRNAs, differentially expressed between SCP and MRP, target genes; B: Expression profiles of miRNAs target genes; C: GO enrichment of target genes (BP: Biological process; CC: Cellular component; MF: Molecular function). (PDF 429 kb)

Additional file 11:

Figure S4. Target gene enrichment and cluster analysis. A: KEGG enrichment of miRNAs, differentially expressed between SCP and PM, target genes; B: Expression profiles of miRNAs target genes. C: GO enrichment of target genes (BP: Biological process; CC: Cellular component; MF: Molecular function). (PDF 501 kb)

Additional file 12:

Figure S5. Three miRNAs were used for the judgment of anther development stage according to the expression relationship were confirmed by real-time quantitative RT-PCR. The order of expression amount is ghr-miR396a > Dt_chr7_8575 > ghr-miR2948-5p in SCP, Dt_chr7_8575 > ghr-miR396a > ghr-miR2948-5p in MP, ghr-miR2948-5p > Dt_chr7_8575 > ghr-miR396a in MRP, and Dt_chr7_8575 > ghr-miR2948-5p > ghr-miR396a in PM. (PDF 67 kb)

Additional file 13:

Figure S6. The four stages of cotton pollen development. Cotton bud and anther sizes at the sporogenuous cell proliferation (SCP), meiotic phase (MP), microspore release period (MRP) and pollen maturity (PM) four continuous stages. A: Bud; B: Stamen; C: Anther; D: Pollen cellular state under optical microscope (stained with acetocarmine). (PDF 165 kb)

Additional file 14:

Table S8. Primers list in real-time qRT-PCR. (XLSX 11 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark


  • Cotton
  • Transcriptome
  • microRNA
  • Anther development
  • Gene expression regulation