Genome-wide identification of Gramineae histone modification genes and their potential roles in regulating wheat and maize growth and stress responses

Background In plants, histone modification (HM) genes participate in various developmental and defense processes. Gramineae plants (e.g., Triticum aestivum, Hordeum vulgare, Sorghum bicolor, Setaria italica, Setaria viridis, and Zea mays) are important crop species worldwide. However, little information on HM genes is in Gramineae species. Results Here, we identified 245 TaHMs, 72 HvHMs, 84 SbHMs, 93 SvHMs, 90 SiHMs, and 90 ZmHMs in the above six Gramineae species, respectively. Detailed information on their chromosome locations, conserved domains, phylogenetic trees, synteny, promoter elements, and gene structures were determined. Among the HMs, most motifs were conserved, but several unique motifs were also identified. Our results also suggested that gene and genome duplications potentially impacted the evolution and expansion of HMs in wheat. The number of orthologous gene pairs between rice (Oryza sativa) and each Gramineae species was much greater than that between Arabidopsis and each Gramineae species, indicating that the dicotyledons shared common ancestors. Moreover, all identified HM gene pairs likely underwent purifying selection based on to their non-synonymous (Ka)/synonymous (Ks) nucleotide substitutions. Using published transcriptome data, changes in TaHM gene expression in developing wheat grains treated with brassinosteroid, brassinazole, or activated charcoal were investigated. In addition, the transcription models of ZmHMs in developing maize seeds and after gibberellin treatment were also identified. We also examined plant stress responses and found that heat, drought, salt, insect feeding, nitrogen, and cadmium stress influenced many TaHMs, and drought altered the expression of several ZmHMs. Thus, these findings indicate their important functions in plant growth and stress adaptations. Conclusions Based on a comprehensive analysis of Gramineae HMs, we found that TaHMs play potential roles in grain development, brassinosteroid- and brassinazole-mediated root growth, activated charcoal-mediated root and leaf growth, and biotic and abiotic adaptations. Furthermore, ZmHMs likely participate in seed development, gibberellin-mediated leaf growth, and drought adaptation. Supplementary Information The online version contains supplementary material available at 10.1186/s12870-021-03332-8.

Gramineous grain crops, including Triticum aestivum, Hordeum vulgare, Sorghum bicolor, Setaria italica, Setaria viridis, and Zea mays, are widely cultivated and provide important caloric intake for humans [22]. In Gramineae species, growth and development are closely related to grain yield and quality [23,24]. Biotic and abiotic stresses markedly affect crop development and yield [25][26][27][28][29][30]. Although the functions of HMs in plant growth and environmental adaptations have been identified in some plant species [8,18,20], their characteristics and functions in T. aestivum, H. vulgare, S. bicolor, S. viridis, S. italica, and Z. mays remain unclear. The publication of the genomes of these species allows for the systematic characterization of HM genes via bioinformatics analysis.
In this study, 245, 72, 84, 93, 90, and 90 HMs were identified in the T. aestivum, H. vulgare, S. bicolor, S. viridis, S. italica, and Z. mays genomes, respectively. Their location on chromosomes, conserved domains, evolution, synteny, promoter sequences, and gene structures were analyzed. The expression patterns of TaHMs and ZmHMs in developing wheat grain and maize seed were investigated. Moreover, the responses of TaHMs to growth regulators (BR, brassinazole (BRZ), and activated charcoal (AC)) and to biotic and abiotic stresses (heat, drought, salt, insect feeding, nitrogen (N), and cadmium (Cd)) were explored, and changes in the expression profiles of ZmHMs after gibberellin (GA 3 ) and drought treatment were also analyzed.

