Comparative transcriptome analysis reveals molecular regulation of salt tolerance in two contrasting chickpea genotypes

Salinity is a major abiotic stress that causes substantial agricultural losses worldwide. Chickpea (Cicer arietinum L.) is an important legume crop but is salt-sensitive. Previous physiological and genetic studies revealed the contrasting response of two desi chickpea varieties, salt-sensitive Rupali and salt-tolerant Genesis836, to salt stress. To understand the complex molecular regulation of salt tolerance mechanisms in these two chickpea genotypes, we examined the leaf transcriptome repertoire of Rupali and Genesis836 in control and salt-stressed conditions. Using linear models, we identified categories of differentially expressed genes (DEGs) describing the genotypic differences: salt-responsive DEGs in Rupali (1,604) and Genesis836 (1,751) with 907 and 1,054 DEGs unique to Rupali and Genesis836, respectively, salt responsive DEGs (3,376), genotype-dependent DEGs (4,170), and genotype-dependent salt-responsive DEGs (122). Functional DEG annotation revealed that the salt treatment affected genes involved in ion transport, osmotic adjustment, photosynthesis, energy generation, stress and hormone signalling, and regulatory pathways. Our results showed that while Genesis836 and Rupali have similar primary salt response mechanisms (common salt-responsive DEGs), their contrasting salt response is attributed to the differential expression of genes primarily involved in ion transport and photosynthesis. Interestingly, variant calling between the two genotypes identified SNPs/InDels in 768 Genesis836 and 701 Rupali salt-responsive DEGs with 1,741 variants identified in Genesis836 and 1,449 variants identified in Rupali. In addition, the presence of premature stop codons was detected in 35 genes in Rupali. This study provides valuable insights into the molecular regulation underpinning the physiological basis of salt tolerance in two chickpea genotypes and offers potential candidate genes for the improvement of salt tolerance in chickpeas.


Introduction
Chickpea (Cicer arietinum L.) is a nutritious legume crop, but salinity limits its production in areas where it is widely grown (Ryan, 1997;Flowers et al., 2010). Climate change will intensify soil salinity globally (Shokat and Grosskinsky, 2019); therefore, improved salinity tolerance is essential for sustained chickpea productivity in major salt-affected regions worldwide. Chickpea is a salt-sensitive species (Flowers et al., 2010); however, despite its narrow genetic diversity between genotypes, some variation for salt tolerance has been reported (Vadez et al., 2007;Turner et al., 2013;Varshney et al., 2021), which can be exploited for varietal improvement. Developing salt-tolerant crop varieties requires effective genetic variations, selection procedures and insights into salt tolerance mechanisms (Munns and Tester, 2008). The physiological and biochemical changes under salt stress are modulated by numerous genes and mediated via highly complex gene regulatory networks (Vadez et al., 2012). Therefore, identifying candidate genes involved in key physiological processes of salinity tolerance could help direct gene selection in chickpea breeding programs.
Several chickpea studies have investigated transcriptional changes of numerous genes to understand the regulatory mechanisms for tolerance to salt and other abiotic stresses like heat, drought, desiccation, and cold (Mantri et al., 2007;Varshney et al., 2009;Molina et al., 2011;Jain et al., 2013;Garg et al., 2015;Garg et al., 2016;Kaashyap et al., 2018;Kudapa et al., 2018;Kumar et al., 2021;Kaashyap et al., 2022). These studies have explored many chickpea genotypes and tissues under various stress conditions; however, the research was limited to one developmental stage or a single genotype until recently, when RNA-seq was used to examine salt stress tolerance in four genotypes (Kumar et al., 2021). Some studies identified the role of certain genes in salt stress tolerance and adaptation by exploring gene expression profiles using qPCR in tolerant and susceptible chickpea genotypes (Singh et al., 2018;Arefian et al., 2019). In chickpeas, salinity tolerance levels and mechanisms vary among different accessions and cultivars (Vadez et al., 2007;Turner et al., 2013;Sweetman et al., 2020) and the molecular basis of gene regulation for salinity tolerance has not been investigated in all genetic material.
Salt tolerance is a complex trait in plants, typically achieved through osmotic stress tolerance, Na + and Cl − exclusion, Na + and Cl − tissue tolerance, and maintenance of an adequate tissue K + /Na + ratio (Munns and Tester, 2008;Kronzucker and Britto, 2011). Ion exclusion is the ability of salt-stressed plants to keep toxic ions (Na + and/or Cl − ) at relatively low concentrations in shoots by retaining them in roots, unloading them from the xylem stream, or storing them away from leaf photosynthetic tissues (Munns and Tester, 2008;Roy et al., 2014). Tissue tolerance is the ability of cells and tissues to function while accumulating high Na + and Cl − , presumably by compartmentalising ions at the cellular and intracellular levels (Flowers et al., 2015;Munns et al., 2016). In addition, plants tolerate salt stress using a suite of mechanisms such as maintaining Na + /K + homeostasis, growth rate, photosynthesis, cell wall integrity, stress signalling and regulatory pathways, cell redox homeostasis, hormone regulation, carbon partitioning, and translocation (van Zelm et al., 2020). Identifying salt tolerance mechanisms in chickpeas requires dissecting the main components in combination with the associated adaptive responses and then identifying the genic regulation of these mechanisms.
Two desi chickpea cultivars, Genesis836 and Rupali, show contrasting phenotypes when exposed to salt, leading to differences in photosynthesis, growth, and seed yield Khan et al., 2016;Khan et al., 2017;Kotula et al., 2019;Atieno et al., 2021). However, their shoot ion (Na + or Cl − ) concentration did not differ Kotula et al., 2015;Khan et al., 2016), suggesting that ion exclusion (at shoot or leaf level) does not explain the salt tolerance difference between the two genotypes and that it could be due to differences in their tissue tolerance to shoot high Na + Kotula et al., 2015;Khan et al., 2016). Moreover, another study associated salt tolerance in Genesis836 with higher photosynthetic rates and less structural damage to chloroplasts than Rupali, possibly through Na + exclusion from the photosynthetically active mesophyll cells and compartmentalising Na + into the non-photosynthetic epidermal cells (Kotula et al., 2019). Thus, the physiological basis of salt tolerance in chickpeas appears to be driven by a combination of Na + exclusion and tissue tolerance to Na + in different leaf tissues. However, the molecular mechanisms underlying salt tolerance in these genotypes remain unexplored.
This study used RNA-sequencing (RNA-seq) of leaf tissues of the two physiologically well-characterised genotypes grown in control and 60 mM NaCl conditions and a comprehensive transcriptome analysis to garner in-depth and unique information on transcriptional reprogramming, pathways, and regulatory networks associated with salt tolerance. In addition, sequence variants between Genesis836 and Rupali were examined for the presence/absence of functional variants. To our knowledge, this is the first study to explore the chickpea leaf transcriptome in two chickpea genotypes contrasting in response to salt stress.

