Transcriptome profiling reveals the expression and regulation of genes associated with Fusarium wilt resistance in chickpea (Cicer arietinum L.)

Fusarium wilt (FW) is one of the most significant biotic stresses limiting chickpea production worldwide. To dissect the molecular mechanism of FW resistance in chickpea, comparative transcriptome analyses of contrasting resistance sources of chickpea genotypes under control and Fusarium oxysporum f. sp. ciceris (Foc) inoculated conditions were performed. The high‐throughput transcriptome sequencing generated about 1137 million sequencing reads from 24 samples representing two resistant genotypes, two susceptible genotypes, and two near‐isogenic lines under control and stress conditions at two‐time points (7th‐ and 12th‐day post‐inoculation). The analysis identified 5182 differentially expressed genes (DEGs) between different combinations of chickpea genotypes. Functional annotation of these genes indicated their involvement in various biological processes such as defense response, cell wall biogenesis, secondary metabolism, and disease resistance. A significant number (382) of transcription factor encoding genes exhibited differential expression patterns under stress. Further, a considerable number of the identified DEGs (287) co‐localized with previously reported quantitative trait locus for FW resistance. Several resistance/susceptibility‐related genes, such as SERINE/THREONINE PROTEIN KINASE, DIRIGENT, and MLO exhibiting contrasting expression patterns in resistant and susceptible genotypes upon Foc inoculation, were identified. The results presented in the study provide valuable insights into the transcriptional dynamics associated with FW stress response in chickpea and provide candidate genes for the development of disease‐resistant chickpea cultivars.

presented in the study provide valuable insights into the transcriptional dynamics associated with FW stress response in chickpea and provide candidate genes for the development of disease-resistant chickpea cultivars. INTRODUCTION Chickpea (Cicer arietinum L.) is one of the most important grain legumes, with an annual production of 15.87 million tons (FAOSTAT, 2021). It is a self-pollinated, diploid annual crop with a genome size of ∼740 Mb and grown predominantly in arid and semi-arid regions of the world . The importance of chickpea lies in its intrinsic potential for symbiotic nitrogen fixation and its dietary proteins, vitamins, and essential minerals. Chickpea production is essential for food security and improving the nutritional quality of diets for people living mainly in developing countries. Global chickpea production has seen a significant rise in recent years owing to the research efforts that have resulted in the development and release of improved varieties (FAOSTAT, 2021). However, meeting the growing demand calls for higher productivity gains in the chickpea crop. Further improving crop productivity will require sustainable management of devastating diseases like Fusarium wilt (FW) and Ascochyta blight (AB), which put chickpea cultivation at great risk. FW, caused by a soil-borne fungus, Fusarium oxysporum f. sp. ciceris (Foc), is one of the most prevalent diseases of chickpea worldwide. FW is reported to cause yield losses varying from 10 to 100% depending on varietal susceptibility and suitable climatic conditions . Since FW is a soil-borne disease, its management through crop rotation strategies or chemical control is difficult. Therefore, using varieties resistant to FW is the most cost-effective, promising, and environmentally sustainable strategy to manage the disease. In this direction, several quantitative trait loci (QTLs) for FW resistance have been reported to develop FWresistant cultivars via molecular breeding Sabbavarapu et al., 2013;Varshney et al., 2014). However, genetic variability in the pathogen is high, causing varied levels of virulence and leading to the breakdown of resistance in available sources . To expedite the molecular breeding process or to develop resistant varieties via gene editing approaches, it is imperative to have deeper insights into the molecular mechanism of FW resistance in chickpea.
Upon stress, plants undergo a series of physiological, morphological, molecular, and biochemical changes mediated by a network of genetic regulations. These genetic regulations mostly include modifications in gene expression at both transcriptional and post-transcriptional levels. The identification of genes and their expression quantification have been concrete activities in molecular biology ever since the establishment of RNA's role as the fundamental mediator between the genome and proteome. The study of gene expression profiles under stress can help to pinpoint important components of resistance mechanisms. The advent of RNA-seq has facilitated transcript identification and gene expression quantification in a robust manner. Several transcriptomic studies have successfully delineated the molecular mechanisms of drought, salinity, and AB stress in chickpea (Garg et al., 2019;Mashaki et al., 2018). Different studies using the long serial analysis of gene expression, or microarray approach were performed for only two genotypes (Ashraf et al., 2018;Upasani et al., 2017); however, high-throughput and comprehensive analysis of genes involved in FW resistance is not yet available.
To get a deeper understanding of FW resistance in chickpea, we employed RNA-seq to analyze the transcriptomes of roots of six contrasting chickpea genotypes for FW resistance under control and stress conditions at two different time points. We investigated these data to unravel the transcriptome dynamics associated with FW stress in chickpea and identified key differences that determine disease resistance in chickpea. The genes, transcription factors, and biological pathways showing specific and differential expression patterns between resistant and susceptible genotypes were identified. The identified genes, pathways, and molecular interactions regulating disease resistance against F. oxysporum in chickpea provide potential candidate targets for the development/improvement of FW-resistant chickpea cultivars.

