Cotton CMS line 07-113A was bred from a mutant of 07-113 by saturation backcrossing. After several years of backcrossing, 07-113A and 07-113B had a similar nuclear genetic background and formed a pair of near-isogenic lines. The plant materials were cultivated at the agricultural farm of Guangxi University, Nanning, Guangxi, China. Our previous cytological observations showed that pollen abortion occurred roughly at the tetrad stage (4 mm–5 mm in diameter). Floral buds at the tetrad stage were collected and frozen in liquid nitrogen and stored at −80°C in a freezer for further determination.

Total genomic DNA was isolated from the floral buds of cotton with DNeasy Plant Mini Kit (Qiagen, CA, United States). Total RNA from the floral buds (tetrad stage) of 07-113A and 07-113B was extracted using a Quick RNA Isolation Kit (Huayueyang, Beijing, China) following the manufacturer’s protocol. A NanoDrop 2000 spectrophotometer and an Agilent 2100 Bioanalyzer were used to check the DNA and RNA concentration and integrity. After filtering, we obtained the clean data. RNA samples that had a 260/280 ratio above 2.0 were selected for use in subsequent processes.

For normal whole-genome bisulfite sequencing (WGBS) library construction, Genomic DNA was broken into fragments with an average size of 250 bp using a Bioruptor (Diagenode, Belgium), and then, the end of DNA fragments was repaired, and A base was added at the 3′ end to connect the methylation junction. Bisulfite was treated using an EZDNA Methylation-Goldkit (ZYMO). Fractions of different insertion sizes were extracted from the same channel of 2% TAE agarose gel. The qualified library was obtained for sequencing after PCR amplification. Clean data were mapped to the reference genome by BSMAP, after making sure the quantity of clean data were sufficient after filtering the raw data. These uniquely mapped data were then used to obtain methylation information of cytosine through the whole genome, after which duplication reads were merged and the mapping results were removed according to each library. Methylation level was determined by dividing the number of reads covering each mC by the total reads (more than four reads effectively) covering that cytosine, which was also equal the mC/C ratio at each reference cytosine (Lister et al., 2009; Xiang et al., 2010). The different methylated region (DMR) from different samples was calculated by tDMR (Li et al., 2010), which was developed by BGI (Beijing Genomics Institute) as mentioned already. Putative DMRs were identified by comparison of the sample1 and sample2 methylomes using windows that contained at least five CpG (CHG or CHH) sites with a 2-fold change in Methylation level and Fisher test value of p ≤ 0.05. Otherwise, the two DMRs were viewed as independent. After iteratively merging interdependent DMRs, the final dataset of DMRs was made up of those that were independent from each other. The gene element annotation of methylated region or DMR from different samples was carried out as well. A comparison of windows containing at least five CGs (CHG or CHH) at the same location in the genomes of the two samples was made to obtain putative differentially methylated regions (DMRs) with a 2-fold change. Methylation level and Fisher’s test value of p ≤ 0.05. We used CIRCOS to compare the difference in methylation levels of DMR between samples to calculate the difference in methylation levels between the two samples and the difference in methylation levels at a site between the two samples can be calculated using the following equation:

Rm₁ and Rm₂ represent the methylation levels of mC of samples 1 and 2, respectively. If the value of Rm₁ or Rm₂ is 0 then use 0.001 instead.

Subsequently, the genes located in DMRs, called differential methylated genes (DMGs), were characterized.

Total RNA was processed by mRNA enrichment or rRNA division:magnetic beads with oligo-dT were used to enrich mRNA with poly-A tails; RNase H was used to selectively digest the DNA/RNA hybridization strand and then digest the DNA probe with DNase I and subsequently purify it to obtain the required RNA. The obtained RNA was segmented using interrupting buffer, the random N6 primers were reverse-transcribed; and then, the two strands of cDNA were synthesized to form double-stranded DNA. The end of the synthesized double-stranded DNA was smoothed and phosphorylated at the 5′ end. The 3′ end formed a sticky end with a protruding “A,” and then a bubble-like connector with a protruding “T” at the 3′ end was connected. The ligation products were amplified by PCR with specific primers. The PCR product was thermally denatured into single strands; then, a bridge primer was used to loop them into a circular DNA library for sequencing.