Materials and methods
Plant growth, stress treatment, and tissue sampling Two genotypes of desi chickpea (salt-tolerant Genesis836 and salt-sensitive Rupali) were selected based on their physiological and seed yield data from previous experiments Kotula et al., 2015;Khan et al., 2016;Atieno et al., 2021). The study was performed in a glasshouse (Waite Campus, The University of Adelaide, Adelaide, SA, Australia) with natural irradiation and photoperiod (22 ± 3°C). Plants were grown in plastic pots (5 L) with a continuously aerated nutrient solution . Fourteen-day-old seedlings were subjected to salt stress (60 mM NaCl, added in four 15 mM increments 24 h apart). The nutrient solution in all pots was renewed weekly across the 40 days of plant growth, both before and after onset of salt treatment. The pH (~7.0) was adjusted twice weekly using potassium hydroxide (KOH); however, the degree of change in pH was similar among the genotypes. Pots were arranged in a completely randomised design, with six biological replicates of each genotype and treatment maintained under controlled conditions. Tissues were sampled for RNA-seq and ion concentrations as soon as the first symptom of leaf damage appeared on older leaves of salt-treated plants (20-day-old plants or 6 days after treatment). The second youngest fully expanded leaf (2nd YFEL) was collected from six different plants/pots for each treatment condition. Leaf tissues were harvested for RNA-seq, quickly frozen in liquid nitrogen, and stored at −80°C until RNA extraction. A similar leaf (2nd YFEL) was harvested at the same time from another set of plants (grown in the same pot) to measure leaf ion concentrations (see below). The plants were grown for another 20 days (40-day-old plants) to observe the growth response, after which shoots and roots were separated and oven-dried at 65°C for 48 h to measure dry mass. The phenotypic data were analysed by two-way analysis of variance (ANOVA) using Genstat (VSN International Ltd. Hemel Hempstead, UK), with the means compared for significant differences using LSD (least significant difference) at the 5% significance level.

Tissue ion analysis
Leaf samples (2nd YFEL oven-dried at 65°C for 48 h) were ground to a fine powder and sub-sampled for analysis of Na + , K + , and Cl − following as previously described . Tissue samples were extracted in 0.5 M HNO 3 by shaking for 48 h at room temperature. Diluted samples of the extracts were analysed for Na + and K + using a flame photometer (Flame Photometer Model 420, Sherwood Scientific Ltd., Cambridge, UK) and Cl − with a chloridometer (Chloride Analyzer Model 926S, Sherwood Scientific Ltd., Cambridge, UK).

RNA preparation and Illumina sequencing
Total RNA was extracted from the 2nd YFEL of 24 samples using the RNeasy Plant Mini Kit (Qiagen), according to the manufacturer's instructions. The quantity and quality of the RNA were assessed using Nanodrop Spectrophotometer (NanoDrop Technologies, Wilmington, USA) and Agilent 2100 BioAnalyzer (Agilent Technologies Inc., Santa Clara, CA, USA). Stranded Illumina TruSeq libraries were prepared using high-quality RNA (RNA integrity number ≥ 8) and run on a HiSeq 2500 (Australian Genome Research Facility Ltd., Melbourne, Vic, Australia) to generate paired-end reads with a length of 100 base pairs (bp). All 24 RNA libraries (2 genotypes × 2 treatments × 6 biological replicates) were spread across four lanes (four technical replicates), totalling 96 RNA-seq datasets produced in FASTQ format. The data discussed in this publication have been deposited in NCBI SRA as PRJNA798198 (https://www.ncbi.nlm.nih.gov/sra/PRJNA798198).

