To identify TaNHX genes in wheat, we employed a bioinformatics approach. Firstly, we obtained protein sequences of NHX from Arabidopsis, cotton, sorghum, and barley, which served as query sequences. These query sequences were then searched against the wheat genome available at the Ensembl Genomes database (fp://fp.ensemblgenomes.org/pub/plants/release-51/fasta/triticumaestivum/pep/) using the BLASTP algorithm. We selected all homologous protein sequences of NHX candidates that met the specified criteria, including an e-value threshold of 1e-10 and a bit score value higher than 100 percent. To further validate the presence of the Na⁺_H⁺_Exchanger domain (PF00999) characteristic of NHX transporters, the obtained sequences were scanned against this domain using the HMMER 3.1b2 online software (https://www.ebi.ac.uk/Tools/hmmer/) (Potter et al., 2018).

For the subcellular localization, Wolf Psort (https://wolfpsort.hgc.jp/) was used. The physio-chemical properties were predicted using the online programme ProtParam (http://web.expasy.org/protparam/).To study the evolutionary linkage, Clustal W software was used to analyze NHX protein sequences from Arabidopsis, barley, sorghum, cotton, and wheat (Larkin et al., 2007). Following that, an unrooted phylogenetic tree was built in MEGA-X with 1000 bootstrap repetitions using Maximum likelihood (ML) with default parameters the JTT model with uniform rates and 4 number of threads (Kumar et al., 2018). NHX members in wheat were classified into subfamilies based on their Arabidopsis, barley, sorghum, and cotton homologs.

Functional annotation of the identified genes was performed using Blast2GO (https://www.blast2go.com/). This analysis enabled the elucidation of molecular functions, biological processes, and cellular components intricately linked to the genes under investigation. Additionally, we harnessed the KEGG (Kyoto Encyclopaedia of Genes and Genomes) database (http://www.genome.jp/kegg/) to annotate the metabolic pathways, thus furnishing invaluable insights into the precise pathways where these genes actively participate. The combination of Blast2GO and KEGG analysis allowed for a comprehensive understanding of the functional properties and metabolic roles of the identified genes.

To unravel the presence of conserved motifs within the identified genes, we used the MEME webserver (Multiple Em for Motif Elicitation). Employing default parameters, which encompassed the exploration of 2-20 motif sites, the extraction of up to 10 motifs, and a motif width range extending from 6 to 50 (Bailey et al., 2009), this analysis provided an intricate glimpse into the genetic signatures. Furthermore, we delved into the gene structure using the Gene Structure Display Server (GSDS) tool (Hu et al., 2015), unravelling the intricacies of these genes’ architectural features. This comprehensive approach not only unveiled conserved motifs but also shed light on the structural nuances inherent to the genes under study. The promoter sequences, of up to 2 kb in size, were used for cis-acting element analysis and submitted to the PLANTCARE webserver. The cis-acting elements that were produced were categorized according to their functional class.

Based on the position sites, the chromosomal locations of the wheat TaNHX genes were examined. The Ensembl database is used to retrieve information about the TaNHX gene’s locations. Physical mapping of the transporters was constructed using Map-Chart software to construct the location on the wheat genome (Voorrips, 2002).

For the identification of functional protein-protein interactions, the STRING v1054 databases were employed (Szklarczyk et al., 2015). However, Arabidopsis was used as the reference species to search the interactive network in the database. Blast the sequences with set parameters to an e-value of 1e-10, and the A. thaliana genome was searched against all known interaction partners. The best-hit gene for each gene was chosen using Cytoscape 58 to create a PPI network. Using Cytoscape’s cyto Hubba plugin, the top hub gene from the interaction network was finally determined. Furthermore, we conducted an analysis to pinpoint significant Gene Ontology (GO) terms (FDR ≤ 0.01) related to the molecular functions and biological processes of the interaction network nodes using the iDEP webtool (Ge et al., 2018). The transcript sequences of the TaNHX gene were obtained using the Ensembl plant database. Now, the transcript sequences of TaNHX and the mature miRNA sequences from miRbase (Kozomara and Griffiths-Jones, 2014) were analyzed by using the psRNATarget service’s default settings (Dai and Zhao, 2011) to find the targets of miRNAs.

We accessed the coding sequences (CDS) of T. urartu (A genome), A. tauschii (D genome), and T. dicoccoides (AB genome) from Ensembl plants (fp://fp.ensemblgenomes.org/- pub/plants/release-42/fasta/triticum_urartu,fp://fp.ensemblgenomes.org/pub/release-42/plants/fasta/aegilops_tauschii/cds,fp://fp.ensemblgenomes.org/pub/release-42/plants/fasta/triticum_dicoccoides/cds/). These genomes served as references for the A, D, and AB genomes, respectively. We conducted BLASTn searches using these sequences against the TaNHX CDS sequences in T. aestivum to identify orthologous genes. Top hits were designated as orthologs within each species, adhering to stringent criteria, including an e-value cut-off of 1e-10 and a 150-bit-score cut-off. This approach was similarly applied to search for orthologous genes in other monocots and dicots, including O. sativa, H. vulgare, Z. mays, and A. thaliana. To examine the synteny relationships of the NHX gene family within and across various species, we employed TBtools for a comprehensive analysis.

To determine in-silico gene expression in different tissues under salinity stress, the log₂ value of the FPKM value was calculated using the SRA data (SRP304900) from NCBI. Heatmaps were created with Clustvis (https://biit.cs.ut.ee/clustvis/).

To investigate the expression patterns of TaNHX genes in response to salinity stress, we sourced wheat genotypes, specifically KRL213 and HD2009, from the Germplasm Unit of the ICAR-Indian Institute of Wheat and Barley Research, located in Karnal, India. In a controlled environment, we initiated the germination process by placing the seeds in Petri dishes at a temperature of 22°C. Prior to germination, the seeds were sterilization using a 1% sodium hypochlorite solution for duration of 10 minutes, followed by thorough rinsing with distilled water, ensuring the removal of any residual sterilization agents. After five days of germination, seedlings were transferred to full-strength Hoagland’s solution phytojars and incubated for 14 days in a BOD incubator with two sets of three biological replicates of each genotype. For salt stress, two contrasting wheat genotypes, HD2009 (salt sensitive) and KRL213 (salt tolerant), were used. Both genotypes at the two-leaf seedling stage were stressed with 150 mM NaCl. After the treatment, the leaf samples were taken at 0, 3, 24, and 48 hours. For total RNA isolation, all acquired samples were immediately wrapped in foil and frozen in liquid nitrogen at -80°C.

To isolate RNA, we employed TRIzol reagent following the manufacturer’s guidelines. The extracted RNA was further processed to eliminate any residual DNA through DNase I treatment (NEB, USA). Utilizing Superscript-III reverse transcriptase (Invitrogen, USA), 1 µg of total RNA was converted into the first strand of cDNA. This resulting cDNA was subsequently diluted at 1:2 ratio, and 1 µl of the diluted cDNA was employed as a template within a 10 µl reaction volume for real-time qRT-PCR analysis. We conducted SYBR Green-based real-time quantitative RT-PCR analysis using the BIO-RAD CFX96 system (Bio-Rad). For normalization purposes, we employed wheat actin as an endogenous control (Muthusamy et al., 2016). The expression levels were quantified as relative fold changes through the application of the 2^△△-Ct method (Livak and Schmittgen, 2001). All the samples were set up in 3 biological replicates. Experimental data was statistically analyzed by a one-way Analysis of Variance (ANOVA) and Tukey’s multiple range test was used for mean value separation by PAST software (Hammer et al., 2001). P-value of <0.05 was considered statistically significant.

The three-dimensional structure for TaNHX was generated using AlphaFold2 (Skolnick et al., 2021), as it employs deep learning techniques, specifically deep neural networks, to predict protein structures with remarkable accuracy. It leverages vast amounts of protein sequence and structural data to generate predictions of protein structures, even for proteins with no known structures.

The molecular dynamics simulations (MDS) were performed using GROMCAS 5.0 (Abraham et al., 2015). The simulations aimed to investigate the conformational changes of the TaNHX protein structures in the presence of a solvent system. To set up the simulations, the OPLS_2005 force field (Dubay et al., 2012) was employed to describe the interactions within the proteins. The system was solvated in a cubic water box using the SPC water model. The protein atoms were kept at a minimum distance of 10 Å from the edges of the box. To achieve a neutral charge for the systems, counter ions were added, and a 0.15 M ionic concentration was maintained by including Na⁺ and Cl ions. First, each system underwent 50,000 steps of steepest descent energy minimization, which helped eliminate steric overlap and stabilize the initial configurations. Following the energy minimization, a two-step equilibration phase was performed. The first step was NVT equilibration, where the number of particles, volume, and temperature were kept constant. This phase ran for 100 picoseconds (ps) to stabilize the system’s temperature. The V-rescale temperature-coupling method was used for NVT, with a constant coupling time of 1 ps and a target temperature of 303.15 K. The second step was NPT equilibration, where the number of particles, pressure, and temperature were maintained constant. This phase also ran for 100 ps, and it involved using the Nose-Hoover pressure coupling method with a constant coupling time of 1 ps and a target temperature of 303.15 K. During this phase, the system was relaxed, and the protein was restrained using position restraints (h-bonds). To account for electrostatic interactions, the Particle Mesh Ewald method was employed for both NVT and NPT simulations. After the equilibration phases, each system underwent a full production run of 30 nanoseconds (ns) without any restraints. The integration time step used was 0.002 ps, and coordinates were recorded every 10 ps using an xtc collection interval of 5,000 steps.

In order to identify members of the TaNHX gene family within the wheat genome, we employed a two-step approach. Firstly, utilized known NHX sequences from other species as queries and performed a blast search against the wheat proteome database ( Supplementary Table S1 ). This enabled us to retrieve putative TaNHX genes in wheat. Subsequently, we eliminated redundant sequences and confirmed the gene identifications based on the presence of the Na^+_H⁺_Exchanger domain (PF00999). This stringent filtering process resulted in the identification of 30 TaNHX genes in T. aestivum.

To understand the physio-chemical properties of the wheat TaNHX genes, we conducted a comprehensive analysis, the results of which are summarized in Table 1 . Hence, to determine the subcellular localization of NHX proteins, prediction tools have been employed, which were found to predominantly reside on the plasma membrane, endosomes, and vacuoles. Furthermore, we assessed several parameters, including the number of amino acids, which ranged from 374 to 1191, reflecting the diversity in NHX protein lengths. Additionally, the molecular weight of the NHX proteins revealed a wide range of values. For instance, TaNHX16 exhibited a molecular weight of 40.84 kDa, whereas TaNHX11 displayed the highest molecular weight of 131.42 kDa. The findings presented in Table 1 shed light on the distinct physiochemical characteristics of the identified wheat NHX genes.

To analyze the evolutionary lineages of the NHX gene family across different species, a dataset comprising 63 NHX protein sequences was utilized. The sequences were obtained from five species, namely A. thaliana, H. vulgare, G. hirsutum, S. bicolor, and T. aestivum ( Supplementary Table S2 ). These sequences were employed to construct a phylogenetic tree, as depicted in Figure 1 . The analysis involved the utilization of previously reported genes, resulting in the classification of the 63 genes into three distinct categories: class I (endo-class), class II (PM-class), and class III (vac-class). Among the five species, the largest group was class III (vac-class), which encompassed 39 proteins. Specifically, this group comprised 4 NHX proteins in A. thaliana, 2 in H. vulgare, 15 in G. hirsutum, 6 in S. bicolor, and 12 in T. aestivum. It is worth noting that 9 TaNHX proteins were distributed evenly between classes I (endo-class) and II (PM-class), as shown in Figure 1 .

GO enrichment analysis was conducted to uncover the potential functions of the NHX genes. Among the twenty TaNHX genes, significant enrichment was observed in biological processes (BP), molecular functions (MF), and cellular components (CC). In BP, the enriched terms included “sodium ion transmembrane transport,” “response to salt stress,” “potassium ion homeostasis,” and “regulation of pH,” among others ( Supplementary Table S3 ). The enriched MF terms encompassed “binding,” “antiporter activity,” and “potassium: proton antiporter activity.” Furthermore, the enriched CC terms were associated with “plasma membrane” and “vacuolar membrane.” Notably, individual TaNHX genes exhibited enrichment in specific GO terms, indicating their involvement in distinct processes. For comprehensive information on the enriched GO terms and their functions see Supplementary Table S3 .

The MEME web-server with default settings was employed to predict conserved motifs, shedding light on the evolutionary preservation of specific functional amino acids within the NHX protein family. Based on the protein sequences of all NHXs, we discovered a total of 10 potential motifs ( Figure 2A ). The length of the TaNHXs’ anticipated motifs ranged from 6 to 50 amino acids. Members of the same subfamily shared a similar motif arrangement. Logos of the motifs were presented in Supplementary Figure S1 . The descriptive details of domains were given in Supplementary Table S4 .

Using Gene Structure Display Server v.2.0, the intron/exon architecture of the genes was examined to determine the structural properties of the wheat TaNHXs. The number of introns and exons differed significantly depending on the intron/exon patterns studied, which further explains the diversity in gene length. Non-coding sequences are frequently found in the genome, which is thought to be a sign of genomic complexity. Thus, examining these intron configurations reveals important details about the development, regulation and functionality of the NHXs. The examination of the wheat TaNHX gene architectures revealed significant variations across the three classes in terms of the number of introns and exons. The bulk of the Vac-class family among the 30 wheat TaNHXs included UTR (untranslated region) sequences at both the 5’ and 3’ ends. However, 6 out of the 12 Vac-class TaNHXs had 14 exons, 3 had more than 20 exons, and the rest were less than 10. The endo-class TaNHX had more than 10 exons except TaNHX26 ( Figure 2B ). However, the PM-class TaNHX (TaNHX9 and TaNHX10) had the highest share of intron-exons, with 27, 37 exons and 14, 20 introns, respectively. It was found that the distribution of introns and exons in the genes belonging to the same clade was very comparable. The exon lengths and intron regions of genes in the same class were mostly conserved. The analysis of amino acid sequence identity further confirmed the wheat TaNHXs’ sequence conservation ( Figure 2B ).

In order to better understand transcriptional control and gene expression, the promoter sequences of wheat TaNHX were identified by PlantCARE software ( Figure 3 ). The major focus of the analysis was on cis-acting components associated with stress and hormones. Twenty two of the cis-acting components that were discovered had a hormonal connection, including ABA, salicylic acid (SA), gibberellin (GA), and jasmonate. The TaNHX comprised of phytohormones were TaNHX1, TaNHX3, TaNHX5, TaNHX7, TaNHX12, TaNHX13, TaNHX19, TaNHX22, TaNHX23, and TaNHX29 ( Figure 3 ). Eleven cis-acting components, including anaerobic environments and low temperatures, were responsible for stress. The number of light-responsive elements was higher among these components. These cis-acting elements’ findings imply that TaNHXs may be crucial for the control of hormones and the response to stress in wheat.

To better comprehend how the TaNHX gene family evolved, we further examined the gene duplication occurrences. Thirty TaNHX genes were unevenly mapped onto six (A, B, D) chromosomes of the 21 T. aestivum chromosomes ( Figure 4 ). Two chromosomes (Chr2 and Chr5) contained six TaNHX genes and three chromosomes (Chr1, Chr3, and Chr4) contained only one TaNHX gene on each A, B, and D sub-genome ( Figure 4 ). No genes were present on chromosome 6. Gene duplication events, which are a primary mechanism of gene family growth, gave chances for the synthesis of additional genes and their functional divergence.

The response to salt stress is one of the major biological processes in which NHXs play a significant role. The PPI network was created by the STRING database to further examine the potential role of TaNHXs during potential interactions with other proteins ( Figure 5 ). TaNHX proteins are not expected to have any direct interactions with one another. Out of thirty TaNHX proteins, six were participating in the protein-protein interaction network. They shared a similar kind of putative interactive protein. TaNHX1, TaNHX2, TaNHX7, and TaNHX20 were involved in proton exchange and belong to monovalent and proton antiporters ( Figure 5 ). Numerous physiological processes, including vesicle trafficking, pH control, K⁺ homeostasis, protein transport, and growth and development, are regulated by the antiporters. In contrast, TaNHX6 interacted with SOS1, which is crucial for the expulsion of Na⁺ ions from cells.

This study employed network analysis, using TaNHX gene expression data, to uncover key insights into sodium ion transport mechanisms in plants. Significant Gene Ontology (GO) terms (FDR < 0.001) were identified, revealing TaNHX genes’ involvement in cellular processes, particularly transmembrane transport. Intriguingly, our analysis revealed intricate interactions between sodium ion transport and the regulation of biological quality, chemical homeostasis, and metal ion homeostasis ( Supplementary Figure S2 ). Subcellular localization indicated that TaNHX proteins primarily reside in plasma membranes, endosomes, and vacuoles, emphasizing their role in Na⁺/-H antiporter activity and ion homeostasis. The findings also imply potential gene interactions affecting activity levels, paving the way for further investigations into these regulatory networks in plant sodium ion transport ( Supplementary Figure S2 ). This research contributes valuable insights to plant biology and ion homeostasis regulation.

The psRNATarget service has been used to look at how miRNAs regulate the expression of TaNHX genes. For twenty distinct miRNAs, we identified 32 TaNHX genes as potential targets ( Table 2 ). Tae-miR395 is implicated in the regulation of eight TaNHX genes (TaNHX6, TaNHX9, TaNHX10, TaNHX11, TaNHX12, TaNHX24, TaNHX25, and TaNHX28). Tae-miR9666 accounted for regulating the expression of four TaNHX genes (TaNHX6, TaNHX9, TaNHX24, and TaNHX25). The expression of TaNHX genes (TaNHX1, TaNHX3, and TaNHX29) may be influenced by Tae-miR9656 ( Table 2 ). Tae-miR9772, Tae-miR9773, and Tae-miR9774 were predicted to regulate the expression of TaNHX22, TaNHX13, and TaNHX9, respectively.

Comparative analysis was utilized to evaluate the orthologous of TaNHX gene duplication in wheat and A. tauschii. We discovered 28 orthologous gene pairs among all the NHX genes from wheat and A. tauschii. Contrarily, there were 27 NHX orthologous gene pairs discovered between wheat and rice, 29 between wheat and T. dicoccoides, and 30 between wheat and Arabidopsis ( Supplementary Table S4 ). Orthologous associations between several species can be found through the analysis of collinearity correlations. NHX gene pairings between the genomes of T. aestivum and A. thaliana were syntenized. The findings revealed that 12 T. aestivum NHX genes shared syntenic relationships with AtNHX genes ( Figure 6 ), indicating that these genes may have aided in the development of the TaNHX gene family.

We conducted an analysis of synonymous substitutions (K_s) and non-synonymous substitutions (K_a) values to investigate the selective pressures underlying the duplication events of TaNHX genes, considering all nucleotide sequences. Our findings unveiled intriguing patterns, notably identifying six gene pairs (TaNHX7/TaNHX5, TaNHX15/TaNHX14, TaNHX22/TaNHX27, TaNHX1/TaNHX3, TaNHX20/TaNHX18, and TaNHX26/TaNHX23) characterized by K_a/K_s ratios lower than 1. These results suggest that these gene pairs experienced purifying selection, indicative of their evolutionary conservation ( Supplementary Table 5 ). Our analysis detected nine distinct segmental duplication events occurring across various chromosomes, along with one notable tandem duplication event within the same chromosome. These findings underscore the pivotal role of segmental duplications in driving the expansion of TaNHX genes within the wheat genome. Furthermore, our observations suggest that certain TaNHX genes may have originated as a consequence of gene duplication events, shedding light on the evolutionary mechanisms responsible for the diversification of this gene family. But it may be due to selection pressure that the number of TaNHX genes was not expended in comparison to other transcription factors in wheat. Genes from three subfamilies (PM, endo and vac.) participated in the tandem and segmental duplication. We also investigated how frequently tandem duplications occur. This region contained 10 TaNHX gene pairs, all of which were closely linked. Notably, the sequence identities among these genes exceeded 80%, strongly suggesting their involvement in tandem duplication events. Given the profound impact of gene duplication on the emergence of novel functionalities and gene families, we conducted a comprehensive exploration of TaNHX gene duplication events within the wheat genome. The paralogous gene pairs were employed to construct a circos plot ( Figure 7 ), enabling us to visualize and elucidate the intricate patterns of duplication events and their implications.

TaNHX expression levels were analyzed in wheat roots and leaf tissues exposed to salt stress. The expression analysis was done using RNA-seq data (SRP304900) from NCBI. Under normal conditions, TaNHX1, TaNHX2, TaNHX3, TaNHX5, TaNHX10, TaNHX11, and TaNHX20 were found to be upregulated in all tissues except in grain ( Figure 8 ). However, TaNHXs displayed differential gene expression levels among the tissue types in different genes. It was observed that 23 out of 30 TaNHX genes were expressed very highly in both tissues. Whereas, the genes e.g., TaNHX4, TaNHX17, TaNHX21, TaNHX22, TaNHX27, TaNHX28, and TaNHX29 expression were downregulated in all tissues under salt stress conditions ( Figure 8 ). However, it was found that TaNHX3, TaNHX5, TaNHX7, TaNHX11, TaNHX16, TaNHX18, and TaNHX20 expressed upregulation in all tissues under heat, drought, and cold ( Supplementary Figure S3 ) stresses.

Using qRT-PCR analysis, we quantitatively evaluated the expression patterns of ten TaNHXs genes under control and salt-stress conditions to observe the potential functions that TaNHXs may play in wheat in responses to salinity stress. The findings demonstrated that salt tolerance (KRL213) and salt-sensitive (HD2009) wheat cultivars both stimulated all TaNHX genes in leaf and root tissues in a concentration-dependent manner ( Figure 9 ). The tissues with the highest induction were the roots, followed by the leaves. TaNHX gene overexpression caused by salinity was greater in KRL213 than HD2009. Out of ten TaNHX genes, five genes [TaNHX2 (~14 folds), TaNHX12 (~9 folds), TaNHX16 (~7 folds), TaNHX20 (~10 folds), TaNHX23 (~3 folds)] showed up-regulation, in the roots of KRL213 at 48 h under 150 mM NaCl stress condition ( Figure 9A ). On the other hand, expression levels of TaNHX2, TaNHX16, and TaNHX20 were expressed more than twofold in roots treated with 150 mM NaCl at different time intervals. It was observed that TaNHX2, TaNHX12, and TaNHX20 expression levels were higher in the leaves of KRL213 (~10, ~5, and ~16 times, respectively) than HD2009 were down-regulated at 48 h of 150 mM NaCl treatment ( Figure 9B ). Salinity stress clearly shows up-regulation of TaNHXs in the leaves of KRL213, whereas in HD2009 two genes TaNHX2 and TaNHX20, were expressed higher in leaf tissue at 0 h and 3 h of stress. In KRL213 leaf tissue, all genes used for the validation were up-regulated at 3 h of salt stress ( Figure 9B ). The gene expression profiles of individual TaNHX members were evaluated in both leaf and root tissues under salinity stress conditions, yielding results that consistently mirrored the salt tolerance characteristics of the two wheat genotypes under study.

The Root Mean Square Deviation (RMSD) serves a crucial statistic for assessing the average changes of a set of atoms concerning a reference frame. Analyzing the RMSD of the protein backbone atoms relative to their initial positions can yield valuable insights into the protein conformational dynamics. In our investigation, the RMSD values for all three proteins were below 2.5 nm. Figure 10A illustrates a stabilizing trend in the RMSD values around 10 ns mark, suggesting the attainment of equilibrium within the simulation. This stabilization milestone serves as an optimal foundation for subsequent analysis.

Moreover, the Root Mean Square Fluctuation (RMSF) analysis delves into the localized fluctuations occurring along the protein chain during the molecular dynamics (MD) simulations RMSF graph uncovered fluctuations primarily at the C-terminals of the proteins, while the core of the TaNHX proteins displayed no substantial variations ( Figure 10B ). Additionally, the radius of gyration (Rg) value for the TaNHX structure exhibited a decreasing trend as the simulation time progressed, indicating a progressive compaction of the modeled structure. This observation aligns with the results from the initial run of the simulations. The radius of gyration is a measure of the overall size and compactness of a protein structure. A lower Rg value suggests a more tightly folded and compact conformation, while a higher Rg value indicates a more extended or flexible structure. Therefore, the lower Rg values obtained for the TaNHX protein structure suggest that they adopt a more compact and stable conformation ( Figure 10C ).