To obtain wheat sequences of dMTases, we constructed a local protein database, which contains all of the wheat protein sequences accessible on the Ensemble database available at http://plants.ensembl.org/index.html (Howe et al., 2020). Then, the database was examined with known Arabidopsis dMTase proteins (retrieved from the TAIR database, https://www.arabidopsis.org) utilizing the BLAST_P program with an e-value of le-5 and an identity of 50%. The dMTase genes were cross-checked and redundant sequences were eradicated. Using the SMART online software program (http://smart.embl-heidelberg.de/) (Schultz et al., 2000) and the Pfam database (http://pfam.xfam.org/) (El-Gebali et al., 2018), all of the detected wheat dMTase protein sequences were also validated for the presence of conserved domains (HhH-GPD glycosylase domain; PF00730 and the RNA-recognition motif in Demeter; PF15628). The ExPASy server (http://www.expasy.org/) (Artimo et al., 2012) was used to determine the relative molecular weight (MW) and isoelectric point (PI) of TadMTase proteins.

The TadMTase genes were assigned to individual chromosomes using information available in the Ensemble database (http://plants.ensembl.org) (Howe et al., 2020). The chromosomal location of TadMTase genes was drawn using TBtools software (Chen et al., 2020). The gene structures were also analyzed using the TBtools software. The motifs of TadMTase were analyzed using the MEME suite (Bailey et al., 2009), and their positions were displayed using TBtools software.

The dMTase protein sequences of wheat (12), Aegilops tauschii (4), Triticum urartu (3), Triticum turgidum (7), Oryza sativa (4), Arabidopsis thaliana (4), Gossypium raimondii (4), Gossypium arboretum (5), Gossypium barbadense (6), and Gossypium hirsutum (10) were aligned using the MUSCLE package available at MEGA X software (Stecher et al., 2020). The phylogenetic tree was created using the neighbor-joining method with the Poisson model, pair-wise deletion, and 1,000 bootstrap values in MEGA X software. The phylogenetic tree was visualized using the iTOL online tool (Letunic and Bork, 2021). Dated phylogeny trees for 31 plant species were retrieved from TimeTree (Kumar et al., 2017).

To analyze the cis-regulatory elements (CREs) of TadMTase, the upstream sequences (1,500 bp) of the start codon were fetched from the Ensemble database (Howe et al., 2020). The CREs were identified using the PlantCARE online server (Lescot et al., 2002).

The NLS in TadMTases was predicted using the cNLS Mapper online tool (Kosugi et al., 2009). Using the CELLO online resource (Yu et al., 2014), the subcellular localization of TadMTases was predicted. The corresponding Ensemble database IDs of TadMTase genes were subjected to the ShinyGO v0.61 database (Ge et al., 2020) to obtain gene ontology (GO) annotation.

The WebSat online tool (Martins et al., 2009) was used to mine simple sequence repeats (SSRs) in the TadMTase gene sequences. The minimum length criteria of six for di-nucleotide repeats and three for tri-nucleotide repeats, tetra-nucleotide repeats, penta-nucleotide repeats, and hexa-nucleotide repeats. The Primer 3 software (available on the WebSat online tool) was used to design primers. The SSR primers were developed using the following criteria: melting temperature (Tm) of 55°C–60°C, primer length of 20–25 bp, and product size ranging from 110 to 400 bp.

The expression levels of TadMTase genes in different tissues, developmental stages, and responses to heat stress were analyzed using the Genevestigator tool (Hruz et al., 2008). The Genevestigator has a broad collection of public microarrays and RNA-Seq study data.

Expression profiling of TadMTase genes was examined in two wheat cvs. HD2329 (heat-sensitive) and HD2985 (heat-tolerant) under control and heat stress conditions. Seeds of two cultivars were surface-sterilized with 1% hydrogen peroxide, washed gently with distilled water, and germinated for 2 days in Petri plates on water-soaked filter paper at 22°C. Then the seedlings were transferred and cultivated in half-strength Hoagland nutrient solution at 22°C with a 16:8 h light/dark photoperiod. HS was induced by exposing one-week-old seedlings to a temperature of 42°C for 2 h (Meena et al., 2022). The environment for the control sample was set to 22°C. Each treatment (control and HS treated) included three biological replicates. All samples were promptly frozen using liquid nitrogen and preserved at −80°C for RNA extraction.

Total RNA was extracted from 100 mg tissues using TRizol reagent (Ambion), and DNA was removed by DNaseI enzyme (TaKaRa, United States). cDNA was generated using the RevertAid First Strand cDNA Synthesis Kit (Thermo Scientific). Quantitative expression was performed with a QuantStudio Real-Time PCR system (Applied Biosystems, United States). Three technical replicates for each sample were used. The details of gene-specific primers are provided in Supplementary Table S1. The wheat GAPDH gene was used as an internal reference gene and the relative gene expression was estimated using the 2^−ΔΔCT method (Livak and Schmittgen, 2001).

After examining the wheat reference genome, a total of 12 full-length dMTase genes were identified. The nomenclature of 12 wheat dMTase genes was based on the corresponding genes reported for rice and Arabidopsis. The 12 TadMTase genes were mainly divided into two groups (ROS1 and DML). The TaROS1a (5A/5B/5D) genes were homologous to AtROS1 (AT2G36490) and OsROS1a (Os01t0218032) in Arabidopsis and rice respectively; the TaROS1c1 (1A/1B/1D) and TaROS1c2 (1A/1B/1D) genes were homologous to AtROS1 (AT2G36490) and OsROS1c (Os05t0445900) in Arabidopsis and rice respectively; the TaDML3a (3A/3B/3D) genes were homologous to AtDML3 (AT4G34060) and OsDML3a (Os02t0496500) in Arabidopsis and rice respectively (see Supplementary Table S2). The length of the transcript of TadMTase genes ranged from 3,003 bp (TaDML3a-3D) to 7,285 bp (TaROS1a-5D), and the average length was found to be 5,377 bp. The protein length ranged from 989 aa (TaDML3a-3A) to 1982 aa (TaROS1a-5B), and the average length was 1,581 aa. The estimated molecular weights of each TadMTase protein varied from 109.6 to 218.4 kDa (Table 1). The dMTase gene names, gene IDs, chromosomal location and other features are summarized in Table 1.

The 12 wheat dMTases were located on nine (1A, 1B, 1D, 3A, 3B, 3D, 5A, 5B, and 5D) of the 21 wheat chromosomes (Supplementary Figure S1). Furthermore, a maximum of six genes were located on group 1 chromosomes (1A/1B/1D), whereas a minimum of three genes were mapped onto group 5 and 3 chromosomes. Since bread wheat is a hexaploid species (i.e., it contains three sub-genomes). Three homologues for each gene were found corresponding to chromosomes A, B, and D (Table 1; Supplementary Figure S1).

To determine the phylogenetic relatedness of dMTases in plants, 59 dMTase protein sequences from five monocots (Triticum aestivum, Aegilops tauschii, Triticum urartu, Triticum turgidum, and O. sativa) and five dicotyledons (Arabidopsis thaliana, Gossypium hirsutum, Gossypium arboretum, Gossypium raimondii, and Gossypium barbadense) were used for multiple sequence alignments and to build a phylogenetic tree (Figure 1A). All the members of dMTase are grouped into three subfamilies/clades: ROS1, DME, and DML. The ROS1 clade comprised 22 monocots and nine dicots, the DML clade comprised eight members of monocots and 13 dicots, while the DME comprised seven members of dicots only (Figure 1A). Furthermore, we also constructed an evolutionary tree of life for 31 plant species belonging to monocot (8) and dicot (23) and showed the total number of dMTase proteins for each species (Figure 1B). Our analysis revealed that the DME subfamily of dMTases was absent in monocots and existed only in dicots.

To gain further insights into the TadMTase genes, we surveyed the gene structure, conserved domain, and motif components of each TadMTase gene. The number of exons in the TadMTase genes ranged from 17 to 22 (Figure 2A; Table 1). We found that the number of exons in dMTase genes from the same homoeologous group did not differ much (maximum up to three exons in TaROS1c2). Among these, the dMTase gene with the maximum number (22) of exons was found to be TaROSc2-1B, while TaDML3a had the least number of exons (17).

The distribution of conserved motifs in the TadMTase proteins was predicted using the MEME online suite. Twenty different types of motifs, which were designated 1 to 20, were predicted in each TadMTase protein. The amino acid length of motifs varied from 21 to 50. The number of motifs in each dMTase varied from 13 to 20. Among them, motifs 3 and 9 were part of Domain A, motifs 1, 2, 6, 8, and 13 were part of the HhH-GPD glycosylase domain, and motifs 4, 5, 7, 10, 11, and 19 were part of the RRM_DME domain (Figure 2B). Furthermore, the Motifs 1, 2, 3, 4, 5, 10, 11, 13, 16, 18, and 19 were highly conserved and found in every dMTase protein. However, some motifs (9, 15, and 20) were specific to the ROS gene group (Figure 2B). It could be inferred that motifs found in all dMTase proteins are possibly associated with conserved functions, but those specific to the ROS gene group may be involved in gene-specific functions.

To further determine the structural characteristics of the wheat dMTase family members, conserved domains were also identified. All the wheat dMTases contained three different domains, including 1) Domain A, 2) HhH-GPD glycosylase domain (PF00730) and 3) Domain B or RRM_DME (RNA-recognition motif in Demeter; PF15628) domain (Figure 3). The domain A contains the stretch of basic amino acids (K and R) that is required for DNA binding activity. The HhH-GPD glycosylase domain contains a helix–hairpin–helix (HhH) motif, a glycine/proline-rich loop with a conserved aspartic acid (GPD), and four cysteine residues in the Fe-S cluster, these are involved in 5-methylcytosine excision activity. The RRM_DME domain facilitates the interaction of the catalytic domain with ssDNA or regulatory RNA. The multiple sequence alignments of wheat dMTases showed that these three domains are conserved (Supplementary Figure S2).

