Identification and characterization of histone modification gene family reveal their critical responses to flower induction in apple

Histone methylation and acetylation regulate biological processes in plants through various histone modifications (HMs) gene families. However, knowledge of HMs genes is limited in horticultural deciduous trees, including apple (Malus domestica). Here, a comprehensive study of identifying and investigating HMs genes was performed using the recently published apple genome. In total, 198 MdHMs were identified, including 71 histone methyltransferases, 44 histone demethylases, 57 histone acetylases, and 26 histone deacetylases. Detailed analysis of the MdHMs, including chromosomes locations, gene structures, protein motif and protein-protein interactions were performed, and their orthologous genes were also predicted against nine plant species. Meanwhile, a syntenic analysis revealed that tandem, segmental, and whole genome duplications were involved in the evolution and expansion of the MdHMs gene family. Most MdHMs underwent purifying selection. The expression profiles of 198 MdHMs were investigated in response to 6-BA treatment and different flowering varieties (easy-flowering ‘Yanfu No.6’ and difficult-flowering ‘Nagafu No.2’) using transcriptome sequencing data, and most MdHMs were involved in flower induction processes. Subsequent quantitative real-time PCR was then performed to confirm the expression levels of candidate MdHMs under different flowering-related circumstances. MdHMs were involved in, and responsive to, flower induction in apple. This study established an MdHMs platform that provided valuable information and presented enriched biological theories on flower induction in apple. The data could also be used to study the evolutionary history and functional prospects of MdHMs genes, as well as other trees.

Apart from their different structures, the number of HMs genes was also different in plants. A total of 136 HMs (47 HMT S , 23 HDMs, 50 HAT, and 16 HDACs) have been identified in sweet orange, and they played important roles in fruit development [6]. Additionally, 125 HMs (32 HATs, 15 HDACs, 52 HMTs and 26 HDMs) were also identified from tomato genome [7]. In total, 35 SDGs members have been identified in the grape genome and some were up-regulated during grape softening [3]. Meanwhile, HMs gene functions were partially characterized, especially in the model plant Arabidopsis. They played important roles in plant growth and development, including in photomorphogenesis, seed germination and dormancy, embryo development, flowering-related processes, fruit development, stress and defense, and hormonal signaling [8][9][10][11][12][13][14][15]. HMs can directly function in regulating flowering through their over expression or down expression. They can also affect the expression of flowering related genes. For example, an Arabidopsis thaliana HDA family member, AtHDA9 (AT3G44680), repressed flowering by affecting the acetylation of AGAMOUS-LIKE 19 (AGL19) [16]. Additionally, AtHDA19 (AT4G38130) influenced flower development together with A-class organ identity gene AP2 (APETALA2), similar as AtHDA6 (AT5G63110), which showed late flowering in the HDA6-RNAi plants [17,18]. Other genes, such as AtHAM1 (At5g64610), AtHAM2 (At5g09740), AtHAC1 (At1g79000) were also responsible for flowering time [19][20][21]. For example, the artificial microRNA AtHAM1 and AtHAM2 showed earlier flowering time, while overexpression AtHAM1 flowered later and had more rosette leaves [20]. In tomato (Solanum lycopersicum), SlHAG22, SlHAG8 and SlHAG18 were involved in vegetative or reproductive development, and SlSRT2 participated in flowering [7]. Additionally, HM genes can also regulate the expression level of flowering-related genes, such as FLOWERING LOCUS C (FLC), LEAFY (LFY), MADS AFFECTING FLOWERING4 (MAF4) and (MAF5) [20,22,23]. For example, the over-expression of HAM1 resulted in a higher H4 hyperacetylation and H4K5ac at FLC in Arabidopsis [20]. In Arabidopsis, an enriched level of histone H3 acetylation and H3K4 trimethylation at FLC and MAFs occurred in the histone deacetylase6 mutant (had6) [23,24]. Meanwhile, FLOWERING LOCUS T (FT) was also influenced by HMs. The Arabidopsis JmjC family protein T-DNA insertion mutant lines (atjmj4, AT4G2040), showed earlier flowering, which might enrich FT mRNA and H3K4me3 levels within FT chromatin [25]. Among various flowering related genes, FLC and FT were the main well researched genes that associated with HMs [25][26][27]. These indicated that HMs affect or interact with their downstream or upstream flowering genes to control flowering.
Apple (Malus domestica) is an economically important fruit tree in temperate regions worldwide, and flower induction is an important issue, which restricted fruit yield and economic incomes [28][29][30]. Hormones mediated flower induction, with GA (gibberellin) inhibiting flowering and 6BA (6-benzylaminopurine) or sugar promoting flowering, has been characterized and researched in apple [31][32][33]. Additionally, some important gene families, including INDETERMINATE DOMAIN (IDD), SQUAMOSA PROMOTER BINDING PROTEIN-LIKE (SPL), MADs-box, and GIBBERELLIC ACID STIMU-LATED ARABIDOPSIS (GASA), have also been well identified and reported to regulate flower induction in apple [33][34][35][36][37]. However, less is known about of HMs and their potential involvement in apple flower induction. In 2017, with the publication of a new apple genome [38], it is able for us to systematically identify HMs gene family in apple and help us to make a comprehensive investigation about their characterizations and potential response to flower induction.
In this study, we identified 198 HMs gene members in the apple genome. They were 71 MdHMTs (64 MdSDGs and 7 MdPRMTs), 44 MdHDMs (16 MdHDMAs and 28 MdJMJs), 57 HATs (50 MdHAGs, 2 MdHAMs, 4 MdHACs, and 1 MdHAF) and 26 MdHDACs (16 MdHDAs, 3MdSRTs and 7 MdHDTs). Additionally, their chromosomes locations, gene and protein structures, gene phylogeny, synteny analysis and protein-protein interaction network were also performed. Meanwhile, transcriptomic sequencing of 6BA treated trees and different flowering varieties (Nagafu No.2 and Yanfu No.6) were performed to investigate their potential involvement during apple flower induction. Furthermore, quantitative real-time PCR (qRT-PCR) was employed to investigate the expression levels of candidate MdHMs in different tissues (stem, leaf, flower, fruit and bud) and different flowering circumstances (alternate bearing and sugar-treated trees), and various hormones (GA 3 , SA, ABA and MeJA) stress treatment. The results revealed valuable information of HMs genes in apple, which might be applicable to other fruit trees.

