Dynamic profiles of DNA methylation and the interaction with histone acetylation during fiber cell initiation of Gossypium hirsutum

Abstract

Background

Fiber, as the main product of cotton, provides main raw material for the textile industry. Many key factors have been revealed a significant role in fiber cell development including Myb proteins, phytohormones, fatty acid metabolites, and epigenetic modifications. DNA methylation is one of the important epigenetic modifications to regulate plant development and responses to abiotic or biotic stimuli. In general, DNA methylation consisting of 5mC and 6mA regulates the chromatin structure and gene transcription to affect plant development, however, the detailed role and underlying mechanism of DNA methylation in the fiber development of cotton are yet vague.

Results

Here, systematical study of the 5mC and 6mA DNA methylation profiles during the fiber initiation period of Xu142 and its glabrous mutant Xu142fl represented a clear alteration of global DNA methylation associated with fiber cell initiation. Then, the genome-wide identification of genes responsible for methylation regulation at the fifth carbon of cytosine and the sixth carbon of adenine of DNA was operated in Gossypium hirsutum. As a result, 13, 10, 6, and 17 genes were identified for 5mC methylation, 5mC demethylation, 6mA methylation, and 6mA demethylation, respectively. We then investigated the tissue expression pattern of all these genes, and some genes showed higher expression levels in fiber initiation, among which some displayed a significant change in transcription between Xu142 and Xu142fl. The possible interaction between histone acetylation and DNA methylation in fiber initiation through in vitro culture was studied by dot blot, and the results showed that repressed histone deacetylation by Trichostatin A (TSA) inhibited the global DNA methylation, and some causal genes (e. g., GhDMT13, GhDAMT2, GhALKBH12, GhDM7) were also identified.

Conclusions

In this study, all the findings indicated the interplay between histone acetylation and DNA methylation, supporting their important roles and providing precious clues for the epigenetic modifications associated with DNA methylation in the fiber development of cotton.

Background

As one of the most important cash crops, cotton has been planted worldwide and was paid much effort to improve its yield and fiber quality. Fiber, as the main product of cotton needs about 50 days to mature from fiber cell initiation with four different yet overlapping stages: initiation, elongation, secondary wall deposition, and maturation (Qin and Zhu 2011). Much progress has been made to explore the underlying mechanisms of fiber cell development, and many key genes and associated pathways have been identified, among which epigenetic mechanism play an important role (Du et al. 2018; Li et al. 2015; Qin and Zhu 2011; Sun et al. 2017).

Epigenetic modification characterizing the modulation of heritable phenotypic changes without any alteration in DNA sequence, has been proved to be a principal regulatory mechanism involved in diverse physiological or developmental processes. The research area of epigenetics covers various topics, typically including DNA methylation, chromatin remodeling, histone modification, non-coding RNAs (ncRNAs), etc. (Katada et al. 2012). All of these have been studied extensively so far and showed their crucial roles in plant development including cotton fiber development (Feng et al. 2018; Guan et al. 2014; Kumar et al. 2018; Wang et al. 2019b). Even so, the detailed key factors or mechanisms associated with DNA methylation, a well-known and crucial epigenetic modification is unclear in regulating fiber development. With the development of detection methodologies and high-throughput sequencing, the genome-wide features of 6mA, another type of DNA methylation are being uncovered and studied to reveal its role and mechanisms in eukaryotes. Consequently, recent studies have illuminated the potential regulatory role of 6mA DNA modification in eukaryotes (Fu et al. 2015; Greer et al. 2015; Liang et al. 2018; Wu et al. 2016; Zhang et al. 2015). In Arabidopsis, 6mA sites are broadly distributed across the genome at different developmental stages, and further study revealed that 6mA DNA modification positively correlates with actively expressed genes with the different context of 6mA sites from other organisms (Liang et al. 2018). However, whether 6mA DNA modification occurs in the cotton genome is still unknown, although some genome-wide DNA methylation assays indicated the distribution and role of DNA methylation (e.g., 5mC at CHH and CHG sites) in fiber development (Song et al. 2015; Wang et al. 2016).

