Two chickpea genotypes, namely KWR 108 and GL 13001, were selected for the study as these two varieties were reported to be highly contrasting in response to Fusarium wilt (Kumar et al., 2019), particularly to Foc4. The experiment was conducted in the Fusarium wilt sick plot of TCA, Dholi farm (a designated plot for chickpea Fusarium wilt pathological trials from the All India Co-ordinated Research Projects on Chickpea located at TCA-Dholi under the Department of Agricultural Biotechnology and Molecular Biology, CBS&H, Dr. Rajendra Prasad Central Agricultural University, Pusa (Samastipur) 848,125, Bihar, India). F. oxysporum isolates were collected from the field of TCA, Dholi and identified as Foc4 (Sharma et al., 2014). Among these chickpea genotypes, KWR 108 was reported as highly resistant and GL 13001 was reported as a highly susceptible genotype against Fusarium wilt (Kumar et al., 2019). The seeds of KWR 108 and GL 13001 genotypes were procured from the ICAR-Indian Institute of Pulse Research, Kanpur, Uttar Pradesh, India. Seeds of both genotypes were sown in a grow bag filled with a mixture of autoclaved soil and compost at a ratio of 6:4. For each genotype, two seeds were sown in each grow bag. After germination, only one plant bag⁻¹ was maintained for 25 days after sowing. Then plants of both genotypes were divided into two groups, namely control and treatment, which had three biological replicates. All the field experiments were carried out under the net house condition during rabi season 2021–2022, and the molecular works were carried out at the Functional Genomics Laboratory of the Department of Agricultural Biotechnology and Molecular Biology (25.98°N latitude and 85.67°E longitude and 52 m altitude from mean sea level), College of Basic Sciences and Humanities, Dr. Rajendra Prasad Central Agricultural University, Pusa 848,125, Bihar, India during the year 2021–2022.

For infecting the plants, broth culture was prepared from the pure culture plates of Foc4 (pure culture was derived from the wilt-infected chickpea plant grown at the wilt sick plot of TCA, Dholi). Small-size inoculums (0.5 cm²) from the pure culture plate were inoculated in broth culture media (Potato Dextrose Broth) with the help of a forcep. The broth culture was incubated at 26°C and 150 rpm in a rotary agitator for 7 days. For inducing wilt infection in the treatment group of plants, the broth culture was thoroughly mixed using a pipette, and 2 mL of the solution was added to each plant at the root zone with the help of a pipette. To further ensure the infection of plants in treatment groups, soil collected from the Fusarium wilt sick plot of TCA Dholi was filled on the top of the grow bags (treatment), followed by light irrigation. Symptoms appeared in the susceptible genotypes within 8–10 days of the treatment (Figures 1A,B). The plant leaf samples were collected after 10 days of treatment for transcriptome sequencing.

The leaf samples were collected from control and treated plants of both the genotypes in three biological replicates and immediately dipped in the 2 mL micro-centrifuge tubes filled with RNA stabilizer solution (G-Biosciences) and outsourced for RNA isolation and transcriptome sequencing at Unigenome, Ahmedabad, Gujarat, India-380015. The leaf samples from three biological replicates were bulked together for RNA isolation. RNA was isolated by using the Alexgen Total RNA kit. RNA quality and quantity were analyzed using 1% agarose gel and Qubit^® 4.0 fluorometer, respectively. The NEBNext^® Ultra^™ RNA Library Prep Kit for Illumina NEB #E 7770 was used for paired-end sequencing library preparation. The mRNA enrichment was performed as per the user manual and mRNA was first fragmented, followed by cDNA synthesis and repairing of ends and adenylation of 3′ end. Subsequently, the adapter was ligated at the fragmented ends and PCR amplification was performed for selective enrichment of adapter ligated DNA fragments. The adapter-based selectively enriched amplified fragments/libraries were analyzed on TapeStation 4,150 (Agilent using RNA ScreenTape^® as per manufacturer’s instructions). The average size of libraries was observed to be 331 bp, 363 bp for GL 13001 (control and treatment) and 353 bp and 344 bp for KWR 108 (Control and Treatment), respectively. The final library was pooled with other samples, denatured, and loaded onto the flow cell. On the flow cell, cluster generation and sequencing were performed using the Illumina Novaseq 6,000 platform to generate 2 × 150 bp paired-end (PE) reads. Data were filtered to remove the adapter and low-quality reads. The raw transcriptome sequencing data obtained from both genotypes under control and F. oxysporum f. sp. ciceris-treated plants with two technical replicates were submitted to the SRA (Sequence Read Archive) database of NCBI (accessed on 20 October 2023).¹

