Investigating NAC Transcription Factor Role in Redox Homeostasis in Solanum lycopersicum L.: Bioinformatics, Physiological and Expression Analysis under Drought Stress

NAC transcription factors regulate stress-defence pathways and developmental processes in crop plants. However, their detailed functional characterization in tomatoes needs to be investigated comprehensively. In the present study, tomato hybrids subjected to 60 and 80 days of drought stress conditions showed a significant increase in membrane damage and reduced relative water, chlorophyll and proline content. However, hybrids viz., VRTH-16-3 and VRTH-17-68 showed superior growth under drought stress, as they were marked with low electrolytic leakage, enhanced relative water content, proline content and an enhanced activity of enzymatic antioxidants, along with the upregulation of NAC and other stress-defence pathway genes. Candidate gene(s) exhibiting maximum expression in all the hybrids under drought stress were subjected to detailed in silico characterization to provide significant insight into its structural and functional classification. The homology modelling and superimposition analysis of predicted tomato NAC protein showed that similar amino acid residues were involved in forming the conserved WKAT domain. DNA docking discovered that the SlNAC1 protein becomes activated and exerts a stress-defence response after the possible interaction of conserved DNA elements using Pro72, Asn73, Trp81, Lys82, Ala83, Thr84, Gly85, Thr86 and Asp87 residues. A protein–protein interaction analysis identified ten functional partners involved in the induction of stress-defence tolerance.


Introduction
In the 21st century, agriculture, crop growth and productivity are being redundantly affected by reiterative environmental pollution, prompting severe threats to global food security for ever-growing populations [1]. In Asia, where most farmers routinely practise rainfed agriculture, drought stress has become a daunting challenge evoking crop failure, yield losses and livestock death [2]. Plants activate many cellular, physiological, biochemical and molecular responses to confront drought stress-induced oxidative damage by regulating the expression of several stress-responsive genes, such as proline/glycine betaine, thus, improving photosynthesis and water-use efficiency [3]. The coordinated regulations of defence-related genes and the accumulation of osmolytes under drought stress are provoked by the expression of functional regulatory proteins that instigate an extensive reprogramming of transcriptome and metabolome. These functional regulatory proteins, i.e., transcription factors (TFs), play a central role in regulating the complex gene regulatory network of signal transduction and stress-defence pathways [4].
No Apical Meristem (NAM), Arabidopsis Transcription Activation Factor (ATAF), and Cup-Shaped Cotyledon (CUC) transcription factors, collectively called NAC TFs, are one of the most prominent families of TFs, which play a significant role in the transcriptional Further, SlNAC1 was distinguished from other NAC members by forming an individual clade with a total bootstrap value of 94%, thus, indicating the recent expansion of SlNAC1 from other members of the solanaceous family. Inspection of the functional domain of SlNAC1 protein using a multiple sequence alignment analysis revealed maximum conservation of core amino acid residues around the WKATGAD region ( Figure 2). Multiple sequence alignment analysis confirmed that the slightest discrepancies occurred in the WKAT motif during evolution and might exert a functional role in the interaction of NAC proteins with cis-regulatory elements. Further, SlNAC1 was distinguished from other NAC members by forming an individual clade with a total bootstrap value of 94%, thus, indicating the recent expansion of SlNAC1 from other members of the solanaceous family. Inspection of the functional domain of SlNAC1 protein using a multiple sequence alignment analysis revealed maximum conservation of core amino acid residues around the WKATGAD region ( Figure 2).

Gene Structure and Conserved Motifs Analysis
Plants 2022, 11, x FOR PEER REVIEW 3 of 23 using MEGA 7.0 employing standard methods. The tree topologies generated by both algorithms revealed the monophyletic origin of SlNAC1 with Solanum pennelli and its close resemblance with Solanum tuberosum ( Figure 1A,B). Further, SlNAC1 was distinguished from other NAC members by forming an individual clade with a total bootstrap value of 94%, thus, indicating the recent expansion of SlNAC1 from other members of the solanaceous family. Inspection of the functional domain of SlNAC1 protein using a multiple sequence alignment analysis revealed maximum conservation of core amino acid residues around the WKATGAD region ( Figure 2). Multiple sequence alignment analysis confirmed that the slightest discrepancies occurred in the WKAT motif during evolution and might exert a functional role in the interaction of NAC proteins with cis-regulatory elements. Multiple sequence alignment analysis confirmed that the slightest discrepancies occurred in the WKAT motif during evolution and might exert a functional role in the interaction of NAC proteins with cisregulatory elements.

Gene Structure and Conserved Motifs Analysis
The conserved motifs of NAC proteins in tomatoes were analysed by scanning the SlNAC1 proteins for the presence of putative motifs through the MEME analysis tool. Twenty-one motifs showed vital conservation around the WKAT motif across all the members of the NAC protein family. The result also illustrated that the phylogenetically related members had similar motif distribution patterns, thus, demonstrating structural and functional conservation among NAC proteins ( Figure 3A). SlNAC1 proteins for the presence of putative motifs through the MEME analysis tool. Twenty-one motifs showed vital conservation around the WKAT motif across all the members of the NAC protein family. The result also illustrated that the phylogenetically related members had similar motif distribution patterns, thus, demonstrating structural and functional conservation among NAC proteins. ( Figure 3A). The reproducibility and reliability of the identified motifs were confirmed using the E-value that corresponds to the frequency of the occurrences of the motif and the p-value which calculates the likelihood of their existence, and in the present study, the N-terminal of SlNAC1 protein ( Figure 3A) was denoted by motif 1 (WYFFSPRDRKYPNGSRP-NRAAGTGYWKATGADKPVGKPKTLGIKKALVFY; p-value 2.1e-64) The multiple sequence alignment results of NAC proteins corroborate the motif analysis result substantiating the conservation of core residues, as reflected in our motif distribution analysis. Interestingly, the presence of some specific motifs was also observed in other NAC protein family members, for instance, in N tabacum motif 2 ( Figure 3A(i)) (DEVEEE; p-value 1.3e-7) and motif 3 ( Figure 3A(ii)) (EEKPKV; p-value 4.7e-7), whereas in A thaliana motif 4 (Figure 3A(iii)) (QSNELI; p-value 2.1e-7) and motif 5 ( Figure 3A(iv)) (QLAPAP ; p-value 2.3e-6) they were observed at the C-terminal end ( Figure 3A). These additional motifs within the same group may be due to sequence divergence during the evolutionary processes. The logo diagram of motif 1 representing the WKAT motif showing maximum conservation among all the solanaceous NAC members is depicted in Figure 3B.

