The genotype was developed at ICAR-Central Soil Salinity Research Institute, Karnal, India, through conventional breeding method (bulk selection) involving the variety HD 1982 (released in the year 1975 for north eastern plain zone, pedigree: E5557/HD H45) as a female and Kharchia 65 (released in the year 1970 for salt affected lands in all the wheat growing zones, pedigree Kharchia local/EG 953) as a male parent. The selected genotype is tall (140–155 cm) in nature and exceptionally unique in salt tolerance, having pale green foliage and higher early seedling vigor similar to their parent Kharchia 65. It possesses high yielding ability along with higher levels of tolerance to salinity, sodicity, and water-logging compared with their parent Kharchia 65. This conclusion was substantiated by 4 years of observations of the All India Salinity/Alkalinity Tolerance Screening Nursery (SATSN) conducted in different salt-affected agro-ecosystems during 2004–2005 to 2008–2009 (Supplementary Table S1). In the previous literature, Dr. Abdul Mujeeb-Kazi (distinguished scientist and ex program leader of CIMMYT) had already mentioned that Kharchia 65 is the most salt-tolerant bread wheat in the world (Wang et al., 2003).

This germplasm already has been registered (IC408331; INGR09087), in ICAR-National Bureau of Plant Genetic Resources (NBPGR), New Delhi, salt and waterlogging tolerance. The salt tolerance of KRL 3–4 is characterized by better regulation over lower uptake of sodium and highest uptake of magnesium salt-affected soils (Singh et al., 2010). Additionally, it possesses higher tolerance to boron, iron, and aluminum concentrations in the soil solution, and their full potential reflected on highly sodic condition (pH >9.3), where it performs as a halophyte, and it competes with the parent Kharchia 65 in yield as well as in other physiological traits. KRL 3–4 also had the highest grain Ca, Mg, Fe, Zn, and S concentrations compared with most of the Indian cultivars, as it substantiates in an experiment conducted in hostile soil conditions (pH: 4.5–9.5) across the six different agro-ecosystems for 2 years (Khokhar et al., 2018).

In our study, a wheat genotype KRL 3–4 was selected to understand the transcriptional changes in response to sodicity stress. The genotype was grown in two concrete lysimeters (dimension: 2 m × 2 m × 1 m) situated at ICAR-Central Soil Salinity Research Institute, experimental farm (29°42′31″N, 76°57′6″E) Karnal, India. The first lysimeter, filled homogeneously with sandy loam soil (Inceptisols, pH: 7.4, ESP:6.3, ECe: 0.9 dS m⁻¹, organic carbon: 2.3 g kg⁻¹ soil) was considered as control. The other lysimeter was filled homogeneously with incubated synthetic sodic soils formulated by mixing the carbonate and bicarbonate salts in the soils (pH: 9.5, ESP: 74.6, ECe: 1.6 dS m⁻¹, organic carbon: 2.1 g kg⁻¹ soil). The sowing seeds were surface sterilized in 1% sodium hypochlorite (NaClO⁻) for 10 min and then washed with distilled water several times and planted in the lysimeters. Plants that were grown in normal soil were considered control plants, and other plants grown in the sodic lysimeters (stress intensity: pH-9.5) were considered treated plants. At the anthesis stage, fully matured flag leaf samples were taken from both control and treated plants with two biological replicates. The leaf samples were immediately frozen in liquid nitrogen and stored at −80°C until analysis.

Total RNA was isolated from the two biological replicates of flag leaves from normal and treated plants (after 12-h exposure to sodicity stress) using RNeasy Plant Mini Kit (Qiagen) according to the instructions of the manufacturer. Equal amounts of total RNA of each two biological replicates were used for the RNA sequencing. Initially, the concentration of RNA was checked by nanodrop and later purity through agarose gel electrophoresis. Furthermore, RNA integrity was estimated by Agilent Bioanalyzer 2100 system (Agilent Technologies Co. Ltd., Beijing, China). Accordingly, the samples with RIN value higher than 9 were used for transcriptome sequencing. After quality control (QC) of the RNA samples, poly (A) enrichment, RNA fragmentation, random hexamer-primed cDNA synthesis, linker ligation, size selection, and PCR amplification reactions were performed to prepare cDNA libraries for each sample. Finally, the qualified libraries fed into HiSeq sequencer according to its effective concentration and expected data volume. Illumina HiSeq 2500 platform was used to sequence libraries, and 150-bp paired end reads were generated.