The filter software SOAPnuke, independently developed by BGI (Shenzhen, China), was used for statistics, and the filter software Trimmomatic was used for filtering. Reads with low quality joint contamination and unknown base N content >5% were filtered to ensure the reliability of the results. We then used HISAT to compare clean reads to the cotton reference genome,¹ and Bowtie2 (Langmead and Salzberg, 2012) to compare clean reads to the reference gene sequence. We then used RSEM (Dewey and Li, 2011) to calculate the expression levels of genes and transcripts. Finally, we added new transcripts with protein-coding potential to the reference gene sequence to form a complete reference sequence and then calculated the gene expression level.

Differentially expressed genes (DEGs) were detected using the Deseq2 method (fold change ≥2.00 and adjusted p value ≤0.05) based on the principle of negative binomial distribution and DEG detection using the previously described protocol (Anders and Huber, 2010; Love et al., 2014). According to the results of differential gene detection, the union of differential genes was analyzed by hierarchical clustering using the PHEATMAP function in R software. According to Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) annotation results and official classification, functional classification of the differential genes was performed. Meanwhile, phyper function in R software was used for enrichment analysis to calculate a p value, and then, FDR correction was performed on the p value; usually a Q-value ≤0.05 was considered as significant enrichment.

A PrimeScript RT Reagent Kit (Takara, Dalian, China) was used to reverse RNA into cDNA to verify the relative expression of the DEGs by using quantitative reverse transcription polymerase reaction (qRT-PCR). The 18 s gene of cotton served as the endogenous control. All the primers were designed by Premier 5.0 and synthesized by BGI (Shenzhen, China). Gene expression levels were determined in triple biological repeats as previously described (Kong et al., 2017b).

Differentially expressed genes associated with DMRs (DMEGs) were screened out by GO analysis and KEGG enrichment analysis, using the GO database,² and the KEGG database for enrichment analysis.³ GO terms or KEGG pathways with a corrected p value ≤0.05 were considered significantly enriched.

Floral organs play an important role in plant reproduction. A stamen abnormality, in which pollen is absent or the anthers do not undergo dehiscence usually results in CMS. In the process of microspore development, the floral organ phenotypes of CMS line 07-113A and its maintainer line 07-113B were basically the same. However, there were obvious differences in the anthers. The anthers of 07-113B were normally full and powdery, and pollen grains were released after dehiscence while the anthers of 07-113A were shriveled (Figure 1).

To explore the potential role of DNA methylation in cotton pollen abortion, the floral buds (tetrad stage) of the CMS line 07-113A and its maintainer line 07-113B were collected for constructing genomic DNA libraries. We sequenced two samples for Gossypium_barbadense species using WGBS technology, generating on average 84.660 Gb clean bases after filtering low quality reads, N reads, and adaptor sequences. Table 1 summarizes sequencing data for each sample. According to the sequence context of cytosine, cytosine in the DNA sequence can be classified into three types: CG, CHG, and CHH (H = A, G, or T). Next, these mapped data were used to retrieve the methylation level (ML) of each cytosine site in the CG, CHG, and CHH contexts. The percentages of mCG, mCHG, and mCHH in the whole genome of 07-113B and 07-113A were calculated in the corresponding cytosine environment. In 07-113B, the number of methylated cytosine sites in the three contexts was 34,500,648 mCG sites (33.3.0% of all mC), 33,334,348 mCHG cites (32.20% of all mC), and 35,672,033 mCHH sites (34.5% of all mC), respectively. Similarly, there were 34,496,729 mCG sites (33.9% of all mC), 33,038,776 mCHG sites (32.5%of all mC), and 34,112,831 mCHH sites (33.6% of all mC) in 07-113A (Figure 2A). Meanwhile, we counted the coverage of cytosine sites in the whole genome in the samples and found that 07-113A had 86.42, 87.69%, and 86.22 coverage on mCG, mCHG, and mCHH, respectively, and its genomic mC coverage was 86.41%; the corresponding values in 07-113B were 85.58, 86.95, and 84.67%, respectively, and its genomic mC coverage was 85.04%. Their genomic mC coverage was 85.04% (Figure 2B). In order to further compare the two samples, we also analyzed the whole-genome methylation level within the scope of various genomic elements, including CDS, Down2k, Up2k, mRNA, repeat sequences, and CpG island. Obviously, both 07-113A and 07-113B have relatively high methylation levels in the CDS region, repeat region, and CpG-island region, indicating that these are epigenetic regulation regions that may be related to the regulation of gene expression (Figure 2C).