Prediction of Conserved Functional Domains
Imposition of SlNAC1 protein to ExPASy PROSITE and InterProScan online tool (http://prosite.expasy.org) (accessed on 2 January 2022) revealed the presence of available signature sequence at the N-terminal end constituting the NAC domain region, thereby confirming its belonging to the NAC gene superfamily. Additionally, an analysis of SlNAC1 protein through InterProScan also clearly indicated that the putative SlNAC1 protein belongs to the NAC superfamily containing DNA-binding domain ranging from 1 to 301 amino acid residues. The results of the PROSITE and InterProScan analysis were The reproducibility and reliability of the identified motifs were confirmed using the E-value that corresponds to the frequency of the occurrences of the motif and the p-value which calculates the likelihood of their existence, and in the present study, the N-terminal of SlNAC1 protein ( Figure 3A) was denoted by motif 1 (WYFFSPRDRKYP-NGSRPNRAAGTGYWKATGADKPVGKPKTLGIKKALVFY; p-value 2.1e-64) The multiple sequence alignment results of NAC proteins corroborate the motif analysis result substantiating the conservation of core residues, as reflected in our motif distribution analysis. Interestingly, the presence of some specific motifs was also observed in other NAC protein family members, for instance, in N tabacum motif 2 ( Figure 3A(i)) (DEVEEE; p-value 1.3e-7) and motif 3 ( Figure 3A(ii)) (EEKPKV; p-value 4.7e-7), whereas in A thaliana motif 4 ( Figure 3A(iii)) (QSNELI; p-value 2.1e-7) and motif 5 ( Figure 3A(iv)) (QLAPAP; p-value 2.3e-6) they were observed at the C-terminal end ( Figure 3A). These additional motifs within the same group may be due to sequence divergence during the evolutionary processes. The logo diagram of motif 1 representing the WKAT motif showing maximum conservation among all the solanaceous NAC members is depicted in Figure 3B.

Prediction of Conserved Functional Domains
Imposition of SlNAC1 protein to ExPASy PROSITE and InterProScan online tool (http://prosite.expasy.org) (assessed on 15 July 2022) revealed the presence of available signature sequence at the N-terminal end constituting the NAC domain region, thereby confirming its belonging to the NAC gene superfamily. Additionally, an analysis of SlNAC1 protein through InterProScan also clearly indicated that the putative SlNAC1 protein belongs to the NAC superfamily containing DNA-binding domain ranging from 1 to 301 amino acid residues. The results of the PROSITE and InterProScan analysis were also corroborated by the multiple sequence analysis performed for core signature sequences revealing vital conservation across all the members (Supplementary Materials, Figure S1).

Modelling and Superimposition of the NAC Domain
In the present study, a total of five models (SlNAC1) were generated and the model (Figure 4) exhibiting the lowest RMSD (Cα trace) and DOPE score compared to the respective templates was used for the interaction study ( Table 1). The modelled SlNAC1 protein displaying the best score values was submitted to PMDB (PM0082340) (Supplementary Materials, Figure S2). Often, proteins of divergent evolution and distinct sequences exhibit structural resemblance;  Figure S3A(I-III)). The distinguishing topology for the modelled protein by the 3ULX template may be due to substitutions of some amino acid residues during evolution.

Modelling and Superimposition of the NAC Domain
In the present study, a total of five models (SlNAC1) were generated and the model (Figure 4) exhibiting the lowest RMSD (Cα trace) and DOPE score compared to the respective templates was used for the interaction study ( Table 1). The modelled SlNAC1 protein displaying the best score values was submitted to PMDB (PM0082340) (Supplementary Materials, Figure S2). Often, proteins of divergent evolution and distinct sequences exhibit structural resemblance; therefore, to test this hypothesis, the predicted SlNAC1 protein was superimposed on the respective templates using the SALIGN web server. The superimposition results indicated a high similarity of the modelled SlNAC1 protein with 1 UT7, 4 DUL and 3 ULX, as confirmed by their RMSD values viz., 1 Figure S3A I-III). The distinguishing topology for the modelled protein by the 3ULX template may be due to substitutions of some amino acid residues during evolution.

Statistical Validation of the SlNAC1 Model
The statistical validation of the predicted SlNAC1 model protein generated through the DS Modeller was performed in terms of probability density functions (PDFs) and discrete optimized protein energy (DOPE) scores. The geometry of the resultant SlNAC1 model was assessed using RAMPAGE and PROCHECK, which showed that 98% of the amino acid residues were observed in the most favoured regions (Table 2). Further, a PROCHECK analysis also predicted that 91.4% (expected 90%) of the residues occurred in the active region (Supplementary Materials, Figure S3BI,II). The ProSA analysis verified