Plant material
A total of four chickpea genotypes and two near-isogenic lines with contrasting phenotypes for FW resistance were selected for the study. The genotype, WR 315, is a desi landrace resistant to FW (Foc race 1A, race 2, race 4, and race 5) (Gaur & Chaturvedi, 2004), while the genotype ICCV 05528 is a breeding line with a moderate level of disease resistance . The C 214 and JG 62 are elite desi cultivars susceptible to FW (early wilter) . The two near-isogenic lines (NILs)-NIL 01 (NIL 10057) and NIL 02 (NIL 10058) are obtained from introgression of foc1 locus conferring resistance to race 1 of FW from WR 315 into the genetic background of C 214 (Sabbavarapu et al., 2013;Varshney et al., 2014).

Stress treatment
For the imposition of stress, seedling raising and inoculum preparation were carried out as described earlier (Sharma et al., 2010). In brief, all the six lines selected were raised in trays filled with sterilized river sand in a greenhouse maintained at 25 ± 1˚C for 10 days. Foc race 1 inoculum was prepared from a single conidial culture of F. oxysporum f. sp. ciceris isolated from wilt-infected plants collected from an International Crops Research Institute for the Semi-Arid Tropics (ICRISAT) wilt sick plot. To prepare mass inoculum, a 7-mm disc of actively growing F. oxysporum f. sp. ciceris culture was put into a 250-mL conical flask containing 100 mL of sterilized potato dextrose broth and incubated for 7 days in an incubator shaker (25 ± 1˚C and 125 rpm). The culture was homogenized in sterilized distilled water and adjusted to 6 × 10 5 conidia mL -1 using a haemocytometer for use as an inoculum (as per Sharma et al., 2010). Ten-day-old seedlings of each test line grown in sterilized river sand were uprooted, cleaned with tap water, and inoculated by dipping in inoculum suspension for 1-2 min to enable conidia to adhere to the roots. Inoculated seedlings were transplanted in preirrigated sterile vertisol and sand (3:1) in pots and incubated in a greenhouse at 25 ± 3˚C.

Sample collection and RNA isolation
The root tissues from both uninoculated/control and inoculated/stressed plants were harvested at 7th day postinoculation (dpi) and 12th dpi in sets of three biological replicates. Total RNA was isolated using the "NucleoSpin RNA Plant" kit (Macherey-Nagel, Düren). Further, the quality of isolated RNA was assessed using Agilent 2100 Bioanalyzer (Agilent Technologies) and Nanodrop 8000 Spectrophotometer (Thermo Fisher Scientific).

Transcriptome sequencing and reference-guided assembly
High-quality total RNA (RIN ≥ 8) was subjected to library construction using TruSeq RNA Sample Prep Kit v2 (Illumina Inc.) as described previously (Garg et al., 2019). The libraries were sequenced on Illumina HiSeq 2500 platform using a 150 bp-paired end strategy. In total, 24 samples representing six different genotypes (C 214, JG 62, ICCV 05528, WR 315, NIL 01, and NIL 02) under control and F. oxysporum inoculated conditions at two-time points (7th and 12th dpi) were sequenced. The raw reads generated from transcriptome sequencing were processed to remove primer/adaptor contam-

Core Ideas
• Transcriptome sequencing of Fusarium wilt (FW) responsive chickpea genotypes elucidates a set of 5182 differentially expressed genes (DEGs). • Identified DEGs were involved in defense response, cell wall biogenesis, secondary metabolism, and disease resistance. • A total of 287 DEGs co-localized with previously reported quantitative trait loci for FW resistance. • Near isogenic lines reduce genetic background in transcriptome studies. • Transcriptome dataset with candidate genes is an excellent functional genomics resource for FW resistance breeding.