The interactions among different epigenetic modifications is common and complex. Studies have disclosed that 5mC and 6mA DNA methylations are the most common types of DNA modifications in eukaryotes and prokaryotes (Johnson et al. 2002; Orouji et al. 2020; Tamaru et al. 2003; Xiao et al. 2018; Zhang et al. 2015). Furthermore, 5mC is generated by DNA methyltransferase (DNMT) and removed by DNA glycosylase (DML), respectively (Jones and Takai 2001). While, AMT and ALKBH have been demonstrated as the methyltransferase and demethylase responsible for 6mA modification, respectively (Xiao et al. 2018). Recently, an RNA methylome presented the crosstalk between RNA and DNA methylation in fruit ripening. The expression of demethylase SIALKBH2 which causes the m⁶A reduction on fruit-ripening gene transcripts, is regulated by the 5mC demethylase SIDML2. In turn, SIALKBH2-mediated m⁶A RNA demethylation modulates the stability of SlDML2, and mutation of SIALKBH2 reduces the SlDML2 transcription levels and delays fruit ripening (Zhao et al. 2021; Zhou et al. 2019), showing a feedback loop between 5mC DNA methylation and m⁶A RNA methylation.

In Neurospora, DIM-5 encodes a specific histone H3 methylase responsible for lysine 9 tri-methylation (K9Me3) and DIM-5 mutation showed complete loss of DNA methylation suggesting that H3K9Me3 is an epigenetic mark that involved in the maintenance of DNA methylation in Neurospora (Tamaru and Selker 2001; Tamaru et al. 2003). However, in Arabidopsis, the mutation in KRYPTONITE (KYP) which is the first reported H3K9-specific methyltransferase in plant, results in normal phenotype and unaffected level of CG methylation, revealing that H3K9Me is not necessary for gene silencing maintenance (Jackson et al. 2002; Johnson et al. 2002). The different correlations between DNA methylation and histone methylation in Neurospora and plants indicate the functional diversification of epigenetic modification in the evolution of organisms.

Studies also implied that the DNA demethylation complex identifies acetylated H3K18 and H3K23 by REPRESSOR OF SILENCING 4 (ROS4)/INCREASED DNA METHYLATION 1 (IDM1) for DNA demethylation at some loci, in which the histone acetyltransferase IDM1 functions on the same genetic pathway of ROS1 in DNA demethylation regulation suggesting the prerequisite of histone acetylation for active DNA demethylation (Qian et al. 2012; Wang et al. 2015; Zhao et al. 2014). Peroxisomal acyl-CoA oxidase 4 (ACX4), which is responsible for the fatty acid β-oxidation, influences the nuclear histone acetylation and DNA methylation at some endogenous genomic loci and is also a target of the demethylation enzyme ROS1 (Wang et al. 2019a). This further supports the implicated interaction between DNA methylation and histone acetylation. Up to now, both DNA methylation and histone acetylation show an important regulatory role in fiber development (Kumar et al. 2018; Wang et al. 2016), however, the underlying mechanisms and the interplay between them are almost blank in the fiber initiation. Here, using the known glabrous cotton Xu142fl and its counterpart Xu142, we systematically investigated the dynamics of DNA methylation level including 5mC and 6mA, the interaction between DNA methylation and histone acetylation during fiber initiation, and then screened the potential responsible genes for DNA methylation regulation to provide clues for the epigenetic mechanisms underlying fiber initiation.

Results

Global DNA methylation assay during fiber initiation of Xu142fl and its counterpart plant

A previous study indicated the important role of DNA methylation in fiber development (Wang et al. 2016), however, the underlying mechanism is obscure. As the known glabrous cotton, Xu142fl has been studied widely to uncover the fiber initiation mechanisms (Hu et al. 2018; Singh et al. 2020; Wan et al. 2016; Wu et al. 2018). Here, Xu142fl and Xu142 were used to investigate the potential role of DNA methylation in fiber initiation. First, with the specific 5mC and 6mA antibodies and dot blot assays, the dynamics of DNA methylation level in Xu142fl and Xu142 during the fiber initiation period was detected. According to the results, the 5mC level decreased weakly in 0 days post anthesis (DPA) ovules, while increased significantly in 5 DPA ovules of Xu142fl compared with the wild type (WT, Xu142). The further study indicated that little 5mC was detected in the 5 DPA fibers. In contrast, the 6mA level decreased weakly in −2 and 0 DPA ovules, while increased significantly in 5 DPA ovules of Xu142fl compared with WT. Same as 5mC, faint dot of 6mA was determined in the 5 DPA fibers of Xu142 (Fig. 1).