Out of 12 wheat dMTase gene sequences screened for SSRs, ten genes have 17 SSR motifs. The tri-nucleotide repeats motif were found most abundant, comprising about nine (53%) followed by the tetra-nucleotide repeat motifs comprising four (23.5%), whereas, di- and penta-nucleotides SSR formed a small share having 2 (11%) each (Supplementary Table S3). The number of repetitions of a motif ranged from 3 to 19. For each of the 17 SSR motifs, primer pairs were also designed. The following details, such as nucleotide sequence, melting temperature, and product size related to SSR motif primers, are provided in Supplementary Table S3.

All TadMTase members contain both monopartite and bipartite NLSs (except TaROS1a-5A). In total, 12 dMTases were predicted with 35 NLS (21 mono and 14 bipartite NLSs) (Supplementary Table S4). Subcellular localization prediction indicated that all the 12 wheat dMTase proteins are localized in the nucleus (Table 1).

To further understand the possible regulatory function of dMTases in plant growth and stress response, we searched for CREs in the 1,500 bp promoter region of wheat dMTase genes. A total of 158, representing 29 types of CREs, were predicted (Figure 4). Most of the identified CREs were hormone response factors (37.34%), followed by light-responsive elements (25.94%), stress response elements (22.15%), and plant development-related elements (14.55%). The CREs involved in hormonal signaling comprised auxin-responsive (TGA-element), abscisic acid-responsive (ABRE), methyl jasmonate responsive (CGTCA and TGACG motif), gibberellins responsive (TATC-motif) and salicylic acid-responsive (TCA-element) elements. Among these CREs, the CGTCA, ABRE, TGACG, and TGA motifs appeared to be the most prevalent and were detected in all TadMTase genes. In TadMTase promoters, all 13 kinds of light-responsive CREs were also found. Out of those, G-box CREs appeared to be the most prominent. Each TadMTase gene had one to three G-box copies (Figure 4A).