Statistical Validation of the SlNAC1 Model
The statistical validation of the predicted SlNAC1 model protein generated through the DS Modeller was performed in terms of probability density functions (PDFs) and discrete optimized protein energy (DOPE) scores. The geometry of the resultant SlNAC1 model was assessed using RAMPAGE and PROCHECK, which showed that 98% of the amino acid residues were observed in the most favoured regions (Table 2). Further, a PROCHECK analysis also predicted that 91.4% (expected 90%) of the residues occurred in the active region (Supplementary Materials, Figure S3B(I,II)). The ProSA analysis verified the superiority of the predicted SlNAC1 model over individual templates in terms of minimum Z score, thus, indicating that the expected model was very close to the template structures (Table 2; Supplementary Materials, Figure S3B(I,II)). In addition, the ERRAT analysis of the predicted model displayed a maximum score (60.17) compared to the respective templates, affirming superior atomic interactions in the resultant model ( Table 2). The quality of SlNAC1 was also verified by a Resolution by Proxy (RESPROX) and Qualitative Model Energy Analysis (Q-MEAN), which significantly attested that the modelled protein had a better atomic resolution that allows more feasible long-range interaction in the expected model compared to the respective templates ( Table 2). The quantitative validation of the predicted SlNAC1 model showed helical type structure 76 (50%), beta 58 (38%) and 18  A good quality protein has ReSProx values ranging between 0 and 1.5, ERRAT score of more than 90, minimum Z and Q-mean scores occupying more than 98% of the residues in the most favoured region, less than 2% residues in additionally allowed region and 0% residues in outlier region in RAMPAGE analysis.

DNA-Protein Docking
The molecular docking analysis was performed through the HexDoc server using the most reliable SlNAC1 model ( Figure 5A(I)), sufficing all the crucial energy parameters with the 3D-modelled DNA element ( Figure 5A(II)). The result indicated that the most stable and favourable docked complex has the binding energy (E total = −1416.42 Kcal/mol; Supplementary Materials, Figure S4), which confirms that the NAC proteins are activated upon binding with the specific DNA elements through conserved WKAT motifs ( Figure 5A(I)). The docking results further identified the key amino acid residues involved in the interaction as Arg 63 , Pro 73 , Asn 74 , Trp 82 , Lys 83 , Ala 84 , Thr 85 , Gly 86 , Ala 87 and Asp 88 ( Figure 5A(I,II)), as compared to the template structure (1UT7), in which the interacting residues were Pro 71 , Asn 72 , Trp 81 , Lys 82 , Ala 83 , Thr 84 , Gly 85 , Ala 86 and Asp 87 ( Figure 5A), thus, indicating the vital condensed nature of the core amino acid residues of WKAT motifs.

Protein-Protein Interaction
The functional protein-protein interaction was assessed using the STRING server at a medium confidence level from the values of the shell of interactors ( Figure 5B). At a medium confidence score, SlNAC1 was found to have a strong interaction with Solyc04g009440.2.1 (score value 0.780) protein, which is an uncharacterized protein with a possible role in regulating NAC protein expression. The SlNAC1 protein also interacted with the ethyleneresponsive transcription factor 1 protein, showing its significant role in regulating gene expression during fruit ripening and abiotic stresses. It is also involved in the modulation of stress signal transduction pathways with a score value of 0.777 from the first shell on interactors. Further, the SlNAC1 protein also showed interaction with the AP2 transcription factor family (101256610; score value 0.780), which is involved in the regulation of a wide array of plant physiological and molecular responses under various abiotic stresses, and the p450 protein involved in multiple biosynthetic and xenobiotic pathways under stress  Table S1).

Protein-Protein Interaction
The functional protein-protein interaction was assessed using the STRING server at a medium confidence level from the values of the shell of interactors ( Figure 5B). At a medium confidence score, SlNAC1 was found to have a strong interaction with Solyc04g009440.2.1 (score value 0.780) protein, which is an uncharacterized protein with a possible role in regulating NAC protein expression. The SlNAC1 protein also interacted with the ethylene-responsive transcription factor 1 protein, showing its significant role in regulating gene expression during fruit ripening and abiotic stresses. It is also involved in the modulation of stress signal transduction pathways with a score value of 0.777 from the first shell on interactors. Further, the SlNAC1 protein also showed interaction with the AP2 transcription factor family (101256610; score value 0.780), which is involved in the regulation of a wide array of plant physiological and molecular responses under various abiotic stresses, and the p450 protein involved in multiple biosynthetic and xenobiotic pathways under stress conditions (Solyc06g076160.2.1; score value 0.776). Apart from these, the SlNAC1 protein also interacted with three uncharacterized proteins viz.,  Table S2).

Gene Ontology Analysis and Subcellular Localization
To assess the biological, molecular function and subcellular localization, the predicted SlNAC1 protein was analysed by the CATH and Gene3D server, which uses a clus-

Gene Ontology Analysis and Subcellular Localization
To assess the biological, molecular function and subcellular localization, the predicted SlNAC1 protein was analysed by the CATH and Gene3D server, which uses a clustering algorithm to identify gene ontology (GO) terms. The GO terms under biological function were regulation of transcription (GO: 0006355); regulation of hyperosmotic and salinity response (GO:0042538); response to salicylic acid and water deprivation (GO: 0009751; GO: 0009414); and proline biosynthesis (GO: 0006561) under stress condition (Supplementary Materials, Figure S5A(I)), whereas the GO terms associated with molecular functions (Supplementary Materials, Figure S5A(II)) were DNA binding (GO: 0003677) and transcription factor activity (GO: 0003700). The Cello2GO results predicted that the SlNAC1 protein is mainly localized in the nucleus (37.4%; score value 1.869) and mainly intricated in regulating DNA-binding and transcription factor activity, whereas 22.5% of SlNAC1 proteins are located in the cytoplasmic region (score value 1.124) involved in various biosynthetic and metabolic processes (Supplementary Materials, Figure S5B). On the other hand, a cellular component analysis revealed that the predicted proteins have a wide range of functions, as inferred by the results of protein-protein interaction study.