Fig. 1

Dynamic profile of DNA methylation during fiber initiation in Xu142fl and Xu142. Dot blot assays using specific anti-5mC (upper) and anti-6mA (lower) antibodies showed the global 5mC and 6mA signals in gDNA extracted from intact ovules of −2, 0, and 5 DPA. Furthermore, the fibers of 5 DPA ovules were also stripped for gDNA extraction and dot blot. The amounts of gDNA loaded are indicated on the left. The experiment was carried out in three replications and showed similar results. One representative was displayed here

Genome-wide identification of genes responsible for the DNA methylation in Gossypium

To investigate the potential genes responsible for DNA (de)methylation in fiber initiation, the genome-wide identification of DNA (de)methylase encoding genes were operated. The worm and human protein sequences were used as queries to identify the DNA 6mA methyltransferase (GhDAMT) and demethylase (GhALKBH) in Gossypium hirsutum and other plant species, respectively (Deniz et al. 2019; Greer et al. 2015). Furthermore, the previous reports and public genome data were also referred for the DNA 5mC methyltransferase (GhDMT) and demethylase (GhDM) identification (Jin et al. 2013; Sun et al. 2018; Yang et al. 2019a).

In total 13, 10, 6, and 17 potential genes were identified for 5mC methylation, 5mC demethylation, 6mA methylation, and 6mA demethylation, respectively (Additional file 2: Table S1). After that, the names were given according to the positions of methyltransferase and demethylase encoding genes on their chromosomes. The genome-wide identification and evolutionary analysis of 5mC methyltransferase and demethylase have been operated as described previously (Jin et al. 2013; Sun et al. 2018; Yang et al. 2019a). The genes responsible for 6mA methylation and demethylation were identified from other species. Consequently, every three genes encoding GhDAMT were identified in Arabidopsis thaliana, Glycine max, Theobroma cacao, Vitis vinifera, Oryza sativa, and Populus trichocarpa. For the GhALKBH encoding genes, 13, 11, 9, 10, 9, and 9 genes were identified in A. thaliana, G. max, T. cacao, V. vinifera, O. sativa, and P. trichocarpa, respectively (Additional file 2: Table S1).

Furthermore, the evolutionary relationships among 24 DNA 6mA methyltransferase genes of all referred to plant species were investigated. The evolutionary analysis divided the 6mA methyltransferase genes into three classes as DAMT-A to DAMT-C (Fig. 2A). Interestingly, the number of genes in each group is equal and two genes from G. hirsutum and each one gene from the other six species constitute a group, which indicated that DAMT genes may play conserved and general roles during plant evolution. Moreover, the 6mA demethylase genes were sorted into six classes as ALKBH-A to ALKBH-F in the phylogenetic tree analysis (Fig. 2B). Among which, groups ALKBH-B (21 genes) and ALKBH-E (7 genes) are the largest and smallest groups, respectively. Other groups contained different gene numbers from 8 to 16, and each group harbored the genes from all the other species is in agreement with the evolution analysis of the DAMT family in part. The phylogenetic analysis further showed that both ALKBH and DAMT genes in G. hirsutum always displayed a close distance with that in T. cacao in the phylogenetic tree, implying their close relationship and potential common progenitor as before (Li et al. 2014). Furthermore, the more predicted number of DNA methyltransferase and demethylase genes in G. hirsutum compared with other species suggested the gene duplication events and expansion due to polyploidization (Additional file 2: Table S1).

Fig. 2