From the NCBI database (accessed on 20 October 2023),² the reference genome of Cicer arietinum (GCA_000331145.1_ASM33114v1) and reference genome annotation GTF file were downloaded. Subsequently, by using STAR (v2.7.10a) aligner, reference genome-based assembly of transcripts was performed first by separately mapping high-quality (HQ) clean reads of control and treatment samples of GL 13001 and KWR 108 on the reference genome, respectively. The STAR Genome-Generate mode option was used to carry out indexing of the reference genome. Then the input reads, in FASTQ format, along with the indexed reference genome generated in the previous step, were given to the STAR aligner. To identify the positions of the origin of the reads, the STAR aligner created the alignment in BAM format for each sample through reference genome-guided mapping of the HQ reads followed by assembly of transcripts using StringTie (v2.2.1). After the assembly of transcripts by StringTie, the quantification of assembled transcripts was performed. The quantification of the assembled transcripts was done in order of most highly abundant transcripts to less abundant by using the network flow algorithm. The annotation of assembled transcripts was performed using the GFF annotation file of reference genome having “known” gene exon structures followed by quantification of the expression of known genes. Simultaneously, to merge the gene structure found in all four samples, the assembled transcripts of all four samples were passed to StringTie’s merge function to avoid assembly of only a partial version of any of the transcripts in the initial StringTie run.

The reference GTF file and the string-tie merged GTF file were taken and the GFF compare utility was run to check or identify if there were any novel transcripts or novel isoforms of transcripts were present in the merge sample. Novel transcripts were then identified using a method described by Gleeson et al. (2022). The transcripts that did not have any available information in the annotation file of the reference genome were considered novel transcripts.

Using the SringTie merge function, a consensus set of transcripts was derived by merging the structures of all the genes found in all four samples and simultaneously estimating the abundances of the merged transcripts in the four samples. The read count information that was obtained from the files generated by StringTie was extracted by using prepDE.py, a Python program (Gill and Dhillon, 2022). The edgeR package (Robinson et al., 2010) was then used to calculate differential gene expression (DEGs) by taking read count results of a Python program prepDE.py as input.

Gene ontology (GO) analysis was done to identify the distribution of significant differentially expressed genes (DEGs) into three major domains, i.e., cellular component, molecular function, and biological process. GO term was assigned to each significant differentially expressed transcript using Blas2GOcli 1.4. KAAS (KEGG automatic annotation server; Moriya et al., 2007) was used for ortholog identification and mapping of significant differentially expressed transcripts to the biological pathways. The criteria used to identify significantly down and up-regulated transcripts were log2FC > 0 and q value<0.05 (means significantly up-regulated), and log2FC < 0 and q value <0.05 (means significantly down-regulated). Using BLASTX with a threshold bit-score value of 60 (default), significant differentially expressed transcripts were compared against the KEGG database. A Heatmap of DEGs was generated for 25 highly significant down and up-regulated genes (a total top 50 transcripts) based on q value (i.e., least q value).

The reference of the genome-based transcript identification revealed that many of the transcripts showed differential expression under Fusarium wilt infections, which are the transcript fragments encoding the same proteins. Further analysis was conducted using the NCBI gene database. However, the protein IDs of the transcripts revealed that several proteins were the products/isoforms of the same genes. Using the sorting option available in MS Excel, the transcripts encoded by the same genes were sorted together to finally get the exact numbers of genes showing differential expression under Fusarium wilt stress in chickpea.