To understand the potential role of DNA methylation in the development of cotton anthers, we first performed GO analysis of DMGs. We found 448 genes that were both hypermethylated and highly expressed, and these genes were classified into 1,605 functional categories, including 469 biological processes (BPs), 645 cellular components (CCs), and 491 molecular functions (MFs; Supplementary Figure S1). Here, 441 hypomethylated DEGs (hypo-DEGs) were annotated to 1,923 functional categories, including 480 BPs, 581 CCs, and 485 MFs (Supplementary Figure S2).

Differentially methylated regions (DMRs) are certain DNA segments in different samples where the genome exhibits different methylation patterns. DMRs are associated with genetic imprinting and exhibit a methylation status in individuals that is consistent with that of the parent or dam. DMRs play an important role in epigenetic variation by modifying biological processes through regulation of gene expression. Abnormally increased methylation levels can trigger abnormalities in a variety of biological processes. In the comparison between 07-113A and 07-113B, a total of 20,774 DMRs were identified, including 11,405 differential regions in CG mode, 9,328 differential regions in CHG mode, and 41 differential regions in CHH mode (Supplementary Table S1). It can be seen that DMRs are mainly concentrated in CG and CHG conditions. The situation on each chromosome is consistent with this result (Figure 3). We divided DMRs into two categories according to the sequence conditions: (1) in CG mode and (2) in CHG mode. Functional annotation of these two types of DMRs using GO and KEGG databases, respectively, showed that the number of total DMRs, on the other hand, was the least in biological processes and the most in functions, which indicates that there are multiple genes involved in different biological processes, but DMRs are more involved in molecular functions (Supplementary Figures S3, S4).

In this study, we also classified DMGs into two subgroups: DMR-associated genes and DMR-associated promoters. The results showed that there were 5,118 DMRs associated with genes and 2,586 DMRs associated with promoters, and the number of hypermethylated DMGs was higher than that of hypomethylated DMGs (Figure 4A; Supplementary Table S2). Data showed that 07-113A was more prone to hypermethylation than 07-113B. To further explore the possible link between DNA methylation and gene expression, we examined the overlap between DMGs and DEGs. In total, 3,025 upregulated genes and 3,655 downregulated genes were found, of which 2,653 were upregulated in the hypermethylated genes and 2,465 were downregulated in the hypomethylated genes (Figure 4B).

To understand the potential role of DNA methylation changes in cotton anther development, we first performed GO analysis of DMEGs. We found that there were 847 DEGs in DMR-associated genes (Supplementary Table S3), and these genes were annotated into a total of 35 functional categories, namely, 13 BPs, 12 CCs, and 10 MFs (Figure 5A). Among these biological activities, DMGs are mainly involved in cellular processes, metabolic processes, membranes, membrane fractions, catalytic activity, and binding. Similarly, in the DMR-associated promoters group, 310 DMEGs were annotated to 36 functional categories, namely, 14 BP, 12 CC, and 10 MF (Figure 5B). Surprisingly, we found that both subgroups were involved in similar biological processes. Similarly, KEGG enrichment analysis was performed in the hyper-DEGs of 07-113B and 07-113A, and the top two enriched pathways in the DMR-associated promoters group were the “plant-pathogen interaction” and “starch and sucrose metabolism” pathways (Figure 6A; Supplementary Table S4). In the DMR-associated genes subgroup, the “starch and sucrose metabolism” pathway was the most enriched pathway, followed by the “galactose metabolism” pathway and the “glycolysis and glycolytic synthesis” pathway (Figure 6B; Supplementary Table S5), with a corrected value of p ≤ 0.05.