Phylogenetic relationship analysis of DNA 6mA methyltransferase (DAMT) and demethylase (ALKBH) in different species. A A phylogenetic tree to analyze the evolutionary relationship among 24 6mA methyltransferase (DAMT) family members from 7 plant species. B The evolutionary relationship among 78 demethylases (ALKBH) family genes from 7 plant species by constructing the phylogenetic tree. A, B The phylogenetic trees were constructed by using NJ method with 1 000 bootstrap while curves with different colors indicate different groups of genes and prefixes At, Os, Gm, Pt, Vv, Tc, and Gh were used to recognize the corresponding genes from Arabidopsis thaliana, Glycine max, Theobroma cacao, Vitis vinifera, Oryza sativa, Populus trichocarpa, and Gossypium hirsutum, respectively

Chromosomal distribution and gene duplication analysis

The distribution of the predicted DNA 6mA methyltransferase and demethylase genes on the chromosomes of G. hirsutum genome (At and Dt sub-genome) was analyzed by TBtools. Chromosomal location analysis revealed that all the genes were unevenly distributed on 18 chromosomes (Additional file 1: Fig. S1). A total of 6 DNA methyltransferase genes were evenly distributed on chromosomes 6, 7, and 12 in At and Dt sub-genomes, while a total of 17 DNA demethylase genes were unevenly distributed on chromosomes 4, 6, 9, 10, 11, and 12 among At and Dt subgenomes. Among them, 6, 9, and 11 chromosomes from At and 6, 9, 11, and 12 chromosomes from Dt sub-genome harbored only one gene, while 4 and 12 chromosomes of At as well as 4 and 10 of Dt sub-genome harbored two genes, respectively. The imbalanced distribution of methyltransferase and demethylase genes for 6mA among chromosomes indicated the complicated gene expansion due to different chromosome sizes and structures in G. hirsutum genome.

Intron/exon structure and protein motifs analysis of GhDAMT and GhALKBH gene families

Gene structure analysis is also an important strategy to study genetic evolution and gene function. The phylogenetic tree was built to analyze the gene structure and conserved motifs distribution of DNA 6mA methyltransferase and DNA demethylase family members in G. hirsutum. The phylogenetic tree classified the DNA methyltransferase GhDAMT genes into three groups although gene structure revealed that 6∼9 exons exist in all GhDAMT genes. The conserved protein motif analysis showed that all 8 protein motifs were randomly distributed between GhDAMT1 and GhDAMT2. Motif 8 was not found in GhDAMT3 or GhDAMT4 however motif 6 and 8 were not identified in GhDAMT5 and GhDAMT6. Interestingly, it is noticed that a big exon spans almost half of the entire gene in GhDAMT3 and GhDAMT4 which did not display any potential functional motif (Fig. 3A).

Fig. 3

Gene structure and protein motif analysis of GhDAMT and GhALKBH protein families. A Evolutionary tree, intron/exon structure, and protein motif of six GhDAMT family proteins. B Evolutionary tree, intron/exon structure, and protein motif of 17 GhALKBH family proteins. A, B The phylogenetic tree was constructed using MEGA while the scale at the bottom presented the length of the protein. Gray boxes indicated exons whereas black lines indicate introns. The different colored boxes at the right side indicated 1–8 protein motifs in 6 GhDAMT and 17 GhALKBH proteins. a, family phylogenetic tree analysis; b, gene structure; c, protein domain

The phylogenetic tree also classified the DNA demethylase (GhALKBH) genes into three groups. The DNA demethylase genes bear a similar gene structure with 4 to 7 exons. The conserved motif distribution analysis showed that most GhALKBH genes exhibit similar motif structures, of which motif 1, 6, 7 are common to most genes and motif 3 and 5 are unique to GhALKBH1, GhALKBH10, GhALKBH16, and GhALKBH7. While, motif 8 only existed in GhALKBH3, GhALKBH12, GhALKBH6, and GhALKBH15 (Fig. 3B).

Collectively, the gene structure of the GhDAMT and GhALKBH gene families was found to be conserved, however, the differences in gene structure and conserved motifs in different groups indicated functional diversity of the GhDAMT and GhALKBH family genes in G. hirsutum.

Tissues specific expression of DNA methyltransferase and demethylase genes