Conserved domain and phylogenetic analyses of HM genes
Conserved HM domains were investigated, with various domains identified in the different HMs (Fig. S2) To clarify the evolutionary relationships among HM genes, unrooted phylogenetic trees were constructed (Fig. S3). All AtHMTs (except AtSDG41), OsHMTs, and TaHMTs were classed into groups A-E, which were further subdivided ( Fig. S3-1). All SDGs were clustered together in classes A-D and F-H, and PRMTs were divided into group E subgroups e1 and e2 ( Fig. S3-2). In Fig. S3-3, AtPRMTs, OsPRMTs, and SbPRMTs were clustered together in group A, which could be divided into subgroups a1 and a2, and all SDGs (except for AtSDG41) were identified in groups B-G. AtPRMTs, OsPRMTs, and SvRMTs were grouped in either class A or B, and the other HMTs were in classes C-H ( Fig.  S3-4). All Arabidopsis, rice, and S. italica SDGs were clustered together in groups B-E, with the exception of OsSDG738, SiSDG17, and SiSDG34, and PRMTs were all clustered in group A (Fig. S3-5). AtPRMTs-OsPRMTs-ZmPRMTs and AtSDGs-OsSDGs-ZmSDGs were clustered together in groups A and B-G ( Fig.  S3-6). In classes A, B, and D, the model and Gramineae JMJs were closely clustered, and HDMAs were divided into subclasses c1 and c2 in group C (Fig. S3-7). The evolutionary relationships among HATs were investigated ( Fig. S3-8). HAGs were classified into groups B, C, and E, and HAFs, HAMs, and HACs were separately divided into groups A, D, and F. HDTs and SRTs were clustered into groups A and B, while HDAs were found in groups C and D ( Fig. S3-9).
We investigated the syntenic relationships among Gramineae and Arabidopsis HMs ( Fig. 2c and   To evaluate selection pressure during duplication of the above gene pairs, their non-synonymous (Ka), synonymous (Ks), and Ka/Ks values were calculated. Data showed that the Ka/Ks values were all less than or generally equal to 1 (Tables S2, S3, S4). However, several gene pairs, such as SiJMJ5-SiJMJ19, AtJMJ13-TaJMJ3, and AtJMJ13-TaJMJ7, shared no non-synonymous mutations based on their Ks values.
We next identified HM gene structures. In general, homologous HM genes, especially those in the same pair, shared similar structures, although gene lengths differed (Fig. S8). For example, most homologous TaHMT genes contained more than one CDS and were more than 3000 bp in length ( Fig. S8-1). All HvHMTs, except for HvSDG4, shared a short non-coding sequence, and most were 2000-5000 bp in length (

Expression patterns of TaHMs in developing wheat grain in response to BR and AC
To investigate the potential roles of HMs in wheat grain growth and development, we examined their expression profiles in the endosperm, inner pericarp, and outer pericarp (Fig. 3). Based on these expression patterns, TaHMs were divided into various clusters ( Fig. 3a-d). In cluster 1, TaSDG53, TaSDG29, TaSDG56, and TaSDG61 were highly expressed in all tissues, especially in the inner pericarp. In cluster 2, several TaSDGs, such as TaSDG15, TaSDG103, and TaSDG21, were also highly expressed in the inner pericarp. In cluster 3, seven TaSDGs were found at relatively low levels in the outer pericarp. In clusters 4 and 5, TaSDGs showed lower expression levels than in the other clusters. In cluster 6, most genes were highly expressed in the inner and outer pericarps (Fig. 3a). The expression levels of TaHDMAs and TaJMJs were generally low compared with other genes (Fig. 3b). TaHATs were classified into two classes according to their expression patterns. Genes in cluster 1 were more highly expressed in all tissues than genes in cluster 2 (Fig. 3c). In clusters 3 and 4, TaHDAs, TaSRT3, and TaSRT5 were highly expressed in the pericarps (Fig. 3d).
In Arabidopsis, rice, wheat, and maize, BR plays an important role in root growth, including lateral root initiation and hair formation [33][34][35][36]. BR treatment significantly increases the number of lateral roots in wheat, but inhibits root length and diameter, whereas the BR synthesis inhibitor BRZ shows the opposite roles on lateral root number and root diameter [33]. Although HM genes are known to regulate various developmental processes, information on their roles in regulating wheat roots is scarce. In this study, we analyzed HM gene expression profiles during BR-and BRZ-mediated root growth (Fig. 4). In cluster 2, TaSDG4, TaSDG23, TaSDG55, and TaSDG112 showed a 2-fold increase after BRZ treatment, whereas, in cluster 1, BRZ treatment repressed TaSDG26, TaSDG68, TaSDG89, TaSDG92, TaSDG95, TaSDG103, and TaJMJ5 expression (Fig. 4a). BR treatment increased the expression of more than 10 TaHMs (especially TaJMJ5 and TaSDG28), while BR1 and BR2 exposure inhibited the expression of several other TaHMs ( Fig. 4b and c). For example, TaSDG26, TaSDG28, and TaJMJ5 were induced by BR1 and BR2 treatment; TaSDG92 and TaSDG101 were up-regulated by BR2 treatment; and TaJMJ21, TaSDG53, and TaHDA18 were repressed by both BR1 and BR2 treatment.
In plant culture, AC is widely used to promote seedling growth [37]. Notably, AC treatment promotes wheat seedling growth, accompanied by an increase in soluble protein, root activity, and total phenol and sugar content [37]. Here, we found that 26 and 31 TaHMs were differentially expressed in roots and leaves after AC treatment, respectively (Fig. 5), with an almost equal number downregulated and up-regulated by AC treatment. For example, after AC treatment, TaSDG68 and TaSDG84 showed a 4-and 8-fold decrease in the roots and leaves, respectively; TaSDG55 was increased in the roots; and TaJMJ21 was up-regulated in the leaves ( Fig. 5a and b).

Responses of TaHMs to abiotic and biotic stresses
To explore whether TaHMs respond to abiotic stresses, we analyzed their expression levels after heat stress (HS), drought stress (DS), and heat stress (HD) treatment using previously published RNA-seq data [38]. In total, 86 TaHMT genes (83 TaSDGs and three TaPRMTs) were differentially expressed at 1 or 6 h after the different treatments (Fig. 6a). These TaHMT genes could be divided into six clusters based on their transcription patterns. In cluster 1, almost all TaSDGs were induced and repressed by DS at 1 and 6 h, respectively, and were up-regulated by HD at 6 h. In cluster 2, 20 TaSDGs were obviously increased at 6 h after HS and HD treatment, but were decreased at 1 h, and several TaSDGs were clearly induced or inhibited by DS. In cluster 3, TaS-DGs were generally induced by both HS and HD at 6 h but were suppressed at 1 h in the HS and HD groups, and increased at 1 and 6 h in the DS group. All TaHDMs were divided into four clusters (Fig. 6b). In cluster 1, TaJMJ21 was highly expressed after HS and HD treatment, and increased at 1 h following DS treatment. In cluster 2, TaJMJ7, TaJMJ11, and TaJMJ3 were generally up-regulated by DS and HD at 1 h, whereas other genes were generally up-regulated at 6 h after HS and HD treatment. TaHATs were clustered into two classes (Fig. 6c). In cluster 1, TaHAG1, TaHAG2, and TaHAG5 were increased after HS and HD treatment. In cluster 2, TaHAM2 and TaHAM3 were obviously up-regulated at 6 h in the HS and HD groups, and other genes were induced by HS, DS, or HD at at least one time point. As shown in Fig. 6d, TaHDA4, TaHDA17, and TaSRT2 were induced by DS in cluster 1, and TaHDA4 was also increased in the HS group. DS treatment increased the expression levels of 10 TaHDACs in cluster 2, and these genes were also affected by HS and HD at several time points. Nine TaHDACs were distinctly induced by HS, DS, or HD treatment in cluster 3. In cluster 4, TaHDACs (TaHDA10, TaHDA12, TaHDA16, TaHDA19, and TaHDA21) were primarily expressed at 6 h in the HS and HD groups.
To investigate the responses of TaHMs to SS, we identified the expression profiles in the salt-sensitive wheat cultivar Chinese Spring (CS) and salt-insensitive wheat cultivar Qingmai 6 (QM). After SS treatment, almost all TaHMs were up-regulated at at least one time point (Fig. 7). The TaHMTs were grouped into three clusters (Fig. 7a). In cluster 1, TaSDG17 and TaSDG21 were significantly induced in QM at 12 h after SS treatment. In cluster 2, TaSDGs and TaPRMTs were up-regulated at several time points in both the CS and QM groups after SS treatment. Most TaHDM genes were up-regulated by SS from 6 to 24 h in cluster 1. Genes in cluster 2 were induced by SS treatment at most time points. The transcripts of TaJMJ18, TaJMJ27, TaJMJ23, TaJMJ25, TaJMJ38, TaJMJ44, and TaJMJ48 increased from 12 to 48 h after SS treatment in cluster 3. In both CS and QM, six TaJMJs were obviously induced by SS in cluster 4 (Fig. 7b). TaHATs were divided into two clusters (Fig. 7c). In cluster 1, TaHAC8 and TaHAC10 were mainly regulated by SS in CS, and TaHAF4, TaHAG3, and TaHAC6 were induced by SS in both the CS and QM cultivars. In cluster 2, SS treatment increased the expression levels of TaHACs and TaHAGs, especially TaHAC1 and TaHAC2, at every time point (Fig. 7c). The expression levels of TaHDACs, especially genes in cluster 3, were also markedly increased at at least one time point after SS treatment (Fig. 7d).
Wheat pests Sitobion avenae and Schizaphis graminum can increase yield losses [39]. Compared with the non-phytotoxic aphid S. avenae, feeding by phytotoxic aphid S. graminum causes more severe damage in wheat leaves [40]. N is an essential macronutrient for plant growth and development, and low N stress can repress wheat leaf and root growth [40]. In addition, Cd can inhibit leaf photosynthesis, carbon and N metabolism, and wheat growth and yield [41]. To clarify how TaHMs respond to biotic, nutrition, and heavy metal stress, their expression patterns were obtained from previous transcriptome research [39][40][41]. In our study, TaJMJ7 increased 3.4-and 4-fold after S. avenae and S. graminum feeding, respectively; TaJMJ11 was induced by S. avenae infection; TaJMJ40 and TaJMJ42 increased in the S. graminum feeding group compared with the control; and TaHDA17, TaSDG73, TaHDA20, TaSDG81, TaHDA22, and TaSDG89 were distinctly controlled following S. graminum feeding (Table 1). In addition, N stress suppressed TaSDG73 and TaHDA20 expression in the leaves, but up-regulated TaJMJ11 and TaJMJ3 in the roots (Table 1). Furthermore, Cd treatment induced a 2.2-to 6.4-fold increase in the expression levels of 12 TaHMs (e.g., TaSDG13, TaJMJ28, and TaHDT1) in the roots but decreased the expression of TaSDG102 (Table 1).

Diverse responses of TaHMs to growth and stress signaling
To investigate the multiple functions of TaHMs in wheat growth and stress adaptations, a Venn diagram was constructed with the above identified DEGs (Fig. 8). The DEGs were clustered into six sets, including the BR or BRZ (BR-BRZ) class, AC class, heat or drought (heatdrought) class, salt class, S. avenae or S. graminum (Sa-Sg) class, and N or Cd (N-Cd) class ( Fig. 8 and Table  S5). Some TaHMs were simultaneously respond to various signals, while several ones were only regulated by single clue. For example, two TaHMs (TaSDG68 and TaJMJ5) were concurrently in response to BR-BRZ and AC; the expression patterns of TaSDG95 and TaSDG103 were altered by both BR-BRZ and salt treatment; and 72 TaHMs simultaneously responded to heat, drought, and salt stress. We found that TaSDG13 and TaJMJ28 were common DEGs after AC, heat-drought, salt, and N-Cd treatment. TaJMJ34 was commonly induced or repressed by BR-BRZ, AC, heat-drought, salt, and N-Cd treatment. In total, 55, 23, five, and two TaHMs responded to heatdrought, salt, AC, and BR-BRZ, stress, respectively ( Fig. 8 and Table S5).

Expression analysis of ZmHMs in developing seed and response to GA treatment
To investigate the functions of ZmHMs in maize growth and development, we analyzed the expression profiles of ZmHMs in different seed growth stages of B73 and SWL01 cultivars ( Fig. 9a and b). The SWL01 cultivar is a mutant of B73 and shows higher viscosity [42]. From 0 to 24 days (d) after pollination (DAP), 80 ZmHM genes were clustered into five classes (Fig. 9a). During the whole experimental period (especially at 2 DAP), ZmSDG36 in cluster 1 showed higher expression than genes in other classes. In cluster 2, ZmHMs showed higher expression at the early stages (from 0 to 8 DAP) than during the later periods (from 16 to 24 DAP). In cluster 3, 11 ZmHMs were highly expressed at all stages, especially from 0 to 4 DAP. There were 81 ZmHM genes detected during SWL01 seed development (Fig. 9b). Like genes in B73, these ZmHMs were distributed into five clusters in SWL01. For example, ZmHMs in clusters 1 and 2 (especially cluster 2) were mainly expressed at 0, 2, and 4 DAP. ZmHMs, such as ZmSDG29, ZmSDG36, ZmSDG40, and ZmHDA1, showed higher expression levels in cluster 3 than genes in other clusters. A total of 79 ZmHMs were commonly expressed in both B73 and SWL01 seeds, but most showed different expression patterns between the two cultivars. For example, ZmSDG41 expression was higher in B73 than in SWL01; ZmHAF1 gradually decreased over time in B73 but showed almost no change in SWL01 (Fig. 9c). GA 3 application significantly promoted leaf sheath growth of D11 [43]. Seven ZmHM genes were differentially expressed between the GA and control groups (Fig. 9d). In cluster 1, ZmH-DMA3, ZmHDA10, ZmJMJ10, and ZmSDG10 were down-regulated by GA, whereas ZmHDA12, ZmHDA3, and ZmSDG33 were up-regulated.

Expression analysis of ZmHMs in response to drought stress
To identify the potential roles of ZmHMs in drought adaptation, their expression patterns were analyzed in drought-tolerant cultivars (ND476 and H082183), drought-sensitive cultivars (ZX978 and Lv28), and C7-2 ( Table 2). In total, 10 ZmHMs were identified as DEGs in response to drought stress. After drought treatment, the transcription level of ZmJMJ2 showed a 6-fold increase in ND476 compared with ZX978, whereas ZmHDA11 was repressed in ND476. ZmSDG5, ZmJMJ4, and ZmSDG24 were induced by drought treatment in C7-2, but ZmSDG33 and ZmJMJ17 were controlled. In Lv28 and H082183, ZmJMJ5 was up-regulated under both moderate and severe drought treatment. The expression level of ZmSDG1 increased in H082183 after moderate drought treatment, whereas ZmHDA2 was markedly down-regulated after severe drought treatment.

Discussion
Although HM genes are known to play essential roles in plant growth and biotic and abiotic stress in model plants [8,18,20], little information has been reported for Gramineae species. Here, we systematically characterized TaHMs, HvHMs, SbHMs, SvHMs, SiHMs, and ZmHMs, including information on their gene location, conserved domains, gene phylogeny, gene expansion, synteny, promoter cis-elements, and gene structure. Moreover, we analyzed their expression levels in wheat and maize in regard to growth and stress adaptations. These findings will provide a basis for further functional analyses of HM genes. 10 ZmHATs, and 18 ZmHDACs) (Fig. 1). In terms of gene number, we found 2.4-and 2.6-fold greater number of TaHMs than AtHMs and OsHMs, respectively. TaS-DGs, TaHDMAs, TaJMJs, TaHAGs, TaHAMs, TaHACs, TaHAFs, TaHDAs, and TaSRTs were increased 1.5-3fold. However, the number HM genes in other species varied slightly compared with those in the model plants ( Fig. 1a and b). In wheat, a total of 144 gene pairs were identified in 10 kinds of HM genes, but there were 4-14 HM gene pairs among S. bicolor, S. viridis, S. italica, and Z. mays, and no gene duplication in H. vulgare. Genome duplication occurs during species evolution [45], and the wheat genome contains three homologous subgenomes [22]. Therefore, the expansions in wheat HM genes may be associated with gene and genome duplications during evolution.
To better understand Gramineae HMs, duplicated blocks between model plants and Gramineae were determined. In this study, 13 orthologous genes were identified between Arabidopsis and the six Gramineae species (Fig.  S5 and Table S3)), and 389 rice-Gramineae gene pairs were found (Fig. S6 and Table S4), indicating that these gene pairs shared common ancestors. Gene pairs showed considerable differences between Arabidopsis-Gramineae and rice-Gramineae in terms of number, which may be due to the diversity in evolutionary history between monocotyledons and dicotyledons. Several AtHMs and OsHMs are involved in plant growth and stress responses [9, 10, 12, 15-17, 20, 21, 46-48]. Although many unknown Gramineae HMs could be inferred from the orthologous genes of model plants, these predictions must be confirmed in future experiments. Gene evolution mode can be determined through Ka/Ks values. Here, the Ka/Ks ratios of all gene pairs were less than 1, indicating purifying selection [49].