Expression profiling and identification of DEGs
The transcript abundance was estimated based on fragments per kilobase of transcript per million mapped reads (FPKM) using Cufflinks. Transcripts with FPKM ≥ 1 in any of the samples and quantification status as "OK" were considered to be expressed and were used for downstream analysis. The differential gene expression between different samples was estimated using Cuffdiff (Trapnell et al., 2013). Differentially expressed genes (DEGs) were identified between different samples under control (noninoculated), and F. oxysporum inoculated conditions at different time points. The different pairwise combinations included stress versus control samples of each genotype at 7th and 12th dpi; resistant versus susceptible genotypes under control conditions at 7th and 12th dpi; and resistant versus susceptible genotypes under FW stress conditions at both 7th and 12th dpi. The DEGs with log 2 fold change ≥ 2 (upregulated) and/or ≤ −2 (downregulated) and an FPKM ≥ 2 for either of the samples in each pair-wise comparison were considered to be significantly differentially expressed. The log 2 transformed FPKM values of the DEGs were further subjected to k-means clustering using pheatmap package in R/Bioconductor.

Annotation and GO enrichment of DEGs
To predict putative function, the identified DEGs were subjected to similarity search using BLASTX (E-value cut-off of ≤ 1E-5) against the National Center for Biotechnology Information nonredundant Viridiplantae protein database and were further assigned gene ontology (GO) terms using Blast2GO v4.1.9 (Conesa et al., 2005). To pinpoint the overrepresented functional categories in DEGs across different sample comparisons, an R package GOseq (Young et al., 2010) was used. GO terms displaying a corrected p-value of ≤ 0.05 for a given set of genes were considered to be significantly enriched. The enriched GO terms were summarized and visualized using REVIGO (Supek et al., 2011).

Pathway and TF analysis
The identified DEGs were studied for their involvement in different pathways using the Kyoto Encyclopedia of Genes and Genomes (KEGG) Automatic Annotation Server (Moriya et al., 2007). The pathway enrichment analysis was performed based on a hypergeometric model (Boyle et al., 2004) with a significance threshold of p-value of 0.05. The enriched pathways were plotted as a scatter plot using data visualization R package ggplot2 (Wickham, 2016). Further, the transcription factor encoding genes were predicted using a similarity search against the PlantTFDB 5.0 database (Tian et al., 2020) with an E-value cut-off of 1E-10. For co-localization studies, information about the FW resistance QTLs was obtained from previous studies Patil et al., 2014;Sabbavarapu et al., 2013). The sequence similarity search using BLASTN was performed to obtain the physical locations of the markers associated with FW resistance QTLs in the chickpea reference genome. Only hits with 100% query and subject coverage were retained.

Transcriptome sequencing and reference-guided assembly
In the present study, a total of 24 samples representing four chickpea genotypes (WR 315-resistant, ICCV 05528-moderately resistant, C 214-susceptible, and JG 62susceptible) and two near-isogenic lines (NIL 10057 and NIL 10058, designated as NIL 01 and NIL 02, respectively) under control and stress (F. oxysporum inoculated) conditions at two (7th dpi and 12th dpi) time points were subjected to transcriptome sequencing. The two-time points identified for this study were selected based on the earlier reports in legumes that the initial symptoms begin at 7th dpi and express to the maximum at 12th dpi (Biswas et al., 2020). The NIL 01 and NIL 02 used in the present study were obtained by introgressing FW-resistant QTLs from WR 315 into C 214 (Varshney et al., 2014). In the study, the samples are named genotype-condition-dpi for ease of understanding. For instance, C 214-C-7d represents genotype C 214 under control conditions at 7th dpi. A total of 1137 million reads were generated from the paired-end sequencing of these 24 samples. After trimming and filtering low-quality reads, a set of 1110 million high-quality reads were retained and mapped on the recently available chickpea genome (Garg et al., 2022). Overall, about 95% of the high-quality reads were mapped on the reference genome (Table 1). The sequencing data and alignment statistics reflect very high-quality transcriptome sequencing. On average, 93.27% of the reads were mapped to the exonic region, 0.83% to the intronic region, and 5.90% to intergenic regions of the chickpea genome (Table S1). The mapped reads from all 24 samples were used to generate reference-guided assembly using the Cufflinks-Cuffmerge pipeline. This assembly generated 76,382 transcripts representing 27,736 gene loci, including 25,695 known and 2041 novel gene loci.