From the RNA sequencing data, 15 genes that showed significant differential expression in the resistant genotype, i.e., KWR 108, were selected for real-time PCR verification of the RNA sequencing results. Primer Blast Program, available at NCBI database (accessed on 20 October 2023),³ was used for designing primers with CDS of 15 selected genes and CDS of Glucose-6-Phosphate Dehydrogenase (G6PD) gene as internal control (Reddy et al., 2016). RNA was extracted by using RNA isolation kit (Qiagen) from control and treatment plants of KWR 108 and cDNA was synthesized from RNA using the cDNA synthesis kit following the manufacturer’s protocol (Himedia, HiGenoMB). The gene expression analysis was done by using the cDNA template and 2× AB HS SYBR Green qPCR Mix (Thermo Scientific) in a qRT-PCR (Stratagene Mx3000P, Agilent). The details of real-time PCR primers specific to the ERF transcription factor family used in the experiment are listed in Table 1.

Sample comparison was made for GL 13001 (control vs. treatment) and KWR 108 (control vs. treatment), respectively, for differential expression analysis. The analyses were performed in the edge R package to identify differentially expressed genes. Reads per million mapped reads or Counts per million mapped reads for both control and treatment samples, LogCPM (Log10 of CPM value), Log2Fold Change [the logarithm (base 2) of the fold change], pval (value of p for the statistical significance of this change), FDR [FDR adjusted value of p (q-value)], and the criteria used to identify significantly down and up-regulated transcripts were log2FC > 0 (means up-regulated), log2FC > 0 and q value<0.05 (means significantly up-regulated), log2FC < 0 (means down-regulated), log2FC < 0 and q value <0.05 (means significantly down-regulated). In real-time PCR analysis, the mean and standard deviation were calculated based on three pieces of biological replicated data and, using the Livak method (Schmittgen and Livak, 2008), the relative gene expression was calculated with G6PD as the internal reference gene. The rest of the data statistics and graphing were performed using MS Excel. The reliability of the RNAseq data was confirmed by comparing the expression data of selected genes using a linear regression model in Excel.

In the present investigation, the total reads obtained from four samples, i.e., GL 13001 (control and treatment) and KWR 108 (control and treatment), ranged from 33.19 million reads in GL 13001-Control to 42.43 million reads in KWR 108-Treatment (Table 2).

The QC passed reads were then subjected to reference-based transcript assembly. The genome size for Cicer arietinum is 531 Mb. There are a total of 30,236 genes and 220,289 CDS as per the GFF file of the NCBI chickpea reference genome. The transcript assembly was performed by mapping clean reads of control and treatment samples of GL 13001 and KWR 108 on the reference genome using STAR (v 2.7.10a). The mapping statistics of HQ reads against the reference genome are presented in Table 3.

The number of assembled transcripts varied from 36,081 to 39,186 with an average transcript size of 2003.40 bp and 1992.60 bp in control and treated GL 1300, respectively. Similarly, the number of assembled transcripts ranged from 37,496 and 37,304 with an average transcript size of 2005.90 bp and 1984.60 bp in control and treated KWR 108 leaves samples, respectively. Simultaneously, StringTie’s merge resulted in several transcripts assembled as 51,353 with an average transcript size of 1898 bp in merge data (Table 4).

StringTie’s merge function was used because individually in some of the samples, transcripts might only be partially covered by reads, leading to the assembly of only a partial version of those transcripts in the initial StringTie run. Thus, the set of transcripts that are consistent in all four samples was created to be compared in the subsequent steps using the merge step.

A total of 4,111 novel transcript isoforms corresponding to or encoded by 2,918 genes were identified as potential novel isoforms. Interestingly, 167 novel transcripts in KWR 108 and 170 novel transcripts in GL 13001, which were not available in the annotation file of the reference genome, were identified. The detailed list of novel transcripts and novel isoform transcripts identified are enlisted in Supplementary raw data file S1.