Gene spatiotemporal expression is linked with the corresponding biological function. To investigate the expression profiles of different methyltransferase and demethylase genes, RNA-seq data were downloaded from NCBI to generate a heatmap. Then, transcription of each gene was analyzed at different developmental stages. The resulting heatmaps demonstrated that about half of the GhDMTs (i.e. GhDMT1, GhDMT4, GhDMT7, GhDMT9, GhDMT11, and GhDMT13) and two genes GhDM2 and GhDM7 encoding 5mC demethylase exhibited higher expression levels in ovules before and after anthesis. Additionally, GhDM2 showed higher expression levels in fiber at 15 and 25 DPA (Additional file 1: Fig. S2A, B). The 6mA methyltransferase encoding genes GhDAMT1-GhDAMT4 showed higher expression levels in 10 and 20 DPA ovules, besides, only two genes (GhALKBH3 and GhALKBH12) responsible for 6mA demethylation displayed higher expression level in the ovules and in fiber initiation stage and development stage (Additional file 1: Fig. S3A, B).

Moreover, the expression profiles of the above-listed potential genes were also analyzed in Xu142fl and Xu142 during the fiber initiation period. The results showed that GhDMT1 exhibited up-regulation while GhDMT4 and GhDMT7 exhibited down-regulation in Xu142 compared with Xu142f, respectively (Fig. 4A). Interestingly, GhDM9 displayed higher expression levels in all tested ovules and fibers of Xu142 than that of Xu142fl (Fig. 4B). Considering the 6mA methylation regulation, the GhDAMT4 showed higher transcript level although GhALKBH1, GhALKBH16 and GhALKBH17 displayed lower transcript levels during fibre initiation in Xu142 compared with Xu142fl (Fig. 5A, B). The aforementioned results indicate some potential causal factors associated with DNA methylation regulation and provide some clues for epigenetic modification mechanisms in fiber initiation.

Fig. 4

Quantitative PCR of GhDMT and GhDM during fiber initiation of Xu142 and Xu142fl. A Expression pattern of seven GhDMT genes responsible for 5mC DNA methylation at −2, 0, 5 DPA ovules, and 5 DPA fiber of Xu142 and Xu142fl. B qRT-PCR analysis of GhDMT genes responsible for 5mC DNA methylation at −2, 0, 5 DPA ovules and 5 DPA fiber of Xu142 and Xu142fl. A, B All expression levels were normalized to UBQ7, however error bars indicated SD (standard deviation) of three biological replicates. Primers are listed in Additional file 3: Table S2

Fig. 5

qRT-PCR of GhDAMT and GhALKBH during fiber initiation of Xu142 and Xu142fl. A, B qRT-PCR analysis was performed for four GhDAMT and seven GhALKBH genes responsible for 6mA DNA methylation. Total RNA was extracted from ovules of −2, 0, 5 DPA, and 5 DPA fiber of Xu142 and Xu142fl. All expression levels were normalized to UBQ7 while error bars indicated SD (standard deviation) of three biological replicates. Primers are listed in Additional file 3: Table S2

To explore the conceivable interaction between histone acetylation and DNA methylation for fiber development, the in vitro TSA (Trichostatin A)—histone deacetylase inhibitor-treated ovules were used to test the alteration of DNA methylation level and transcription of DNA methylation related genes. The results suggested that TSA treatment decreased the global 5mC level weakly, while the 6mA level decreased significantly after 4 and 7 days in vitro (Additional file 1: Fig. S4). These results indicated the interplay between histone acetylation and DNA methylation during fiber development (Fig. 6). So, we also tested expression profiles of some DNA methyltransferases and demethylases after TSA treatment. The results showed that GhDMT1, GhDMT13, and GhDAMT3 were up-regulated, while GhDMT9, GhDMT11, and GhDAMT2 were down-regulated after TSA treatment. The DNA demethylase GhALKBH1 and GhALKBH12 were up-regulated, nevertheless, GhDM4, GhDM9, and GhALKBH14 were down-regulated after TSA treatment (Fig. 7), indicating that histone acetylation impacts DNA methylation through different genes and pathways.