Previous studies reported that changes in mitochondrial genome structure and sequence in higher plants are closely related to plant CMS (Hanson and Bentolila, 2004) and that mitochondria are the main sites of eukaryotic energy metabolism, whose main function is to perform the TCA cycle and oxidative phosphorylation for the synthesis of ATP to provide direct energy for life activities. Therefore, we then focused on the analysis of DMEGs involved in the oxidative phosphorylation signaling pathway and TCA pathway identified in two comparative assemblies of 07-113A and 07-113B. As a result, six genes were found to be enriched in 07-113A and 07-113B, respectively (Supplementary Tables S6, S7). Transcriptome sequencing data showed that a total of four genes were significantly expressed in 07-113B, including GOBAR_AA36705, GOBAR_AA26633, GOBAR_AA39007, and GOBAR_AA30876, of which GOBAR_AA30876 is a v-type proton ATPase subunit A, mainly responsible for ATP hydrolysis, driving cationic transport through the membrane, thus promoting the absorption of nutrients such as amino acids. GOBAR_AA36705 encodes aconitate hydratase. Aconitase exists in the inner membrane of mitochondria and easily binds with lipophilic NO, leading to the inactivation of various enzymes in mitochondria, causing inhibition of mitochondrial respiration and cell death.

In order to investigate the effect of methylation on male sterility in cotton, genes involved in “starch and sucrose metabolism,” “oxidative phosphorylation,” and “tricarboxylic acid cycle” in 07-113B and 07-113A were subsequently selected for further analysis. Notably, 12 genes showed the largest fold change (>2-fold; Supplementary Tables S6–S8). Interestingly, 11 of these genes, namely, GOBAR_AA11600, GOBAR_AA29066, GOBAR_AA35700, GOBAR_AA38817, GOBAR_AA29872, GOBAR_AA28360, GOBAR_AA12958, GOBAR_AA36705, GOBAR_AA26633, GOBAR_AA39007, and GOBAR_AA30876 were downregulated in CMS line 07-113A while only one gene was upregulated, namely, GOBAR_AA09793. To investigate the effect of DNA methylation on gene expression, we examined the transcript levels of these 12 genes in 07-113B and 07-113A using qRT-PCR. Overall, the qRT-PCR expression trends of these 12 genes were consistent with the results of RNA-Seq (Figure 7). This indicates the reliability of the transcriptome data. And among these 12 genes, except the genes GOBAR_AA30876 and GOBAR_AA09793, the methylation levels were consistent with the trend of gene expression levels, the remaining 10 genes showed opposite methylation level trends. This indicated a negative correlation between methylation and gene expression. Among them, we identified five genes most likely to be associated with CMS, namely, GOBAR_AA11600, a gene belonging to the hexokinase family; GOBAR_AA29872, a gene homologous to the Arabidopsis hexokinase superfamily protein KORRIGAN 2 (KOR2); GOBAR_AA28360, a gene involved in pollen outer wall formation; GOBAR_AA36705, a gene belonging to the aconitase/IPM isomerase family; and a polypeptidase gene encoding V-type proton ATPase subunit a3, namely, GOBAR_AA30876 (Table 2).