Global gene expression analysis
The normalized expression level (FPKM) of each gene was estimated in all the samples. The genes with very low expression values (see Section 2) in all 24 samples were filtered out. After filtering, a set of 21,609 genes exhibiting expression in at least one of the samples was obtained. The number of expressed genes varied from a minimum of 18,041 genes in WR 315-S-12d to a maximum of 18,872 genes in ICCV 05528-C-12d (Table 2). Further, based on the expression level, the genes were classified into highly expressed (FPKM ≥10) and low/moderately expressed genes (FPKM ≥1 and <10). The number of low/moderately expressed genes was less than highly expressed genes across all the samples in the study. The number of low/moderately expressed genes ranged from 7563 in JG 62-S-12d to 8390 in ICCV 05528-C-12d, whereas highly expressed genes ranged from 10,257 in ICCV 05528-S-12d to 10,671 in NIL 01-C-12d (Table 2).

Identification of differentially expressed genes in response to FW stress
The fold change of each gene was calculated across different pairwise combinations using Cuffdiff. A gene was considered differentially expressed if it exhibited a log 2 fold change of ≥ 2 T A B L E 1 Statistics of transcriptome sequencing data generated for 24 samples obtained from six different chickpea genotypes.

Raw reads Filtered reads Filtered reads (%) Mapped reads Mapped reads (%)
C (upregulated) or ≤ −2 (downregulated). Using this criterion, a total of 5182 genes were differentially expressed in at least one sample combination (Table S2). The differential expression patterns between control and stress samples of each genotype were studied at both stages of infection (7th and 12th dpi). At 7th dpi, the maximum DEGs were identified in the JG 62 genotype (577 upregulated; 357 downregulated), followed by WR 315 (327 upregulated; 353 downregulated) (Figure 1a-c). At 12th dpi, the maximum transcriptional changes were observed in WR 315 genotype with differential expression of 1078 genes (656 upregulated; 422 downregulated). The minimum number of DEGs was identified between C 214-C-7d and WR 315-C-7d (63), and the maximum was identified between ICCV 05528-S-12d and JG 62-S-12d (1452). Under stress, the pattern of upregulation was more profound at the 12th dpi than the 7th dpi across all genotypes (Figure 1c). Further, an overlap of DEGs between different genotypes under control and stress conditions at both stages of infection was studied. A significant number of DEGs were unique to each combination. In total, 1230 genes showed genotypespecific expression patterns under stress at 7th dpi. At 12th dpi, 1549 DEGs were genotype-specific. At 7th dpi, the maximum number of genotype-specific DEGs were found in JG 62 (611) and the minimum in NIL 02 (49) (Figure 1d). Similarly, at 12th dpi, WR 315 (724) and NIL 02 (22) had the maximum and minimum genotype-specific DEGs, respectively (Figure 1e). A total of 1915 genes exhibited differential expression patterns at both 7th and 12th dpi upon infection in any genotypes. These findings indicate an extensive genotypeand stress stage-specific response of chickpea to Foc infection.