Fig. 6

DNA methylation levels of ovules after TSA treatment in vitro culture. Ovules of −1 DPA (days post anthesis) were treated with 10 μmol·L ⁻¹ TSA for 4 and 7 days in vitro. Dot blot assay was performed using specific anti-5mC (upper) and anti-6mA (lower) antibodies with the genomic DNA from ovules treated with TSA or Mock. The amounts of gDNA loaded are indicated on the left. Images of the experiment were captured by Electron Microscope. Three replications were carried out and all showed similar results. One representative was displayed here

Fig. 7

qRT-PCR of some DNA methyltransferase and demethylase encoding genes after TSA treatment. A–D The expression pattern of 7 GhDMT, 6 GhDM, 4 GhDAMT, and 7 GhALKBH (D) genes after TSA treatment was analyzed by qPCR. Total RNA was extracted from the ovules of Xu142 treated with (10 μmol·L ⁻¹ TSA) and without (Control) for 4 days. All expression levels were normalized to UBQ7. Error bar indicated SD (standard deviation) of three biological replicates. Primers are listed in Additional file 3: Table S2

Discussion

Heritable DNA methylation is well documented and known to be crucial for diverse development stages in many organisms. Although, a large body of literature concerns the role and regulation of DNA methylation in eukaryotes including the correlation with gene expression, roles in various development regulations and stress tolerance (Hao et al. 2020; He et al. 2020; Liu and He 2020; Sun et al. 2021), but little is known about the functions and the underlying mechanisms in fiber development of cotton.

Until now, methylations at the fifth carbon of cytosine (5mC) and sixth nitrogen of adenine (6mA) dominate the DNA methylation in eukaryote cells. However, studies about 5mC and 6mA DNA methylation in cotton are scarce, and whether the 5mC and 6mA DNA methylation play roles in fiber development is yet ambiguous. In this study, a classic glabrous cotton cultivar Xu142fl and its counterpart line Xu142 were applied to investigate the profiles of 5mC and 6mA DNA methylation during fiber initiation and the potential function. Dot blot assays showed the abundant 5mC and 6mA methylation in the genomic DNA of the Xu142 and Xu142fl. Furthermore, both 5mC and 6mA DNA showed higher and lower levels in the 0 and 5 DPA ovules of Xu142 compared with Xu142fl, respectively, implying the potential and complex roles of DNA methylation in fiber initiation and primary elongation.

To reveal the possible causal factors resulting in the DNA methylation regulation in fiber development, the genome-wide identification of methyltransferase and demethylase encoding genes for 5mC and 6mA methylation were performed. Based on recent PacBio SMRT sequencing results of the G. hirsutum genome (Yang et al. 2019b), the genome assemblies of two important G. hirsutum lines TM-1 and ZM24 were used together in this study. Thirteen methyltransferase and ten demethylase encoding genes for 5mC were identified, which are not in agreement with the previous results (Sun et al. 2018; Yang et al. 2019a), the reason may be the difference between genome assemblies from Next-Generation sequencing and PacBio SMRT Sequencing technologies. Because PacBio SMRT Sequencing provides a long read, which can deal with the repetitive regions and areas with high GC content that can not be coverd and solved by Next-Generation Sequencing, as well as eliminate the amplification errors caused by PCR. In contrast, the Next-Generation sequencing provides a highly dispersed assembly with short reads, which leads to the different genome assemblies from Next-Generation sequencing and PacBio SMRT sequencing (Jayakumar and Sakakibara 2019; Midha et al. 2019; Ono et al. 2013). Moreover, 6 methyltransferase and 17 demethylase encoding genes for 6mA were identified. Almost all the associated genes were identified in these two genome assemblies accordingly except the GhALKBHs, two of which failed to be identified in TM-1 (Additional file 2: Table S1). A previous work demonstrated that a lot of genetic variations and some distinct structural variations located on chromosome A08 exist between TM-1 and ZM24, which may cause a different number of identified genes among them (Yang et al. 2019b). Furthermore, the expression profile analysis of the genes responsible for DNA methylation revealed that some genes showed special and high expression levels during ovule/fiber primary development and some tissues development (e.g., reproductive organs, ovules from −3 DPA to 20 DPA) such as GhDM2, GhDM7, GhDAMT1-GhDAMT4, and GhALKBH3/12; some genes displayed higher expression levels in Xu142 than in Xu142fl such as GhDMT1, GhDAMT4. Together, the expression analysis provided valuable information for the underlying function of genes in DNA methylation regulation to mediate fiber development.

