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).
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).
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).
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.
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.