Arabidopsis MT, FIP and ALKBH cDNA and protein sequences were downloaded from the ensemble (https://plants.ensembl.org/index.html) using the gene ID provided in the literature (Ougland et al., 2015). The protein sequence was then used for blast search in the Legume Information Database (https://www.legumeinfo.org/) to see the orthologous proteins in pigeon pea. E-value threshold was kept at zero for blast search with 98-100% coverage.

The different Expasy (https://www.expasy.org/) tools like ProtParam, compute pI/Mw etc was used to have a basic understanding of the identified genes in terms amino acid length, molecular weight, iso-electric point, GRAVY (Kyte and Doolittle, 1982), instability index and aliphatic index.

The chromosomal position was identified from the LIS database (https://www.legumeinfo.org/) and subsequently, the locus ID was also identified from the NCBI database (https://www.ncbi.nlm.nih.gov/). Chromosome map was constructed using MapGene2Chrom web v2 (http://mg2c.iask.in/mg2c%5Fv2.1/).

A phylogenetic tree of the identified proteins was constructed to see their relative closeness. MEGA11 software was used for the phylogenetic tree construction (Tamura et al., 2021). First of all, after identification and retrieval of all the sequences Clustal omega (https://www.ebi.ac.uk/Tools/msa/clustalo/) was used to check for similarity among sequences. Further of full-length amino acid sequences of Arabidopsis thaliana, Oryza sativa, Glycine max and C. cajan were fed to MEGA and there again multiple sequence alignment was performed with ClustalW tool. The IDs for Oryza sativa and glycine max is provided in Table 1 .While aligning the sequences of four species of crops in MEGA, alignments were made selecting “with gap option” and during construction of the phylogenetic tree, gap parameters were selected as ‘Use all site’. Phylogenetic tree was constructed using the Maximum Likelihood method and JTT matrix-based model taking bootstrap value 1000 (Tamura et al., 2021). For visualization of the phylogenetic tree an ‘Interactive Tree Of Life’ (iTOL) v6 (https://itol.embl.de/) was used.

For identification of gene structure GSDS 2.0 was used (http://gsds.gao-lab.org/). This gave the idea of exon-intron structure in MTs, FIPs and ALKBHs genes. The conserved motifs of the protein were examined using the MEME online software tool (https://meme-suite.org/meme/). The motif number was kept as 20. The motif width range was kept as minimum 6 and maximum of 50 (6-50) and in site distribution zero or one occurrence per sequence was selected.

The conserved domain of the genes was predicted using an online ‘CD Search tool’ (https://www.ncbi.nlm.nih.gov/Structure/bwrpsb/bwrpsb.cgi). The sub-cellular localization of pigeon pea MTs, FIPs, and ALKBHs was predicted using the WoLF PSORT web tool (https://wolfpsort.hgc.jp/).

Upstream genomic sequences of 2 kb from transcription start site (including 5’ UTR) were retrieved from the LIS database for MTs, FPIs and ALKBHs. Cis-regulatory elements were identified using the Plant Pan v3.0 (http://PlantPAN.itps.ncku.edu.tw/). Data obtained from the web tool was analyzed in MS Excel V.2013 and visualized in the TB tool (https://bio.tools/tbtools). Upstream sequences for MTs, FPIs and ALKBHs were analyzed to predict the m⁶A- methylation using the EpiSemble R-package v.0.1.1 (http://cabgrid.res.in:5799/). MethSemble 6mA tool was used in EpiSemble R-package to predict the methylation site. This package uses three models viz gradient boosting, random forest and Support vector machine.

Pigeon pea genotype, Pusa 992, was grown in the net house in a 4-inch pot (loamy soil) under natural day length (14hr light and 10 hr. dark) and temperature(30-32 °C) in July 2023 at ICAR-NIPB, New Delhi. Initially watering was done on every alternate day up to three leaved stage. After that when soil used to dry based on that watering was done. Plants were grown in triplicates. For heat stress, 20 days old seedlings are kept in a heat chamber at 42 °C and 60% relative humidity. Plants were kept for 6 hrs. (from 11:00 am to 5:00 pm) for two days ( Supplementary Figure 1 ). After the second day leaf samples were taken. For drought stress, 20% PEG 6000 was prepared by adding 200gm of PEG 6000 in 1000ml of autoclaved water and 100ml of PEG was given per pot which contained a single plant. ( Supplementary Figure 1 ). For salt stress 150 mM of NaCl was prepared by adding 8.766 gm in 1000 ml autoclaved water. and 100 ml of the solution was given to 20 days old seedlings (single plant in one pot) (Mi et al., 2024; Dokka et al., 2024) ( Supplementary Figure 1 ). Helicoverpa armigera was used for biotic stress. Larvae were obtained from an in-house culture facility The larvae of H. armigera larvae were raised on an artificial diet with a 16 h light and 8 h dark photoperiod, at a temperature of 26± 1°C and 70-80% relative humidity was maintained. The second instar larvae of the polyphagous insect pest, Helicoverpa armigera, was reared on leaves of 20 days old seedlings of pigeon pea. Larvae were given 7 hrs. starvation and then two larvae were released per pot (1 Plant in each pot). Plants were covered with the perforated polythene which was secured with rubber band on the pots to prevent escape of larvae. After four days of infestation leaf samples were collected.

For tissue-specific expression studies different tissues were collected at different stages, but for stress related studies leaf tissues were collected from 20 days old seedlings. For abiotic stress-induced plants leaf samples were collected after 48 hours of treatment. Leaf samples were collected after 4 days of H. armigera infestation for gene expression study under biotic stress. Total RNA was isolated from different tissues (leaf, roots, internode, shoot apical meristem, flower apical meristem and immature pod) for tissue-specific qPCR and from leaf samples for stress-specific qPCR using RNA isolation kit (Genes2Me; India) according to the ‘manufacturer’s instruction. Isolated RNA was treated with DNaseI (RNase free) to remove any genomic DNA contamination. The quality of total RNA was checked using a Nanodrop spectrophotometer (Thermo Scientific). Total RNA was then immediately stored at -80°C. cDNA was prepared using a PrimesScript cDNA synthesis kit (TaKaRa) and stored at -20 °C for further use. MTs, FIPs and ALKBHs-specific primers ( Supplementary Table 2 ) were designed using the IDT web tool (https://www.idtdna.com/). qPCR assay was performed in Light Cycler 96 PCR detection system (Roche, Basel, Switzerland) using TB green master mix (TaKaRa) using the following conditions: initial denaturation at 95 °C for 5 min, 40 cycles of amplification, each cycle of 95 °C for 30 sec, 60 °C for 30 sec and 72 °C for 20 sec. Also, three biological and three technical (cDNA replicates) replicates were taken for each sample. The CcIF4 was used as a reference gene (Bhattacharjee et al., 2023). The Sequence of the internal primer pair for the reference gene is included in the Supplementary File ( Supplementary Table 2 ). The relative abundance of CcMTs, CcFIPs and CcALKBHs was calculated using the 2^–ΔΔCt method (Livak and Schmittgen, 2001).

For the gene expression study, three biological (separate plants grown in separate pots under different abiotic and biotic stress) and three technical replicates were taken. Wherein, an equal amount of each biological replicate was pooled for RNA isolation. Mean values were given with an error bar (standard error of means) for all the parameters. At 5%, the least significant difference (LSD) was calculated to see the significance of different treatment effects. After that, the significance level between and among the treatments in each experiment was checked by performing a range test.

The sequence information of methyl transferase and demethylase was retrieved from the LIS database and cross-checked through blast at NCBI database, which gave the Locus ID of respective genes. Arabidopsis has 2 MT, 1 FIP and 14 ALKBH genes, whereas in pigeon pea, 2 MT, 2 FIP and 10 ALKBH genes were identified ( Table 2 ).

This exercise provided a complete framework of basic information on iso-electric point (pI), molecular weight (MW), instability index (II), aliphatic index (AI) and grand average of hydropathicity (GRAVY) of the proteins identified in pigeon pea ( Table 3 ). When the protein sequence of identified MTs, FIPs and ALKBHs of pigeon pea was analyzed it was found that there was variation in the genes. For instance, in case of predicted protein length, the amino acid sequence varies from 761 aa to 1089 aa for MTs, 337 aa to 338 for FIPs and 205 aa to 515 aa for ALKBHs family of pigeon pea. The iso-electric point for MTs ranged from 6 to 7, and for FIPs, it was between 5 and 6. The iso-electric point ranged from 5.57 to 8.7 for ALKBHs with CcALKBH1C having the highest PI and CcALKBH8B the lowest. All proteins of MTA, FIPs and ALKBH were hydrophilic as confirmed by GRAVY. Also, instability index analysis showed that MTA was more stable than MTB, and FIPA was more stable than FIPB. Among ALKBH proteins CcALKBH8 was the most stable, and CcALKBH8B was the least stable protein. Aliphatic index analysis indicated that MTA and FIPA had more aliphatic amino acids compared to MTB and FIPB, respectively. CcALKBH8 had the highest aliphatic index and CcALKBH9 had the lowest aliphatic index. The higher aliphatic index, the better the thermo-stability of protein.

From the analysis of the chromosomal position of the genes encoding pigeon pea MTs, FIPs and ALKBHs it was found that all the genes were localized within six chromosomes, viz. chr.01, chr.02, chr.03, chr.05, chr.06 and chr. 11. However, the majority of the genes were localized on the chr.03 ( Table 4 ).

A chromosomal map had been constructed showing the distribution of genes on chromosomes. CcMTA and CcMTB were located on chromosomes 01 and 02, respectively. CcFIPA and CcFIPB were located on chr.03 and 11, respectively. Two ALKBH genes, CcALKBH1A and CcALKBH1B, were located on chr.01, one ALKBH gene, CcALKBH10A, was on chr.02, and four genes, viz, CcALKBH1C, CcALKBH2, CcALKBH8 and CcALKBH10B, were onchr.03. CcALKBH8B, and CcALKBH9 were present on chr.05 and CcALKBH8A was found on chr. 06 ( Figure 1 ).

Deduced protein sequences of MTs, FIPs and ALKBHs from Arabidopsis (A. thaliana), rice (O. sativa), soybean (G. max) and pigeon pea (C. cajan) were taken and a phylogenetic tree was constructed using MEGA11 to find out the relationships between the identified genes and to see the evolutionary relics (Tamura et al., 2021) ( Figure 2 ). The tree sub-clades were clubbed into groups to understand their evolutionary relations. One group for MTs and FIPs, viz. MTA/B and FIPA/B, respectively, and four groups for ALKBHs were made, viz, (ALKBH1/2, ALKBH8, ALKBH9, ALKBH10) following earlier nomenclature (Marcinkowski et al., 2020). The number of genes of MTs and FIPs was almost the same in above mentioned species. However, the number of ALKBH genes varied among species. The highest number of ALKBH genes (14) was found in Arabidopsis, while the least number of ALKBHs (10) was found in pigeon peas. The ALKBH6 group was found absent in pigeon pea. The ALKBH1 group had the highest number of genes (4), while ALKBH2 and AlKBH9 had the lowest number of genes (1). ALKBH10 was similar and related to m⁶A RNA demethylation in Arabidopsis.

Gene structure plays a major role in the evolution of gene families. A phylogenetic tree was constructed using the neighbor joining method grouped CcALKBs into 4 paralogous clades. The members of CcALKBH1, CcALKBH2 and CcALKBH8 were 3 distinct clades whereas, the members of CcALKBH9 and CcALKBH10 together grouped as a separate clade ( Figure 3A ). Analysis of the exon-intron structure of MTs revealed that seven and six exons were present in MTA and MTB, respectively, but intronic portion was more in MTB ( Figure 3B ). In case of FIPA and FIPB there were eight and thirteen exons, respectively, and for FIPA, UTR was found only at the 3’ end ( Figure 3B ). Intronic portion was found to be more in FIPB. Further, for ten ALKBHs, it has been found that variations were present among the genes. Seven exons were present in five of the ALKBHs, viz. CcALKBH1A, CcALKBH8A, CcALKBH9, CcALKBH10A and CcALKBH10B; and the rest five ALKBHs, viz. CcALKBH2, CcALKBH1B, CcALKBH1C, CcALKBH8 and CcALKBH8B, had varying numbers of exons (five to one, respectively). CcALKBH10A and CcALKBH10B contained the largest intronic regions. It was noticed that the CcALKBH2 gene was devoid of any UTR region ( Figure 3B ).

The MEME suite was used for conserved motif analysis, and 20 motifs were identified. According to the phylogenetic analysis ( Figure 4A ), motif distribution was found conserved for closely related genes as shown in Figure 4B . All 20 motifs were found to be present in both the MTs, but their distance varied. FIPs had all the motifs conserved and at the same distance. For ALKBHs, CcALKBH10A and CcALKBH10B had the greatest number of genes conserved at the same distance ( Figure 4B ). Motif 1 was found to be conserved in all the genes except for CcALKBH2 and CcALKBH8. Motif 20 was specific to CcALKBH8A and CcALKBH8B. Similarly, motif 18 was specific to CcALKBH1B and CcALKBH1C, and motif 19 was specific to CcALKBH2 and CcALKBH8 ( Figure 4B ). CcALKBH2 and CcALKBH8 had shown the least conservation of motifs which might be an indication of performing different functions.

Domains are the self-stabilizing polypeptide chain that works independently. It was observed that MTs (CcMTA and CcMTB) contained domain MTA-70 which is a major domain involved in methylation ( Figure 5 ). Microbial surface components recognizing adhesive matrix molecules domain was found in FIPs, viz. CcFIPA and CcFIPB ( Figure 5 ). Most of the ALKBHs contained 2OG-FeII_Oxy_2 domain which is a characteristic feature of the ALKBH gene family. Although the 2OG-FeII_Oxy_2 domain was absent in CcALKBH2, it had a completely different domain of DUF4057 superfamily, which is yet to be characterized ( Figure 5 ). Apart from this, CcALKBH8A was found to have an RRM (RNA recognition motif) domain, which is an essential domain involved in tRNA modification. In addition, CcALKBH8 also had a methyl transferase domain. So, it might have a role in both demethylation and methylation.

Regarding subcellular location, it was found that MTs, FIPs and a majority of ALKBHs had nuclear localization signal. However, a few ALKBHs, viz. CcALKBH1B, and CcALKBH10B, had chloroplast targeting signals, CcALKBH8 had signal peptide sequence targeting the plasma membrane ( Supplementary Table 1 ), and CcALKBH9 had both nuclear and cytoplasm localization signals.

The cis-regulatory elements in the promoter region of Pigeonpea MTs, FIPs and ALKBHs genes were predicted using the 2kb upstream sequence retrieved from the available database for pigeon pea (LIS database). The identified cis-elements were then selected based on their role in growth and development, hormone response and stress. Growth and development regulatory elements like ARID, AT-Hook, Dof, NAC, LOB, SBP, HD-ZIP, PLATZ and FAR1 were selected. AT-Hook, B3, BBR-BPC, BES1 were selected for hormone response and AP2, bHLH, MADS Box, GATA, WOX, WRKY, C3H-Zinc finger, Dehydrin and VOZ were selected for stress response. Cis-elements varied between genes based on the presence and absence and also on the frequency by which they appear. For instance, both the MTs, viz. MTA and MTB, had almost equal number of cis-elements ( Figure 6 ). MTs had the highest number of AP2 binding sequences followed by Dof ( Figure 7 ). Between the two FIPs, FIPA had a greater number of elements as compared to FIPB. CG, FAR1 and HD-ZIP binding sequences were absent in FIPB but present in FIPA ( Figure 6 ). In case of ALKBHs, CcALKBH9 had the highest number of cis-regulatory elements followed by CcALKBH8 and CcALKBH1B ( Figure 6 ). MYB cis-regulatory binding elements were found to be the highest among ALKBHs followed by cis-regulatory binding elements for GATA, bZIP and Dof. ( Figure 7 ). A table of cis-elements with their numbers for all three genes, viz, MTs, FIPs and ALKBHs, were provided in the Supplementary Files ( Supplementary Table 3 ).

EpiSemble R-package was used to predict the m⁶A methylation pattern in MTs, FIPs and ALKBHs. This exercise was carried out to understand the epigenetic regulation of the genes. In case of MTs five and four methylation sites were found in the upstream 2 kb region ( Figure 8 ). For FIPs, three and two sites were found for FIPA and FIPB, respectively ( Figure 8 ). Further in case of ALKBHs, CcALKBH9 and CcALKBH10A had the highest number of methylation site (six), but the lowest methylation site was found for CcALKBH1B and CcALKBH10B (two) ( Figure 8 ).

Quantitative polymerase chain reaction (q-PCR) was performed to understand the expression pattern of the identified MTs, FIPs and ALKBHs in pigeon pea. Different tissues (leaf, root, internode, shoot apical meristem, flower apical meristem and immature pod) were checked for the relative abundance of the transcripts ( Figure 9 ). In case of MTs, it was found that overall expression of CcMTA was higher in the six selected tissues compared to that of CcMTB. The highest expression for MTA was observed in leaf tissues (~4.3 fold), while the highest expression for MTB was detected in FAM tissues (~3.7 fold). A similar kind of expression pattern was observed for these two genes in other tissues (root, internode, SAM and immature pod), but with varied expression levels i.e., CcMTA (~2.5 fold) had significantly higher expression compared to that of CcMTB (~1.0 fold) in root tissues. But, CcMTB (~3.7 fold) had more expression in SAM tissues compared to that of CcMTA (~3.1 fold) ( Figure 9 ). In case of FIPs, both the genes, viz. CcFIPA and CcFIPB, showed the highest expression in leaf and internode. However, higher expression of CcFIPB was detected in the leaf (~6.0 fold) and root (~4.0 fold), and relatively more expression of CcFIPA was detected in FAM (~4.4 fold) tissues ( Figure 9 ).

Majority of the genes encoding ALKBHs displayed similar kind of expression patterns with the highest level of expression in leaf tissues, except for CcALKBH8 (~5.9 fold), CcALKBH9 (~4.1 fold) and CcALKBH10A (~4.9 fold), which showed the highest expression in internode tissue. The second highest expression of seven ALKBH genes was also detected in internode tissue, but three genes, viz. CcALKBH1A (~4.5 fold), CcALKBH1C (~4.3 fold) and CcALKBH8A (~4.5 fold) showed the second highest expression in FAM tissue. Among the ten ALKBH genes, CcALKBH10B had the highest expression in all the six tissues analyzed and CcALKBH2 had the lowest expression ( Figure 9 ). Overall, the highest level of expression of genes encoding MTs, FIPs and ALKBHs was detected in leaf and the lowest expression in root tissues ( Figure 9 ).

We wanted to see the expression level changes in the identified genes under various abiotic and biotic stresses. So, we subjected pigeon pea seedlings under various stresses and the morphological changes which was found is provided in Supplementary Figure 1 Further relative expression of MTs, FIPs and ALKBHs genes of pigeon pea during biotic and abiotic stresses was studied by qPCR analysis. During heat stress the highest induction in expression was observed in CcALKBH8 (~9.5 fold) followed by CcALKBH10B (~8.0 fold), but no induction in expression was found in CcALKBH2 (~1 fold) compared to that of control ( Figure 10 ). Among the methyl transferase genes induction in expression was not so prominent, and about two-fold induction in expression was detected for CcMTA (~2.1fold) and CcMTB (~1.9 fold), whereas very low induction was observed for CcFIPB (~1.2 fold) and CcFIPA (~1.0 fold) during heat stress ( Figure 10 ). Under drought stress, CcALKBH10B showed nine-fold more expression, followed by CcALKBH9 (~7.5 fold) and CcALKBH10A (~7.3 fold). CcALKBH2 (~1.0 fold) showed negligible induction ( Figure 10 ). In case of methyl transferase genes, about four-fold induction in expression was observed in CcMTA (~4.3 fold) and CcMTB (~4.0 fold), but induction was not so prominent in CcFIPB (~1.1 fold) and CcFIPA (~1.0 fold) ( Figure 10 ). The highest level of induction in gene expression of 13-fold was detected in CcALKBH10B (~13.3 fold) during salt stress. Two other genes, CcALKBH10A (~7.6 fold) and CcALKBH9 (~5.7 fold) showed 8- and 6-fold induction, respectively, during salt stress. Whereas, CcALKBH2 (~1) showed negligible induction ( Figure 10 ). Two methyl transferase genes, CcMTB (~5.5 fold) and CcMTA (~5.3 fold), showed about five-fold more expression during salt stress. Again, CcFIPB (~1.6 fold) and CcFIPA (~1.0 fold) showed very little induction in expression during salt stress ( Figure 10 ).

A higher level of induction in gene expression, ranging from 13 to 9-fold, was detected in ALKBH genes in pigeon pea upon H. armigera infestation. The highest induction was observed in CcALKBH10B (~13.5 fold), followed by CcALKBH10A (~12.2 fold), CcALKBH9 (~9.4 fold) and CcALKBH1C (~8.6 fold). Again, CcALKBH2 (~1.0 fold) showed negligible induction during biotic stress ( Figure 10 ). Less pronounced induction was observed for MTs genes with about four-fold induction in CcMTB (~4.0 fold) and CcMTA (~3.8 fold) followed by 2-fold induction in CcFIPB (~1.9 fold). However, induction in CcFIPA (~1.0 fold) was not significant ( Figure 10 ).