Thirty-five stress-responsive CREs were found in TadMTase genes and out of those, ARE (anaerobic responsive elements) elements comprised 28% of them (10). AREs were most commonly detected at the TadMTase gene promoters. The CAAT-box, CCAAT-box, GC-motif, LTR (low temperature responsive) and MBS were among the other stress-related CREs detected. The GC-motif and CAAT-box have been linked to a variety of abiotic stressors. MYB binding sites (MBS) are known to control drought stress, ARE are known to regulate anaerobic induction under oxidative stress, and LTRs regulate the temperature response.

Furthermore, 24 CREs were discovered that were involved in the growth and development of plants. These included A-Box (5), CAT-box (14), motif I (3), and O2 site (1) (Figure 4A). The CAT-box element appeared to be the most frequent and was detected in ROS and DML dMTase genes. We also found Motif-I, a root-specific regulatory element only in the ROS sub-group genes. Overall, these findings reveal that TadMTase genes may be involved in a variety of abiotic stress responses as well as plant growth regulation.

To further characterize the function of the wheat dMTase genes, we performed GO enrichment analysis using the ShinyGO tool (Ge et al., 2020). GO enrichment analysis showed that the TadMTase genes were involved in DNA demethylation, the BER pathway, DNA modification, developmental processes, and stress responses in GO biological processes (Figure 5A). The GO molecular process showed that TadMTase genes were involved in DNA demethylase activity, metal/Iron-sulphur cluster binding and catalytic activity (Figure 5B). This annotation indicates that the TadMTase genes may be involved in stress regulation in plants via DNA demethylation.