Differential expression analysis
Read counting for the genes was performed and a count matrix was created using the featureCounts() function of the Rsubread package in R (Liao et al., 2019). A metadata file containing the complete information for bam files including genotype and treatment for all 24 samples was created. Four experimental groups were created according to the genotype and treatment conditions: Genesis_control, Genesis_treated, Rupali_control, and Rupali_treated. Next, the DGEList() function from the edgeR package was used to calculate the counts per million (CPM) for each experimental group and calcNormFactors() function was used to calculate normalisation factors (Nikolayeva and Robinson, 2014). Genes with no aligned reads in any sample were filtered out.
Counts of aligned reads were normalised to CPM and fitted with a linear model using Limma (Linear Models for Microarray and RNA-Seq Data) (Ritchie et al., 2015) to identify high-and lowexpressed genes. Linear models allow one to assess differential expression in the context of multi-factor designed experiments, resulting in meaningful comparisons using interaction models. A design matrix was created with a group-mean parametrisation approach for multi-level comparison as explained in the manual (https://www.bioconductor.org/packages/devel/bioc/vignettes/ limma/inst/doc/usersguide.pdf ). A contrast matrix was created to identify differentially expressed genes (DEGs) between control and treated samples in (a) Rupali and (b) Genesis836, and in response to a (c) genotype effect, (d) salt treatment effect, and (e) genotype and treatment interaction effect. False discovery rate (FDR) corrections of p-values were carried out using the Benjamini and Hochberg method (Benjamini and Hochberg, 1995). A gene was considered differentially expressed if it showed a corrected p-value ≤ 0.01, AveExpr ≥ 0, and B ≥ 1. AveExpr is the average log 2 -expression level for that gene across all the samples in the experiment and the Bstatistic (lods or B) is the log-odds that the gene is differentially expressed. The lists of DEGs in all five comparisons were generated along with their respective log 2 FC (fold change) values.
The categories of DEGs are described below, which explain how the differential expression analysis was performed and constitutes the structure of how we present the results:

Salt-responsive DEGs in Rupali
The genes exhibiting significant differences in experimental groups of Rupali (Rupali_control vs. Rupali_treated)

Salt-responsive DEGs in Genesis836
The genes exhibiting significant differences in experimental groups of Genesis836 (Genesis836_control vs. Genesis836_treated)

Salt-responsive DEGs
The genes exhibiting significantly higher expression in one salt treatment compared to the other, independent of plant genotype (both control samples vs. both treated samples)

Genotype-dependent DEGs
The genes exhibiting significantly higher expression in one genotype compared to the other, independent of salt treatment (both Rupali samples vs. both Genesis836 samples)

Genotype-dependent salt-responsive DEGs
The genes with significant genotype and treatment interaction effect (the genes responding differently to salt treatment and among the two genotypes)

Functional annotation
Chickpea gene sequences were retrieved from the reference sequence fasta file using the start and end coordinate information for genes given in the gff3 file. BLASTX searches of chickpea gene sequences were performed against Arabidopsis thaliana and Medicago truncatula Mt4.0v1 protein sequences with an e-value cutoff of 10 −15 . Protein sequences for Arabidopsis and Medicago were downloaded from https://www.arabidopsis.org/download_files/ Proteins/TAIR10_protein_lists/TAIR10_pep_20101214 and https:// phytozome.jgi.doe.gov/pz/portal.html#!info?alias=Org_Mtruncatula, respectively. BLASTX results were inspected for their top hit using an in-house perl script, and subsequently, putative annotations were added to the chickpea genes. Gene Ontology enrichment analysis and KEGG pathways were determined using the DAVID functional annotation tool (https://david.ncifcrf.gov/home.jsp ).

Identifying variants in two genotypes
The bam files generated by the STAR aligner were sorted by coordinates and indexed using SAMtools version 1.2. Next, bcftools mpileup and call was used to identify variants that were annotated with SnpEff (version 5.1). Only variants with quality ≥ 20, depth ≥ 10, and annotated as high and moderate categories were considered.

Phenotypic and physiological responses to salt stress
In the control, Rupali had 25.2% and 16.6% higher shoot and root dry masses, respectively, than Genesis836 (Table 1). Salt stress severely reduced plant growth in both genotypes; however, Genesis836 had a higher shoot and root dry mass (34.3% and 45.8% of controls, respectively) than Rupali (8% and 10.8% of controls, respectively). Thus, the tolerant genotype showed better growth under salt stress after 26 days of treatment.
The salt treatment increased leaf Na + and Cl − concentrations (tissue dry mass basis) compared to controls, which did not differ between the two genotypes (Table 1). Moreover, the salt treatment did not affect leaf K + concentration in either genotype. However, the salt treatment significantly decreased the leaf K + /Na + ratio in both genotypes due to the increased leaf Na + concentration ( Table 1). Genesis836 had a higher leaf K + /Na + ratio for control plants than Rupali, but both genotypes had similar ratios in the salt treatment.

RNA-sequencing and identification of DEGs
RNA-seq of 96 leaf samples (2 genotypes × 2 treatments × 6 biological replicates × 4 technical replicates) generated more than 22 million paired-end 100-bp reads per sample. After trimming the reads to remove adapter contamination and low-quality bases, 95% −97% of high-quality reads were retained for each sample, with >87.3% uniquely aligned to the chickpea reference genome (CDC Frontier kabuli version 2.6.3) containing 33,351 genes (Edwards, 2016). Only 4.9% of reads were multi-mapped at various positions while the remaining 7.8% did not align with the reference genome. Supplementary Table S1 summarises the generated sequence data, trimmed reads, and aligned reads.

Salt-responsive DEGs in Rupali and Genesis836
Salt stress significantly altered gene expression between control and treated samples in both genotypes. Rupali had 1,604 saltresponsive genes ( Supplementary Table S2  Values are means ± SE (n = 6). The least significant differences (LSD) for genotype, treatment, and treatment × genotype interaction are given at the bottom of each column of data (p = 0.05). ***p < 0.001, **p < 0.01, *p < 0.05, n.s. = non-significant. The four experimental groups (Rupali_Control, Rupali_treated, Genesis836_Control, and Genesis836_treated) were compared by LIMMA in five ways (shown as five colours) to test each contrast resulting in a list of DEGs.   Figure 2A). Rupali and Genesis836 had 697 common salt-responsive genes; Rupali had 907 and Genesis836 had 1,054 unique salt-responsive genes ( Figure 2B). Functional annotation and GO (gene ontology) enrichment analysis (Supplementary Table S4) revealed that the biological processes of Rupali's downregulated genes were involved in the oxidation-reduction process, protein phosphorylation, and photosynthesis. In contrast, the highly represented GO terms for most of the upregulated DEGs were translation, response to salt stress, response to cytokinin, ribosome biogenesis, cell redox homeostasis, and leaf morphogenesis. For Genesis836, the enriched GO terms of downregulated genes included protein phosphorylation, protein ubiquitination, response to light, phosphorylation, response to ethylene, response to jasmonic acid, and photosynthesis. In contrast, most of the upregulated genes were associated with response to salt stress, response to heat, cell-redox homeostasis, response to high-intensity light, and transport.
Interestingly, the 697 common salt-responsive genes in the two genotypes ( Figure 2B) had the same direction of expression change in both genotypes, suggesting a similar primary salt response mechanism in both the tolerant and susceptible genotypes, except for three transcription factor (TF) genes (Ca10694, Ca24486 and Ca07280) that exhibited the opposite expression patterns [i.e., genes induced in one genotype were repressed in the other (discussed later in the genotype-dependent salt-responsive genes section)]. Other common salt-responsive genes in the two chickpea genotypes included ABC transporter family proteins, alternative oxidase family proteins, vacuolar H + -pumping ATPase, CIPKs, chloride channel proteins, dehydration-responsive genes, electron transport family proteins, heat shock factors, cyclic nucleotide-gated cation channel proteins, and voltage-dependent anion channel proteins, with rRNA processing, ribosome biogenesis, embryo development ending in seed dormancy, and transport being the enriched GO categories associated with these groups of genes. This suggests that these common salt-responsive genes provide cellular membrane stability, signal transduction, stress response, and transporter roles and are associated with rendering key functions under salt stress (Supplementary Table S4).
In contrast, the salt-responsive genes unique to Rupali and Genesis836 had striking differences in expression pattern as delineated by MapMan pathway views ( Figure 3). Upregulated salt-responsive DEGs unique to Genesis836 were mainly involved in stress (biotic and abiotic), development, hormones, RNA synthesis and processing, regulation, and redox. In contrast, saltresponsive DEGs unique to Rupali were downregulated in the same functional categories. Detailed pathway views showed that heat shock proteins, DEGs encoding TFs (ERF, bZIP, WRKY, MYB, and DOF) and secondary metabolites were mainly upregulated in Genesis836 and downregulated in Rupali. In Genesis836, unique downregulated DEGs were associated with proteolysis and signalling, while unique upregulated DEGs in Rupali were involved in protein synthesis and amino acid activation, vesicle transport and DNA synthesis. GO analysis revealed that unique salt-responsive genes in Rupali were associated with translation, embryo development ending in seed dormancy, response to cold, response to cytokinin, and mRNA processing. In contrast, the unique genes in Genesis836 were involved in the regulation of transcription, DNA-templated, response to salt stress, protein phosphorylation, metabolic process, response to light stimulus, response to cold, protein folding, flower development, circadian rhythm, and transmembrane transport (Supplementary Table S4).

Salt-responsive DEGs
The salinity treatment affected the expression of 3,376 genes (Supplementary Table S6). Mapman analysis showed that the downregulated genes in this category were mainly involved in cell division, stress (heat, light), IAA, ethylene signalling, and calcium regulation. In contrast, the upregulated genes were implicated in cell cycle, stress (drought, salt), jasmonate, and ABA signalling. The GO enrichment analysis revealed that genes with treatment-dependent expression differences were mainly associated with response to light, glycolytic processes, stomatal movement, translation initiation, lignin biosynthesis, and chloroplast organisation (Supplementary Table S4). We did not explore this DEGs category further as investigating genes that represent common or general salt responses was not of interest.

Genotype-dependent DEGs
Comparative transcriptome analysis of the two chickpea genotypes revealed divergent gene expression of 4,170 DEGs, representing the genes expressed in both genotypes ( Figure 4; Supplementary Table S5). Of these, 2,322 genes had higher expression in control and treated Genesis836 (↑Genesis836 expression) than the Rupali control and treated samples. Meanwhile, 1,848 genes had higher expression in Rupali (↑Rupali expression) than the Genesis836 samples. MapMan analysis of these DEGs indicated that a large fraction of genes with ↑Genesis836 expression were involved in DNA synthesis, RNA processing, protein synthesis and activation, and cell redox (thioredoxin), while most of the DEGs with ↑Rupali expression were lightresponsive, receptor kinases, and involved in cell redox (glutaredoxin), cytokinin signalling, and calcium regulation. In this category, 767 DEGs exhibit at least a twofold change in expression with 401 and 366 genes highly expressed in Genesis836 and Rupali samples (control and treated), respectively.
DEGs contributing to the difference in salt response mechanism in two genotypes GO and pathway analysis has highlighted some important functional categories for the unique Genesis836 and Rupali saltresponsive DEGs (Table 2) and DEGs exhibiting genotypedependent expression (Table 3), which are discussed below. The categories below indicate the main differences between the genotypes in their response mechanism to salt stress.
Chloride influx, transport, and regulatory mechanisms are other important aspects of plant salt tolerance studies. The salt treatment upregulated a gene encoding chloride conductance regulatory protein (ICln-Ca04418) in Genesis836, which also presented ↑Genesis836 expression. Similarly, chloride channel C (ClC1-Ca05659) and chloride channel B (ClC-Ca24944) also delineated ↑Genesis836 expression, while only one chloride channel C (ClC1-Ca01155) exhibited ↑Rupali expression.
Cyclic nucleotide-gated ion channel (CNGC) is an interesting category of proteins potentially involved in sodium transport. Interestingly, two CNGCs (Ca17327 and Ca17275) were downregulated in Rupali under salt stress but showed ↑Rupali expression along with Ca32080. One CNGC (Ca26987) was upregulated in control Rupali, while the other CNGC (Ca18209) exhibited ↑Genesis836 expression. Other ion transport DEGs with ↑Genesis836 expression included a sodium/calcium exchanger (CHX-Ca29649). Ca21958 and Ca09655 were among the unique salt-responsive DEGs in Rupali.

FIGURE 5
Heatmap of differentially expressed genes in two contrasting chickpea varieties, salt-sensitive Rupali and salt-tolerant Genesis836 involved in (A) ion transport and (B) photosynthesis.
Thus, the unique salt-responsive DEGs in Rupali and Genesis836 and DEGs exhibiting genotype-dependent expression changes indicate that the two chickpea genotypes exhibit different tolerance mechanisms to combat salt stress.

Genotype-dependent salt-responsive DEGs (interaction)
Genes with genotype-dependent expression change due to salt treatment are of particular interest since they directly represent the difference in salt tolerance mechanisms between genotypes. Genesis836 and Rupali had 122 genes that responded to salt stress (interaction effect in the linear model). Of these, the expression of 88 genes changed direction between control and treated samples in both genotypes. The remaining 34 maintained the same direction but had manifold different expressions between genotypes. Of the 88 genes, 66 were upregulated in Genesis836 but downregulated in Rupali and 22 were downregulated in Genesis836 and upregulated in Rupali due to salinity. Supplementary Table S7 lists all the genes in this category.
Thus, the study highlights many potential candidate genes and their roles in salt response and tolerance mechanisms in the two chickpea genotypes, which can be further potential candidates for breeding salt-tolerant chickpea cultivars.

Validation of differentially expressed genes
Quantitative real-time PCR (qRT-PCR) was performed on RNA samples for eight genes (Ca30477, Ca13456, Ca02100, Ca14863, Ca10383, Ca29966, Ca19227, and Ca01215) randomly selected to validate the RNA-seq results. Table 4 lists the gene-specific primer sequences. A significant correlation occurred between the relative expression levels obtained from RNA-seq and qRT-PCR analysis (r 2 = 0.62 with P < 0.001; Figure 7), validating the biological significance of RNA-seq data.

Variants in the two genotypes
Uniquely mapped reads were used to identify variants in Rupali and Genesis836. We identified 13,829 variants in 6,484 genes in Genesis839 and 14,002 variants in 6,564 genes in Rupali compared to the reference genome. Of these, 1,449 variants were present in 701 Rupali salt-responsive DEGs and 1,741 variants were present in 768 Genesis836 salt-responsive DEGs (Figure 4; Supplementary Table S8).

FIGURE 6
Heatmap of genotype-dependent salt-responsive DEGs and MapMan pathway views delineating the possible role of these genes.
Furthermore, stop_gained mutations were particularly interesting as they indicate premature stop codons and knockouts. Forty variants in 35 genes of the salt-sensitive variety Rupali were stop_gained mutations (Supplementary Table S9). In other words, 35 genes were knocked out due to premature stop codons, suggesting that salt intolerance in the sensitive genotype Rupali is a consequence of loss of function. Of the 35 genes, Ca20374 and Ca17320 were related to salt stress tolerance: Ca20374 is a predicted B-box zinc finger protein, which encodes a transcription factor BBX20 (also named STH7/ salt tolerance homolog 7), involved in the brassinosteroidmediated signalling pathway and photomorphogenesis and Ca17320 is a phosphoinositide kinase involved in protein autophosphorylation and response to salt stress.
Other DEGs with stop codons included a peptidyl-prolyl cis-trans isomerase FKBP62-like protein (Ca16597) involved in response to heat and osmotic stress, a phosphoenolpyruvate carboxylase kinase 1 (PPCK1 Ca10979) whose expression is induced by light and involved in signal transduction and phosphorylation, an RNA polymerase II transcription subunit (Ca06996) involved in regulating hormonemediated signalling pathways and photoperiodism, and a gene involved in morphogenesis, growth, and meristem development (Ca03638) ( Supplementary Table S9). Moreover, Ca16824, Ca02694, and Ca16823 are involved in vesicle transport. TF genes with  Validation of differential expression of genes obtained by RNA-seq. Correlation of gene expression results obtained from RNA-seq and quantitative real-time PCR for eight randomly selected genes. Each value is an individual replication from different treatment combinations. *** Significant with p-value < 0.001. knockout mutations included Ca02188 (RING-H2 finger protein), Ca07014 (AP2/B3 transcription factor family protein), Ca10947 (zinc finger CCCH domain-containing protein), and Ca15012 (auxin response factor protein). Ten genes containing stop_gained mutations were uncharacterised and may serve as potential candidates for further evaluation (Supplementary Table S9).

DEGs associated with salinity tolerance QTLs
A recently published study developed a Recombinant Inbred Line (RIL) population from a cross between Rupali (salt-sensitive) and Genesis836 (salt-tolerant) chickpea (Atieno et al., 2021) and reported three specific salinity tolerance QTLs, i.e., Ca4 (6.8-7.5 Mb),  Mb) from the coordinates of flanking markers based on CDC Frontier v.2.6.3. Thus, we investigated our DEGs categories for genes positioned within the three QTL regions (Supplementary Tables S2, S3, S5, S7). Among the DEGs positioned within these QTLs, 15 genes were identified with a potential role in chickpea salt tolerance (Table 5), which included Ca20374 at Ca5 with a stop_gained mutation in Rupali, a thylakoid lumenal 17.9-kDa protein (Ca21063 at Ca5) and a glutathione S-transferase-related protein (Ca20449 Ca5). Moreover, some dehydration-responsive proteins, heat shock proteins, and MYB-like transcription factor family proteins were among the DEGs within the salinity tolerance QTLs. The gene names highlighted in bold letters are common between different comparisons.