Raw reads of control and treated samples were evaluated for their quality using FastQC (version 0.11.8) (http://www.bioinformatics.bbsrc.ac.uk/projects/fastqc/). The adaptor sequences and low-quality sequence reads were removed from the data sets using trim galore (version 0.6.2) considering the threshold parameters, including phred score 33, quality score 20, and length 20 (https://github.com/FelixKrueger/TrimGalore/releases). Raw sequences were transformed into clean reads after data processing. The quality of cleaned reads was also checked through base quality score distribution, sequence quality score distribution, average base content per read, and GC distribution in the reads. The obtained high-quality reads were mapped to the reference genome of Triticum aestivum, i.e., IWGSC v.1.0 (ftp://ftp.ensemblgenomes.org/pub/plants/release-42) using bwa (version 0.7.5, http://bio-bwa.sourceforge.net/) with default parameters. These mapped reads were spliced using Cufflinks software based on the reference genome sequence. Quantification of the gene expression levels was estimated as fragments per kilo base of transcript per million fragments mapped (FKPM). Differential expression analysis between control and treated samples was performed using the DESeq R package (1.10.1) (Anders and Huber, 2010). DESeq provides reliability for determining differential expression in digital gene expression data using a model based on the negative binomial distribution. The resulting p-values were adjusted using Benjamin and Hochberg’s approach for controlling the false discovery rate. Genes with an adjusted p-value <0.05 were found using DESeq and were assigned as differentially expressed. A fold change ≥2 and FDR <0.01 were considered as the thresholds for determining the differential expression of a gene.

Gene ontology (GO) enrichment analysis of the DEGs was implemented using topGO R package (https://bioconductor.org/packages/release/bioc/html/topGO.html). The GO terms, i.e., biological process (BP), cellular component (CC), or molecular function (MF) were assigned to the DEGs. Kyoto Encyclopedia of genes and Genomes (KEGG) is a database resource for understanding high-level functions and utilities of the biological systems, such as the cell, organism and ecosystem from molecular-level information, especially large-scale molecular datasets generated by genome sequencing and other high-throughput experimental technologies (http://www.genome.jp/kegg/). KOBAS software (http://kobas.cbi.pku.edu.cn/) was used to analyze the statistical enrichment of differential expression of genes in KEGG pathways. Furthermore, MapMan (version 3.5.1; http://mapman.gabipd.org/web/guest) was used for pathway analysis of DEGs as p-value <0.05.

Alternative splicing (AS) events for identified DEGs under sodicity stress were determined using an AStalavista web tool (version 3; http://genome.crg.es/astalavista/) with default parameters. The results obtained for all AS events were further analyzed to unravel the molecular mechanisms of AS during exposure of wheat crop to the sodicity stress.

The freebayes tool (Garrison and Marth, 2012) and SNP effect predictor (Cingolani et al., 2012) were used to identify and annotate the SNPs, respectively. SNPs and INDEL variants were filtered through bcftools. Based on genomic locations, the variants were annotated as intronic, exonic, 5′UTR, and 3′UTR.

The candidate genes for sodicity stress tolerance were identified on the basis of similarity (blastx with similarity >80% and e-value <0.001) with published or validated genes for their role in sodicity stress tolerance (Singh et al., 2018). The heatmap was made using ggplot2 R package (https://cran.r-project.org/web/packages/ggplot2/index.html). Furthermore, the chromosomal localization of the identified candidate genes was performed using MapChart (Voorrips, 2002), and the 2-kb upstream region was used for the identification of cis-acting regulatory elements using PlantCARE database (Rombauts et al., 1999).

In order to understand the mechanism involved at the transcriptional level for sodicity stress tolerance, transcriptome analysis was carried out for control and treated samples of wheat flag leaves. The libraries were sequenced to obtain 150-bp paired-end reads using the Illumina high-throughput sequencing platform. In total, approximately 37.24 Gb of raw data was generated from the two biological replicates of the wheat flag leaf. After trimming, the volume of the total clean reads for all the samples was 36.54 Gb. The clean reads of each sample was more than 8.75 Gb, and the percentage of Q20 bases was ≥97.40%. Sequencing quality distribution was examined over the complete sequence length to detect the sites (base positions) with an unusually low sequencing quality, where incorrect bases may be incorporated at abnormally high levels. Approximately, 97.60%–99.06% of reads are mapped to the wheat reference genome (IWGSC RefSeq v1.0, http://www.wheatgenome.org/), whereas GC content varied from 51% to 53% for the control and sodicity stress treated samples (Table 1).

In total, we identified 97,514 differentially expressed genes (DEGs) between control and treated flag leaf wheat samples. The genes with expression level greater than two fold (adjusted p-value < 0.05) were considered as significant DEGs under sodicity stress condition with reference to the controlled conditions. At this threshold, there were only 1,980 significant DEGs between control and sodicity stress conditions (Supplementary Table S2). Among these genes, 872 DEGs showed upregulation, and 1,108 showed downregulation under sodicity stress condition (Figure 1A).

Furthermore, GO enrichment analysis was carried out to characterize the main biological functions of DEGs in wheat flag leaves under sodicity stress condition. All the DEGs were annotated into three classes, i.e., BP, MF, and CC. Annotation of DEGs revealed that a total of 1,384 genes were assigned 2,267 GO terms corresponding to 502 BP, 638 CC, and 1,127 MF (Supplementary Table S3). The dominant terms in CC category were “membrane,” “membrane part,” “cell,” “cell part,” and “organelle,” whereas in MF, “catalytic activity,” “binding,” “transporter activity,” “molecular function regulator,” and “molecular transducer activity” were the most dominant GO terms. However, in case of BP, most of the DEGs were classified in “response to stimulus,” “cellular process,” “metabolic process” followed by “biogenesis,” and “localization” (Figure 1B).

As expected, the most enriched BP terms for over-presented DEGs were associated with metal ion transport, response to auxin, water and oxidative stress, transmembrane transport, oligopeptide transport, lipid transport, carbohydrate metabolic process, plant-type cell wall organization, which acted as indicators of significant biological processes underlying the specific sodicity stress responses of plants. The most enriched MF terms for over-presented DEGs were metal ion binding, ATP binding, calcium ion binding, catalytic activity, channel activity, iron and magnesium ion binding, oxidoreductase activity, transferase activity, voltage-gated anion channel activity, etc. Meanwhile, with regard to the over-representation of CC terms among the DEGs, we found extracellular region, membrane, intrinsic to membrane, integral to the membrane, and intrinsic to plasma membrane to be the most enriched.

Furthermore, a large number of TF genes were found to be differentially expressed under sodicity stress condition. In total, 985 DEGs were encoded for 53 TF families (Figure 2A). The important TF families that were represented by DEGs included bHLH, ERF, bZIP, NAC, MYB-RELATED, MYB, and WRKY. Among the TF families, the highest number of genes belonged to the bHLH family (109) followed by NAC (99), ERF (89), C2H2 (59), WRKY (58), MYB_RELATED (56), and GRAS (37), and so on. We also checked the expression of these TF in control as well as in sodicity treated samples and visualized using heat map. Moreover, we have also identified 37 DEGs encoding for serine/threonine kinase family proteins, out of which 22 showed upregulation, and the remaining 15 showed downregulation under sodicity stress. Similarly, genes encoding for other stress tolerance and signaling-associated proteins were also identified, including 22 genes of calmodulin ion binding and calcium signaling, 11 genes for metal ion transport, 4 genes of ABC transporter, 4 genes of SWEET sugar transporter, and 1 gene encoding potassium transporter. In this study, 43 DEGs appear to protect plant cells from oxidative damage by ROS scavenging, such as peroxidase genes (26) and glutathione S-transferase genes (17). For peroxidase genes, 6 of the 26 DEGs were upregulated, and the other 20 genes were observed to be downregulated. Besides, the peroxidase gene family, a number of glutathione S-transferase genes were found differentially expressed. Most of these genes (12/17) were noticed to be upregulated during sodicity stress (Figure 2B; Supplementary Table S4).

The KEGG pathway analysis was performed to unveil the active biological pathways involved in DEGs under sodicity stress in KRL3-4 flag leaf. KEGG is a database resource for interpretation of large-scale molecular datasets generated by genome sequencing for understanding high-level functions and utilities of the biological system (Kanehisa et al., 2017). The pathway analysis showed that out of 127 KEGG pathways, 40 pathways were observed to be significantly enriched at a p-value <0.05 covering the five main KEGG categories of metabolism, environmental information processing, genetic information processing, organismal systems, and cellular processes (Figure 3). The KEGG enrichment pathways indicated that the specific DEGs, found under sodicity stress in flag leaf, were widely enriched in the pathways of metabolic, biosynthesis of secondary metabolites, phenylpropanoid biosynthesis, starch and sucrose metabolism, photosynthesis, and plant hormone signal transduction. According to previous studies, these pathways were found related to abiotic stress tolerance. The list of DEGs involved in 40 KEGG enriched pathway is given in Supplementary Table S5.

The differentially expressed sodicity stress-responsive genes were also visualized using MapMan in order to find out alterations in metabolic pathways during sodicity stress (Supplementary Table S6). The MapMan analysis revealed that lipid metabolism and protein homeostasis, nutrient uptake, etc., were the most enriched pathways under sodicity stress. In lipid metabolism, particularly, the expression of 12 genes involved in lipid degradation, i.e., 9 phospholipases A1 (PC-PLA1) and 3 phospholipase A2 (pPLA2-II) were highly upregulated in response to sodicity stress, whereas the genes related to fatty acid biosynthesis and elongation (11 genes), i.e., 3-ketoacyl-CoA reductase (KCR) and 3-ketoacyl-CoA synthase (KCS) showed mixed expression. Glycerol-3-phosphate acyltransferase (GPAT1-3) and acyl-CoA:diacylglycerol acyltransferase (DGAT2), involved in glycerolipid biosynthesis, were found to be up and downregulated, respectively. Similarly, in the case of protein homeostasis, the genes encoding Bowman-–Birk protease inhibitor, pepsin-type protease, amino acid hydrolase was found to be upregulated under sodicity stress. On the other hand, we have also observed the upregulation of phosphate signaling regulatory protein (SPX), phosphate transporter 1 (PHT1), and nitrogen limitation adaptation (NLA) genes in nutrient uptake pathway. We identified one gene for PHT1 (Traescs2a02g047300), one for NLA (Traescs5b02g373000) and seven for SPX (Traescs2a02g169600, Traescs2b02g195900, Traescs2d02g177100, Traescs7a02g376200, Traescs7a02g554100, Traescs7b02g478000, and Traescs7d02g372600) ( Figure 4). The cellular overview pathway exposed that the genes encoding for cell wall organization showed mixed expression under sodicity stress. For example, the genes encoding glutaredoxins (GRXs) and arabinogalactan protein (xylogen) were found highly upregulated while alpha class expansins and glucuronosyltransferase (GUX) showed downregulation. GRXs are small redox proteins, which use glutathione to catalyze the reduction of disulfide bonds of substrate proteins to maintain cellular redox homeostasis. Furthermore, diverse functions such as transcriptional regulation of defense responses, flower development, oxidative stress response, redox signaling, hormonal regulation, iron homeostasis, and environmental adaptation have been reported for various plant GRXs.

Alternative splicing (AS) is an important post-transcriptional mechanism that increases the complexity and diversity of transcriptome and proteome in higher eukaryotes. In plants, AS is found to be involved in a range of functions including plant growth and development, and responses to biotic and abiotic stresses. We have categorized the DEGs into four primary classes of AS events that are retained intron (RI), skipping exon (SE), alternative 5′ splice site (A5SS), and alternative 3′ splice site (A3SS) (Figure 5B). Our study revealed, 271 AS events in the DEGs under sodicity stress. In total, 237 (11.96%) genes out of 1,980 have shown alternative splicing of which the highest number belonged to alternate 3′ acceptor (36.16%) followed by intron retention (26.19%), alternate 5′ donor (16.60%) and exon skipping (13.28%) (Figure 5A). Functional categorization of these alternatively spliced DEGs revealed their involvement in different stress-responsive pathways, including the plant hormone signal transduction, calcium signaling pathway, biosynthesis of secondary metabolites, glutathione metabolism, AMPK signaling pathway, etc. In terms of GO biological processes, these were found to be involved in metal ion transport, response to water, transmembrane transport, lipid and carbohydrate metabolic process, cell wall organization, and abscisic acid-activated signaling pathway, while, the major molecular functions represented by these alternatively spliced DEGs involve abscisic acid binding, metal ion binding, iron ion binding, phospholipase A1 activity, transmembrane transporter activity, etc.

Trait variations are primarily due to sequence level changes in the coding region. The wheat genotype KRL3-4 included in our study represents one of the most sodicity tolerant genotypes of the world. Therefore, one might expect that alignment of the flag leaf transcripts from the control and susceptible conditions against a reference genome of the cv. Chinese spring (a relatively susceptible variety) may unravel novel SNPs and INDELs that may have a role in sodicity stress tolerance. A total of 36.64 Gb clean reads was obtained after RNA sequencing of both control and sodicity treated wheat flag leaf samples to discover SNP and INDELs. On an average, one variant was called within every 17 kb in the treating sample and 32 kb in the control sample. In total, 2,35,985 and 1,63,420 polymorphisms (SNPs and INDELs) were identified in the control and treated samples, respectively. In the case of control samples, 39,295,381, and 39,200 SNPs were found as missense, nonsense, and silent mutations, respectively, whereas in sodicity treated wheat flag leaf samples, 26,575, 234, and 26,377 were observed as missense, nonsense, and silent mutations, respectively.

SNPs and INDELs were mapped to exonic regions that included 8.02% of the total reads (63,239 read counts) in treated and 12.95% (94,262 read counts) in control samples. Moreover, 4.10% of the total reads (32,339) were mapped on 5′UTR and 3′UTR regions in treated samples (Figure 6). In addition, polymorphism due to transitions (Ts) and transversions (Tv) was also determined in each sample. In control, the Ts and Tv were estimated at 7,17,799 and 3,68,641, while in the sodicity treated sample, 4,94,965 and 2,54,817 transition and transversion were detected, respectively. In addition, amino acid substitutions were also identified for SNPs mapped to IWGSC wheat reference genome in control and sodicity treated samples. In control sample, the major amino acid substitutions observed were alanine to valine (1,584), alanine to threonine (1,386), and valine to isoleucine (1,004), whereas in treating samples, alanine to valine (1,031), alanine to threonine (902), and valine to isoleucine (709) were noticed to be the major amino acid substitutions.

We have identified 19 sodicity stress responsive genes in our dataset. These candidate genes were identified on the basis of similarity with published validated salt-responsive genes (Table 2). The identified candidate genes encoded for a diverse group of proteins including LEA, salt-responsive protein, potassium/sodium cation transporter (HKT), calcium-dependent signaling calcium sensor, glutathione-S-transferase, serine/threonine protein kinase, etc. These candidate genes were found to be upregulated under sodicity stress. In terms of GO annotation, these genes were observed to be involved in metal ion binding, hydrolase activity, protein dimerization activity, diacylglycerol O-acyltrasferase activity, DNA binding, polygalacturonase activity, etc. Some candidate genes corresponded to different TF families including C2H2 (TraesCS5D02G026300), NAC (TraesCS2D02G061500), WRKY (TraesCS5D02G145800), AP2/ERF (TraesCS5D02G245300), ARF (TraesCS3D02G276600), and bHLH (TraesCS5B02G054800) and showed upregulation in case of sodicity stress compared with control.

Chromosomal distribution analysis showed that 19 identified candidate genes were unevenly mapped in 10 out of 21 wheat chromosomes. On chromosome 5, a total of seven candidate genes were located that represent 36.84% of the entire candidate genes identified. Besides, four genes were located on chromosome 3, whereas on 3B, 3D, and 6D chromosomes, one gene each was located. On the other hand, 2D, 4D, 6A and 7A chromosomes were found to contain two genes each. In contrast, no genes were found in the remaining 11 chromosomes. The sodicity responsive candidate genes were observed to be unevenly distributed among the sub-genomes A, B, and D, with 1, 4, and 10 members representing 5.26%, 21.05%, and 52.63% of the total candidate genes, respectively (Figure 7). The chromosomal distribution of the 19 candidate genes identified in our study is similar to a previous study of bread wheat where the maximum candidate salt-responsive genes were found on D sub-genome than on A and B sub-genomes (Singh et al., 2018) suggesting that sodicity stress controlling genes is preferentially located on the D sub-genome of wheat. The difference in expression patterns of the 19 candidate genes between control and sodicity stress was represented using heatmap (Figure 8). Heat map analysis indicated that nine genes showed contrasting expression between control and sodicity stress, while the rest of the genes were having a mixed expression level. Therefore, this analysis suggests that these genes might have a more certain functional role in sodicity stress tolerance.

The cis-acting regulatory elements (CARE) have been shown to play a major role in the regulation of gene expression at the transcriptional level in plants. Indeed, these elements are the central component in plant abiotic stress responses since transcription regulation involves association between TFs and, in particular, cis-acting regulatory elements (CAREs) of a specific gene (Kaur and Pati, 2016). From the 19 candidate genes identified in this study, the majority belonged to TF categories, which have cis-binding motifs that binds to cis-elements of the downstream stress-responsive genes and modulate their expression. Identified 19 putative candidate genes were analyzed for their cis-regulatory elements using PlantCARE with default parameters (Lescot et al., 2002; Doi et al., 2008). A total of 17 CARE were A-box, ABRE3a, CAAT-box, CCGTCC-motif, CCGTCC-box, CGTCA-motif, DRE-core, G-Box, MTB, STRE, Sp1, TGACG-motif, as-1, AAGAA-motif, MBS, MYC, and TATA-box were observed among identified sodicity responsive candidate genes that have a significant role in stress tolerance. The cis-regulatory elements located in promoter regions of putative 19 sodicity responsive genes with their functional description are given in Table 3.

In our study, 604 DEGs were mapped using the MapMan tool to analyze the impact of sodicity stress on wheat flag leaves. Our results suggested that in response to sodicity stress, a large number of genes were found to be associated with lipid metabolism, protein homeostasis, cell wall synthesis, and modification. In addition, transcriptional regulation of defense responses, flower development, oxidative stress response, redox signaling, hormonal regulation, iron homeostasis, and environmental adaptation were also reported for various plant GRXs.

Lipids are structural, storage, and signaling molecule components in plant cells. They are the principal constituents of plant membranes and the cuticle (Harwood, 1980). In this study, Mapman analysis showed upregulation of various lipid metabolism genes including 12 genes encoding for phospholipase, 9 genes encoding for phosphocholine phosphatase, and 1 gene encoding for glycerol-3-phosphate acyltransferase under sodicity stress. All three together are involved in the synthesis of triacylglycerol (TAG) that is an important component of lipids in cell membranes. Prior studies have shown that different stress conditions induce TAG production in vegetative tissues leading to upregulation of several key genes involved during TAG biosynthesis. Abscisic acid (ABA), a plant hormone, accumulates under stress and functions as a signaling molecule to regulate plant development and metabolic pathways. A transcription factor, ABSCISIC ACID INSENSITIVE 4 (ABI4), that is a key component in the ABA signaling pathway was upregulated under stress conditions and was observed to bind the DGAT1 promoter, which induces its expression in Arabidopsis (Lu et al., 2020). Therefore, this important observation supports the fact that salt stress generally disturbs the integrity of the plasma membrane and the high level of TAG might be contributing to salt stress tolerance in KRL3-4 by maintaining integrity of plant membranes. Another component involved in lipid metabolism is 3-ketoacyl-CoA synthase (11), which was found to be downregulated under sodicity condition. Earlier research showed that cuticular wax was found to be associated with plant resistance in response to abiotic stress. 3-Ketoacyl-CoA synthase (KCS) catalyzes the biosynthesis of very-long-chain fatty acid (VLCFA) wax precursors. 3-Ketoacyl-CoA (KCS) synthase encoded by three salt-responsive genes (LOC_Os02g11070, LOC_Os05g49900, and LOC_Os02g56860) have a role in wax biosynthesis (Yuenyong et al., 2018). Furthermore, a recent study in navel orange showed that a novel KCS family gene named CsKCS6 regulate the production of fatty acid precursors involved in wax synthesis and delivers tolerance of transgenic Arabidopsis plants in abiotic stress. Their experiment demonstrates that the expression of CsKCS6 significantly increased the amount of VLCFAs in the cuticular wax in transgenic Arabidopsis plants. The mechanism of KCS6-like genes to influence cuticle permeability in transgenic plants involves ectopic expression of CsKCS6, which decreased the amount of primary alcohols and aldehydes suggesting that wax biosynthetic genes might be suppressed by ectopic expression of CsKCS6 in transgenic plants. Hence, a significant decrease in leaf chlorophyll permeability, water loss rate, ion leakage, increased root length, and survival rate were observed in the transgenic plants under drought and high salt stress indicating that ectopic expression of CsKCS6 can reduce the cuticle permeability and enhance plant tolerance to drought and high salt (Guo et al., 2020). Furthermore, two peroxygenases/caleosin encoded by CLO/PXG-like genes were also found, showing downregulation during sodicity-induced condition in our analysis. Caleosins are a calcium sensor carrying one EF-hand calcium-binding motif. A transcriptomic study on fungus has shown that the CLO/PXG-like genes have complex patterns of developmental and tissue-specific expression and exhibited response toward biotic and abiotic stresses and play key role in lipid metabolism, signaling, reproduction, and pathogenesis (Rahman et al., 2018). In another study on rice, caleosin family protein 5 (OsClo5) or peroxygenase (PXG) modulates salt stress tolerance via calcium-dependent membrane associated hemoprotein that catalyzes the PXG reaction (Blée et al., 2012) indicating its potential roles in enzymatic activities and Ca²⁺ signaling in plants. OsClo5 interacts with OsDi19-5 and inhibits transcription of OsUSP and OsMST and that activity leads to the enhanced sensitivity of rice seedlings in response to salt stress compared with the wild-type, where OsClo5 with T-DNA knockdown mutation distorted heterodimer formation with OsDi19-5 and promotes the weak transcriptional repression effect on the two target genes, i.e., OsUSP and OsMST. This confers the reduced sensitivity to salinity or the enhanced tolerance to salt stress in rice seedlings in comparison with wild-type. Thus, OsClo5 negatively affected salt stress tolerance in rice in cooperation with the transaction repression control of OsDi19-5 to two target genes OsUSP and OsMST in rice (Jing et al., 2021).

In addition to lipid metabolism, many DEGs encoding the Bowman–Birk protease inhibitor (BBI), such as TraesCS3D02G034200, TraesCS3A02G046000, TraesCS3A02G046100, TraesCS1D02G022000, TraesCS5A02G484800, TraesCS5B02G498100, and TraesCS5D02G498300, showed upregulation. The Bowman–Birk inhibitor (BBI) is one of the subfamilies of serine protease inhibitors. It has been reported to function in controlling abiotic stresses, such as salinity and drought stresses (Prakash et al., 1996). In a previous study on wheat, a salt-responsive gene WRSI5 that contains a Bowman–Birk domain was found responsible for tolerance to salt stress. It was characterized from the salt-tolerant cultivar Shangrong No.3 (SR3) and salt-susceptible cultivar JN177. The WRIS5 gene functions as trypsin inhibitor as it exhibits a high level of trypsin-inhibiting activity, which signifies that BBI suppressed the trypsin activity, thus, promoting elevated WRSI5 expression level in SR3 roots in salt, drought, or oxidative stress. Furthermore, the SR3 allowed limited transport of K⁺ than Na⁺ from root to shoot leading to low Na⁺ concentration in SR3 leaves resulting in increased leaf growth rate. In addition, overexpression of WRSI5 in Arabidopsis improves seedling growth at 150 mM NaCl condition, suggesting that WRSI5 controls plant growth rate or long-distance Na⁺ transport in SR3 plants under salt stress (Shan et al., 2008). Similarly, BBI family genes expressed in our study might function in sodicity tolerance.

Furthermore, our study showed upregulation of 13 glutaredoxin (GRX) genes in response to sodicity stress. The GRX family proteins are involved in various cellular functions, such as redox regulation and protection under abiotic and biotic stresses (Lillig et al., 2008). A recent report showed overexpression of glutaredoxin (OsGrx_C7) ubiquitously expressed in rice including root and shoot, thus indicating positive response in salt induced stress. However, silencing of this gene exhibited increased sensitivity of rice plants to salt stress. Furthermore, OsGrx_C7 was suggested to be a positive regulator of salt tolerance by reinforcing the expression of transporters employed in Na⁺ homeostasis in overexpressing rice plants (Verma et al., 2021). Therefore, GRX genes found in our study can be considered to play a role in sodicity tolerance through Na⁺ homeostasis. Additionally, eight downregulated DEGs encoding for α-expansins were also identified in our study. Expansins are key regulators of cell wall formation and required for cell enlargement under various environmental stresses (Wu and Cosgrove, 2000). Several studies have shown that the expression of expansin-encoding genes are induced under high salt and tends to promote salt tolerance in plants (Qiu et al., 2014). Expansin is a superfamily of plant proteins and comprised of four families, such as α-expansin, β-expansin, expansin-like A, and expansin-like B. The α and β-expansin proteins entail cell wall loosening activity and are involved in cell expansion and plant developmental processes. Recently, transgenic tobacco plants, including α-expansin 4 (EXPA4) overexpression and RNA interference (RNAi) mutants were assayed for their tolerance to abiotic stress. The study revealed that transgenic plants with altered EXPA4 demonstrated that RNAi mutants have increased hypersensitivity to salt and drought stress and vice-versa in overexpressed tobacco lines, consequently, leading to lower cell damage, higher fresh weight, soluble sugar, proline accumulation, and induces several stress responsive gene expression. Therefore, EXPA4 indicated its involvement in salt and drought stress response via changes in the expression levels of some stress-responsive genes (Chen et al., 2018). Moreover, downregulation of arabinogalactan proteins (AGPs) was also seen in our study. AGPs were seen involved in abiotic stress response such as salinity regulating cell wall expansion. Several studies have reported downregulation of AGPs in response to salt treatments. For example, in tomato, different AGPs were strongly repressed in salt-treated tomato roots (Ouyang et al., 2007). Similarly, two rice OsFLA (OsFLA10/18) and two wheat TaFLA3/4 genes were significantly downregulated under salinity (Faik et al., 2006; Ma and Zhao, 2010). AGPs also act as a sodium carrier through the mechanism of vesicle trafficking (Garcia de la Garma et al., 2015). Moreover, the tobacco (Nicotiana tabacum L.) cells cultured in saline medium, which were deficient in AGPs on the plasma membrane, demonstrates a reduced rate of cell enlargement and a decreased cell wall extensibility. However, non-saline tobacco cells showed higher upregulation of AGPs in salt-stress indicating an important increase rate in AGP diffusion through a much more highly porous pectic network leading to improved cell wall formation (Zhu et al., 1993; Lamport et al., 2006). This mechanism may contribute to sodium homeostasis during salt adaptation to high saline concentrations (Olmos et al., 2017). Overall, our study revealed diverse classes of genes as key mediator of salinity stress response and could be further validated using functional genomics tools to uncover their exact role in salt tolerance.

Sodic soil is known to decrease the availability of nutrients such as phosphorus required for plant growth and development (Dotaniya and Meena, 2015). In plants, the primary source of phosphorus is inorganic orthophosphate (Pi), but in sodic soil, it forms an insoluble compound with calcium ions and becomes inaccessible (Péret et al., 2011). As a result, a deficiency in Pi stimulates a set of plant adaptive responses with the purpose of decreasing Pi usage and increasing Pi uptake and recycling via upregulation of regulators and transporters involved in P-homeostasis (Ticconi and Abel, 2004; Desnos, 2008; Lin et al., 2009; Rouached et al., 2010). Various genes related to low phosphorus stress have been identified in different plant species (Rausch and Bucher, 2002; Jain et al., 2007; Wang et al., 2011). For example, phosphate transporter 1 (PHT1) gene family participates in the uptake of Pi from the soil (Muchhal et al., 1996; Mudge et al., 2002; Shin et al., 2004; Jia et al., 2011). In our study, we have observed one PHT1 (Traescs2a02g047300) gene, which showed upregulation under sodicity stress. PHT1 transporters are found to be involved in the uptake of phosphate from the soil and transport it to the shoot and mycorrhizal symbiotic interface. In Arabidopsis, rice, and maize, PHT1 expression is strongly upregulated during phosphorus starvation (Smith and Read, 2008; Huang et al., 2011). Similarly, in wheat also, two PHT1 members, i.e., TAPHT1.1 and TaPHT2.1, overexpressed in Pi-deficient condition (Shukla et al., 2016). Our study also revealed the upregulation of seven phosphate signaling regulatory protein (SPX), mapped to the 2A, 2B, 2D, 7A, 7B, and 7D chromosomes (Traescs2a02g169600, Traescs2b02g195900, Traescs2d02g177100, Traescs7a02g376200, Traescs7a02g554100, Traescs7b02g478000, Traescs7d02g372600). These genes perform multiple functions in plant tolerance to phosphorus starvation and plays a vital role in sensing the concentration of phosphorus in cytosol (Duan et al., 2008; Wild et al., 2016). During Pi starvation, SPX1 dissociates with the phosphate response 1 (PHR1) TF. PHR1 controls many Pi related including miRNA399, which targets phosphate 2 (PHO2). Furthermore, PHO2 reduction leads to PHT1 and PHO2 accumulation and, last, an increase in the plant capacity to uptake Pi and translocate it to shoots. Various genes participate in the proper functioning of this signaling pathway, such as SUMO E3 ligase SIZ1, phosphate transporter traffic facilitator (PHF1) and nitrogen limitation adaptation (NLA) (Lin et al., 2013). According to our MapMan results, we have also found one gene encoding NLA (Traescs5b02g373000), highly upregulated under sodicity stress. NLA functions of the plasma membrane to direct the degradation of PHT1s is a process necessary for Pi uptake ability of plants (Lin et al., 2013). Hence, we can conclude that although we have observed the expression of many transporter proteins and their respective genes under sodicity stress, the physiological role of Pi transporters in response to sodicity stress is still limited. Understanding the molecular mechanism involved in the regulation of phosphate transporters is very important, due to the constant lowering availability of phosphorus reserves in nature and will facilitate the development of more efficient Pi-utilizing plants.

Plants undergo adaptation during stress conditions by triggering a network of signaling cascade that leads to activation of a series of target gene expression. Studies have shown that several TFs, when binding to specific cis-elements, act as molecular regulators of gene expression during stress responses (Gujjar et al., 2014). PlantCARE analysis revealed recognition of many cis-elements in genes expressed in response to sodicity stress (Table 3). For instance, A-box, CAAT-box, CGTCA-motif, G-Box, STRE, TGACG-motif, as-1, TATA-box were present in all 19 genes followed by MYB and MYC observed in 18 genes, whereas CCGTCC motif, CCGTCC-box, DRE core, and Sp1 were seen in 17 genes, and ABRE3a, AAGAA-motif, MBS were located in 13 genes. Previous studies have reported that these cis-elements in many genes are expressed in response to abiotic stresses. For example, “CGTCA-motif” was reported in drought stress-responsive genes of trifoliate orange and suggested to play crucial role in drought response via jasmonic acid-mediated signaling (Xiong et al., 2020). Similarly, “G-box” cis-acting DNA regulatory element are recognized by the bZIP and bHLHs family TFs to mediate responses toward abiotic and biotic stresses (Qian et al., 2021). Prior studies have shown that bZIP and bHLHs family TFs bind to “G-box” in the promoter of cold, drought, and salinity-responsive genes in Arabidopsis and other plant species (Liu et al., 2008; Shah et al., 2021). In wheat, a TabHLH13 gene that responds via “G-box” elements were upregulated under salt stress (Fu et al., 2014). Another important cis-element, i.e., “STRE” is recognized by zinc-finger DNA binding TF (Estruch, 2000). “STRE” is functional in both orientations (“CCCCT” or “AGGGG”). A study conducted on yeast showed that the “STRE” cis-element located in the promoter of CTT1 and DDR2, a DNA damage-responsive gene mediates transcriptional induction in response to different stress conditions (Kobayashi and McEntee, 1993). In contrast, “TGACG-motif,” with consensus sequence “TGACG” is targeted by MeJA family TFs as it acts as JA-responsive element (Liao et al., 2017). In rice, two “TGACG-motif” (JA-responsive element) cis-elements were found in the promoter region of OsTMP14 that is related to stress tolerance. Subsequently, activation sequence-1 (“as-1”) with consensus sequence “TGACG” situated in all identified genes showed the role in abiotic stress and salicylic acid responsiveness (Banerjee et al., 2015). In the current study, two cis-acting elements (“MYB” and “MYC”) were located in 18 candidate genes that might be involved in sodicity stress tolerance. The MYB cis-element with consensus sequence “CAACAG” mostly interacts with MYB TFs. In past reports, overexpression of the wheat MYB TFs, such as TaMYB30-B, TaMYB33, TaMYB56-B, and TaSIM, import tolerance to salt and drought stresses in transgenic Arabidopsis (Yu et al., 2017), whereas MYC cis-element with consensus sequence “CATGTG” functions in salt and drought stress response through interaction with bHLH and NAC family TFs. The role of “MYC” has also been seen in drought tolerance through the dehydration-inducible expression of EARLY RESPONSIVE TO DEHYDRATION STRESS 1 (ERD1) gene in Arabidopsis. Since ERD1 promoter contains the “CATGTG” motif that interacts with NAC, which specifically binds to the MYC-like sequence “CATGTG,” the expression of NAC TFs genes, such as ANAC019, ANAC055, and ANAC072, was induced by drought, high salinity, and abscisic acid (Tran et al., 2004).

Dehydration-responsive element (DRE) core having “GCCGAC” consensus sequence is found in promoters of various environmental stress-responsive genes. DREB TFs binds at “DRE-core” (Liu et al., 1998). In Arabidopsis, overexpression of DREB1A promotes enhanced expression of downstream drought stress-responsive genes (Liu et al., 1998). In addition, DREB proteins also have a role in phytohormone signaling pathway such as abscisic acid, salicylic acid, jasmonate acid, ethylene, and gibberellic acid (Choi et al., 2002). “Sp1” is another cis-element present in 17 identified candidate genes. Homeodomain leucine zipper (HD-ZIP) TFs targets the consensus sequence “GGGCGG” of “Sp1” for abiotic stress response. The presence of “Sp1” in the promoter region of banana HD-ZIP genes have shown its participation in various functions and were found to be upregulated by drought and salinity stress (Pandey et al., 2016).

SNP identification may have a critical role in conferring sodicity stress tolerance in wheat flag leaf. We noticed 39 SNPs mapped to the important sodicity stress-responsive genes associated with various pathways, such as ROS scavenging, serine/threonine protein kinase, calcium signaling, and metal ion transporters) that are known to have a role in abiotic stress responses. From these 39, a total of 29 SNPs were mapped to ROS-scavenging genes that were differentially expressed during sodicity stress condition, thus, stimulating stress tolerance through glutathione metabolism and the biosynthesis of secondary metabolite pathways in wheat. In addition, SNPs were also mapped to highly upregulated metal ion transport, AMPK signaling, and ABC transport activity (Table 4). These results suggest that besides variation in expression pattern, sequence level polymorphism in key abiotic stress genes may also have a role in conferring enhanced level of sodicity stress tolerance in KRL 3–4.