Genome-wide distribution of DEGs and localization with FW QTLs
The identified DEGs were studied for their distribution in the chickpea genome. Out of 5182 DEGs, 5108 were mapped to the pseudomolecules of the chickpea genome and the remaining on scaffolds (Figure 2). The highest number of DEGs were present on pseudomolecule Ca6_v2.0 (803), followed by pseudomolecule Ca4_v2.0 (741), and the lowest was present on pseudomolecule Ca8_v2.0 (431). Several QTLs have been reported for FW resistance to date Patil et al., 2014;Sabbavarapu et al., 2013). The identified DEGs were further co-localized with previously reported FW resistance QTLs. By deploying the physical position of the simple T A B L E 2 Distribution of expressed genes identified from different chickpea genotypes. sequence repeats (SSRs) markers on the chickpea pseudomolecules, these QTLs were mapped on pseudomolecules Ca2_v2.0, Ca4_v2.0, and Ca6_v2.0 of the chickpea reference genome (Garg et al., 2022). Further, a total of 287 out of the 5182 DEGs, were co-localized with seven FW resistance QTLs. Of these 287 co-localized DEGs, 39 were novel. The maximum DEGs were present in two QTLs, namely, FW-Q-APR-2-1 (87, Garg et al., 2018) and Wilt 2 (76, Patil et al., 2014) located on pseudomolecule Ca2_v2.0 ( Figure 2). Further, functional annotations were assigned to the identified DEGs using BLASTX. From the total, 4767 DEGs were annotated, and 4213 of these DEGs were assigned GO terms. A total of 13,823 GO terms were retrieved for these DEGs, and these GO terms were evenly assigned to biological process (4261), molecular function (5813), and cellular component (3749) categories. GO enrichment analysis of the identified DEGs suggested the enrichment of GO terms such as defense response (GO:0006952), response to reactive oxygen species (GO:0000302), ion homeostasis (GO:0050801), signal transduction (GO:0007165), defense response to fungus (GO:0050832), cell wall biogenesis (GO:0042546), response to chitin (GO:0010200), brassinosteroid biosynthetic process (GO:0016132), abscisic acid-activated signaling pathway (GO:0009738), and response to auxin (GO:0009733) in the biological process category (Figure 3a; Table S3). Among the enriched GO terms related to biological processes, response to chitin, defense response, lipid transport, abscisic acid-activated signaling pathway, MAPK cascade (GO:0000165), and flavonoid biosynthesis (GO:0009813) were mainly associated with genes showing increased expression upon stress, whereas genes showing downregulation were related to circadian rhythm (GO:0007623) and response to ethylene (GO:0009723). Among the cellular component category, extracellular region (GO:0005576), plant-type cell wall (GO:0009505), apoplast (GO:0048046), photosystem (GO:0009521), and anchored component of plasma membrane (GO:0046658) were the most overrepresented GO terms (Table S4). In the molecular function category, oxidoreductase activity (GO:0016491), transcription regulator activity (GO:0140110), protein serine/threonine kinase activity (GO:0004674), peroxidase activity (GO:0004601), hydrolase activity (GO:0016798), and monooxygenase activity (GO:0004497) were the most significantly enriched GO terms (Table S5).

F I G U R E 1 Statistics of differentially expressed genes (DEGs) identified between different combinations of chickpea genotypes in response to
Further, pathway analysis of the DEGs was performed to identify pathways enriched in response to FW stress in chickpea. We identified a total of 130 diverse pathways using the KEGG database. The enrichment analysis (p-value < 0.05) indicated that genes involved in the biosynthesis of secondary metabolites (ko01110), MAPK signaling pathway (ko04016), plant-pathogen interaction (ko04626), plant hormone signal transduction (ko04075), flavonoid biosynthesis (ko00941), and fatty acid elongation (ko00062) were significantly differentially expressed (Figure 3b).
Transcription factors (TF), being the upstream regulatory proteins, play significant roles under both abiotic and biotic stresses in plants (Bhoghireddy et al., 2020;Nuruzzaman et al., 2013). The reference-guided transcriptome assembly identified members of 58 TF families represented by 1564 genes using the PlantTFDB database (Tian et al., 2020).
Among these, 382 TF-related genes from 39 families were found to show differential expression patterns under FW stress (Table S6). The maximum members of the ERF family (56) were differentially expressed, followed by members of bHLH (43), WRKY (35), MYB (32), and C2H2 (23) TF families. In addition, NAC (19), GRAS (17), Dof (17), and bZIP (16) TF families also exhibited differential expression patterns under stress conditions (Figure 3c). Taken together, these analyses provided an important resource for identifying specific processes, functions, and pathways associated with FW resistance in chickpea.