The complex interactions between DNA methylation and other epigenetic modifications play an important role in diverse eukaryote development (Garbo et al. 2021; Gui and Yuan 2021; Jackson et al. 2002; Johnson et al. 2002; Wang et al. 2019b; Zhao et al. 2021). In plants, histone acetylation and DNA methylation regulate reciprocally in gene expression and development (Wan et al. 2016; Wang et al. 2015, 2019a). In this study, the interplay between histone acetylation and DNA methylation during fiber development was also analyzed. After the TSA-histone deacetylase inhibitor treatment, the DNA methylation level including 5mC and 6mA showed clear alteration during the fiber initiation with dot blot assay (Fig. 6), indicating that histone acetylation influences DNA methylation in fiber cell initiation; and some potential genes (e.g., GhDMT1/3/13, GhALKBH1/12/14, GhDAMT2/3) involved in the interaction between them were identified, however, the explicit relationship between histone acetylation and these genes associated with DNA methylation still need much more researches. Recent studies have applied gene-editing technology to modify the methylated cytosine to mediate plant development successfully (Liu et al. 2021), supporting the great potential of DNA methylation research in plant development and crop molecular breeding.

Conclusions

In this study, we detected the 5mC and 6mA DNA methylation pattern during the fiber initiation period of Xu142 and Xu142fl, and found a fluctuation of global DNA methylation level along with the fiber cell initiation and primary elongation. Then, the responsible genes for methylation regulation were identified in G. hirsutum. Among which, some genes showed higher expression levels in fiber initiation and significant difference of transcription between Xu142 and Xu142fl. The positive correlation between histone acetylation and DNA methylation in fiber initiation was also determined through in vitro approach, some causal genes (e. g., GhDMT13, GhDAMT2, GhALKBH12, GhDM7) were also suggested. These work provide potential key genes in DNA methylation and interesting clues about epigenetic modifications involved in fiber development.

Methods

Plant materials

The cultivar Xu142 and its natural isogenic glabrous mutant (Xu142fl) were planted in the experimental field in Institute of Cotton Research of Chinese Academy of Agricultural Sciences (Zhengzhou research base, Henan). The flower buds and flowers in the full-bloom stage were labelled on consecutive days. The ovules of −2 DPA (day post-anthesis), 0 DPA, and 5 DPA, in which the ovules and fibers were separated and used, were collected from the labelled flowers with a sterile knife. The collected ovules/fibers were frozen in liquid nitrogen and stored at −80 °C for further experiments.

Dot blot analysis of DNA 5mC and 6mA levels

The genomic DNA of each sample was extracted using the Vazyme kit (DC104) from the ovules cultured in BT medium with or without TSA in vitro. Purified genomic DNA was treated with RNase and denatured at 95 °C for 7 min, followed by chilling on ice directly. Then, the denatured DNA was spotted on a Hybond-N⁺ membrane (Amersham RPN203B) optimized for nucleic acid transfer. After that, Hybond-N⁺ membrane harboring DNA was UV-crosslinked (UPV Crosslinker CX-2000) four times, blocked with 5% of non-fat milk in PBST, and incubated overnight with anti-6mA (1:500; Epigentek A-1014) and anti-5mC (1:2 000; Synaptic Systems 202 003) antibodies at 4 °C, respectively. After incubating with horseradish peroxidase (HRP)-conjugated anti-rabbit IgG secondary antibody (EpiZyme LF101/LF102), the membrane was visualized by ECL Western Blotting Detection Kit (EpiZyme SQ201).

Identification of cotton DNA demethylase and methyltransferase family members