Large numbers of transcripts were shown to be differentially expressed as revealed by differential expression analysis using edgeR software. Simultaneously, 7,504 and 6,698 transcripts were found to be significantly differentially expressed in GL 13001 and KWR 108, respectively (Table 5).

By transforming the data onto A (mean average) and M (log ratio) scales, the differences between measurements taken in control and treatment samples were visualized in MA plot. Simultaneously, a volcano plot that arranges differentially expressed transcripts along dimensions of statistical as well as biological significance was obtained from the edgeR software (Figures 2A–D). Similarly, the heatmap of 25 highly significant down and up-regulated genes (total 50 transcripts) based on q-value (i.e., least q value) and transcripts with proper annotations generated for both the genotypes are represented in Figures 2G,H.

The reliability of sequencing data was confirmed by real-time-based expression analysis of 15 selected genes that showed significant up-regulation in the resistant genotype under Fusarium wilt infection as revealed by RNAseq data. Similar to RNA sequencing results, the real-time PCR data of 15 selected genes showed significant differential expression of the genotype KWR 108 (Figure 2E). The linear regression analysis of RNAseq and real-time PCR-based expression data also showed significant similarity (R² = 0.871), which confirmed that the RNAseq data were reliable for further analysis (Figure 2F).

Identification of genes encoding differentially expressed transcripts revealed that, in GL 13001, 427 genes were significantly down-regulated and 794 were significantly up-regulated (Supplementary raw data file S2). Similarly, in KWR 108, 495 genes were significantly down-regulated and 608 genes were significantly up-regulated (Supplementary raw data file S2). When the genes that showed significant differential expression were compared between two genotypes, i.e., GL 13001 and KWR 108, it was revealed that 81 genes were down-regulated and 211 genes were up-regulated in both the genotypes. Ninety-one and 172 genes showed significant down-regulation and up-regulation only in GL 13001 and 93 and 97 genes showed significant down-regulation and up-regulation only in KWR 108, respectively. Similarly, 111 genes were significantly up-regulated in GL 13001 but significantly down-regulated in KWR 108 and 50 genes were significantly down-regulated in GL 13001 but significantly up-regulated in KWR 108.

The GO enrichment analysis of significant differentially expressed genes (DEGs) summarized the genes into three main GO categories, namely cellular component, molecular function, and biological process (Supplementary raw data file S2). The GO enrichment results of the top 50 significantly up-regulated and down-regulated genes in two genotypes are represented as pie charts (Figure 3).

The GO enrichment results of the top 50 significantly up and down-regulated genes revealed that the proteins encoded by most of these genes are the components of the membrane and nucleus. The products of most of the top 50 genes down-regulated in GL 13001 have nucleotide binding, catalytic, metal ion binding, transferase, and hydrolase activity involved in biological processes like metabolic processes, mRNA processing, DNA repair, oxidation–reduction process, and macroautophagy. Apart from membrane and nucleus, many of the proteins encoded by the top 50 genes up-regulated in GL 13001 are the components of proteasome and ribosome having nucleotide binding, catalytic, metal ion binding, oxidoreductase, and transferase activity involved in biological processes like metabolic processes, translation, DNA repair, and oxidation–reduction process.

Similar to the GO terms of the top 50 genes significantly up and down-regulated in GL 13001, the GO terms of the top 50 genes down-regulated in KWR 108 showed the highest cellular component as membrane and nucleus and molecular functions as nucleotide binding, metal ion binding, transferase, catalytic, and nuclear pore activities. However, the GO terms of the top 50 genes up-regulated in KWR 108 revealed that, apart from membrane and nucleus, the products of the genes are the components of ribosome and cytosol having molecular functions such as kinase activity, oxidoreductase activity, nucleotide binding, transferase, and metal ion binding activities involved in biological processes such as phosphorylation, oxidation–reduction process, metabolic and biosynthetic process, and cell redox homeostasis (Figure 3).