Potential functions of TaHMs and ZmHMs in plant growth and stress responses
Like transcription factors, HMs are important regulators of many biological processes, including plant growth and development [1][2][3]. We proposed that TaHMs and ZmHMs share similar roles with known HMs. Candidate TaHMs involved in wheat grain development and ZmHMs involved in maize seed development were characterized in this study. Expression patterns showed that almost all TaHMs (especially TaSDGs in cluster 1 (Fig. 3a), TaHDMs in cluster 3 (Fig. 3b), TaHATs in cluster 1 (Fig. 3c) and TaHDACs in cluster 4 (Fig. 3d)) were expressed in developing wheat grains, and many genes were highly expressed in specific grain tissue layers (Fig. 3). About 80% ZmHMs showed different expression   Table 2 Expression analysis of ZmHMs during different drought stresses TD tolerant cultivar ND476 drought treatment, SD sensitive cultivar ZX978 drought treatment, CD C7-2 drought treatment, CC C7-2 control, LMD Lv28 moderate drought treatment, LMC Lv28 control, LSD Lv28 severe drought treatment, LSC Lv28 control, HMD H082183 moderate drought treatment, HMC H082183 control, HSD H082183 severe drought treatment, HSC H082183 control

Gene_id TD/SD CD/CC LMD/LMC LSD/LSC HMD/HMC HSD/HSD
patterns in developing maize seeds ( Fig. 9a and b). Several ZmHM genes specifically expressed in B73 (ZmSDG23) or SWL01 (ZmSDG14 and ZmJMJ4) were found (Fig. 9c). In addition, several commonly expressed ZmHMs between B73 and SWL01 were found but showed varied expression patterns (Fig. 9c). Moreover, seed-specific motifs of ZmHMs were identified. These findings suggest that TaHM genes affect grain growth and development, most ZmHM genes play roles in wax and regular maize seed development, and several ZmHM genes specifically participate in regulating seed viscosity. BR is an essential plant hormone and stimulates wheat root hair formation and lateral root initiation [33]. However, responses of TaHMs to BR and BRZ are not known. In this study, four TaSDGs were induced by BRZ, but six TaSDGs as well as TaJMJ5 and TaHDA14 were repressed (Fig. 4a). In addition, BR respectively increased or decreased the expression of 11 TaHMs (Fig. 4b). We also found that GA treatment stimulated leaf sheath elongation of maize seedlings and altered the expression of seven ZmHMs (Fig. 9d). The above results indicate that these TaSDGs and ZmSDGs are likely involved in BRmediated root growth and GA-mediated leaf development. AC is a positive growth regulator in wheat culture [37]. However, the relationship between TaHMs and AC is unclear. Here, 26 TaHMs were differentially expressed between the control and AC-treated roots, with about half repressed or induced by AC, respectively (Fig. 5a). In leaves, 16 TaHMs were regulated by AC, with 15 found to be highly expressed (Fig. 5b). Thus, these up-and downregulated TaHMs are speculated to play important roles in AC-promoted wheat seedling growth.
In addition to their important functions in growth, HM genes also play essential roles in plant defenses [9,17,21,46]. Here, TaHM-mediated stress responses were explored (Figs. 6-7 and Table 1). In total, 86 TaHMTs were differentially expressed after HS, DS, or HD treatment (Fig. 6a), and 45 TaHDMs, 20 TaHATs, and 27 TaH-DACs were induced by stress treatment (Fig. 6b-d). In response to SS, almost all TaHMs were increased, especially TaSDGs in cluster 3 (Fig. 7a), TaJMJs in cluster 4 ( Fig. 7b), TaHATs in cluster 2 (Fig. 7c), and TaHDACs in cluster 3 (Fig. 7d). The expression patterns of 10 TaHMs, including TaSDG73, TaSDG81, TaSDG89, TaJMJ7, TaJMJ11, TaJMJ40, TaJMJ42, TaHDA17, TaHDA20, and TaHDA22, were affected by S. avenae or S. graminum feeding (Table 1). Furthermore, N stress regulated the expression of four TaHMs (TaSDG73, TaJMJ3, TaJMJ11, and TaHDA20) ( Table 1). Transcriptions of 13 TaHMs were influenced by Cd treatment, with most found to be increased (Table 1). Several ZmHMs were up-regulated or down-regulated by drought treatment (Table 2). A number of stress-related elements were identified in TaHMs and ZmHMs, which may partly explain their responses to stress. The above findings suggest the occurrence of methylation when wheat and maize experience biotic or abiotic stresses. The multiple functions of TaHMs are discussed in Fig. 8 and Table S5. In total, 85 TaHMs were simultaneously regulated by two signals; 25 TaHMs were simultaneously regulated by three treatments; nine TaHMs were up-regulated or down-regulated by four signals; and one wheat gene was simultaneously regulated by five treatments. The diverse functions of these TaHMs indicate that they are essential for wheat growth and stress adaptations, and thus warrant further study. Moreover, all ZmHMs, except for ZmJMJ2 and ZmJMJ4, that responded to drought stress were also identified in developing seeds, indicating their roles in maize growth and stress adaptations.