Epigenetics are a phenomenon in which the DNA sequence remains unchanged, but gene expression is altered. DNA methylation is one of the most deeply studied and important epigenetic modifications; many studies have reported its very close association with genome stability, gene transcriptional regulation, and disease development. Several genome-wide DNA methylome studies have shown that epigenetics plays an important role in gene regulation during seed germination in Arabidopsis, maize, and soybean (Hsieh et al., 2009; Song, 2013; Li et al., 2015b). These results suggest that gene expression may be influenced by epigenetic modifications. Whole-genome sulfite sequencing can determine methylation patterns at single-base resolution (Harris et al., 2010). Recently, WGBS has been applied to decrypt increasing numbers of plant methyl bodies from Arabidopsis to rice and cotton (Cokus et al., 2008; Lister et al., 2008; Zhang et al., 2015; Feng et al., 2016; Song et al., 2017). In this study, genome-wide DNA methylation patterns were described in the CMS line 07-113A and its near-isogenic line 07-113B using WGBS to investigate the potential relationships between DNA methylation and cotton pollen abortion.

Previous studies have shown that aberrant DNA methylation genes impair the development of pollen and flowers and are associated with CMS in plants (Li et al., 2017; Han et al., 2018). Studies have shown that epigenetic regulatory mechanisms play an important role in regulating the expression of relevant flowering genes (Khan et al., 2014). As early as 1993, Burn et al. treated Arabidopsis with the DNA methylation inhibitor 5-A07-113Ac (a cytosine analogue), which reduced genomic methylation levels and led to a reduction in FLC expression, resulting in earlier flowering in Arabidopsis (Burn et al., 1993). Previous studies have shown that methylation occurs in promoter and intron regions, while DNA methylation in the mulberry genome occurs mainly in intergenic and exonic regions (Li et al., 2020); in addition, many intergenic regions also tend to be highly methylated, while most promoter regions have low methylation levels. In our study, the number of hypermethylated genes occurring in the methylation region was more in 07-113A than in 07-113B, both in the gene and promoter regions, which indicated that 07-113A was more susceptible to hypermethylation than 07-113 B.

These results suggest that the DNA methylation pattern of cotton is similar to that of other plants and that there may be some degree of conservation in the DNA methylation pattern of different plants. It is generally believed that DNA methylation, a major transcriptional silencing pathway, is negatively correlated with gene expression levels in plants (Polansky et al., 2008). However, some studies have also shown that DNA methylation is positively correlated with gene expression. Transcriptome analysis of mulberry at different developmental stages has shown that some genes have higher DNA methylation levels along with higher expression (Li et al., 2020). In our study, a total of 847 DMEGs were identified by associated analysis of methylomics and transcriptomics. Among them, 198 genes were upregulated and 250 genes were downregulated in the hyper-DEGs. In the hypo-DEGs, 204 genes were upregulated and 206 genes were downregulated (Supplementary Tables S9, S10). This result suggests that the occurrence of methylation suppresses gene expression. KEGG analysis of this DMEG showed that the hyper-DEGs of 07-113B and 07-113A were significantly enriched in the “starch and sucrose metabolic” pathway (Figure 6A and Supplementary Table S11), followed by the “galactose metabolic” pathway.