KEGG enrichment analysis of the significant DEGs identified in both the genotypes of chickpea was performed to estimate the numbers of significant DEGs involved in different KEGG pathway levels. It revealed that, in the susceptible genotype, i.e., GL 13001, the number of significantly up-regulated genes involved in different pathways such as DNA replication and DNA repair, cell motility, lipid metabolism, carbohydrate metabolism, and DNA degradation is significantly higher than in KWR 108. Whereas the number of DEGs specifically involved in energy metabolism and environmental adaptation pathways are significantly higher in KWR 108 than in GL 13001 (Figure 4A). The same genes that were up-regulated in GL 13001 but down-regulated in KWR 108 were observed to be higher in most of the pathways; more specifically, the genes involved in cell motility and membrane transport were significantly up-regulated only in GL 13001 (Figure 4B). However, the same genes involved in carbohydrate metabolism were up-regulated in resistant genotype KWR 108 but down-regulated in susceptible genotype GL 13001 (Figure 4B).

In our present study, a comparison between the shoot transcriptome of KWR 108 and GL 13001 showed a significantly higher number of genes involved in DNA repair mechanisms were up-regulated in GL 13001, a susceptible genotype (Figure 4A). Although, under normal circumstances, DNA damage can be thought of as a frequent occurrence, it is likely to occur even more in response to several stress conditions (Nisa et al., 2019) either directly or indirectly due to increased level of reactive oxygen species (ROS; Szurman-Zubrzycka et al., 2023). The resistant or tolerant plants protect the DNA from damage by activating enzymatic and non-enzymatic antioxidant machinery (Moradbeygi et al., 2020). In the present study, the significantly up-regulated genes related to DNA repair mechanism in susceptible genotypes are crossover junction endonuclease EME1, DNA excision repair protein ERCC-1 & 4, DNA mismatch repair protein MLH3, DNA repair and recombination protein RAD54 and RAD54-like protein, DNA repair protein REV1, DNA-3-methyladenine glycosylase I, flap endonuclease-1, telomere length regulation protein, UV excision repair protein RAD23, and WD repeat-containing protein 48. This suggests that DNA replication and other similar processes of DNA synthesis are more stable in resistant genotypes as compared to susceptible genotypes under Fusarium wilt stress condition.