Discussion
Several studies have reported that salt stress induces complex regulatory mechanisms and major transcriptional reorganisation in chickpeas. However, the molecular mechanisms for salt tolerance and the different adaptive responses in salt-tolerant Genesis836 and salt-sensitive Rupali have not been explored. Some studies on Genesis836 and Rupali have demonstrated that Na + exclusion mainly determines salt sensitivity (Khan et al., 2016) and the tissue tolerance to Na + differed between the two genotypes (Khan et al., 2017), with Genesis836 moving ions from photosynthetically active mesophyll tissues into epidermal cells (Kotula et al., 2019). Hence, identifying DEGs in Genesis836 and Rupali can shed light on their tolerance mechanisms and determine what contributes to the enhanced salinity tolerance in Genesis836.
This study compared the transcriptome profiles of Genesis836 and Rupali under salt stress. Both genotypes had a similar number of salt-responsive DEGs, indicating that salinity highly disturbed these genotypes at the transcriptional level. Detailed investigations of Genesis836 and Rupali gene expression profiles grouped the DEGs as salt-responsive, genotype-dependent salt-responsive, and genotype-dependent and treatment-dependent, providing novel insights into the molecular processes involved in salt tolerance. Comparing the control vs. treated samples revealed unique saltresponsive genes in both genotypes. Such unique genes have also been identified in other susceptible and tolerant chickpea genotypes under salt stress (Kumar et al., 2021). The unique salt-responsive genes in each genotype, genotype-dependent expression differences, and genotype-dependent salt-responsive genes collectively indicate that tolerant Genesis836 and sensitive Rupali have distinct salinity responses and tolerance mechanisms. Our study agrees with previous physiological results ( Kotula et al., 2019) that tissue tolerance via leaf Na + regulation and photosynthetic maintenance are the key mechanisms that differ in Genesis836 and Rupali.
Genesis836 maintained higher shoot and root growth than Rupali with similar leaf Na + concentrations, consistent with our previous findings (Khan et al., 2016;Khan et al., 2017) and suggests higher tissue tolerance to leaf Na + ion in Genesis836. Plants maintain adequate growth by keeping low Na + levels in the cytoplasm (Munns and Tester, 2008) or storing Na + away from photosynthetic tissues. Principally, this is achieved by excluding Na + from the cell cytosol via plasma membrane Na + /H + antiporters (SOS1) or sequestration of Na + /K + into vacuoles via tonoplast Na + / H + antiporters (NHXs) while using energy from H + -ATPase in the plasma membrane or V-ATPase (VHA and AVA) and H + -PPase (VP) in intracellular compartments (Bassil and Blumwald, 2014). Our results delineated the upregulation of a tonoplast-localised Na + / H + antiporter (NHX2) under salt-stressed Genesis836 in combination with Na + /H + antiporter (NHD1) and genes encoding vacuolar proton ATPase and assembly proteins (VHA-A2, VHA-B-like and VMA21like), all exhibiting ↑Genesis836 expression, suggesting improved Na + sequestration into vacuoles. Nevertheless, Rupali uniquely upregulated vacuolar H + -translocating inorganic pyrophosphatase (VP1), along with two genes encoding plasma membrane H + -ATPase (HA1) with ↑Rupali expression. Coordination between Na + antiporters and vacuolar H + -ATPase and H + -PPase activities is critical for Na + compartmentalisation in vacuoles (Hasegawa, 2013). In addition, the electrochemical gradient of H + generated by vacuolar H + -ATPase and H + -PPase might have provided additional energy for NHX activity to transport Na + against high vacuolar concentrations in Genesis836 (Peleg et al., 2011).
High-affinity K + transporter proteins (HKTs) also reduce shoot Na + by removing Na + from the xylem and keeping it to roots (Davenport et al., 2007). However, in Arabidopsis leaves, AtHKT1 is expressed in the plasma membrane in xylem parenchyma cells, selectively unloading Na + directly from xylem vessels to xylem parenchyma cells and protecting plant leaves from salinity stress (Sunarpi et al., 2005). Moreover, time-and tissue-dependent expression of AtHKT1 determines Na + distribution in plant organs/ tissues (Hamamoto et al., 2015). We observed downregulation of HKT1 in Rupali under stress, but the gene showed ↑Rupali expression compared to Genesis836. Similarly, a tolerant rice genotype had lower HKT1 expression in shoots and roots than a sensitive genotype (Kader and Lindberg, 2005;Kader et al., 2006;Zhou et al., 2016). In wheat, the Nax1 and Nax2 loci (HKTs) prevented higher Na + concentrations in leaf blades by unloading more Na + from the xylem into the leaf sheath (James et al., 2006;James et al., 2011). Higher expression of HKT1 in Rupali tissues indicates a higher Na + influx into the cytosol of mesophyll cells, making it more salt-sensitive. Similarly, HKT1 expression may also differ between chickpea tissue types, e.g., petioles and leaflets, affecting the unloading and storage of more Na + into petioles and avoiding high Na + concentrations in leaf blades. Therefore, investigating the cell-specific expression of HKT1 will help understand how its expression influences Na + distribution in Genesis836 and Rupali leaf tissues.
Salt overly sensitive (SOS) is another important signalling pathway regulating Na + efflux from root cytosol (Ji et al., 2013). The pathway comprises a Ca 2+ -binding protein SOS3 (calcineurin B-like protein-CBL), interacting with SOS2 (CIPK24) to form a protein kinase complex that upregulates the expression of SOS1, a putative plasma membrane Na + /H + antiporter responsible for Na + efflux from the cytoplasm (Shi et al., 2002;Ji et al., 2013). Moreover, SOS1 also requires an H + gradient generated by the plasma membrane H + -ATPase. Although our study found increased expression of SOS3, SOS3-like (CIPK10 and CIPK11), and H + -ATPase (HA1) in Rupali leaves and SOS2/CIPK24 in Genesis836, SOS1 was not detected as differentially expressed. SOS2/CIPK24 forms a protein complex with CBL10 and then interacts with NHXs to regulate Na + transport across the tonoplast in shoots (Qiu et al., 2004;Kim et al., 2007;Quan et al., 2007;Luan et al., 2009). Both SOS2/CIPK24 and CBL10 had ↑Genesis836 expression, which suggests facilitated NHX activity to sequestrate more Na + into the vacuoles of Genesis836 leaves.
Most DEGs related to Cl − transport exhibited ↑Genesis836 expression, such as genes encoding Cl − conductance regulatory protein ICln protein and Cl − channel proteins (CLC-B and CLC-C), whereas only one gene encoding Cl − channel (CLC-C) delineated ↑Rupali expression. CLC-C is a vacuolar Cl − transporter linked to salt tolerance in several studies (Jossier et al., 2010;Wei et al., 2013;Henderson et al., 2014). Hence, upregulation under salinity and higher expression of CLC-C in Genesis836 might have contributed to enhanced tissue tolerance to Cl − . These DEGs indicate that the ion transport mechanisms differ between the two genotypes and appear to function better in the tolerant chickpea genotype. All the potential candidate genes involved in ion transport with known physiological and/or molecular regulation are described in Figure 8.
Salt stress also reduced photosynthesis in salt-sensitive Rupali through photosystem II damage (Khan et al., 2016;Khan et al., 2017) due to structural damage in the chloroplasts from increased Na + accumulation in mesophyll cells (Kotula et al., 2019). Rupali under salt stress downregulated several DEGs involving the photosystem II reaction centre, thylakoid lumen, and plastid proteins, while genes encoding light-harvesting complex I chlorophyll A/B-binding proteins exhibited ↑Rupali expression. In contrast, salt-stressed Genesis836 downregulated most of the lightharvesting chlorophyll A/B-binding proteins and photosystem Iassociated and chloroplast-associated proteins, and upregulated genes related to plastid development. In addition, several genes exhibiting ↑Genesis836 expression, mainly related to plastid development, RuBisCo binding and function, and chlorophyllbinding proteins, may contribute to salt tolerance in Genesis836. Various studies have reported the effect of salt stress on photosynthesis and photosystem damage in diverse crops (reviewed by Zahra et al., 2022); however, the combination of physiological and molecular insights remains limited and our study assessed both aspects of salt tolerance in two chickpea genotypes exhibiting contrasting phenotypes.
Ascorbate-glutathione is a crucial ROS-scavenging pathway that responds to oxidative stress and correlates with improved photosynthesis and salt tolerance in plants (Tang et al., 2015). The main components of the ascorbate-glutathione pathway are genes encoding glutathione S-transferase (GST), glutathione reductase (GR), glutathione peroxidase (GPX), and ascorbate peroxidase 3 (APX3).
DEGs encoding GST were upregulated in Genesis836 under salt stress. DEGs encoding GR and glutaredoxin family proteins (GRX) were highly expressed in Genesis836 tissues, whereas Fe superoxide dismutase 2 (FSD2) and GPX5 and GPX8 had ↑Rupali expression. GR mainly operates in photosynthetic tissues (Gill et al., 2013), playing an essential role in the glutathione redox state to protect rice plants against salt stress . Thus, these results support that Genesis836 improves photosynthesis via maintaining better functional photosystem against ion toxicity and ROS damage than Rupali.
Our results primarily presented three groups of protein kinases among the DEG list: calcium-dependent protein kinase (CDPK/ The schematic representation of some potential candidate genes involved in various biological mechanisms and differentially expressed between the two genotypes. CPK), CBL-interacting protein kinase (CIPK), and mitogen-activated protein kinase (MAPK). In our study, most of the DEGs corresponding to CIPKs and CPKs (CPK6, CPK30, CPK8, CIPK3, CIPK12, CIPK25, CPK33, and CIPK9) exhibited ↑Rupali expression. CPK30 is a positive regulator of stress response in barley (Sheen, 1996). CDPKs are generally positive regulators of abiotic stress response and enhanced expression of CDPKs increased salt tolerance (Boudsocq and Sheen, 2013;Thoday-Kennedy et al., 2015). MAPK6 exhibited ↑Genesis836 expression and its overexpression correlated with increased salt tolerance in Arabidopsis (Ichimura et al., 2000). Lower expression of MAPK3 and MAPKKK5 in Genesis836 is supported by earlier findings where the activation of MAPK3 mediated salt sensitivity in rice via signalling interactions and ultimately increased plant survival (Li et al., 2014). Overall, the protein kinases identified in our study were either downregulated under salinity or had comparatively lower expression in Genesis836 than Rupali, except CPK13 and MAPK6, indicating a possible role of protein kinases in genotypic variation under salt stress.
In addition, the genotype-dependent salt-responsive genes represent another category of DEGs that responded differently to salt stress in Genesis836 and Rupali. These DEGs mainly comprised heat shock, transport, and regulatory proteins, which were induced in Genesis836 and repressed in Rupali. DEGs within the salinity tolerance QTLs were involved in various processes for salt tolerance, e.g., photosystem activity, antioxidative pathway, and stress response such as MYB TFs, dehydration-responsive proteins, and heat shock proteins. Furthermore, the salt-sensitive Rupali contained many stop_gained mutations, including two genes (Ca20374 and Ca17320) with a known role in salt stress tolerance, which will be useful for future investigations.
Overall, we found that genotypic variability in salt tolerance is attributed mainly to differences in the expression of genes involved in Na + transport and photosynthetic maintenance. Moreover, we identified uncharacterised DEGs without known functions in abiotic stress tolerance, serving as ideal candidates for further experimental validations. Hence, this study highlighted the role of potential candidate DEGs and their regulatory networks in salt tolerance, as described in Figure 8 based on their biological role, which can be used to breed salt-tolerant chickpea cultivars.

Conclusion
This study presents a comprehensive analysis of the leaf transcriptome of two chickpea genotypes contrasting in response to salt stress. The comparative differential expression of genes between salt-sensitive Rupali and salt-tolerant Genesis836, together with previously identified physiological mechanisms, extends our understanding of the cultivar-specific molecular mechanisms contributing to salt tolerance in chickpeas. Most of the DEGs identified in the tolerant and susceptible phenotypes were related to sodium transport, photosynthesis, stress signalling, cell redox homeostasis, and heat shock proteins, which serve as a valuable repository of candidate genes for developing salt-tolerant chickpea varieties.

Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.