Expression Analysis of Defence-Related Genes
To gain functional insight into the mechanism that induces drought stress tolerance in tomato hybrids, we examine the expression pattern of essential genes related to the stress-defence pathway. The genes viz., NAC 1, NAC1.1, NAC 2, heat shock protein (HSP), late embryogenesis abundant protein (LEA), zinc finger protein (ZFP) and elongation factor protein (EF) were selected for their possible role in the regulation of transcription, protein folding and activation of signal transduction pathways under various abiotic stresses ( Figure 7B). The qRT analysis of the above genes was performed among all the tomato hybrids/parents at 80 days of water deficit condition. In the present study, the relative abundance of transcript encoding NAC1 (45.2-66.6%), NAC1.1 (39.9-54.6%) and NAC 2 (48.8-74.3%) was significantly increased more prominently in VRTH-16-3 and VRTH 17-68 compared to VRTH-17-2 and VRTH-17-81, as well as respective parents at 80 days of irrigation deficit condition. Similarly, the fold change of HSP (38.7-51.4%) and ZFP (22.4-30.9%) was increased in VRTH-16-3 and VRTH 17-68 and down-regulated in VRTH-17-2 and VRTH-17-81, as well as respective parents at 80 days of irrigation deficit condition ( Figure 7B). Furthermore, the expression of ZFP and EF also showed a differential response in all the hybrids, where maximum upregulation was observed in VRTH-16-3 and VRTH 17-68. However, both LEA (20.3-28.9%) and EF (18.4-25.4%) were also significantly induced in VRTH-17-2 and VRTH-17-81 hybrids, as well as in tolerant parent Punjab Kesari. Altogether, the result indicated that the expression of the critical genes of the signal transduction pathway was more prominent in VRTH-16-3 and VRTH 17-68 hybrids compared to VRTH-17-2 and VRTH-17-81, which could be the possible reason behind the improved physiological growth, increased productivity and enhanced antioxidant capacity of VRTH 16-3 and VRTH 17-68 hybrids, thus, labelling them as drought-tolerant hybrids compared to their respective counterparts.