The genes specifically belonging to environmental adaptation pathways are significantly up-regulated in the resistant genotype. The environmental adaptation genes basically encode for the proteins or enzymes involved in plant-pathogen interaction upon infection. In the present study, the significantly up-regulated environmental adaptation genes in KWR 108 were calcium-dependent protein kinase (CDPK), pathogenesis-related genes transcriptional activator (Pti5 & Pti6), calmodulin (CaM), calmodulin-like proteins (CML), cyclic nucleotide-gated channel (CNGC), mitochondrial translation elongation factor Tu (TUFM), chitin elicitor receptor kinase 1 (CERK1), enhanced disease susceptibility 1 protein (EDS1), disease resistance protein (RPM1/ RPS3), mitogen-activated protein kinase 3 (MPK3), mitogen-activated protein kinase 4 (MPK4), mitogen-activated protein kinase kinase 4/5 (MKK4_5), pathogenesis-related protein 1 (PR1), and WRKY transcription factor 33 (WRKY33). The important roles and involvement of environmental adaptation genes like CDPK (Asano et al., 2012; Baba et al., 2019; Wu et al., 2021), Pti5 & Pti6 (Yang et al., 2020; Wang Y. et al., 2021; Sun et al., 2022), CaM, and CML (Zeng et al., 2015; Ranty et al., 2016) in several biotic and abiotic stress responses have been extensively studied and reported. The plant accumulates Ca⁺ upon experiencing stress as a secondary messenger to modulate the cellular processes (Luan, 2009). CNGC, which showed significant up-regulation under Fusarium wilt in the present study, is mostly localized to the plasma membrane and acts as ligand-gated cation channels (Chin et al., 2009) for the accumulation of Ca⁺ in the cell during stress. Several reports suggested the potential role of CNGC in plant disease resistance through hypersensitive response (Clough et al., 2000; Ahn, 2007; Keisa et al., 2011; Ma and Berkowitz, 2011; Guo et al., 2018; Peng et al., 2019; Dietrich et al., 2020; Chakraborty et al., 2021; Zhao et al., 2021). Mitochondrial translation elongation factor Tu (TUFM) is a mitochondrial protein encoded by the nuclear genome that potently inhibits RLR (RIG-I-like receptors) signaling and promotes autophagy during a pathogen infection, particularly viruses (Lei, 2011; Wu et al., 2022). In the study, TUFM showed significant up-regulation in Fusarium wilt-resistant genotype KWR 108, suggesting its potential role in imparting resistance against fungal pathogens and other stresses as reported in crops like maize (Rao et al., 2004) and soybean (Yin et al., 2009). Chitin, which is an essential unit of fungal cell walls including chickpea wilt, causes F. oxysporum to act as an elicitor of plant immunity (Gong et al., 2020). Chitin elicitor receptor kinase 1 (CERK1), which showed significant up-regulation in resistant genotypes in the present study, is a heterotrimeric complex (Gigli-Bisceglia and Testerink, 2021) required for the perception of chitin and subsequent activation of the chitin-responsive cellular processes as well as for triggering the innate immunity of the plants against abiotic stress (Espinoza et al., 2017). CERK1 mutant (cerk1) plants showed increased susceptibility to fungal pathogens like Erisiphe cichoracearum, Alternaria brassicicola, and Blumeria graminis (Gigli-Bisceglia and Testerink, 2021), indicating its vital role in perceiving pathogen signals and subsequent induction of defense mechanism in resistant plants. The enhanced disease susceptibility 1 (EDS1) gene, which showed significant upregulation in KWR 108 in the present study, encodes a lipase-like protein that has been reported to be essential for resistance against biotrophic pathogens (Ochsenbein et al., 2006). In plants, it is basically required by Toll/Interleukin-1 receptor (TIR) domain of NLR (nucleotide-binding leucine-rich repeat) immune receptors (Johanndrees et al., 2023) and for inducing expression of PR (pathogen-related) proteins through the salicylic acid (SA) pathway against several pathogens (Yan et al., 2016; Bhandari et al., 2019; Neubauer et al., 2020). Likewise, RPM1/RPS3 (disease resistance protein), which in the present study showed significant up-regulation in the resistant genotype, is basically a CC-NBS-LRR protein that was first identified as a protein in Arabidopsis essential for resistance mechanism against P syringae pv. Maculicola. It was later reported to have roles in resistance against diseases caused by other pathogens like Puccinia striiformis f. sp. tritici through hypersensitive response and Ca⁺ concentration dynamics (Grant et al., 2000; Wang et al., 2020). In the present study, mitogen-activated protein kinases (MAPK) such as MPK3, MPK4, and MKK4_5 showed significant up-regulation in resistant genotypes that are basically involved in intracellular responses (Sinha et al., 2011). In plants, MAPKs are reported to induce stress tolerance (Lin et al., 2021) including biotic stresses (Zhang and Zhang, 2022) through JA (Jasmonic acid) and SA signaling pathways (Jagodzik et al., 2018). The SA signaling in plants under stress including fungal infection induces the expression of defense proteins like PR1 (pathogenesis-related protein 1; Almeida-Silva and Venancio, 2022; Anuradha et al., 2022). Antifungal and several other antimicrobial activities of PR1 have been well-reviewed and reported (Kattupalli et al., 2021; Almeida-Silva and Venancio, 2022; Anuradha et al., 2022; Pečenková et al., 2022). Similar to all these reports, in the present study also, PR1 genes showed significant up-regulation in chickpea under Fusarium wilt stress, confirming its role as a defense mechanism against Fusarium wilt disease. Likewise, in plant defense responses, an important oxidoreductase enzyme is peroxidase, which is involved in cell wall modifications and the generation of reactive oxygen species (ROS; Chang et al., 2021). Xue et al. (2017) reported that, in common bean, the peroxidase gene PvPOX1 enhanced resistance against Fusarium wilt. In the present study, a large number of peroxidase genes up-regulated in KWR 108, but none in GL 13001, suggesting its role in resistance against Fusarium wilt in chickpea also. Similarly, WRKY33, an important member of WRKY transcription factor family in plants, showed significant up-regulation in the resistant genotype. The role of WRKY transcription factor has been studied extensively and it plays a role in responses to several stimuli irrespective of physiological, biotic, or abiotic stress (Phukan et al., 2016; Cheng et al., 2021).

