The whole-genome protein sequence dataset for G. arboreum, G. barbaense, and G. hirsutum was downloaded from the CottonGen database (https://www.cottongen.org/), and the dataset for G. raimondii was obtained from Phytozome, version 12 (https://phytozome.jgi.doe.gov/pz/portal.html) (Nordberg et al., 2014). Total protein sequences for other plant species from different taxonomic groups were downloaded from the website of the National Center for Biotechnology Information (NCBI) (http://www.ncbi.nlm.nih.gov/). A total of 58 HD2 protein sequences from 42 plant species were obtained from the NCBI and utilized to construct Hidden Markov model (HMM) profiles. This profile of the HD2 domains (the HD2 label and the catalytic, regulatory, and zinc- finger domains) was employed as a query to identify HD2 gene family members using HMMER (V3.0) (Finn et al., 2015). Protein sequences and CDSs for G. arboreum, G. raimondii, G. hirsutum, and G. barnadense were also downloaded from the CottonGen database (https://www.cottongen.org/) (Yu et al., 2014). All hits were queried in the Pfam (http://pfam.xfam.org/) and InterProScan (http://www.ebi.ac.uk/interpro/search/sequence-search/) databases to verify the presence of conserved domains. The ProtParam (http://web.expasy.org/protparam/) tool offered by Expasy was used to estimate the physicochemical parameters of Gossypium HD2 proteins. The ProtParam tool was also used to estimate biophysical and biochemical properties, such as number of amino acids, molecular weight, grand average hydropathy (GRAVY), theoretical isoelectric point (pI), aliphatic index, and instability index. The cotton HD2 gene subfamilies were named as per the orthologous HD2 members in the A. thaliana genome.

The CELLO server, v.2.5 (http://cello.life.nctu.edu.tw/), was used to predict the probable subcellular locations of all the cotton HD2 proteins identified. The Pfam database was used to analyze the protein domain structures, and the Illustrator for Biological Sequences software package (http://ibs.biocuckoo.org/) (Liu et al., 2015a) was subsequently utilized to construct a schematic diagram of protein functional domains. Nuclear localization signals were examined using Motif Scan (https://myhits.isb-sib.ch/cgi-bin/motif_scan).

HD2 protein sequences were extracted from the complete genome sequences of 39 plant species. Multiple sequence alignments (MSAs) of full-length HD2 protein sequences with identified cotton HD2 proteins were performed using the Clustal X program (http://www.clustal.org/) (Larkin et al., 2007) with default parameters. Subsequently, these aligned sequences were used to construct a phylogenetic tree. The MEGA 7 software package (http://www.megasoftware.net/) (Kumar et al., 2016) was utilized to build an unrooted phylogenetic tree using the maximum likelihood (ML) method with the following parameters: pairwise gap deletion, JTT matrix-based model, and 1000 bootstrap values (Felsenstein, 1985).

The physical (genomic) locations on chromosomes of all Gossypium HD2 genes were obtained by cross-referencing the results of searching for HD32 CDS sequences using the BLASTN search program against the resources of the CottonGen cotton database (https://www.cottongen.org/) and the Phytozome database, version 12 (https://phytozome.jgi.doe.gov/pz/portal.html). All HD2 genes of Gossypium species were mapped on the chromosomes using the MapInspect software (http://mapinspect.software.informer.com/).

Paralogous HD2 genes were identified through reciprocal BLAST analysis with e-value <10^-5 in order to deduce the evolutionary mechanism of the HD2 gene family in Gossypium species (G. arboreum, G. raimondii, G. hirsutum, and G. barbadense). As per the reciprocal BLAST output, duplication events in Gossypium HD2 genes were determined using the MCScanX toolkit (Wang et al., 2012b). Additionally, Ka/Ks analysis of ortholog and paralog sequences was carried out using the PAL2NAL program. Subsequently, Ks values were used to compute the approximate dates of duplication and speciation events via the formula T = Ks/2λ, where λ is the assumed value of the clock-like rate. In addition, the Ka/Ks ratio was used as an indicator of selection pressure for the duplicated HD2 genes. Specifically, Ka/Ks ratios greater than 1, less than 1, and equal to 1 were taken to suggest positive, negative, and neutral (purifying selection) evolution, respectively (Zhang et al., 2006).

In G. arboreum, G. raimondii, G. hirsutum, and G. barbadense species, the gene structures of identified HD2 genes were examined via comparison of the predicted HD2 coding sequences with their corresponding genomic sequences using the GSDS online tool (Gene Structure Display Server; http://gsds.cbi.pku.edu.cn/) (Hu et al., 2015).

A MEME (Multiple Expectation maximization for Motif Elicitation) tool (http://meme.nbcr.net/meme/) (Bailey et al., 2009) was used for the identification of conserved protein motifs in identified cotton HD2 proteins. The following parameters were used in this analysis: optimum width = 6 to 300; one per sequence; maximum number of motifs to find = 20. Furthermore, these motifs were annotated using the InterProScan program (Quevillon et al., 2005).

To investigate the expression pattern of cotton HD2 members under salt and drought stress conditions, raw Illumina RNA-Seq data for G. hirsutum (accession number: PRJNA532694) were obtained from the Sequence Read Archive of the National Center for Biotechnology Information (the NCBI-SRA database). The raw reads were filtered using the Trimmomatic tool (Bolger et al., 2014), and high-quality reads were subsequently mapped on the G. hirsutum genome (Wang et al., 2019) using the STAR aligner (Dobin et al., 2013) with default parameters. Transcript abundance was estimated using the StringTie software (Pertea et al., 2015), and the differential expression of genes was evaluated using the edgeR Bioconductor package (Robinson et al., 2010).

The regulatory mRNA sequences of HD2 family members were validated via qRT-PCR. Stress conditions were employed on eight-week-old G. hirsutum plants. The plants were raised in pots (30 cm diameter, 25 cm height, 12 liter capacity) at the CSIR-NBRI garden, with soil made up of 40% silt, 20% clay, and 40% sand. They were exposed to a minimum temperature of 27°C and a maximum of 34°C; there were 14 hours of light and 10 hours of darkness per day. To create conditions of drought and salt stress, plants (in triplicate) were administered a single treatment of 20% PEG8000 solution (Shafiq et al., 2015) and 300-mM NaCl (Wei et al., 2017). Untreated plants were taken as controls. During the experiment, leaves were taken from the topmost juvenile section of each plant at the 0h, 12h, 24h, 48h, and 72 h time points after administration of the stress treatment. RNA was extracted from these samples using a Spectrum™ Plant Total RNA Isolation kit (Sigma, USA). cDNA was prepared with 1µg RNA using the Verso cDNA synthesis kit (Thermo Scientific) following the provided protocol. 2× diluted cDNA was utilized to examine the expression of nine selected putative genes via a fluorescent quantitative detection system using qRT-PCR (HiMedia Insta Q48 M4); a 20 µl reaction was prepared with 10 µl SYBER™ green master mix (Applied Biosystems), along with 5 pmol of primer concentration. Primers were designed with the help of Primer-BLAST, and ubiquitin was taken as a normalizing control. The qRT-PCR program was set to apply a temperature of 95°C for 2 min, followed by 40 cycles of denaturation at 95°C for 15 s and annealing at 60°C for 1 min. Gene expression levels were calculated with reference to ubiquitin for normalization, and the relative expression levels of the selected genes at each time point were compared with their expression levels at 0h. CT values were calculated with mean ± SD by the2-∆CT method (Livak and Schmittgen, 2001).

Based on the higher expression of GhHDT3D.2 observed in RNA-Seq and in qRT-PCR values, the co-expression network of GhHDT3D.2 was constructed for both stress conditions (drought and salt). This co-expression network was built using FPKM values using the “expression correlation networks” module of Cytoscape version 3.8.2 (Otasek et al., 2019). Specifically, this module was used to identify positive and negative Pearson correlations (r ≥ 0.95 and r ≤ - 0.95) between interacting members of the network. Co-expressed genes and networks were visualized using Cytoscape in a circular force-directed layout. The key metabolic pathways of genes positively and negatively co-expressed with GhHDT3D.2 (PCoEGs and NCoEGs) were evaluated using the PageMan software, version 3.5.1 (Thimm et al., 2004). To determine their functional categories, the average statistical test Benjamini Hochberg test were employed. Gene ontology (GO) analysis for three categories (biological processes, molecular functions, and cellular components) was carried out using agriGO, v2.0 (http://systemsbiology.cau.edu.cn/agriGOv2/). Singular enrichment analysis (SEA), with statistical testing at a p-value threshold of < 0.05, was used to retrieve GO annotations.

To analyze the cis-regulatory elements in the upstream sequences of HD2 genes, 1.5 kb upstream region sequences of 27 cotton HD2 genes were obtained from the CottonGen (https://www.cottongen.org/) and Phytozome v12 (https://phytozome.jgi.doe.gov/pz/portal.html) databases. Cis-regulatory elements were identified using the PlantCARE database (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) (Lescot et al., 2002).

Four cotton species, viz., G. arboreum, G. raimondii, G. hirsutum, and G. barbadense, were used to identify HD2 members. HD2 domain sequences were obtained from the NCBI and used to construct Hidden Markov model (HMM) profiles. The HMMER search program was used to identify HD2 orthologs by querying against these four Gossypium species. The HD2 gene was later analyzed for similarity and conserved domains by comparing data from the Pfam and InterproScan databases. A total of 27 HD2s were identified ( Table 1 ): 4 GaHD2s (G. arboreum), 5 GrHD2s (G. raimondii), 9 GbHD2s (G. barbadense), and 9 GhHD2s (G. hirsutum). Number of amino acids, instability index, grand average hydropathy (GRAVY), theoretical isoelectric point (pI), molecular weight, and aliphatic index are important physiological parameters in determining the primary structure of proteins. The properties of the identified cotton HD2s, including protein length (aa), gene name, locus ID, isoelectric point (pI), molecular weight (Da), chromosome location, number of introns, sub-cellular localization, aliphatic index, instability index, and GRAVY, were analyzed using the ProtParam tool. Additionally, the subcellular localizations of predicted HD2 members were identified using the CELLO v2.5. There was wide variation in the length of cotton HD2 proteins, ranging from 201 to 299 amino acids in G. arboreum, 272 to 304 amino acids in G. raimondii, 240 to 304 amino acids in G. hirsutum, and 285 to 376 amino acids in G. barbadense. Isoelectric points ranged from 4.52 to 6.02, 4.19 to 5.06, 4.50 to 4.86, and 4.86 to 5.29 in G. arboreum, G. raimondii, G. barbadense, and G. hirsutum, respectively. Polypeptide molecular weight ranged from 21.588 kDa to 32.484 kDa in G. arboreum, 29.604 kDa to 32.838 kDa in G. raimondii, 26.563 kDa to 32.875 kDa in G. hirsutum, and 30.424 kDa to 40.694 kDa in G. barbadense. Interestingly, Gossypium HD2s contained a large number of introns; this count varied from 5 to 8, 7 to 9, 5 to 9, and 7 to 11 in G. arboreum, G. raimondii, G. hirsutum, and G. barbadense, respectively. Cotton HD2 members were found to be located in the nucleus ( Table 1 ).

Proteins with an instability index < 40 are considered stable, whereas those with a value > 40 are considered unstable. The volume occupied by aliphatic side chains relative to the overall volume of the protein is termed the “aliphatic index”. In G. arboreum, G. raimondii, G. hirsutum, and G. barbadense, instability index values for HD2 genes ranged from 42.44 to 63.31, 37.93 to 62.94, 35.89 to 62.35, and 42.13 to 62.35, respectively; and aliphatic index values varied from 39.35 to 56.65, 45.48 to 60.62, 44.75 to 57.75, and 44.41 to 65.80, respectively. Thus, the majority of the predicted HD2 members were hydrophilic in nature and unstable, a finding consistent with the instability index measure. However, some HD2s (such as GhHDT3D.1, GhHDT3A.1, and GrHDT3.1) were stable. Moreover, HD2 proteins had negative GRAVY scores, indicating their hydrophilic behavior; however, the degree of hydrophilicity demonstrated a greater variability ( Table 2 ).

In order to perform multiple sequence alignment (MSA) of cotton HD2 protein sequences, ClustalX 2.1 was used. Comparison of cotton HD2s via the Pfam and InterPro protein domain databases revealed that they contained highly conserved domains relative to the structure of a typical HD2 protein ( Figure 1 ). These results indicated that GaHDT1.1, GaHDT1.2, GaHDT3.1, and GaHDT3.2 contained the HD2 label, a deacetyl/catalytic domain, phosphorylation sites, mono- and bipartite NLS motifs, and a zinc finger domain from N-terminus to C-terminus. Similarly, GrHDT1.1, GrHDT 1.2, GrHDT 3.2, GhHDT1A.1, GhHDT1A.2, GhHDT1D.2, GhHDT3A.1, GhHDT3A.2, GhHDT3D.1, GhHDT3D.2, GbHDT1A.1, GbHDT1A.2, GbHDT1D.1, GbHDT1D.2, GbHDT3A.1, GbHDT3A.2, GbHDT3D.1, and GbHDT3D.2 contained all the conserved domains. In contrast, GrHDT3.1, GrHDT4, GhHDT2A, and GbHDT2A did not contain a bipartite NLS motif or zinc finger domain, and GhHDT1D.1 did not contain a zinc finger domain at the C-terminus end ( Figures 2 , 3 ). A monopartite NLS motif was predicted to be present in all the cotton HD2 proteins, suggesting that this might generate a signal that drives cotton HD2s to the nucleus. Most HD2 proteins were between 232 and 384 amino acids in length. A highly correlated relationship exists between amino acid sequence variation and the length of the highly variable regulatory domain, which resides in the center (Bourque et al., 2016). These results suggest a high degree of homology between the amino acids in these conserved cotton HD2s and AtHD2s ( Figure 3 ).

To evaluate the classification of HD2 genes in cotton, the sequence features of 44 plant species were analyzed, and an unrooted phylogenetic tree of the HD2 genes was constructed ( Figure 4 ). Subsequently, an investigation of the evolutionary relationships between HD2 proteins in cotton and those in other plant species was carried out. For this purpose, MSA was conducted on 27 cotton HD2 proteins and 198 HD2 proteins from other plant species (basal angiosperm, bryophytes, lower plants, Lycopodiophyta, monocots, and eudicots) ( Table 3 ). The maximum likelihood (ML) method was used for phylogenetic tree construction. The nomenclature of cotton HD2s was done as per A. thaliana HD2s.

Based on the phylogenetic relationships plant HD2 genes were clustered into ten major groups (I to X), and cotton HD2 proteins were classified into four of these groups, with none belonging to Group II, IV, VI, VIII, or IX. Groups I, II, III, VII, and VIII were further divided into two, two, two, three, and two subgroups respectively (labeled a, b, and c). Group III contained the largest number of cotton HD2 genes (13 members) followed by Group I with 11 and Groups V, VIII, and IX, with one HD2 member each ( Figure 4 ).

To understand the evolutionary relationships and gain insight into the structural diversity and similarity of Gossypium HD2s, a separate maximum likelihood tree was constructed with 1000 bootstrap values using the HD2 protein sequences of Gossypium species. The topology of the tree, HD2 duplicating nodes, conserved motifs, and exon/intron distribution enabled the categorization of the cotton HD2s into six sub-families (I - VI, shown in different background colors in Figure 5 ). Genes within the same sub-family with the highest levels of identity (>80%), indicate divergent evolution from a common ancestor, in the orthologous gene pairs ( Figure 5A ).

To understand the structural variations in cotton HD2s, coding sequences and the corresponding genome sequences were compared to determine their exon/intron arrangements. Interestingly, the majority of cotton HD2 gene structures within the same clade were similar, i.e., group-specific exon/intron patterns were observed. The majority of cotton HD2 members had introns ranging from 5 to 1. Specifically, the number of introns varied from 5 to 8 (GaHD2s), 7 to 9 (GrHD2s), 5 to 9 (GhHD2s), or 7 to 11 (GbHD2s) ( Figure 5B ). The gene structure of all cotton HD2s was analyzed and compared in order to examine the stability of the exon/intron members with respect to the structure of the phylogenetic tree. It was found that the majority of cotton HD2s falling under the same clade shared similar arrangements of exons/introns. For example, in subfamily I, the HD2 gene comprised 9 introns with 10 CDS (with the exceptions of GrHDT1.1 and GaHDT1.1), while members of the V group comprised 7 introns with 8 CDS (with the exceptions of GaHDT3.2, GhHDT3A.2, and GrHDT4, which consisted of 5, 8, and 7 introns with 0, 8, and 7 CDS and 6, 0, and 0 exons, respectively).

A motif-based sequence analysis tool (MEME) was utilized to examine the conserved motifs in HD2 proteins. The InterProScan was used to annotate these motifs. A total of 20 conserved motifs were determined in the Gossypium HD2 members. Motifs 1 and 4 were found in the database and annotated as nucleoplasmin-like domains. These have been reported to be important players in governing the dynamic architecture of the chromatin (Singh et al., 2022) ( Figure 5C and Supplementary Table 1 ). Motif 4 (HD2 pentapeptide, shown in violet) was approximately present in all the cotton HD2s and showed a conserved N-terminus pentapeptide sequence. Motif 2, annotated as a zinc finger domain, was found in all the cotton HD2 proteins except GhHDT1D.1, GbHDT2A, GhHDT2A, GrHDT3.1, and GrHDT4 ( Figure 5C ). Additionally, motifs 3 and 5 were annotated as histone deacetylase- like domains and were found to be present in all cotton HD2 proteins except GaHDT1.2, GbHDT2A, GhHDT3A.1, and GhHDT3D.1 members. Most of the cotton HD2s exhibited a similar motif composition to their close evolutionary relatives; therefore, it is speculated that they might have similar functions ( Figure 5C ).

In order to determine the chromosomal distribution within the HD2 gene family in Gossypium, BLASTN search was carried out for all GaHD2, GrHD2, GhHD2, and GbHD2 members in G. arboreum, G. raimondii, G. barbadense, and G. hirsutum. After all the Gossypium genes had been mapped to their corresponding chromosomes, GaHD2 genes were localized on chromosomes 2, 4, 10, and 11 ( Figure 6A ), while the distribution of GrHD2 genes consisted of locations on chromosomes 2, 6, 7, and 9 ( Figure 6B ). Similarly, in G. hirsutum, more HD2 genes were positioned on A_T (1, 5, 7, 9, 11) than on D_T chromosomes (1, 5, 9, 11), with 5 and 4 genes, respectively ( Figures 6C, D ); in G. barbadense, a larger number of HD2 genes were also located on A_T (1, 5, 7, 9, 11) than on D_T chromosomes (1, 5, 9, scaffold3329), with 5 and 4 genes, respectively ( Figures 6E, F ).

Analysis of gene duplication events identified two pairs of paralogous HD2 members in G. arboreum ( Figure 7A ), while no paralogous genes were detected in G. raimondii ( Figure 7B ). Furthermore, four pairs of paralogous HD2 members were identified in G. hirsutum and G. barbadense ( Figures 7C, D ). As shown in Figure 7 , all paralogous gene pairs were found to be located on distinct chromosomes, providing a considerably strong indication that Gossypium HD2 members were expanded through segmental duplication. In G. arborium, two segmental duplications (GaHDT1.1/GaHDT3.1 and GaHDT3.1/GaHDT3.2) occurred between 45.53 and 18.06 MYA. Furthermore, in G. hirsutum, four segmental duplications (GhHDT3A.1/GhHDT3D.1, GhHDT1A.1/GhHDT1D.1, GhHDT1A.2/GhHDT1D.2, and GhHDT3A.2/GhHDT3D.2) occurred 2.51, 1.32, 2.12, and 1.42 MYA; finally, in G. barbadense, four segmental duplications (GbHDT3A.1/GbHDT3D.2, GbHDT1A.2/GbHDT1D.2, GbHDT1A.1/GbHDT1D.1, and GbHDT3A.2/GbHDT3D.1) were detected before 1.76, 1.44, 2.09, and 1.60 MYA ( Table 4 ).

The ratio of non-synonymous to synonymous substitutions (Ka/Ks) provides an important measure of the pressure of evolutionary assortment on the relevant amino acids (Hurst, 2002). This ratio was measured for the duplicated Gossypium HD2 gene pairs (paralogous). The result of this calculation showed that all the duplicated HD2s had a Ka/Ks ratio < 1 ( Table 4 ). Therefore, duplicated gene pairs in cotton HD2 members were under powerful purifying selection pressure.

A pair of homologous genes that have diverged in dissimilar species during a speciation event is known as an ortholog. The orthologous associations among the HD2 members in the four cotton species were determined. Specifically, orthologous genes among the HD2 members were identified via sequence similarity. Orthologs having sequence similarity equal to or greater than 90% in both cDNA and amino acid (aa) composition were selected for further evolutionary study. The selection pressure and probable functional divergence of Gossypium HD2s genes were examined through the calculation of Ka, Ks, and Ka/Ks ratio among the orthologs (A vs. D, A_T vs. A, D_T vs. D, A_T vs. A_T, and D_T vs. D_T) and within the homeologs (A_T vs. D_T). Notably, the Ka values of Gossypium HDT3 orthologs (group IV and V HD2s; GaHDT3.2/GrHDT3.2, GhHDT3A.1/GhHDT3D.1, GbHDT3D.1/GrHDT3.2, GhHDT3A.2/GbHDT3A.2, and GhHDT3D.2/GbHDT3D.1) and HDT1 orthologs (GbHDT1A.1/GaHDT1.2) were higher than those of other HD2 gene pairs, suggesting that these ortholog pairs underwent more rapid evolution of the relevant genes. With a Ka/Ks ratio <1, this analysis revealed that negative selection had been exerted on HD2 orthologous genes ( Supplementary Table 2 ), and some orthologous gene pairs had experienced directional selection. The Ka/Ks ratio was higher for HDT3 orthologous gene pairs, in A vs D, A_T vs D_T, A_T vs A, D_T vs D, A_T vs A_T, and D_T vs D_T. There were three HDT3 orthologous gene pairs (GaHDT3.2/GrHDT3.2, GbHDT3A.1/GaHDT3.1, and GhHDT3D.2/GrHDT3.2) that had a Ka/Ks ratio > 1, meaning that they might have undergone adaptation to certain advantageous alleles that might play an important role in cotton. These outcomes suggest that the diploid cotton HDT3 member was under greater evolutionary pressure, and that the evolution of the A genome might have been faster than that of the D sub-genome ( Figure 8 and Supplementary Table 2 ).

Promoters are regions of DNA that drive the initiation of transcription of a certain gene; these promoters are situated near the transcription start site (TSS) of the corresponding gene. In the present study, we searched for cis-regulatory elements in the 1.5 kb upstream promoter region of the identified GhHD2s. The outcomes of this analysis are listed in Supplementary Tables 3A, B . The results revealed that cis-regulatory elements like ABRE were most abundant, followed by TGACG and CGTCA motifs. Additionally, TCA-element was abundant in the promoter regions of identified GhHD2s ( Figure 9 ). These outcomes indicate that GhHD2 genes are associated with plant resistance to environmental stress conditions, including drought and salt stress.

It is interesting that the GhHD2 genes, namely GhHDT1A.1, GhHDT1A.2, GhHDT1D.1, GhHDT1D.2, GhHDT2A, GhHDT3A.2, GhHDT3D.1, and GhHDT3D.2, were found to comprise several hormone-responsive cis-elements, such as TCA-element (salicylic acid responsiveness), TGA-element (auxin-responsive element), ABRE (abscisic acid responsiveness), CGTCA-motif (MeJA-responsiveness), P-box (gibberellin-responsive), TGACG-motif (MeJA-responsiveness), and GARE-motif (gibberellin-responsive), in their promoter region ( Figure 9 ). MeJA is reported to stimulate the production of defensive compounds that may be used against pathogens, drought stress, heavy metal stress, low temperatures, and salt stress (Yu et al., 2019), while auxin plays a crucial role in plant responses to adverse biotic and abiotic conditions (Sharma et al., 2015). In addition to pathogenesis-related resistance, drought and heat stress also trigger the production of salicylic acid (SA), and this also reduces the concentrations of Na+ and Clˉ ions caused by salt stress conditions (Yuan and Lin, 2008; Emamverdian et al., 2020). Plant growth hormones, such as gibberellins (GAs) are important in increasing resistance under abiotic stresses (Emamverdian et al., 2020), and plants respond to downstream stress by integrating various stress signals via the abscisic acid (ABA) hormone (Tuteja, 2007). Overall, these promoter analyses show that these genes play an important role in providing tolerance to drought and salt stress conditions.

In order to determine the expression pattern of HD2 genes under drought and salt stress, HD2 genes were analyzed via the leaf RNA-Seq data under both types of stress. Nine identified GhHD2 genes showed differential expression, with three genes (GhHDT1A.2, GhHDT3D.1, and GhHDT3D.2) exhibiting higher expression during drought ( Figure 10B ) and salt stress ( Figure 10D ) compared to a control. Furthermore, these nine putative genes (GhHDT2A, GhHDT3A.1, GhHDT3A.2, GhHDT1D.2, GhHDT1A.1, GhHDT1A.2, GhHDT1D.1, GhHDT3D.1, and GhHDT3D.2) were validated through qRT-PCR. The list of primers is provided in Table 5 . In dought, the majority of GhHD2s exhibited stronger responses at 12h ( Figure 10A ), with the exception of GhHDT2A, while the strongest response was that of GhHDT3D.2 with respect to fold change (FC) compared with the control at 12 h (FC >29), 24h (FC >33), 48h (>16), and 72h (>15) ( Figure 10A ). Under salt stress conditions, the entire set of GhHD2 members exhibited higher levels of expression at 24h compared with the control ( Figure 10C ), while GhHDT3D.2 exhibited the highest level of expression in comparison to the control at 12 h (FC >4), 24h (FC >9), 48h (>14), and 72h (> 6). Overall, this analysis shows that these HD2 members might play a crucial role in regulating responses to drought and salt stress conditions.

The GhHDT3D.2 gene was selected for analysis of the co-expression network, pathways, and gene ontology due to its higher levels of expression under both types of stress. Co-expressed genes with GhHDT3D.2 were explored using log2FPKM expression values. A total of 15 genes positively and 14 negatively co-expressed with GhHDT3D.2 (PCoEGs and NCoEGs) were identified under drought stress conditions ( Supplementary Figures 1A, B ; Supplementary Table 4 ), while 80 such PCoEGs and 237 such NCoEGs ( Supplementary Figures 1C, D ) were identified under salt stress conditions ( Supplementary Table 4 ).

In order to understand the functional pathways and molecular importance of the PCoEGs and NCoEGs of GhHDT3D.2, a PageMan pathway analysis was performed. For drought stress conditions, this analysis demonstrated that there was higher expression of fatty acid metabolism and auxin response factors (ARFs) in relation to NCoEGs of GhHDT3D.2 and vacuolar-sorting protein (SNF7) in relation to PCoEGs of GhHDT3D.2 ( Supplementary Figure 1E ). Moreover, for salt stress conditions, diacylglycerol kinase ( Supplementary Figure 1F ) exhibited higher expression at all time intervals, while the highest level of expression among the PCoEGs of GhHDT3D.2 at 12h was that of brassinosteroid hormone metabolism; among the NCoEGs of GhHDT3D.2, the highest level of expression was that of phenypropanoid and flavonoids ( Supplementary Figure 1G ). Transcription factors such as basic leucine zipper (bZIP), WRKY, plant homeo-domain (PHD), and heat shock factors (HSF) exhibited expression with PCoEGs of GhHDT3D.2, and abiotic stress-responsive genes exhibited expression with NCoEGs of GhHDT3D.2 under salt stress conditions ( Supplementary Figure 1H ).

A gene ontology (GO) analysis was carried out to examine all the identified co-expressed genes and identify their functional roles. GO divided into three categories: biological processes (BPs), cellular components (CCs), and molecular functions (MFs). The significant PCoEGs and NCoEGs of GhHDT3D.2 were selected for gene ontology analysis. Under drought stress, enriched MFs were hydrolase activity, carbohydrate derivation binding, and purine nucleotide binding; enriched BPs were establishment and localization; and enriched CCs were the membrane, endomembrane system, and cytoplasmic part ( Supplementary Figure 1I and Supplementary Table 5 ). Likewise, under salt stress conditions, enriched BPs were response to water deprivation, abiotic stimulus, salt stress, ethylene, the auxin-activated signaling pathway, and gibberellins ( Supplementary Figure 1J and Supplementary Table 5 ); enriched MFs were transcription factor activity, transporter activity, calmodulin binding, phosphor-transferase activity, and hydrolase activity; and enriched CCs were the vacuolar membrane, golgi apparatus, plasma membrane, cytoskeletal part, and cell part ( Supplementary Figure 1J and Supplementary Table 5 ).