In the process of microspore development, many previous studies have shown that insufficient energy supply plays a key role in causing pollen abortion. This energy is provided by the mitochondria. Mitochondria are the main center of energy metabolism in eukaryotes, and their main functions are the tricarboxylic acid cycle and oxidative phosphorylation, synthesizing ATP to provide direct energy for life activities. Previous studies have shown that decreased sugar and biosynthetic fluxes in starch and altered redox are likely to result in altered gene expression, ultimately leading to pollen sterility (Datta et al., 2002). Therefore, we focused on developments related to energy metabolism and biological development. Seven genes were significantly downregulated and one gene was significantly upregulated in the enrichment pathway of “starch and sucrose metabolism” (Supplementary Table S8). This is consistent with our qRT-PCR results. By analyzing the KEGG pathway of genes with different levels of DNA methylation, we found that these genes were involved in a variety of metabolic pathways, which may have a complex regulatory effect on the development of anthers. Combined studies related to male sterility and bioinformatics analysis of these DMEGs, three genes associated with CMS were selected. GOBAR_AA11600 was found to be homologous to the Arabidopsis HXK1 (hexokinase-1) gene, which belongs to the hexokinase family and may be involved in the phosphorylation of glucose during its export from the mitochondria to the cell membrane. As a sugar sensor, it can regulate the inhibition or activation of sugar-dependent genes, mediating the effects of sugar on plant growth and development and possibly regulating the execution of programmed cell death in plant cells, which is consistent with kong’s study (Kong et al., 2020), and we therefore speculated that the sterility of 07-113A may be related to aberrant PCD during pollen development. Of which GOBAR_AA29872 is homologous to KORRIGAN 2 (KOR2), a superfamily protein of Arabidopsis hexahairpin glycosidases with hydrolase activity that hydrolyzes o-glycosyl compounds, and is involved in carbohydrate metabolism. GOBAR_AA28360 is homologous to β-1,3-glucanase of Arabidopsis, an enzyme involved in pollen outer wall formation and possibly in pollen development around the microspore tetrads. Callus wall degradation during pollen development. Callus plays multiple roles in male gametophyte development, and it is hypothesized that GOBAR_AA28360 expression is repressed in 07-113A, interfering with callus synthesis during male gametophyte development and leading to CMS. In addition, we also focused on the TCA cycle and the oxidative phosphorylation pathway. Four genes (GOBAR_AA36705, GOBAR_AA26633, GOBAR_AA39007, and GOBAR_AA30876) were significantly enriched and all of them were downregulated (Figure 7). Of which, GOBAR_AA36705 belongs to the aconitase/IPM isomerase family, encoding an aconitase and a member of this family (ACO1—At35830) were shown to specifically bind the 5’ UTR of CSD2 in vitro. This protein is thought to be accumulated in the mitochondria and cytoplasm. This enzyme catalyzes the conversion of citric acid to isocitric acid via a cis-aconitate intermediate and may be involved in major metabolic pathways such as the TCA cycle. It affects CSD2 (At2g28190—a superoxide dismutase) transcript levels and may play a role in the oxidative stress response. The oxidative phosphorylation process occurs in the inner mitochondrial membrane and is closely linked to energy metabolism and the electron transport chain. We found that GOBAR_AA30876 encodes V-type proton ATPase subunit a3, a multimeric enzyme that catalyzes proton transport through the cell membrane and is involved in vesicular nutrient storage (e.g., nitrate accumulation and storage) and tolerance to certain toxic ions (e.g., zinc ion sequestration in vesicles), which leads to the inference that 07-113A anther development. It is speculated that abortion may have occurred during anther development due to blockage of the electron transport chain in 07-113A.

Recent studies have shown that the expression of genes involved in energy metabolism and redox in cotton anthers is significantly regulated by DNA methylation. For example, DNA methylation has a significant effect on the expression of cotton anthers under high-temperature stress (Ma et al., 2018b; Zhang et al., 2019b). Our KEGG analysis revealed that super-DEGs, in 07-113B versus 07-113A, were significantly enriched in the “starch and sucrose metabolism” pathway (Figure 6A). Furthermore, all these hypermethylated genes were significantly upregulated in 07-113B (Figure 7). These results suggest that DNA hypermethylation may play a key role in anther development.

Based on the results of previous studies and our findings, we hypothesized that DNA methylation in 07-113A leads to the disturbance of energy metabolism and disrupts the homeostasis of mitochondria. It suppresses the expression of CMS-related genes, which eventually leads to microspore abortion. However, since gene expression is regulated at multiple levels and by multiple factors, the mechanisms of DNA methylation-mediated gene expression regulation are very complex. It remains to be further analyzed whether these differences in gene expression levels are caused by changes in the level of DNA methylation, or whether they are co-regulated by DNA methylation and other factors.