The important metabolites involved in redox homeostasis of plant cells are nicotinamide adenine dinucleotide (NAD) and its other derivatives, which are mostly nucleotide coenzymes (Kapoor et al., 2015). The redox state of NAD, i.e., NAD (P)+/NAD(P)H (ratio of oxidized to reduced form) in the cell plays an important role as a signal to ameliorate the bridge between cellular metabolic state and gene expression under various conditions, including biotic and abiotic stress condition (Kapoor et al., 2015), and are vital for homeostasis of cellular energy in plant cells under stress condition for optimum growth and development. In the present study, genes related to NAD (P)H, such as NADH-dependent glutamate synthase, NAD(P)H-quinone oxidoreductase subunit 4, NADH dehydrogenase (ubiquinone) Fe-S protein, NADH:quinone reductase, and NADH–ubiquinone oxidoreductase, showed significant up-regulation in resistant genotypes, confirming its role in energy homeostasis in plants under stress conditions. Other than NADH-related genes, shoot transcriptome data of resistant genotypes in the present study also showed significant up-regulation of other genes involved in energy metabolism, such as phosphoenolpyruvate carboxylase (Wang et al., 2016; Waseem and Ahmad, 2019), alcohol dehydrogenase (Su et al., 2020), alanine transaminase (Bashar et al., 2020), ATP citrate (pro-S)-lyase (Liu F. et al., 2022), ATPase (Li J. et al., 2022; Li Y. et al., 2022), carbonic anhydrase (Rudenko et al., 2021), cytochrome c, cytochrome c oxidase (Guerra-Castellano et al., 2018; Analin et al., 2020), fructose-bisphosphate aldolase (Lv et al., 2017; Cai et al., 2022), malate dehydrogenase (Song et al., 2022), 6-phosphofructokinase 1 (Wang H. et al., 2021), photosystem II P680 reaction center D1 protein (Landi and Guidi, 2022), phosphoserine aminotransferase (Wang et al., 2022), pyruvate-orthophosphate dikinase (Yadav et al., 2020), and serine O-acetyltransferase (Mulet et al., 2004; Liu D. et al., 2022), and their molecular mechanism and role in energy homeostasis in different crops is well reviewed. It indicates that the significant up-regulation of the genes involved in energy metabolism in KWR 108 might have helped it to cope with Fusarium wilt stress by maintaining energy homeostasis in cells under stress conditions.

The genes particularly involved in environmental adaptation, energy metabolism, and others which showed significant up-regulation in wilt-resistant genotypes will be a target for designing gene-based markers for marker-trait association studies. The marker-trait association study, especially for the traits attributing to wilt resistance and yield under wilt infection, will be our next target to check whether these up-regulated genes have any roles in imparting resistance against Fusarium wilt through some desirable phenotypic traits or whether there are no roles in phenotypic traits of chickpea under wilt. If these gene-based marker/s are found to be robustly associated with Fusarium wilt-resistant traits, then they could be utilized in marker-assisted breeding to develop a wilt-resistant high-yielding variety of chickpea.