Identification and chromosomes location of HMs gene family in apple
To identify HMs gene family members in the apple genome, a HMM file of each domain was obtained from Pfam database (http://pfam.sanger.ac.uk/), as previous studies [6,7]. These files were then used as a query to search against the apple genome (GDDH13 V1.1) with HMMER3.0 [39]. The detailed accession number of each file was listed as Additional file 1: Table S1. However, there was no available HDT in the Pfam database. Thus, the protein sequences encoded by four Arabidopsis HDT genes, AtHDT1 (At3g44750), AtHDT2 (At5g22650), AtHDT3 (At5g03740), and AtHDT4 (At2g27840), were downloaded from the TAIR database (The Arabidopsis Information Resource; http:// www.arabidopsis.org/) and used as a query to search against the Genome Database for Rosaceae [apple genome (GDDH13 V1.1; https://www.rosaceae.org/] to predict candidate MdHDTs family members. Finally, the putative HMs genes, including HMTs (SDGs and PRMTs), HDMs (HDMAs and JMJs), HATs (HAGs, HAMs, HACs, and HAFs) and HDACs (HDAs, SRTs and HDTs) were manually checked to confirm their highly conserved segments. The relative locations of HMs were obtained from the apple genome [38]. They were then named according to their chromosome orders, as previous study [6].
Phylogenetic tree construction, gene structure, protein motif and domain, and orthologous genes analysis For phylogenetic analysis, MEGA 7.0 [40] was used to investigate the phylogenetic interactions of HMs between apple and Arabidopsis. The Arabidopsis and apple HMs protein sequences were aligned by the ClustalW program with default parameters. The multiple sequence alignment generated files were analyzed and then used to build phylogenetic trees with Maximum Likelihood method, pairwise deletion for sequences analysis, and a bootstrap value of 1000 times. For the gene structural analysis, gene models were downloaded from apple genome database (https://iris.angers.inra.fr/gddh13/) [38]. They were then submitted into Gene Structure Display Server (http://gsds.cbi.pku.edu.cn/) for structural analyses [41]. Additionally, the HMs protein sequences were employed to investigate conserved motifs with MEME Suite platform (http://meme-suite.org/ ), and 10 motifs were found within each gene family. Protein-protein interactions were analyzed with http:// string-db.org/. For orthologous genes identification, each pair of gene was used to BLAST with sequences homology more than 60% and e-value less than 1e-20. Gene Ontology (GO) terms analysis were performed with online database (http://www.geneontology.org).

Tandem duplication and synteny analysis
Circos v. 0.63 (http://circos.ca/) was employed to investigate their tandem duplication and synteny relationships as previous methods [33][34][35][36][37]. Tandem duplication MdHMs genes were identified according to their physical locations on individual chromosomes in the apple genome. Adjacent homologous HMs genes on the same apple chromosome with no more than one intervening gene were considered tandem duplicates. The Plant Genome Duplication Database (http://chibba.agtec.uga.edu/duplication/) was used to identify synteny blocks between Arabidopsis and apple. For orthologous gene pairs, synonymous (Ks) and non-synonymous (Ka) nucleotide substitutions were calculated according to the comparative synteny map between the apple and Arabidopsis thaliana genomes, with ClustalX and PAL2NAL programs for protein and coding sequences alignment. They were calculated with DNASP v5.10 software.

Plant materials and treatment
Trees with different flowering capabilities In total, 18 six-year-old apple trees of contrast flowering varieties (Yanfu No.6 and Nagafu No.2), which had been planted at the Apple Demonstration Nursery of Yangling Modern Agriculture Technology Park (108°70′E, 34°52'N), were used in this study. They were all grafted T337/ Malus robusta Rehd. Additionally, 'Yanfu No.6' , a spur mutation of 'Nagafu No.2' , had a higher flowering rate and greater bud morphological development [34,35]. Trees were divided into three blocks, with three of each, respectively. Terminal buds were collected from current spurs at 30, 50 and 70 days after full bloom (DAFB) [29,35]. They were then stored for further use.

Sugar and hormones treated apple trees
An additional 18 uniform 'Fuji'/ T337/ Malus robusta Rehd were used for sucrose treatment experiments in the same orchard. They were also divided into six blocks. Three of them were sprayed with 15,000 and 20,000 mg L − 1 sucrose at 29 and 36 DAFB [32], and the remaining blocks were sprayed with water as control. For 6BA treatment, 18 similar trees were used, and 300 mg L − 1 6BA was performed at 27 and 30 DAFB. They were all sprayed on the whole trees with a low-pressure hand wand sprayer in a clear morning. Samples were collected at 30, 50 and 70 DAFB and stored for further use. Meanwhile, 100 mM GA 3 , 300 μM ABA, 100 μM SA and 50 μM MeJA were treated on 2-year-old 'Nagafu No.2' trees, as water as control; leaves were collected at 0, 3, 6, and 12 h for each treatment as previous study [36].

Alternate bearing trees
Six14-year-old alternate bearing 'Fuji' trees were used in Tiandu Town, Fufeng, Baoji, Shaanxi (107°57′ E, 34°28′ N) were sampled. Samples were collected from trees in their 'ON' years (with a higher flowering rate) and 'OFF' years (with a low flowering rate) in 2014 at 30, 90, and 150 DAFB in the morning. Terminal buds of current spurs were collected from trees in their 'ON' or 'OFF' years and stored for further use.

Tissues collection
For the tissue-specific expression analyses, various tissues or organs were collected from 'Fuji'/T337 M. robusta Rehd. Flowers were collected on April 9, 2015 during the full-bloom period. Additionally, stems were collected from new shoots with diameters of 2-3 mm, while mature leaves were collected from the adjacent buds. Fruits with diameters of 2-3 cm were also collected. All samples were immediately frozen in liquid nitrogen and stored at − 80°C until used in the gene expression analyses.

RNA extraction and cDNA synthesis
Total RNA was extracted from plant tissue samples using the cetyltrimethyl ammonium bromide method with slight modifications [42]. Briefly, 900 μL extraction buffer (2% cetyltrimethyl ammonium bromide, 2.5% PVP-40, 2 M NaCl, 100 mM Tris-HCl [pH 8.0], 25 mM EDTA [pH 8.0], and 2% b-mercaptoethanol) was preheated at 65°C and added to 2-mL microcentrifuge tubes just before use. Samples containing 200 mg of bud tissue stored at − 80°C were ground to a powder, added to the tubes, and mixed with extraction buffer. After shaking and inverting each tube vigorously for 5 min and incubating at 65°C for 30 min, an equal volume of chloroform:isoamyl alcohol (24:1, v/v) was added. Each tube was shaken and inverted vigorously and then centrifuged at 12,000×g for 10 min at 4°C. For each sample, the supernatant was collected into a new tube and re-extracted with an equal volume of chloroform:isoamyl alcohol (24:1, v/v). The resulting supernatant was then transferred into a new 2-mL tube and LiCl (3 M final concentration) was added. The mixture was incubated at − 20°C for 4 h and the RNA was selectively pelleted by LiCl after centrifugation at 18,000×g for 20 min at 4°C. The pellet was resuspended in 500 μL of SSTE buffer (10mMTris-HCl [pH 8.0], 1 mM EDTA [pH 8.0], 1% SDS, and 1 M NaCl) that had been preheated to 65°C and an equal volume of chloroform:isoamyl alcohol. The mixture was then centrifuged at 12,000×g for 10 min at 4°C. The supernatant was transferred to a new microcentrifuge tube, and the RNA was precipitated with 2.5 volumes of cold ethanol at − 80°C for at least 30 min and centrifuged at 1,2000×g for 20 min at 4°C. Finally, the pellets were washed with 70% ethanol and resuspended in diethylpyrocarbonate-treated water. Total RNA integrity levels were verified by running the samples on 2% agarose gels. First-strand cDNA was synthesized from 1 μg of total RNA using a PrimeScript RT Reagent kit with gDNA Eraser (Takara Bio, Shiga, Japan) following the manufacturer's instructions.

Gene expression analysis
The expression levels of the 12 candidate HM genes were analyzed using quantitative real-time PCR (qRT-PCR). Primers were designed to span an intron-exon junction with Primer Premier 6.0 software. And they were designed with the preferred values to specific amplification with high yield as follows. (1) Length of PCR primers (18-24 bp); (2) Melting temperature (60°C); (3) GC content (40-60%); (4) GC Clamp: more than 3 G's or C's should be avoided in the last 5 bases at the 3′ end of the primer. (5) Avoided hairpins, self and cross dimer, and repeats. (6) Avoid template secondary structure and cross homology. The qRT-PCR mix (20 μL) consisted of 2-μL cDNA samples (diluted 1:8), 10 μL 2× SYBR Premix ExTaq II (Takara Bio), 0.8 μL of each primer (10 μM) (Additional file 2: Table S2), and 6.4 μL distilled deionized H2O. Each PCR assay was run on an iCycler iQ5 Real Time PCR Detection System (Bio-Rad, Plano, TX, USA) with an initial denaturation at 95°C for 3 min, followed by 40 cycles at 94°C for 15 s, 62°C for 20 s, and 72°C for 20 s. The resulting fragments were subjected to melting-curve analysis to verify the presence of gene-specific PCR products. The melting curve analysis was performed directly after real-time PCR and included an initial step of 94°C for 15 s, followed by a constant increase from 60°C to 95°C at a 2% ramp rate. The apple EF-1α gene (GenBank accession No. DQ341381) was used as an internal control to normalize all mRNA expression levels under different treatments and different tissues [34][35][36]. Experiments were performed using three biological replicates with three technical replicates. The 2 −ΔΔCt method was used to calculate the relative amount of template present in each PCR amplification mixture [43].

Statistical analysis
Gene expression data of RT-qPCR were subjected to an analysis of variance (ANOVA) at the 5% level with the SPSS 11.5 software package (SPSS, Chicago, IL, USA). Meanwhile, a detailed GO annotation was provided for all the HMs (Additional file 3: Table S3). These HMs genes were divided into three categories, including biological process, cellular component and molecular function.

Chromosome distributed of different HMs
To clearly identify the HMs, each of the HMs were named based on their chromosomal locations ( Fig. 1,

Phylogenetic and synteny analysis of HMs genes between apple and Arabidopsis
To understand their evolutionary relationship among HMs genes, four rooted phylogenetic trees, including HMTs, HDMs, HATs and HDACs genes, were built with Arabidopsis and apple HMs proteins (Fig. 3). All Arabidopsis and apple HMs genes were classified and clustered into different trends. For HMTs, all the SDGs and PRMTs genes were clustered together, with an exception of AtPRMT16 (Fig. 3a). Additionally, HDMAs and JMJs were also clustered with each other (Fig. 3b). For HATs, three HAFs genes (AtHAF1, AtHAF2, and MdHAF01) were clustered and surrounded by HAGs gene members, and other HAMs and HACs were also tightly grouped with themselves (Fig. 3c). For HDACs, three subfamilies (HDAs, HDTs and SRTs) were also clustered (Fig. 3d).
To characterize the expansion patterns of the HMs in the apple genome, a diagram together with Circos algorithm was performed and generated to investigate the duplicated blocks within the apple genome. A total of 67 pairs of HMs were identified from 18 chromosomes (Fig. 4, Additional file 4: Table S4 Additionally, a syntenic map of MdHMs and AtHMs were also created to help better understand their potential evolutionary and functional relationships. As shown in (Fig. 5, Additional file 5: Table S5), 72 orthologous pairs of MdHMs and AtHMs were found in the apple and Arabidopsis genome, including two pairs of HACs, SRTs and HDTs, one pair of HAFs and HAMs, three pairs of HAGs and PRMTs, 13 pairs of HDMAs, 14 pairs of JMJs and 31 pairs of SDGs. The remaining HMs genes did not have ortholog pairs. In addition, to understand the divergence among the orthologous gene pairs of apple and Arabidopsis, the ratio of the non-synonymous to the synonymous substitution rate (Ka/Ks) was used to evaluate the selection pressure during duplication. The Ka and Ks value was smaller in apple than betwen Arabidopsis and apple. However, the Ka/Ks values between gene pairs in apple less than 1, which was similar to apple and Arabidopsis (Additional file 6: Figure S1). The average Ka/Ks ratio between gene pairs in apple was 0.267, which was larger than between gene pairs in apple and Arabidopsis (0.106).

Structure analysis of MdHMs
As mentioned above and (Additional file 7: Figure S2) shown, different HMs genes had different typical domains.      We then investigated the structures of the 11 kinds of HMs gene families to confirm the present of each domain in apple, and a random gene was selected and to analyze their DOMAIN structure (Fig. 6) Gene structures and motifs also play important roles during gene evolution. Therefore, we performed detailed exon-intron structure and protein motif analysis for t seven candidate gene families. Seven individual phylogenetic trees (MdSDGs, MdPRMTs, MdHDMAs, MdJMJs, MdHAGs, MdHDAs and MdHDTs) were built based on protein sequences (Additional file 8: Figure S3, Additional file 9: Figure S4, Additional file 10: Figure S5, Additional file 11: Figure S6, Additional file 12: Figure S7, Additional file 13: Figure S8, Additional file 14: Figure S9 and Additional file 15: Table S6). As shown in Additional file 8: Figure S3 Figure S4). MdHDMAs also showed conserved motifs. For example, MdHDMA03, 12, 07, 08, 13, 06, 11, 02, and 01 shared only motif 1, 6, and 7; while the remaining MdHDMAs proteins, except MdHDMA09, shared more motif. MdHDMA14 and MdHDMA9 were encoded by genes having similar structures (Additional file 10: Figure S5). Gene structures and protein motifs of the MdJMJs were similar to MdHDMAs. For example, 9 MdJMJs proteins (MdJMJ17, 08, 03, 22, 06, 10, 14, 09, and 20) shared similar motifs. In addition, their gene structures showed less variability, especially among closely connected genes (MdJMJ9 and MdJMJ20, MdJMJ27 and MdJMJ16, MdJMJ23 and MdJMJ7, and MdJMJ28 and MdJMJ11) (Additional file 11: Figure S6, Additional file 15: Table S6).
The structures of one HATs and two HDACs gene families were also determined. Among the MdHAGs gene family members, most of them shared only one CDS (Additional file 12: Figure S7). For example, MdHAG18, 19, 17, 08, 14, 34, 49, 38, 13, 33, and 01 contained a coherent CDS within their gene structures. Their motifs were also conserved among some closely related genes, as seen with other HMs. As for HDACs and HDTs, they also had similar gene structures and encoded proteins with similar motifs, such as MdHDA15 and MdHDA09, MdHDA10 and MdHDA16, MdHDT01 and MdHDT03, and MdHDT02 and MdHDT04 (Additional files 13-14: Figure S8 and S9).

Analysis of HMs orthologous genes against in other species
BLASTP algorithm was employed to identify MdHMs orthologous genes with other sequenced plant species, and they were identified with e-value lower than 1e-20 and sequence homology more than 60%, as previous reported methods [44]. The  Their orthologous relationships were divided into three kinds [44]: a) genes that existed in apple and were absent from a given species; b) apple genes that had one to one orthologs in a given species; and c) apple genes that had homologs in a given species but not orthologs. As shown in Fig. 7

Interactions prediction of MdHMs protein
To further predict their biological interactions, we visualized HM proteins with Cytoscape v3.5.1 [45,46]. As shown in Additional file 16: Figure S10, Additional file 17: Table  S7 members from four HM-related clusters, − HMTs, HDMs, HATs, and HDACs, directly or indirectly interacted with other proteins. Among them, the HMTs interacted with the greatest number of proteins (25), followed by the HDACs (11). The HDMs and HDACs only interacted with five and four proteins, respectively. Some proteins, such as MdSDG29, MdHAM02, MdJMJ01, MdJMJ25, MdHAG28, and MdSDG14, could also directly or indirectly interact with at least three kinds of proteins. Totally, MdHMs regulated downstream genes or were regulated by their up-regulated genes to participate in various processes.

qRT-PCR analysis of candidate MdHMs genes
In total, 12 MdHMs (MdHAG07, 08, 24, and 34, and MdSDG07, 27, 29, 48, and 55, MdHDT03, MdHAM01, and MdJMJ28) were selected and their expression levels assessed using qRT-PCR under different flowering-related circumstances and in different tissues (stems, leaves, flowers, fruits, and buds). These candidate MdHMs showed different expression patterns in the five tissues, and 11 of them showed higher levels in leaves and buds, with the exception of MdSDG55, which was higher in flowers (Fig. 10). The 12 candidate MdHMs' responses to various hormones (GA3, ABA, SA, and MeJA) were also investigated (Additional file 18: Figure S11). MdHAG34 and MdSDG55 were not sensitive to these hormone treatments. The remaining 10 genes were up-or down-regulated at different time points after treatment, indicating that they might have roles in hormone stress responses or apple development.
Exogenous sugar treatments can promote flowering and lead to a higher flowering rate [32]. Here, we analyzed the expression profiles of the genes after sugar treatments. As seen in Additional file 19: Figure S12, these candidate genes were also responsive to sugar-mediated flowering induction during the flower-induction period, especially at 70 DAFB. For example, most candidate genes, such as MdHAG07, MdHAG08, MdSDG29, MdSDG55, MdSDG48, MdHDT03, and MdHAM01, showed different expression patterns at 70 DAFB after the sugar treatment. We further analyzed their expression profiles under different flowering-related circumstances (i.e., sugar treatments and alternate bearing). We investigated their expression levels in alternate-bearing 'Fuji' trees. The 12 MdHMs were expressed during the flowering periods of both the 'ON' and 'OFF' tree buds (Fig. 11). Among them, the MdHAG08 level was higher in the 'OFF' year in all three developmental stages, while those of MdHAG07 and MdSDG29 showed the opposite trend. MdHAG24, MdHAG34, and MdHDT03 were expressed higher at 30 DAFB in 'ON' trees but decreased at 90 and 150 DAFB. The expression patterns of MdSDG07 and MdJMJ28 also differed. At 30 and 90 DAFB, MdSDG07 was higher in 'ON' trees and then decreased, while MdJMJ28 was higher in 'OFF' trees and then decreased. Thus, the 12 MdHMs appeared to be involved in sugar-or hormone-mediated flower induction, as well as in alternate bearing.

Discussion
HMs played important roles during plant growth and development processes. Although, great advances have been made in some model plants, less has been reported in fruit trees except water deficit and fruit development traits [38,47,48]   MdHDTs), were identified in the apple genome. They were further characterized, including gene phylogeny, protein-protein interactions, and expansion and synteny analyses. In addition, we investigated their potential expression levels and roles in response to flower induction. These results will add to the knowledge in this field.
Comparison HMs genes in apple and other sequenced plant species Identification of HMs genes, including HMTs (SDGs and PRMTs), HDMs (HDMAs and JMJs), HATs (HAGs, HAMs, HACs, and HAFs), and HDACs (HDAs, SRTs, and HDTs) have been systematically or partially reported in Citrus, S. lycopersicum, Arabidopsis, Z. mays, and O. sativa [3,6,7,49]. Little is known about HMs gene families in the important economic apple trees. With the republication of apple genome, it is useful for us to explore more information for genomic analysis. In present study, we totally identified 198 putative MdHMs (Table 2). They were divided into four classifications (HMT, HDM, HAT, and HDAC), and they belonged to 11 different subfamilies, which were similar to those of other species [6,7]. Of all the HMs, SDGs were the most conserved among the species. The number of MdSDGs was nearly 1.5-fold greater than the numbers of SlSDGs, CsSDGs, AtSDGs, OsSDGs, and ZmSDGs. Additionally, the number of MdHDMAs genes was nearly two to four times greater than those of other species. The number of HAGs in apple was greatly different from other species, especially Arabidopsis, O. sativa and Z. mays. This great difference was partially diminished when the AT1 domain was used as the query in a BLAST algorithm-based search, which led to 33 HAGs being identified in Arabidopsis [4]. Other genes, such as JMJs, HADs, and HDTs, were also present two times more in apple than in other species (Table 2). We also searched orthologous genes against in other species of HMs genes (Fig. 7), which would be a useful tool for further analysis.
The density of apple HMs was followed by Arabidopsis and citrus (Additional file 20: Table S8). Arabidopsis were the most redundant, while Z. mays were the least dense, which might be the result of its large genome size [50]. The apple genome is nearly 1.7 times larger than that of citrus, but their gene numbers are not positively correlated with genome size. Similar relationships exist in other plants ( Table 2). The complex connections between genome size and gene number are not well characterized, partially because of duplication events in the genomes of different species and/or their complicated species characterizations. Meanshile, 198 MdHMs were not equally distributed on the 18 chromosomes of the apple genome ( Fig. 1), similar to those of citrus [6]. Additionally, this irregular distribution was noticed for the MdGRAS and MdGASA gene families [33,35]. Theoretically, apple chromosome 15 was longer and larger than other chromosomes, it could easily contain more genes (Fig. 2) [38].
The gain or loss of an exon or intron is always associated with the diversification of gene families. These events can be caused by chromosomal rearrangements or fusions, and they can result in distinct functional characterizations [51]. In the present study, MdHMs genes with different structures were always distantly clustered, while genes with similar structures were tightly clustered (Additional file 8: Figure S3, Additional file 9: Figure S4, Additional file 10: Figure S5, Additional file 11: Figure S6, Additional file 12: Figure S7, Additional file 13: Figure S8 and Additional file 14: Figure  S9). This was also observed in the IDD and GASA gene families in the apple genome, indicating potential relationships among phylogeny, gene structures, and protein motifs [33][34][35]. The typical domains of gene clusters in 12 MdHMs were investigated. Generally, these domains were conserved in apple (Fig. 6). For example, the SET domain is conserved in other plants [6,7,52], and in apple, MdSDGs also shared this typical SET domain. Additionally, other dispensable domains were also found in MdSDG family members, as in citrus, which shared 19 different domains among its 40 CsSDGs [6]. All of the MdPRMT shared a typical PRMT5 domain. When compared with citrus and tomato, JmjC and SWIRM were also typical conserved domains in the JMJ and HDMA gene families, respectively (Fig. 6) [6,7]. Similar typical structures are found in other family members among different species. For example, the AT-N domain is in the HAGs, C-terminal MOZ-SAS is in the HAMs, HD is in the HDAs, and SIR2 is in the SRTs (Fig. 6, Additional file 7: Figure S2). Although their sequence characterizations and structures were varied and numerous, their prerequisite domains were conserved, indicating a common characteristic that was conserved among various species [4,6,7].
Evolution and expansion analysis of HM gene family To better understand their phylogenetic interactions, four phylogenetic trees (HMTs, HDMs, HATs, and HDACs) were constructed using the all of the gene members from apple and the model plant Arabidopsis. Interestingly, one PRMT gene, AtPRMT16 was clustered with other SDG genes (Fig. 3a), which might be caused by their partly matching protein sequences. The HDMs, HATs, and HDACs were also clustered. Among the HDMs, the subfamily HDMA clustered separately with the JMJ subfamily (Fig. 3b). The remaining HATs and HDACs were well organized and clustered in a logical fashion, as previous found in other species [6,7]. Totally, we firstly analyzed the subfamilies within the four trees, HMTs (SDGs and PRMTs), HDMs (HDMAs and JMJs), HATs (HAGs, HAMs, HACs, and HAFs), and HDACs (HDAs, SRTs, and HDTs). The emergence of four different trees from HM subfamilies was useful to investigate their complex phylogenetic interactions. Gene duplication contributed to the evolution of species [51]. Additionally, in apple, a recent (more than 50 million years ago) whole-genome duplication event took place, which led to a change of apple chromosomes from 9 (ancestral) to 17 (present) [53]. Using improved sequencing technology, a new version of the apple genome was recently published [38]. In the present reference genome, many identified MdHMs showed duplicated genes according to the Circos diagram (Fig. 4). In tomato, eight SlHAGs gene members, including SlHAG11, 19,20,21,22,24,25, and 26, underwent tandem duplications [7]. Additionally, SlHACs, SlSDGs, and other subfamilies were also analyzed for gene duplications, as in our study. In apple, it was reported that MdSPL, MdGASA, and MdGRAS also experienced tandem, segmental duplications or whole genome duplications, similar as the MdHMs family members [33,35,36]. This gene duplications or gene expansion was associated with the genome duplication [53].
Previously, the Ka/Ks ratio was determined to be a good indicator of positive selection (Ka/Ks > 1), neutral selection (Ka/Ks = 1), or purifying selection (Ka/Ks < 1) [54,55]. Interestingly, our duplicated gene pairs within apple, or between apple and Arabidopsis, were all less than 1 (Fig. 7), which was similar to a previous report for WRKY in Brachypodium distachyon [46], indicating their important relationships during evolution. Totally, these duplications were associated with the expansion of MdHMs, led to their diverse structures and functions.
The synteny between duplicated blocks in Arabidopsis and apple was also determined of the HMs. Because Arabidopsis is a model plant and functions of the AtHMs are better understood. A comparative genomic comparison investigation helped us understand information on AtHMs to MdHMs, and possible functions of MdHMs can be well inferred [56,57]. Here, several orthologous genes were also detected in the syntenic maps, and these orthologous genes were located in different duplicated genomic regions of the Arabidopsis and apple genome (Fig. 6), indicating that these genes were derived from a common ancestor. Previously, AtHMs genes, including, AtSDG8 [9], AtHDA9 [8], AtHDA19 [17], AtHDT1 [58], AtHDT3 [59], AtHDA15 [60], AtHAM1 and 2 [19,20], AtHAF1 [61], and AtSRT1 and 2 [62,63], were shown to be involved in flower induction. Therefore, based on the orthologous genes between apple and Arabidopsis, several MdHMs could be inferred according to their Arabidopsis comprasion. However, these need to be confirmed by further experiments.

MdHMs were putatively involved in apple flower induction
Histone modifications related genes in plants have been reviewed [12,13]. Like transcription factors, the HMs genes were also involved in various biological processes during plant growth and development, especially flower induction [64][65][66]. Various genes and gene families involved in flowering have been well characterized in plants. In apple, the MdMAD-box, MdIDD, MdGASA, and MdGRAS gene families were involved in regulating apple flowering [33][34][35]37]. However, whether MdHMs respond to flower induction was reported. Here, we proposed that MdHMs were also responsible for flower induction in apple. In the model Arabidopsis, several HM genes, such as AtSDG8 [9], AtHDA9 [8], AtHDA19 [17], AtHDT1 [58], AtHDT3 [59], AtHDA15 [60], AtHAM1 and 2 [19,20], AtHAF1 [61], and AtSRT1 and 2 [63,64] have been functionally confirmed, and they are involved in flower development. Thus, we identified candidate apple flowering-related genes by referring to their orthologous genes and their expression patterns. For example, MdHDA13, orthologous to AtHDA9, showed a consistent expression pattern during the flower stages and was expressed higher under higher flowering circumstances ('Yanfu No.6' and 6BA treatment). Similarly, MdHDA16, an orthologs of AtHDT15; MdHAM01, an orthologs of AtHAM1; and MdHAF01, an orthologs of AtHAF1, were expressed higher in 'Yanfu No.6' than in 'Nagafu No.2'. However, MdHDT04, an orthologs of AtHDT1, was more highly expressed in 'Nagafu No.2' (Figs. 8 and 9). Thus, this comparative analysis of HMs genes in apple and Arabidopsis, together with their expression patterns, provided valuable information for the involvement of MdHMs in regulating flower induction.
Leaves and buds are important organs that influence flower development [28,29,67]. Here, 11 of the candidate MdHMs were expressed higher in leaves or buds than in other tested tissues (stems, flowers and fruits), which indicated their involvement in flowering. We analyzed their expression patterns in two varieties Nagafu No.2 and Yanfu No.6. 'Yanfu No.6' is a 'Nagafu No.2' mutant that has a higher proportion of spurs, shorter shoots, larger buds and a higher flowering rate [34,35]. Most of the MdHMs were expressed and showed consistent patterns during the three developmental stages. A majority of the MdSDG genes were higher in 'Yanfu No.6' during the flower development stages (Fig. 8), indicating that methylation is occurring in 'Yanfu No.6' and 'Nagafu No.2'. Similarly, higher acetylation-related activities occurred in 'Nagafu No.2' (Fig. 9). Similar epigenetic interactions were also reported among some somatic mutations [68][69][70]. Therefore, we speculated that the up-or down-regulation of MdHMs contributed to different flowering phenomena, which directly or indirectly affected flowering. The continuous differential expression patterns of MdHMs could partly illustrate their modification processes and affect flowering. We also determined their expression levels in response to sugar treatments and hormonal stresses (Additional files 18-19: Figure S11 and S12). In general, they were also partly involved in sugar-mediated flower induction in apple.
Although crosstalk about hormones or sugar-mediated alternate bearing has been reported in perennial trees, unsloved problems were still remained [29,71,72]. Here, we investigated the expressions of 12 candidate MdHMs in alternate bearing apple trees. With less reported literature about HMs and alternate bearing, we could not make better propose about this. But they were indeed induced and showed different expression patterns in 'ON' and 'OFF' trees at different time points, indicating that they were responsible for different development stages (Fig. 11). Further researches needed to be performed to confirm this.

Conclusions
In this study, we systematically identified HMs genes in the apple genome. Their chromosome locations, gene and protein structures, phylogenetic and synteny relationships, and protein-protein interactions were also characterized. Their expression levels in different flowering ability varieties and 6BA treatment were also investigated using high-throughput RNA sequence data in the apple buds, indicated they were responsible to flower induction. Further some candidate HMs genes were then analyzed by qRT-PCR in different tissues (stems, leaves, flowers, fruits, and buds), in different hormones stresses (GA3, ABA, SA and MeJA), and different flowering related circumstances (sugar treatment and alternate bearing buds). Totally, our identification and characterization of HMs genes in apple provided useful information and enriched biological theories, which could be foundation for further analysis.