Expression trends across resistant and susceptible chickpea genotypes under FW stress
To gain deeper insights into the molecular mechanisms of resistance induced by Foc, the identified DEGs were F I G U R E 2 A circular plot representing genome-wide distribution of Fusarium wilt (FW) resistance quantitative trait loci (QTLs), differentially expressed genes (DEGs), transcription factors (TFs), and expression of genes at 7th and 12th day post-inoculation (dpi). Six different tracks (out to in) of the circular plot shows the following: (a) FW resistance QTLs; (b) DEGs; (c) differentially expressed TFs; (d) expression of genes (log 2 fragments per kilobase of transcript per million mapped reads, FPKM) in different genotypes at 7th dpi (moving from out to in; C 214, JG 62, NIL 01, NIL 02, ICCV 05528, and WR 315); (e) expression of genes (log 2 FPKM) in different genotypes at 12th dpi (moving from out to in; C 214, JG 62, NIL 01, NIL 02, ICCV 05528, and WR 315).
Further, we identified key genes that might be responsible for different phenotypes of chickpea genotypes to FW stress. We investigated those genes that showed similar expression among resistant (WR 315 and ICCV 05528) and near-isogenic lines (NIL 01 and NIL 02) but contrasting expression patterns in susceptible genotypes (C 214 and JG 62). Upon FW infection, genes (Cluster I) involved in oxidation-reduction (Ca_v2.0_07841), fatty acid biosynthesis (Ca_v2.0_24909), and defense response (Ca_v2.0_20495 and Ca_v2.0_20402) exhibited induced expression at both 7th and 12th dpi in resistant genotypes F I G U R E 5 Response specific to resistant and susceptible chickpea genotypes under stress conditions. The genes with similar expression patterns have been grouped together into four clusters and selected representative genes from each cluster are shown in the heatmap. Clusters I-IV represent differentially expressed genes (DEGs) identified between resistant (including WR 315, ICCV 05528, NIL 01, and NIL 02) and susceptible genotypes (C 214 and JG 62) upon Foc infection at 7th and 12th dpi. The color scale at the top shows log 2 fold change. In the heat map, blue signifies upregulation/induced expression and red signifies downregulation/repressed expression in resistant genotypes. The asterisk shows that the gene is present in one of the previously reported Fusarium wilt (FW) quantitative trait loci (QTLs).