To decipher the function of wheat dMTase genes, the expression pattern in ten developmental stages and 21 tissues was analyzed (Figure 6). The expression levels indicate that all the wheat dMTase genes were expressed in different developmental stages/tissues but varied significantly. During various developmental stages, the expression of dMTases gradually increased from the seedling to the booting stage and declined continuously till the ripening stages. However, TaROS1a-5A/5B/5D had the highest expression at anthesis stages, and TaROS1c1/2 1A/1B/1D and TaDML3a 3A/3B/3D had the highest expression at the stem elongation stage (Figure 6A). In the case of tissue-specific expression, TaROS1a 5A/5B/5D, TaROS1c2 1A/1B/1D and TaDML3a 3A/3B/3D exhibited relatively higher expression in all tissues. However, TaROS1c1 1A/1B/1D had tissue-specific expression and was expressed especially in the shoot apex, pericarp, and blade (lamina) tissues (Figure 6B). Furthermore, we also studied the expression of TadMTase genes in hormone-treated [abscisic acid (ABA) and gibberellic acid (GA3)] wheat tissues. TaROS1a-5A/5B/5D were downregulated (up to 1.96 fold) in root tissues at the seedling stage upon ABA treatment. Other genes do not show any significant change in their expression pattern (Supplementary Figure S3).

Heat stress is one of the major stresses that severely affects wheat productivity. Therefore, there has been recent research to understand the mechanism by which plants respond to HS and identify the genes associated with HS tolerance (Samtani et al., 2022). To investigate the putative biological role of TadMTase under HS, the expression pattern of these genes under HS was determined using wheat expression data available in the Genevestigator tool and by performing the qRT-PCR experiments. The wheat RNA-seq data (available on the Genevestigator tool) of the four different experiments when were subjected to ambient (22°C) and heat stress (37°C–40°C) treatment for 5 min to 5 days (as shown in Figure 7A) were evaluated. Results showed that there were significant differential expression patterns and most of the TadMTase genes were suppressed (up to 2.24 folds) in response to HS as compared to controls (Figure 7A). For instance, in response to HS, TaROS1a-5A, -5B, and -5D, and TaROS1c2-1A, -1B, and -1D were downregulated in the seedling stage (leaf and caryopsis tissue; up to ∼2.3 fold) and TaDML3a-3A, -3B, and -3D were downregulated in the anthesis stage (flag leaf tissue; up to ∼2 fold).

Furthermore, we selected five genes (TaROS1a-5A, TaROS1c1-1B, TaROS1c2-1A, TaDML3a-3B, TaDML3a-3D) out of the 12 TadMTase genes for quantitative real-time PCR analysis. For this purpose, contrasting wheat genotypes (sensitive cv. HD2329 and tolerant cv. HD2985) were used. All the five TadMTase genes showed considerable differential expression after being stimulated by HS (Figure 7B). The expression levels of TaROS1a-5A and TaDML3a-3B were repressed (up to 2.7-fold) in both the genotypes. The expression levels of TaROS1c2-1A and TaDML3a-3D were repressed (up to 2.12-fold) in HS tolerant cv. HD2985 and upregulated in HS sensitive cv. HD2329. The expression level of TaROS1c1-1B was upregulated in HS tolerant cv. HD2985 and repressed in HS sensitive cv. HD2329 (Figure 7B). The above results indicate that active demethylation is possibly involved in wheat’s responses to thermal stress. Furthermore, the differentially regulated TadMTase might regulate thermotolerance in crops such as wheat.