Discussion
In the present study, we tried to characterize SlNAC1 TFs using in silico approaches functionally. Then, we validated their role through a gene expression analysis in regulating physiological, biochemical and molecular drought stress responses in four tomato hybrids, along with respective drought-tolerant/susceptible parents. Initially, the presence of the NAC domain in the SlNAC1 protein was envisaged and established by a PROSITE and InterProScan analysis. The target sequence was then subjected to evolutionary computation, which revealed a close phylogenetic relationship among orthologs, monophyletic origin with S. pennelli and close resemblance with S. tuberosum NAC proteins. However, other members of the tomato family, such as S. melongena, C. annum and N. tabacum, form a separate cluster.
Interestingly, when SlNAC1 was compared with A. thaliana, the phylogenetic tree showed three distinct clusters showing S. lycopersicum, S. pennelli and S. tuberosum in one cluster and S. melongena, C. annum and N. tabacum in another cluster, whereas A. thaliana was out-clustered from the group (Figure 1A,B). The high resemblance among the NAC genes of tomato and potato might be due to their common ancestral origin, thus, sharing 97.1-100% homology and identity [13,14]. The results were also corroborated by multiple sequence alignment (Figure 2) done in the present study, which further confirms that the core residues around the WKAT motif in S. lycopersicum, S. pennelli and S. tuberosum are highly conserved compared to those of S. melongena, C. annum and N. tabacum, thus, indicating that the slightest disturbance occurred during the evolution [15,16].
The DNA-protein interaction ability of any transcription factors relies on the presence of functionally conserved motifs, specifically at the N-terminal domain, that are displayed as core residues capable of regulating various signalling pathways other than DNA-protein interaction [17]. Twenty-one motifs were predicted using MEME suit ( Figure 3A), of which five motifs can be classified as A-E type, as they marked the presence of a complete DNA-binding domain and C-terminal activation domain [18]. In contrast, the remaining 16 motifs do not contain an entire DNA-binding domain; however, they exhibit location matching with the above five motifs [19]. We also identified two conserved motifs in all the members of the NAC family. First, the WKATGAD motif ( Figure 3B) at the N-terminal domain is biochemically necessary for DNA-protein interaction and nuclear shuttling [20]. Furthermore, it has been reported that WKAT motifs can also modulate protein-protein interaction, thus, improving plant growth and development against pathogen attack and abiotic stresses [21]. Second, the LVFY motif ( Figure 3B) at the C-terminal domain is essential for the activation/expression of transcription and protein-protein interaction [20].
The motif and phylogenetic analysis revealed that the identified genes are functional homologs of the SlNAC1 protein, thereby confirming the close evolutionary origin of the identified proteins. Several researchers have reported that genes with diverse, active roles might evolve from the same evolutionary process [22][23][24][25]. In recent years, evidence has suggested that structurally different homologs of NAC TFs play a crucial role in improving plant defence response when expressed in a tissue-specific manner [5].
The present study predicted the 3D structure of the SlNAC1 protein by choosing template proteins with the maximum query coverage and sequence identity using BLASTP search against the RCSB protein data bank (http://www.rcsb.org/) (assessed on 15 July 2022) [26]. Three template proteins viz., Template 1-PDB ID: 1UT7 (sequence identity 70%; query coverage 52%); Template 2-PDB ID: 3ULX (sequence identity 74.5%; query coverage 56%); and Template 3-PDB ID: 4DUL (sequence identity 70%; query coverage 52%) were used for modelling SlNAC1 protein. A large body of literature has indicated that the protein model generated using homology modelling usually suffers from poor bond angles, inappropriate bond lengths and unfavourable angle contacts [27]. A good protein model should have minimum energy to maintain appropriate bond angle geometry and, thus, normalize the close angle contacts [27]. In the present study, out of five models generated, the model ( Figure 4) with a low DOPE score (−14,933.81; Table 1) was selected for interaction studies, as a model with such a score has minimum steric collision and strain. Therefore, they are much more stable than models with high DOPE scores [28]. The superimposition of modelled proteins with template proteins is also an alternate way to predict their relative accuracy and sequence homology [29]. Consequently, the superpose program was used for the 3D alignment of the predicted SlNAC1 model with the respective templates. The predicted model showed high structural ( Figure 4) and sequence ( Figure 4) similarity, as indicated by the root mean square deviation (RMSD) values of Cα superimposition (1.74, 2.20 and 2.97), suggesting that the predicted 3D structure of SlNAC1 protein is similar to the templates (Supplementary Materials, Figure S3A) [30].
The SlNAC1 model was then subjected to a quantitative analysis to verify/validate its topology and geometry. The qualitative estimation was performed using RAMPAGE. The PDBsum server indicated that about 98% of the residue lies in the favoured region and 2% in the outlier region from the RAMPAGE server (Table 2). In contrast, from the PDBsum server, 91.4% were in the favoured area, and no residues were present in the disallowed region (Supplementary Materials, Figure S3B(I,II)). Generally, a model is considered reliable if 90% of the residues fall in the favoured area [31]; hence, this putative model can be exploited for future studies in tomatoes. The position of the core residues in the predicted model was determined using a Z-score energy plot analysis of the ProSA webserver (Supplementary Materials, Figure S3B(I,II)). The Z-score analysis results further confirmed the expected model's stability and reliability by measuring the total energy deviated from the total points distributed from random confirmations [32].
To qualify as a satisfactory model, the predicted models should have an ERRAT score of at least 50 and a minimum Q-mean score. In contrast, the ERRAT score of the expected model was high (60.1) compared to the respective templates and the Q-mean score was −1.87, indicating a minimum error in the structure ( Table 2). The ResProx program was also used to predict the atomic resolution of the forecast model, and the expected model exhibited a value of 1.98 compared to the respective templates (Table 2). Usually, the nuclear resolution is sorted out based on the correlation coefficient values between the observed and calculated atomic coordinates of the nuclear model, and a model with a minimum ResProx score (<2.0 Å) is considered a good model [33]. The analysis mentioned above evidently verifies and validates the reliability of the modelled SlNAC1 protein.
In the present study, the DNA-protein docking analysis of the predicted SlNAC1 protein with the DNA element (CACATG) through the HexDoc server indicated that the amino acid residues viz., Arg 63 , Pro 73 , Asn 74 , Trp 82 , Lys 83 , Ala 84 , Thr 85 , Gly 86 , Ala 87 and Asp 88 were exclusively involved in the interaction with the DNA element ( Figure 5A(I,II)). Further, this server reveals the presence of three-way binding sites in the predicted SlNAC1 protein, of which the last binding site was found to be more prominent ( Figure 5A). To validate the effectiveness of the DNA protein interaction of the predicted model, the template protein (1UT7) was docked with the DNA element as a reference [34]. The results of the DNA-protein interaction indicated close resemblance in the interacting residues of both modelled SlNAC1 and the reference template (1UT7), thus, confirming that the predicted SlNAC1 protein ( Figure 5A(II)) was modelled with the right folds and domains [35,36]. Therefore, the expected SlNAC1 model structure was effectively submitted to the PMDB database with PMDB ID PM0082340.
The string analysis reveals ten functional partners that showed a direct interaction with the predicted SlNAC1 protein, where the highest interaction was observed with uncharacterized protein LOC544041 with the identity of 84.1%, which is believed to be a member of the AP2 transcription factor family involved in the regulation of growth and development in plants exposed to climate extremes [37]. Further, the predicted SlNAC1 protein also interacted with ERF1 (Solyc03g093610.1) protein with an interaction score of 0.777, which is involved in the modulation of the expression of several defence-responsive genes, transcriptional activator and the regulation of diseases resistance pathway [38]. Other functional protein partners observed in the analysis scored between 0.551 and 0.699 with unknown functions (Supplementary Materials, Table S1), which could be involved in several other metabolic pathways, as suggested by their KEGG annotation. The present study's interaction results pinpointed the differential regulatory pattern of the SlNAC1 protein, which could be responsible for imparting resistance in plants under drought and other abiotic stress conditions. Therefore, the functional characterization of SlNAC1 protein and its interaction partners can further open a new realm for gaining valuable insights into this protein's molecular and regulatory function in other crop plants.
The gene ontology (GO) analysis was also used to categorize the SlNAC1 protein functionally based on its cellular, biological and molecular function [39]. Following the previous reports [40], the GO terms under biological processes (Supplementary Materials, Figure S5A(I,II)) were the genes that were involved in salicylic acid signalling, proline biosynthesis and salinity stress response, which may have been upregulated after interacting with NAC proteins and henceforth might have resulted in the enhanced tolerance of tomato hybrids under drought stress condition [33]. The GO terms governing molecular function (Supplementary Materials, Figure S5B) were the genes involved in DNA binding and transcription activation that stimulated the expression of NAC and associated proteins, resulting in the increased tolerance of tomato hybrids under drought stress [33]. Furthermore, we also assessed the relative subcellular localization of SlNAC1 protein through the Cello2GO web server, which indicated that 37.4% of the protein is of nuclear origin and mainly involved in DNA binding and transcriptional activation activity, and 22.5% of the protein to cytoplasm imparting biosynthetic function, as suggested by the results of the gene ontology analysis.
Drought stress instigates severe physiological and biochemical changes in plants, which differentially vary in stress intensity, agronomic factors and plant developmental stages [41]. In this study, we observed the differential tolerance response of tomato hybrids to drought stress. Higher relative water content and lower electrolytic leakage were observed in VRTH-16-3 and VRTH-17-68 hybrids compared to their corresponding counterparts/respective parents, demonstrating the competency of these hybrids to withhold water loss (Tables 3 and 4). Leaf RWC and EL are considered reliable indicators of the physiological strength of plants exposed to abiotic stresses. Higher RWC/lower EL values are regarded as an indication of enhanced stress tolerance. The plants exhibiting higher RWC and lower EL are more capable of reducing stress-induced oxidative damage than their counterparts [42]. Based on this notion, the tomato hybrids VRTH-16-3 and VRTH-17-68 were considered tolerant hybrids under drought stress. Compared to their corresponding counterparts, both of these hybrids also showed balanced assimilation/metabolic processes, as evident by their improved reproductive growth measured in terms of fruit length, fruit width, single fruit weight, the total number of fruits, total soluble solids and total yield (Tables 3 and 4; Figure 6A,B).
The significant increase (Table 5) in the yield and related attributes of VRTH-16-3 and VRTH-17-68 hybrids may have been due to their improved ability to maintain efficient transcriptional control of stress-responsive genes' increase in photosynthesis by avoiding transpiration [43]. Furthermore, increased catalase activity may have resulted in the enhanced scavenging of reactive oxygen species, particularly hydrogen peroxide, which might lead to decreased membrane damage in the tomato fruits. Correspondingly, low membrane damage and enhanced proline content could have significantly improved the fruit diameter, size and yield of the tomato hybrids VRTH-16-3 and VRTH-17-68 more efficiently than their respective counterparts and control. Changes in plant metabolites are often associated with impairment in the photosynthetic process driven by accelerated water loss, thus, severely affecting the biosynthesis of essential nutritional components [44]. Our results have shown that in VRTH-16-3 and VRTH-17-68 hybrids, a positive correlation was observed between the RWC and biosynthesis of critical dietary components such as lycopene and ascorbic acid photosynthetic pigment contents, i.e., chlorophyll and carotenoids (Table 4). In contrast, the rest of the hybrids/parents showed a significant decrease in the nutritional components and chlorophyll and carotenoid contents. The enhanced biosynthesis of dietary components and photosynthetic pigment contents in the above hybrids could be due to their increased water-use efficiency, which improved their CO 2 assimilation, thereby maintaining their photosynthesis and water loss under drought stress conditions [45]. Suppressing nutritional factors such as lycopene, ascorbic acid and other related traits under drought stress is a natural phenomenon that varies greatly depending on the plant's idiotype and organ [42]. In this study, lycopene and ascorbic acid were affected more than total soluble solids (TSS), where the effect was observed less for the tomato hybrids VRTH-16-3 and VRTH-17-68. The firm level of lycopene and ascorbic acid content for these tomato hybrids could be due to the enhanced expression of stress-responsive genes that have strengthened the catalase activity, leading to improved drought stress tolerance [44].
A higher proline and lower malondialdehyde accumulation are considered the best biochemical markers for plant drought tolerance [45]. In our research, higher concentrations of proline and lower concentrations of malondialdehyde were observed in VRTH-16-3 and VRTH-17-68 hybrids compared to the rest of the plants (Table 4). Higher proline accumulation under stress conditions has improved stress tolerance in several plant species [46]. Drought-exposed plants suffer from an extensive loss of reducing power due to limited CO 2 assimilation that results in the accumulation of reactive oxygen species; thus, main-taining a higher level of antioxidant enzymes becomes imperative to prevent oxidative damage. We observed a higher activity of catalase enzyme coupled with a low level of H 2 O 2 for VRTH-16-3 and VRTH-17-68 hybrids compared to their corresponding counterparts (Tables 3 and 4). The increased catalase enzyme activity in the above hybrids could have been due to improved physiological activity, such as low electrolytic leakage, decreased lipid peroxidation and enhanced chlorophyll content, and the up-accumulation of proline has rendered the harmful effect of generated H 2 O 2 [46].
The overexpression of specific stress-responsive genes has been implicated in improving drought stress tolerance by increasing plants' antioxidant capacity, thus, enhancing the ROS scavenging and reducing the toxic effect of drought stress on plants [47]. We observed that all three variants of NAC proteins viz., NAC1, NAC1.1, NAC2 and HSP showed high constitutive expression in VRTH-16-3 and VRTH-17-68 hybrids compared to their respective counterparts and were clustered in one group. In contrast, LEA, ZFP and EF genes showed moderate expression ( Figure 7B). The enhanced tolerance of VRTH-16-3 and VRTH-17-68 hybrids could be due to the stimulatory growth role of the NAC genes in drought tolerance that might have enhanced the expression of several defence-responsive genes that have rendered better ion/osmotic homeostasis, induced antioxidant capacity and reduced water loss, thereby minimizing drought stress-induced oxidative damages [48]. However, future studies are required to decipher the detailed functional mechanism of NAC genes that can further our understanding of their molecular roles in drought tolerance in tomato plants.

Database Search, Sequence Retrieval and Comparative Phylogeny
The NCBI database was used to retrieve the sequence of tomato NAC1 protein (Locus: NP_001234482.1), and the NCBI BLAST server http://blast.ncbi.nlm.nih.gov/Blast.cgi (assessed on 15 July 2022) [49] was exploited to identify relevant sequential and functional homologues available for this protein in tomatoes. The NAC1 protein sequence showing maximum percentage identity and query cover was used for its functional characterization and phylogenetic analysis. The Bio-Edit tool [50] was used to analyse the sequence alignment and construct a phylogenetic tree using the MEGA 7 suite http://www.megasoftware.net (assessed on 15 July 2022) [51]. A multiple sequence alignment analysis of the selected protein was performed using a CLC bio workbench. In addition, functional NAC domain regions occupied by NAC proteins were analysed using InterProScan http://www.ebi.ac.uk/Tools/pfa/iprscan (assessed on 15 July 2022) [52], NCBI CDD server http://www.ncbi.nlm.nih.gov/Structure/cdd/cdd.shtml (assessed on 15 July 2022) [53] and ExPASy-Prosite scan http://prosite.expasy.org/scanprosite (assessed on 15 July 2022) [54] analysis.

Identification of Conserved Motifs
MEME (Multiple EM for Motif Elicitation) https://meme-suite.org/meme/ (assessed on 18 July 2022) Suite 5.0.5 was exploited to analyse the functional motifs in NAC1 protein.
A MEME analysis was performed by default selection parameters [55]. The MAST alignment tool was further used for the structural annotation of conserved motifs to identify functional motif patterns from other orthologs.

Homology Modelling of Tomato NAC Protein and Superimposition
The homology modelling of the SlNAC1 protein was carried out by the MODELLER software package (version 7.7) [56]. The template proteins viz., the structure of the conserved domain of Arabidopsis NAC (PDB ID: 1UT7), the crystal-structural of the conserved domain of rice stress-responsive NAC1 (PDB ID: 3ULX) and Arabidopsis NAC019 NAC domain crystal (PDB ID: 4DUL) were used to generate the 3D structure of proteins using the MODELLER module of Discovery Studio 3.0 [56]. The query SlNAC1 protein was refined for Cα traces using the ModRefiner tool http://zhanglab.ccmb.med.umich.edu/ ModRefiner (assessed on 18 July 2022) [57]. The best model structure was then used for superimposition and docking study [58].

DNA-Protein Docking Analysis
The Hex 8.0 molecular docking server was used to analyse DNA-protein interaction [67]. The DNA sequence to structure server http://www.scfbio-iitd.res.in/software/ drugdesign/bdna.jsp was used to model DNA sequences that specifically interact with the modelled SlNAC1 protein (assessed on 19 August 2022). The docking was performed using default parameters, and the docked complex was analysed using Discovery Studio 3.0 to reveal critical amino acid residues involved in DNA-protein interaction.

Active Site Prediction
The amino acid residues involved in creating the active site for mediating DNAprotein interaction were analysed using the metapocket server http://metapocket.eml.org (assessed on 13 September 2022) [73].

Plant Materials and Stress Condition
The tomato plants were subjected to drought stress at the experimental farm of the Indian Institute of Vegetable Research, Varanasi, UP, India, during the post-rainy season (October to March). The physical and chemical properties of soil and weather during the experimental period are represented in Table 5, respectively. During the experiment, the mean day and night temperature (DNT) and mean relative humidity (RH) were 22/32 • C and 75/70% with a 12 h light/dark period. The mean water content in the soil during the entire experimental period was 40-45% of field capacity. The tomato hybrids viz., VRTH 17-2, VRTH 17-68, VRTH 17-81 and VRTH 16-3 and their respective parents viz., Punjab Keshari (drought tolerant, DT) and Suncherry (drought susceptible, DS) were screened for their ability to tolerate drought stress under open field condition. The randomized block design (RBD) was used to design experiments with three treatments and replicated thrice. Each treatment (two treatments (T1 and T2) and control) had six plots, and each plot had 35 plants. The control plots were irrigated as per the requirement of the agronomic condition, whereas in the treatment plots (a) treatment 1 (T1), drought stress was imposed at a vegetative state in which irrigation was withheld for 60 days; (b) treatment 2 (T2), wherein the irrigation was withheld to the crops for 80 days (between 40 and 120 days after transplanting). The first and last rows were used as border rows in each treatment plot, including the control plot. In contrast, four central rows were used for recording the yield data and morphophysiological and biochemical evaluations at the end of the stress imposition. After completion of the experiment, the fully expanded third leaf from the top was stored at −80 • C for physiological, biochemical and molecular analysis. The fruit samples were analysed 14 days after the completion of 80 days of drought stress conditions when the tomato fruits attained the red stage.

Estimation of Leaf Water Potential and Chlorophyll Contents
The relative water content (RWC) was estimated following the standard protocol [74]. Leaf discs of a uniform size devoid of midrib were used to measure the fresh weight (FW), turgid weight (TW) and dry weight (DW) and the leaf RWC was calculated as per the standard equation. Chlorophyll contents were analysed by grinding fresh leaf samples in 80% acetone per the protocol [75]. The chlorophyll content was measured per Arnon's equation [76] and expressed as mg g −1 FW. Total suspended solids (TSS) in tomato fruits were measured using a refractometer. Fruit diameter, single fruit weight, the number of fruits per plant and yield per plant were counted manually. All the measurements were carried out in triplicate.

Estimation of Electrolytic Leakage and Lipid Peroxidation
The electrolytic conductivity of leaf latchets was measured following the standard method devised by Deshmukh et al. [77]. The electrolytic leakage was recorded at 40 • C (C1) and 100 • C (C2) using a conductivity meter (Century Instruments, Chandigarh, India). The measurement of lipid peroxidation was performed using the standard method [78]. Fresh leaf samples were extracted in 10% tri-chloroacetic acid (TCA) containing 0.25% of 2-thiobarbituric acid, and the malondialdehyde level was calculated using the extinction coefficient of 155 mM −1 cm −1 and expressed as µmol g −1 FW.

Estimations of Lycopene and Ascorbic Acid
The lycopene content was estimated as per the standard protocol [79]. The tomato fruits were extracted in petroleum ether and acetone (3:2) simultaneously. After centrifugation, the OD of the supernatant was measured at 503 nm and expressed as mg of lycopene/100 g following the standard equation. The ascorbic acid [80] in the tomato fruit samples was estimated by titrating the fruit samples homogenized in 4% oxalic acid solution against 2, 6-dichlorophenol indophenol dye solution. The end point of the titration was measured with the development of a light pink colour, which persisted for 10 s. The ascorbic acid in tomato fruits was expressed as mg/100 g following the standard equation.

Measurements of Hydrogen Peroxide and Catalase Activity
The hydrogen peroxide content in the fresh samples was measured using the standard method [81]. The leaves were extracted in sodium phosphate buffer (pH 7.2) and 0.1% titanium sulphate, and the absorbance was measured at 410 nm. The hydrogen peroxide generated was measured using the extinction coefficient of 0.28 µM −1 cm −1 and expressed as µmol g −1 FW. As Bates et al. [82] described, the proline level in the leaf samples was estimated. Leaf samples were extracted in 3% sulfosalicylic acid and incubated at 100 • C for 1 h. After incubation, the samples were cool, and toluene was added. The absorbance containing acid ninhydrin was recorded at 520 nm and expressed as µg g −1 FW. The activity of the antioxidative enzyme catalase was assayed per the standard method [83]. The leaf samples were homogenized in 50 mM of Tris NaOH buffer (pH 8.0), and the absorption of the supernatant was recorded at 240 nm.

RNA Isolation and cDNA Preparation
TRIZOL reagent (Invitrogen, Waltham, MA, United States) was used to isolate total RNA following the manufacturer's recommendations. The quality of the total RNA was checked using a nanophotometer (Implen, Westlake Village, CA, USA), and integrity was confirmed on 1.0% agarose gel electrophoresis. The isolated RNA (1.0 µg) was then used to synthesize the first strand of cDNA using an iScriptTM cDNA synthesis kit (Bio-Rad Laboratories, Hercules, CA, USA).

Real-Time (RT-PCR) Gene Expression Analysis
An SsoFast TM EvaGreen ® Supermix detection chemistry (Bio-Rad) iQ5thermocycler (BioRad Laboratories, Hercules, CA, USA) was used for the real-time PCR analysis using gene-specific primers (Supplementary Materials, Table S2). The reactions were formulated in 20 µL, and relative quantification was carried out per the standard 2 −∆∆CT method given by Livak and Shmittgen [84], using ACTIN mRNA level as an internal control. The gene expression analysis was performed as three biological and technical replicates to generate a heat map.

Statistical Analysis
The statistical analysis was performed using two-way analyses of variance (two-way ANOVA) with two fixed factors, i.e., genotype and treatment. All the data were recorded in three biological replicates. The Duncan Multiple Range Test (DMRT) was used to compare means at a 0.05 significance level using SPSS 21.0 software, IBM Corporation, Armonk, NY, USA.

Conclusions
In the present study, the hybrids viz., VRTH-16-3 and VRTH-17-68 ameliorated drought-induced oxidative stress more efficiently by competently modulating their waterholding capacity, membrane damage, proline biosynthesis and up-accumulation of NAC genes, as well as other stress-responsive genes, compared to different hybrids. The present study confirmed that drought stress differentially modulated the expression of various stress-responsive genes, enzymes and metabolites in tomato hybrids. Likewise, the differential expression of NAC genes in the present study attested to their regulatory roles in ameliorating drought stress-induced oxidative damage in tomato hybrids by stimulating their physiological, biochemical and molecular defence response, thereby contriving the destructive effect of drought-induced oxidative stress. However, the molecular mechanisms behind this seminal role of NAC genes in tomatoes need to be addressed in detail.
Further, the pilot study's bioinformatic analysis of SlNAC1 genes also provided valuable insight into their structural and functional classification and evolutionary relatedness with other members of NAC gene families. The docking analysis recognized that Arg 63 , Pro 73 , Asn 74 , Trp 82 , Lys 83 , Ala 84 , Thr 85 , Gly 86 , Ala 87 and Asp 88 residues improve NAC protein plant growth and developmental process under abiotic stress conditions. Furthermore, the generated SlNAC1 model was reliable in quantitative and qualitative assessment and can be exploited for future studies in other crop plants. Overall, the information generated in the present study will unravel the complex regulatory networks associated with NAC TFs in plants exposed to drought stress conditions. In addition, integrating NAC genes in the tomato breeding program will help breeders develop climate-smart crops with enhanced adaptability and survival under adverse conditions, thus, ensuring agricultural sustainability and global food security.

Supplementary Materials:
The following supporting information can be downloaded at: https://www. mdpi.com/article/10.3390/plants11212930/s1, Figure S1. Prediction of the location of NAC domain in SlNAC1 protein using InterProScan scan analysis; Figure S2. Submission of our predicted model to PMDB and assigned PMDB ID (PM0081904) for our submitted structure; Figure S3. Superimposition results represented with their respective global and local RMSD-values showing the structural conservation of SlNAC1 domain. Superimposition results revealed that SlDREB1were found to be structurally closer to 1UT7 template compared to 4DUL, 3ULX templates as evident from their variations for RMSD values. Sequence alignment between SlDREB1 with the respective templates (1UT7, 4DUL, 3ULX) predicts the synonymous (blue) and non-synonymous substitutions (white) (B) Qualitative analysis of predicted model using PROCHECK and ProSA analysis. (I) The stereo chemical spatial arrangement of amino acid residues in predicted model (SlDREB1) were compared with experimentally resolved structures (1UT7, 4DUL, 3ULX) through PROCHECK server. Most favoured regions are coloured red while additionally, generously and disallowed regions are indicated with yellow, light yellow and white fields respectively. (II) Qualitative estimation by ProSA server, that measures structural errors in each amino acid residues; Figure S4. Binding energy of docked complex of SlNAC1 protein with DNA element. Figure S5. Gene ontology analysis using ReviGO web server. The functional and significant GO terms involved in (I) biological processes and (II) molecular function are shown on scattered plot diagram using hypergeometric test distribution of SlNAC1. (B) The subcellular localization and functional gene annotation of SlNAC1 using CELLO2GO web server. The significant terms are represented in terms of their percentage contribution; Table S1. Functional annotation, accession ID and interacting score values for both first and second shell of interactors that form mutual interactive associative network along with SlNAC. The higher value of score indicate the more frequent interaction exist between two associated proteins; Table S2. Primers used for quantification of mRNA levels by qRT-PCR.