DISCUSSION
Understanding the molecular basis of variability in stress response is crucial for all research designed to develop new chickpea cultivars with FW resistance. In this direction, the present study reported a comprehensive transcriptome analysis of 24 samples representing six well-characterized chickpea genotypes (C 214, JG 62, ICCV 05528, WR 315, NIL 01, and NIL 02) under control and stress conditions at two stress stages. The data generated from these samples were used to construct a reference-guided assembly that assembled 27,736 genes and comparison of these genes with the gene models of the chickpea reference genome led to the identification of 2041 novel genes from unannotated loci, thus demonstrating the potential of RNA-seq in the identification of novel genes. Similar results were also reported from other transcriptome studies in chickpea (Garg et al., 2019(Garg et al., , 2023Kudapa et al., 2018). An extensive transcriptional reprogramming was observed in chickpea plants under FW stress at 7th and 12th dpi. A total of 5182 genes showed differential expression patterns between contrasting stress-responsive genotypes under control and stress conditions at different stages of infection. A significant fraction of identified DEGs were present in the QTL regions associated with FW resistance. The functional annotation, GO enrichment, and pathway analysis of these differentially expressed genes throw light on the FW resistance mechanism of chickpea. Several genes involved in the recognition of pathogens, ion homeostasis, cell wall modifications, formation of reactive oxygen species (ROS), signal transduction, and defense response were significantly enriched during FW infection. Pathogen invasion is known to affect the structure and integrity of the plant cell wall. In this study, several genes associated with "cell wall biogenesis," "cell wall organization," and "cell wall modification" were differentially expressed, suggesting the role of structural defense genes in FW response in chickpea similar to the observations reported from other legumes (Chang et al., 2021;Sharma et al., 2016). Several pathways related to specific physiological, metabolic, and biochemical changes were also upregulated in the present study. Transcriptome studies conducted in various plant species have reported the enrichment of similar GO terms and pathways in response to fungal stress (Garg et al., 2019;Nayak et al., 2017). Defense response involves several TFs which regulate the target gene expression by specifically binding to cis-acting elements in the promoter regions (Park et al., 2001;Wang et al., 2009). In the present study, different classes of TFs, such as bHLH, WRKY, ERF, and bZIP exhibited differential expression patterns upon FW stress similar to those observed in previous studies (He et al., 2016;Wang et al., 2009).
When chickpea plants were challenged with Foc, several genes related to defense response, phosphorylation, redox homeostasis, and signal transduction exhibited significantly increased expression levels in all the genotypes. In contrast, genes involved in basic biological processes such as plant growth and photosynthesis were downregulated in all the genotypes. These genes showing similar expression in all the genotypes upon stress might be involved in the basal defense of chickpea against F. oxysporum infection. Under both biotic and abiotic stresses, the plant cells accumulate ROS as a part of the basal defense response. Several antioxidant enzymes such as flavonoids and peroxidases can quench the ROS generated due to infection (Gupta et al., 2019). The present study also demonstrated increased expression of peroxidases and HSPs in all genotypes upon FW infection. Genes encoding for peroxidases were also found to positively regulate defense responses against FW stress in a previous study (Ashraf et al., 2018). Various studies have reported that HSPs play an integral role in developing resistance against biotic stresses by interacting with other defense-related proteins. The increased expression of HSP70 in this study is in concurrence with a recent report on powdery mildew (caused by Golovinomyces orontii) infection in sunflower (Kallamadi et al., 2018).
Several studies have reported an increase in intracellular Ca 2+ levels in plant cells in response to both biotic and abiotic stress (Li et al., 2008;Vadassery & Oelmüller, 2009). The increased Ca 2+ levels induce expression of defense genes via the activation of ion fluxes at the plasma membrane, an oxidative burst, and MAPK activation. In the present study, several genes encoding for calcium-dependent protein kinases (CDPKs), calcium transporting ATPase, and calmodulin-like proteins, the major components of calcium signaling, were differentially expressed in resistant genotypes during the early stages of infection, suggesting the role of calcium signaling in FW response in chickpea. Similar observations were also reported from a previous study on FW in chickpea (Upasani et al., 2017). The increased expression of different stressresponsive genes like serine-threonine kinases, dirigent, and PPR protein was also seen in all resistant genotypes at different time points. Dirigent proteins are known to be implicated in the regulation of lignan biosynthesis and are important for secondary metabolism and pathogen resistance (Li et al., 2017). Their increased expression in resistant genotypes upon stress might be indicative of their involvement in FW resistance. In contrast, several genes encoding for senescenceassociated, accelerated cell death, MLO, DMR6, WRKY25, and RAP2.1 were upregulated in susceptible genotypes upon infection. Genes such as MLO, DMR6, and WRKY25 are classified as susceptibility (S) genes and are known to mediate plant-pathogen compatibility and facilitate pathogen infection (Henningsen et al., 2021). Loss-of-function mutations of MLO genes are reported to reduce susceptibility to fungal pathogens in various crops such as barley, tomato, pea, and wheat (Bai et al., 2008;Büschges et al., 1997;Humphry et al., 2011;Wang et al., 2014). In Arabidopsis, DMR6 is known to positively affect the susceptibility to downy mildew oomycete Hyaloperonospora parasitica (Van Damme et al., 2008). RAP2.1 is a DREB-type, EAR-motif-containing transcriptional repressor known to negatively regulate plant stress responses (Dong & Liu, 2010). The upregulation of these S genes in both the susceptible (C 214 and JG 62) genotypes indicates their involvement in the susceptibility of these genotypes to F. oxysporum infection. In the present study, we identified many genes exhibiting similar expression in two resistant genotypes (WR 315 and ICCV 05528) but different in NILs (NIL 01 and NIL 02). The expression of these genes in NILs was similar to their recurrent parent (C 214), demonstrating the use of NILs in distinguishing the genes that are most likely to be implicated in resistance. Similar observations are reported from other transcriptome studies where NILs helped reduce genetic background noise and identify candidate genes responsible for resistance/tolerance (Kim et al., 2011;Zhang et al., 2018).
In summary, the comparative analysis of the transcriptomes of resistant and susceptible chickpea genotypes during F. oxysporum f. sp. ciceris infection has provided insights into different genes encoding for CDPK, MLO, DMR6, WRKY25, Dirigent, NBS-LRR, etc., and pathways that coherently induce plant defense mechanisms through basal and induced resistance. The in-depth understanding of the molecular basis of FW resistance and the candidate genes identified in the study might contribute to the deployment of either genetic engineering or molecular breeding approaches to develop FW-resistant chickpea varieties.

A C K N O W L E D G M E N T S
The authors are thankful to Food Futures Institute, Murdoch University (Australia), for financial support. HB is thankful to the SERB of DST, Government of India, for providing the SERB Women Excellence Award (SB/WEA-01/2017). This work was supported by resources provided by the Pawsey Supercomputing Research Centre with funding from the Australian Government and the Government of Western Australia.

C O N F L I C T O F I N T E R E S T S T A T E M E N T
The authors declare no conflicts of interest.

D A T A AVA I L A B I L I T Y S T A T E M E N T
The RNA-sequencing data generated in the study are available in NCBI Sequence Read Archive (SRA, http://www.ncbi.nlm. nih.gov/sra) with the BioProject ID PRJNA761666.