Previously, DNA methyltransferase and demethylase for 5mC methylation were identified in G. hirsutum (Sun et al. 2018; Yang et al. 2019a), which were used as queries to generate targets in the updated genome assemblies (Yang et al. 2019b). It was found that the DNA methyltransferase genes CotAD_49037 and CotAD_14980 in JGI correspond to Gh_D01G2041001 in TM-1 in CRI and CotAD_00990, CotAD_00992 and CotAD_41398 correspond to Gh_D08G17330 in TM-1. The protein sequences of 6mA DNA methyltransferases in worms and demethylases in humans (Deniz et al. 2019; Greer et al. 2015) were used as queries to identify the homologous proteins for other plant species. The protein sequences of the candidate genes were analyzed by using SMART (http://smart.embl-heidelberg.de/) to determine the conserved domain, then the amino acid sequences corresponding to the conserved domain were determined by PFAM (http://pfam.xfam.org/) and the genes containing the conserved domain were determined by the amino acid sequence from different species. The methyltransferase and demethylase protein sequences were confirmed by SMART (Letunic et al. 2015), Interproscan 63.0 (http://www.ebi.ac.uk/InterProScan/) (Jones et al. 2014) and hidden Markov model (HMM) as described previously (Ali et al. 2019).

Prediction of the basic structure of DNA demethylase and methyltransferase gene families

The GFF (General Characteristic Format) files for the gene were obtained from the genome annotation file. The genes structure map was drawn by TB tools. Motif analysis was performed by the online tool MEME (http://meme.nbcr.net/meme/). The physical map of the chromosome was established by the software TB tools.

Cotton DNA demethylase and DNA methyltransferase family evolution analysis

The amino acid sequence of G. hirsutum was used as a reference and the E <1× e⁻⁵ was used as a threshold to determine the homologous sequences amino acid sequence in Arabidopsis thaliana, Populus trichocarpa, Theobroma cacao, Vitis vinifera, and Oryza sativa. MEGA7.0 software was used to perform multiple sequence alignment (Clustal W) for DNA demethylase and methyltransferase sequences, and neighbor-joining (NJ) was used to establish trees. The 1 000 bootstrap replicates with 50% cutoff values were used and the tree was then visualized by Tree View 1.6 (http://etetoolkit.org/treeview/).

Cotton ovule culture

Flowers were harvested at −1 DPA, and ovaries were separated and surface sterilized by using 75% ethanol. Ovules were carefully dissected from the ovaries under sterile conditions and immediately floated on liquid media supplemented with optimized phytohormone concentrations (0.5 μmol·L⁻¹ gibberellic acid and 5 μmol·L⁻¹ indole acetic acid) and various concentrations of TSA in culture plates (Beasley and Ting 1974). The ovules were incubated at 32 °C in the dark without agitation. TSA (Millipore, 647925) was dissolved in 95% DMSO to make a 1 mmol· L⁻¹ stock solution. Fiber development for all cultured ovules was analyzed after indication time of incubation by SEM (Hitachi SU3500). The cultured ovules were used for DNA extraction and dot blot assay as well as for RNA extraction for qRT-PCR.

Quantitative reverse transcription-PCR (qRT-PCR)

Total RNA was extracted using the RNAprep Pure Plant Kit (Tiangen, Beijing, China) from ovules of −2, 0, and 5 DPA and then quantified using a NanoDrop 8000 (Thermo Scientific™). Reverse transcription was performed using transScript® First-Strand cDNA Synthesis SuperMix (Vazyme R233-01). cDNA was synthesized using 1 μg of total RNA in 20 μL reactions with HiScript® III All-in-one RT SuperMix according to manufacturer’s instruction (Vazyme R333-01). The qRT-PCR primers were obtained from the qPrimerDB database (http://biodb.swu.edu.cn/qprimerdb). UBQ7 was used as an internal control for the calculation of the relative transcript level of target genes. The primers used are listed in Additional file 3: Table S2. The PCRs were performed for each of the mRNA samples using the SYBR Green mix (AQ101-01, TransGen) and the universal conditions of amplification provided by the LightCycler 480 (Roche). The 2^−∆∆C_t method was used to calculate the relative expression of each gene with three technical and biological repetitions.