Conclusions
TaHMs, HvHMs, SbHMs, SvHMs, SiHMs, and ZmHMs were systematically explored in our study to clarify their chromosome locations, protein structures, gene duplications, promoters, and gene structures. Phylogenetic and synteny comparisons between model plant and Gramineae HMs were performed and the potential roles of Gramineae HMs were posited through their known homologs. The unique characteristics of the HM genes were investigated based on their domains and expansions. Specific domains were identified in several Gramineae species, e.g., SDGs, PRMTs, JMJs, HDAs, and HDTs, which may exhibit unique functions. The expansion patterns of Gramineae HMs were analyzed to elucidate differences in gene number and function among Gramineae species. Using previously published RNA-seq data, we also investigated the potential roles of TaHMs in developing grain, as well as BR-mediated root growth and AC-regulated seedling development and explored the functions of ZmHMs in seed development and GAmediated leaf growth. Candidate wheat HMs involved in high temperature, drought, salt, insect feeding, nutrition and heavy metal stress were analyzed. In addition, the ZmHMs involved in drought response were examined, and their responses to the above-mentioned stresses were inferred through promoter analysis. In summary, based on bioinformatics analysis, we predicted the functions of the Gramineae HMs, with the potential functions of wheat and maize in growth and stress adaptations verified based on expression profile analysis. The results of this study will lay the foundation for future research.
Protein domain composition, gene structure, promotor cis-acting elements, phylogenetic tree, orthologous genes, and heatmap analyses Using protein sequences, protein domain files were generated from the Batch CD-Search database (https:// www. ncbi. nlm. nih. gov/ Struc ture/ bwrpsb/ bwrpsb. cgi). Visualization of conserved motifs was completed with 'Visualize NCBI CDD Domain Pattern' in TBtools [51]. 'Gene Structure View (Advanced)' in TBtools [51] was employed to investigate the HM gene structures. Promoter analysis was performed in the PlantCARE database (http:// bioin forma tics. psb. ugent. be/ webto ols/ plant care/ html/) and visualized using 'Simple Biosequence Viewer' in TBtools [51]. MEGA X was used to construct a phylogenetic tree with the maximum-likelihood method, and protein sequence alignment was completed using the ClustalW method [57]. Synteny blocks were obtained from genome sequences and general feature format v3 (gff3). The 'Fasta stats,  [51] were used to visualize syntenic relationships among homologous HM genes. Using 'Simple Ka/Ks Calculator (NG)' in TBtools, the non-synonymous (Ka) and synonymous (Ks) nucleotide substitutions of orthologous gene pairs were calculated with coding sequences [51]. Heat maps were generated with the 'HeatMap' tool in TBtools [51].

Plant material and treatment
All plant materials mentioned here were used in previous studies, which provide the original sources and formal identification of plant materials.

Seed development
Immature wheat grains were dissected into endosperm and inner and outer seed pericarp tissues at 12 d postanthesis and sampled for library construction and RNA sequencing (RNA-seq). RNA-seq expression data [fragments per kilobase per million (FPKM)] in the above tissues were uploaded to the WheatExp database (https:// wheat. pw. usda. gov/ Wheat Exp/ about. html) [58]. For regular field-grown maize (B73), seeds were sampled at each growth phase. The waxy maize inbred line (SWL01) was grown and sampled in the same way. We downloaded transcriptome data from a previous study [42].

GA 3 treatment
Dwarf D11 is a GA-sensitive maize mutant. Seedlings in the control and GA-treated groups were treated daily with distilled water or10 − 4 M GA 3 (Sigma), respectively. The second leaf sheaths of D11 were collected at the three fully-expanded-leaves (V3) stage for RNA-seq, as described previously [43]. The RNA-seq data were obtained from a previously published study [43].

AC treatment
After sterilization, immature wheat (T. aestivum) embryos were used to investigate the effects of AC (a widely used growth regulator in plant tissue culture) on wheat seedling growth in medium. The roots and leaves were used for RNA isolation and library construction, respectively. The DEGs were obtained from the supplementary data of a previously published paper [37].

Heat, drought, and salt treatment
In the control group, seven-day-old TAM 107 seedlings were grown in hydroponic solution under 16 h day (22 °C) and 8 h night (18 °C) conditions. Seedlings were treated with 40 °C in the HS group, with 20% (m/V) PEG-6000 in the DS group, and with combined heat (40 °C) and drought stress (40 °C and 20% PEG-6000) in the HD group [38]. Seedling leaves were sampled for RNA-seq after 1 and 6 h of stress treatment. DEG analysis was performed by Liu et al. (2015), and DEGs were downloaded from the WheatExp database (https:// wheat. pw. usda. gov/ Wheat Exp/ about. html) [38].
In terms of drought treatment for maize, three published studies were cited [59][60][61]. Firstly, drought-sensitive ZX978 and drought-tolerant ND476 cultivars were used in early study. Maize seedlings planted in soil with 70-80% and 15-20% water content were set as the control and drought-treated groups, respectively. After 12 d of treatment, flag leaves of ND476 and ZX978 under both control and drought conditions were sampled for RNAseq [60]. Secondly, maize ChangC7-2 (C7-2) seedlings were used for identifying drought-tolerant mechanisms. After 7 d of drought treatment, the expanded third leaves of seven-day-old C7-2 seedlings were sampled for RNAseq analysis. The DEGs detected between the control and drought-treated groups were obtained from the additional files of a previously published study [59]. Thirdly, two maize inbred lines (drought-sensitive inbred line Lv28 and drought-tolerant inbred line H082183) were used for maize drought tolerance analysis. In the control group, Lv28 and H082183 seedlings were well-watered. In the moderate drought (MD) and severe drought (SD) treatment groups, maize seedlings were subjected to 27 and 46 d of drought treatment, respectively. The roots of Lv28 and H082183 were sampled for RNA-seq, as described in Zhang et al. (2017) [61]. DEGs between the drought (moderate and severe drought) and control groups were obtained from Zhang et al. (2017) [61].
Salt-tolerant cultivar QM and salt-sensitive cultivar CS were used to detect responses of wheat to salt [62]. The growth conditions of the QM and CS seedlings were the same as for TAM 107. For salt treatment, 150 mmol L − 1 NaCl was added to solution. The roots of QM and CS were collected at 6, 12, 24, and 48 h after salt treatment [62]. Salt-related DEGs were obtained from previous study [62].

Insect feeding and N and cd treatment
Two-leaf stage Zhongmai 175 wheat seedlings were used for adult aphid (non-phytotoxic S. avenae and phytotoxic S. graminum) infestation [40]. Leaf samples (~ 2.5 × 2.5 cm) from aphid feeding sites were used for RNA-seq. Information on DEGs between control and treated samples were obtained from the additional files of a previously published study [40].
Two-leaf stage Wanmai No. 52 seedlings were used for N stress treatment. The roots and leaves were sampled for transcriptome analysis at 10 d after treatment [63]. We downloaded the DEGs related to N stress from a previous study [63].
For cadmium (Cd) stress treatment, 50 μM CdCl 2 was applied to Zhengmai 379 seedlings. The roots of seedlings were harvested at 12 d after Cd treatment and used for transcriptome sequencing [41]. The identified DEGs were obtained from a previous study [41].