1. Introduction {#sec1-ijms-20-05311}
===============

Nucleosomes, the basic unit of chromatin, are composed of 147 bp of DNA wrapped around a histone octamer (two copies of each of H2A, H2B, H3, and H4 histone proteins). Nucleosomes are dynamic in response to developmental and environmental signals, thus altering the DNA accessibility and DNA-template processes to regulate the various processes in plants, including flowering time, root growth, and response to environmental changes \[[@B1-ijms-20-05311],[@B2-ijms-20-05311]\]. Cells use several mechanisms, including post-translational histone modification and DNA methylation, to regulate the gene expression. The N-terminal tails of histone are subjected to various post-translational modifications, including histone acetylation, methylation, etc. \[[@B3-ijms-20-05311]\]. Histone acetylation is carried out by histone acetyltransferases (HATs) in eukaryotes and is associated with transcriptional activation. Histone acetylation can be detected on different lysine residues of histone H3 and H4. For example, in *Arabidopsis*, K9, K14, K18, K23 and K27 of histone H3, and K5, K8, K12, K16, and K20 of H4 are acetylated \[[@B4-ijms-20-05311],[@B5-ijms-20-05311]\]. HATs are well-conserved in yeast, animals, and plants, suggesting the functional conservation of histone acetylation in transcriptional activation. Plant HATs were classified into different subclasses based on their sequence homology to yeast and animal HATs and their mode of action: (1) The Gcn5-related N-acetyltransferase (GNAT)/MYST (Moz, YBF2, Sas2p, Tip) family, (2) CREB-binding Protein (CBP) family, and (3) TBP-associated factor~II~ 250 (TAF~II~250) family \[[@B6-ijms-20-05311]\].

In plants, the genome-wide analysis of HATs has been performed in several species, including *Arabidopsis* \[[@B7-ijms-20-05311]\], rice \[[@B8-ijms-20-05311]\], *Vitis vinifera* \[[@B9-ijms-20-05311]\], and *Citrus sinensis* \[[@B10-ijms-20-05311]\]. HATs have been widely reported to play an important role in various aspects of plant development, including floral development \[[@B11-ijms-20-05311],[@B12-ijms-20-05311],[@B13-ijms-20-05311]\], root growth \[[@B14-ijms-20-05311]\], and gametophyte development \[[@B15-ijms-20-05311]\]. In addition to developmental functions, HATs are also involved in plant adaptation to various environmental fluctuations, such as salt stress \[[@B16-ijms-20-05311]\], cold stress \[[@B17-ijms-20-05311]\], heat stress \[[@B18-ijms-20-05311]\], light signaling \[[@B19-ijms-20-05311]\], abscisic acid (ABA) \[[@B20-ijms-20-05311]\], and other hormone signaling \[[@B21-ijms-20-05311]\]. Therefore, the understanding of HATs functions in field crops may play an important role in sustainable agriculture and food security.

Cotton, as a major source of natural and renewable textile fiber, holds a significant value in the world economy and in daily human life \[[@B22-ijms-20-05311]\]. *Gossypium hirsutum* is a natural allotetraploid (AADD) that arose from the interspecific hybridization between *Gossypium arboretum* (AA) and *Gossypium raimondii* (DD), which occurred approximately 1--2 million years ago \[[@B23-ijms-20-05311]\]. Because allopolyploid cotton produces a superior quality of fiber with high yield compared with their diploid progenitors \[[@B24-ijms-20-05311]\], *G. hirsutum* is widely grown in many parts of the world and contributes to more than 90% of commercial cotton production \[[@B25-ijms-20-05311]\]. On the one hand, world climate is changing quickly, and the abiotic stresses are severely affecting the cotton yield and fiber quality \[[@B26-ijms-20-05311],[@B27-ijms-20-05311]\]. On the other hand, the increasing world population demands the improvement of the cotton yield to meet the requirements of the textile industry. Therefore, identifying the potential genes conferring resistance to different stresses for the molecular breeding of cotton is of the utmost importance \[[@B28-ijms-20-05311]\]. Although the role of HATs has been investigated in the plant development \[[@B29-ijms-20-05311],[@B30-ijms-20-05311]\], response to environmental signals \[[@B16-ijms-20-05311],[@B18-ijms-20-05311],[@B31-ijms-20-05311]\], and hormone signaling \[[@B16-ijms-20-05311]\] in other plants, little is known regarding the function of HATs in cotton. In this scenario, *G. hirsutum*, *G, raimondii*, and *G. arboretum*, having recently available genomic data \[[@B32-ijms-20-05311],[@B33-ijms-20-05311],[@B34-ijms-20-05311],[@B35-ijms-20-05311]\], provide an excellent opportunity to identify the candidate genes involved in fiber development and biotic and abiotic stress tolerance and to expand our understanding of underlying epigenetic mechanisms.

In this study, we identified HAT genes from the whole genome of *G. hirsutum*, *G raimondii*, and *G. arboretum*. The identified HATs were comprehensively analyzed for phylogenetic classifications, gene structures, identification of conserved motifs and domain organization, and the presence of cis-regulatory elements in their promoters. Furthermore, gene expression profiles of *G. hirsutum HATs* were analyzed during the different stages of fiber development and in response to various abiotic stresses. Such a comprehensive analysis of HATs provides a fundamental understanding of their roles in cotton growth and development. Furthermore, this study will be useful for functional genomic studies on the regulations of histone acetylation and will eventually lead to the long-term improvement of stress tolerance in cotton.

2. Results {#sec2-ijms-20-05311}
==========

2.1. Identification of HATs in Cotton {#sec2dot1-ijms-20-05311}
-------------------------------------

A systematic blast search was performed to identify the HATs in the genomes of *G. hirsutum*, *G. raimondii* and *G. arboretum* with the query sequence of *Arabidopsis*, and candidate HAT were identified in the cotton genomes. Then, Pfam and InterProScan databases were used to further verify the candidate HAT, and a total of 18 *G. hirsutum* (GhHATs), 9 *G. raimondii* (GrHATs), and 9 *G. arboretum* (GaHATs) were identified ([Table S1](#app1-ijms-20-05311){ref-type="app"}). The properties of the identified *Gossypium* HATs were analyzed by ExPASy and we found an open reading frame (ORF) length ranging from 1407--6876 bp, which encoded the polypeptides ranging from 468--2291 amino acids with a predicted molecular weight of 53--256 KD. In addition, the theoretical isoelectronic point (PI) values ranged from 0.28--8.84 ([Table S1](#app1-ijms-20-05311){ref-type="app"}).

2.2. Phylogenetic Analysis of the HAT Gene Family {#sec2dot2-ijms-20-05311}
-------------------------------------------------

To investigate the evolutionary relationship of HATs among cotton (*G.raimondii, G. arboretum*, and *G. hirsutum*) and other species, we constructed a phylogenetic tree using MEGA 6.0. We used the newly identified *Gossypium* HATs and previously identified HATs from *Arabidopsis thaliana, Vitis vinifera*, *Oryza sativa*, and *Brachypodium distachyon* to confirm their evolutionary relationship with the un-rooted phylogenetic tree using the neighbor-joining (NJ) method ([Figure 1](#ijms-20-05311-f001){ref-type="fig"}). The phylogenetic tree results showed that similar to *Vitis* and *Arabidopsis*, *Gossypium* HATs could also be grouped into three distinct classes: CPB, GNAT and TAF~II~250. MYST homologs were found absent in the *Gossypium* genomes. To validate the phylogenetic tree constructed using the NJ method, we also used the minimum evolution method to construct a tree and found that HAT genes could be naturally classified into three groups and the members within each clade were stable with little difference between topology, which indicates that the NJ tree method could be used for further analysis.

We found that the number of genes in CPB (30) was greater than that in either GNAT (27) or TAF~II~250 (9). Furthermore, we found that all the three groups comprised mono and dicot species. It is noteworthy that the genes within each group clustered with a dicot- or monocot-specific pattern. The number of HAT from each species was different within each group. Compared to other species, cotton HATs showed a closer relationship with *Vitis* HATs because they always clustered closely to each other in the phylogenetic tree. Nevertheless, their gene numbers were not similar within the group.

2.3. Phylogenetic Tree Construction, Conserved Motif and Gene Structure Analysis of Cotton HATs Genes {#sec2dot3-ijms-20-05311}
-----------------------------------------------------------------------------------------------------

To confirm the subgroup classification and to determine the evolutionary relationship between *Gossypium* HATs, we generated another un-rooted phylogenetic tree using the neighbor-joining method ([Figure 2](#ijms-20-05311-f002){ref-type="fig"}A). We observed one subgroup related to CBP (HACs) and TAF~II~250 (HAFs) in *Gossypium*, whereas three subgroups related to GNAT (HAGs) were observed. The CBP and GNAT subgroups contained four members in *G. raimondii* and *G. arboretum*, while TAF~II~250 contained only one member. Furthermore, each diploid progenitor had its own ortholog in allotetraploid.

To gain insight into the evolutionary origin and putative functional diversification, a Multiple Expectation Maximization for Motif Elicitation (MEME) analysis was performed, which identified a total of 10 conserved motifs in *Gossypium* HATs ([Figure 2](#ijms-20-05311-f002){ref-type="fig"}B and [Figure S1](#app1-ijms-20-05311){ref-type="app"}). All members of the CPB subfamily contained all the conserved motifs, whereas GNAT and TAF~II~250 only contained a few conserved motifs. However, the motif 8 was found conserved in all the *Gossypium* HATs, suggesting that the potential catalytic site of *Gossypium* HATs was conserved. We further dissected the motif analysis of *Gossypium* HATs by investigating their domain organization ([Figure 2](#ijms-20-05311-f002){ref-type="fig"}C). The CBP subfamily of *Gossypium* HATs contained plant homeodomain (PHD), Znf-TAZ, Znf ZZ and a CBP-type HAT domain, while the GNAT subfamily of *Gossypium* HATs contained a conserved GNAT and bromodomain. Moreover, the TAF~II~250 subfamily of *Gossypium* HATs contained a bromodomain, ubiquitin (UBQ), TATA box--binding protein (TBP), and domain of unknown function (DUF) domain. Furthermore, an alignment information was produced to explore the amino acid conservation of the GhHATs domain sequence ([Figure S2](#app1-ijms-20-05311){ref-type="app"}). The multiple sequence alignment analysis revealed that all GhHATs proteins shared regions of conserved polypeptide sequences, which could be involved in their molecular functions ([Figure S2A--C](#app1-ijms-20-05311){ref-type="app"}). In general, *Gossypium* HATs contained a domain organization similar to that of their counterparts in other species.

We then analyzed the gene structure of HATs from *G. hirsutum, G raimondii* and *G. arboretum* ([Figure 2](#ijms-20-05311-f002){ref-type="fig"}D). Our results showed that the number and length of intron/exons were different among CPB, GNAT and TAF~II~250 classes. For example, intron/exon numbers and gene length of TAF~II~250 were greater than the CPB and GNAT. Furthermore, the gene structure of *G. raimondii* and *G. arboretum* was extremely similar to their orthologs in the *G. hirsutum*. However, in general, gene structure in terms of intron/exon was greatly similar within the subgroup, which was consistent with the phylogenetic analysis.

2.4. Chromosomal Distribution and Duplication Analysis of HATs {#sec2dot4-ijms-20-05311}
--------------------------------------------------------------

The *G. hirsutum* HATs were mapped to their corresponding chromosomes ([Figure 3](#ijms-20-05311-f003){ref-type="fig"}A), and all the *HATs* were unevenly distributed on the chromosomes of *G. hirsutum*. For example, all the 18 *G. hirsutum* HATs were assigned to 12 chromosomes out of 26 ([Figure 3](#ijms-20-05311-f003){ref-type="fig"}A). The CBP subfamily of *G. hirsutum* HATs was localized on the chromosomes 5, 6, 8 and 10, while the GNAT subfamily of *G. hirsutum* HATs was localized on the chromosomes 6, 7, 10 and 11. We also mapped the HATs from the diploids *G. raimondii* and *G. arboretum* and found that all the HATs from *G. raimondii* and *G. arboretum* were also unevenly distributed on their chromosomes, similar to *G. hirsutum* ([Figure 3](#ijms-20-05311-f003){ref-type="fig"}B). However, compared with *G. raimondii*, the HATs from *G. arboretum* were more evenly distributed on their chromosomes. We also observed that subfamilies of *Gossypium* HATs were localized to different chromosomes in both diploid progenitors (*G. raimondii* and *G. arboretum*), as well as in the allotetraploid *G. hirsutum*. For example, the TAF~II~250 subfamily member HAF1501 was localized on the chromosome 10 in *G. hirsutum*, and was localized on the chromosome 8 and 11 in *G. arboretum* and *G. raimondii*, respectively.

We also investigated the contribution of gene duplication to the expansion of *Gossypium* HATs ([Figure 3](#ijms-20-05311-f003){ref-type="fig"}A,B). Two segmental duplications were found in *G. hirsutum* (*GhHAC1501A/GhHAC1502A*, *GhHAG1501A/GhHAG1502A*), indicating that duplications occurred in the A genome of *G. hirsutum.* Similarly, two segmental duplications were also found in *G. raimondii* (*GrHAC1501/GrHAC1502*, *GrHAC1502/GrHAC1503*) and in *G. arboretum* (*GaHAC1501/GaHAC1502, GaHAG1501/GaHAG1502*). Furthermore, we also checked the Ka/Ks ratio to explore the selective constraints on each pair of duplicated *Gossypium HAT* ([Table S2](#app1-ijms-20-05311){ref-type="app"}). The Ka/Ks ratio was found to be less than 0.30 in all the duplicated gene pairs of *Gossypium* HATs, suggesting that all the duplicated gene pairs of *Gossypium* HATs experienced strong purifying selection pressure.

2.5. Putative Cis-Elements in the Promoter Regions of GhHATs {#sec2dot5-ijms-20-05311}
------------------------------------------------------------

To gain more insight into the putative functions of HATs, the putative cis-regulatory elements were scanned in the 1000 bp upstream of the transcription start sites of *G. hirsutum* using the Plant CARE database ([Figure S3, Tables S3 and S4](#app1-ijms-20-05311){ref-type="app"}). In addition to TATA- and CAAT-box core cis-elements, phytohormone response elements, stress response elements and development response elements were found in the promoters of *G. hirsutum* HATs. Most of the cis-elements were conserved among the CBP, GNAT and TAF~II~250 subfamilies of *G. hirsutum* HATs. However, some cis-elements were absent in some subfamilies. For example, CAT-box (cis-acting regulatory element related to meristem expression), MRE (MYB binding site involved in light responsiveness), P-box (Gibberellin-responsive element), and O2-sites (cis-acting regulatory element involved in zein metabolism regulation) cis-elements were absent in the CBP subfamily, but were found present in GNAT and TAF~II~250 members. Similarly, circadian (cis-acting regulatory element involved in circadian control), skn-1 motifs (cis-regulatory element required for endosperm expression), Box-W1 (fungal elicitor-responsive element), and 5′ UTR Py-rich stretch (a cis-regulatory element conferring high transcription levels) cis-elements were absent in the TAF~II~250 subfamily. Moreover, cis-elements in the promoters of A and D genomes were largely conserved in *G. hirsutum HATs*. However, TGA cis-elements (a cis-acting regulatory element involved in the MeJA responsiveness) were only present in the D genome. We also observed that cis-elements in the promoters of orthologous gene pairs of *G. hirsutum* were largely similar. However, some exceptions were also found where they differed much for the cis-elements, e.g., GhHAC1504-A/GhHAC1504-D and GhHAG1503-A/GhHAG1503-D.

2.6. Gene Expression Analysis {#sec2dot6-ijms-20-05311}
-----------------------------

### 2.6.1. Expression Analysis of *GhHATs* in Different Tissues, Developmental Stages and Multiple Abiotic Stresses by RNA-sequencing {#sec2dot6dot1-ijms-20-05311}

To better understand the potential physiological functions of GhHATs in allotetraploid cotton, we investigated the expression of *GhHAT* genes. RNA-sequencing data were downloaded from the National Center for Biotechnology Information (NCBI) and analyzed. Their analysis revealed that *GhHATs* were widely expressed in the vegetative (root, stem, and leaf) and reproductive (torus, petal, stamen, pistil, calyx, and −3, −1, 0, 1, 3, 5, 10, 20, 25 and 35 days post-anthesis (DPA) ovule) tissues ([Figure 4](#ijms-20-05311-f004){ref-type="fig"}A), highlighting the diverse biological functions of HATs in different tissues. We also noted that some *GhHATs* did not express in vegetative tissues, but had a very weak expression in the reproductive tissue, e.g., *GhHAG1502-D*. We found that A and D genomes showed preferential expression for only some genes in the leaf, root, and stem. For instance, the expression of *GhHAG1502-A* was higher than *GhHAG1502-D* in all the analyzed tissues. However, the opposite correlation was also observed, e.g., the expression of *GhHAG1504-D* and *GhHAG1503-D* was higher in the root than *GhHAG1504-A* and *GhHAG1503-A*. In addition, we also found that *GhHATs* expression was also differently regulated during the course of development ([Figure 4](#ijms-20-05311-f004){ref-type="fig"}B). For example, *GhHAC1501-D* and *GhHAC1502-A/D* expressions were increased with the development of the root, while the expression of *GhHAC1502-D* decreased with seed cotyledon development. We further investigated the gene expression pattern of *GhHATs* in different abiotic stresses. The expression of *GhHAG1501-A/GhHAG1501-D* and *GhHAG1504-A/GhHAG1504-D* was strongly induced by multiple stresses, indicating their potential involvement in stress responses. However, no clear changes in expression levels were observed for more than half of the *GhHATs* under cold, hot, salt, and polyethylene glycol (PEG) 6000 conditions.

### 2.6.2. Expression Pattern of *GhHATs* Genes by qPCR {#sec2dot6dot2-ijms-20-05311}

#### Tissue Specific Expression Patterns of HATs

We further validated the *GhHATs* gene expression in some vegetative and reproductive tissues by quantitative real-time polymerase chain reaction (qRT-PCR) ([Figure 5](#ijms-20-05311-f005){ref-type="fig"}A). Because the A and D genomes of allotetraploid cotton were extremely similar in mRNA levels, we considered *GhHAT-A* and *GhHAT-D* as one combination (*GhHAT*) and checked the expression levels by qRT-PCR. Among all the nine analyzed *GhHATs* belonging to three different classes, *GhHAG1501, GhHAG1502, GhHAC1503* and *GhHAG1504* showed the most prominent expression levels in the analyzed tissues, indicating their roles in the development of the leaf, root, stem and flowers. Although *GhHAG1502* and *GhHAC1503* were expressed in the root, stem and leaf, the highest expression was found in leaves. However, only *GhHAG1501, GhHAG1504* and *GhHAG1502* showed the highest expression in flowers, implying the specific function of these *GhHAGs* in flower development. It also indicated the dramatic functional divergence among different classes of GhHATs. Furthermore, *GhHAG1501* and *GhHAG1502* formed a paralogous pair and showed the expression conservation in all the analyzed tissues.

#### Expression of HATs at Different Fiber Development Stages

To explore the potential role of GhHATs in fiber development, we investigated their expression at different developmental stages of fiber development (0--25 DPA) by qRT-PCR ([Figure 5](#ijms-20-05311-f005){ref-type="fig"}B). Among all the *GhHATs*, only *GhHAG1501*, *GhHAG1502*, *GhHAC1503* and *GhHAF1501* showed the prominent expression levels during the different stages of fiber development ([Figure 5](#ijms-20-05311-f005){ref-type="fig"}B), indicating that these four genes may play dominant roles in the fiber development. However, among these four, the expression of *GhHAG1501* and *GhHAG1502* was the highest among the *GhHAGs* subgroup members at 0 DPA, and then gradually decreased during the fiber development (0--25 DPA). Similar to *GhHAGs* family members, the expression of *GhHAF1501* was also decreased during the fiber development, except at 20 DPA, where only the expression of *GhHAF1501, GhHAG1501* and *GhHAG1502* slightly increased and then again decreased at 25 DPA. Contrary to *GhHAF* and *GhHAG* family members, the expression of *GhHAC1503* from *GhHACs* family was gradually increased during the fiber development with a maximum at 20 DPA. Afterward, the expression of *GhHAC1503* was decreased at 25 DPA. We further separated the contribution of A and D genomes for the fiber development of *G. hirsutum* from RNA-seq data analysis and found that the A genome preferentially expressed more than D genome partners in the different fiber developmental stages ([Figure 4](#ijms-20-05311-f004){ref-type="fig"}A).

#### Expression of HATs in Response to Heavy Metals and Abiotic Stresses

To better understand the role of GhHATs in abiotic stresses, we investigated the gene expression of *GhHATs* in response to metal stress (Cd and Zn), salt stress (NaCl), cold stress and drought stress by qRT-PCR. All the abiotic stresses differentially regulated the expression of *GhHATs* ([Figure 6](#ijms-20-05311-f006){ref-type="fig"}), indicating the functional specificity of GhHATs in response to a particular stress. In response to Cd stress, the expression of *GhHAC1501, GhHAC1502* and *GhHAG1501* decreased, while the expression of *GhHAC1503, GhHAC1504, GhHAG1502* and *GhHAG1503* increased compared with the control. However, the expression of *GhHAG1501* and *GhHAG1504* did not change compared with the control. In response to Zn stress, only the expression levels of *GhHAC1503*, *GhHAC1504*, *GhHAG1501* and *GhHAG1502* increased, while all other studied genes did not show differential expression compared with the control. However, the expression levels of *GhHAC1503*, *GhHAC1504* and *GhHAG1502* were different in response to Cd and Zn stresses despite their increased expression. Furthermore, the expression of *GhHATs* was investigated in response to cold, salt, and drought stresses. The results showed that in response to cold stress, the expression of *GhHAC1502*, *GhHAC1503*, *GhHAC1504, GhHAG1501, GhHAG1502, GhHAG1503, GhHAG1504* and *GhHAF1501* decreased compared with the control. In response to salt stress, the expression of *GhHAC1501, GhHAC1502, GhHAC1503, GhHAG1501, GhHAG1504* and *GhHAF1501* increased, while the expression of *GhHAC1504* and *GhHAG1503* decreased compared with the control. In response to PEG treatment, the expression of *GhHAG1501* and *GhHAF1501* increased, while the expression of *GhHAC1503* and *GhHAC1504* decreased compared with the control.

#### Expression of HATs in Response to MMS

We also investigated the *Gossypium HATs* gene expression in response to the DNA damage agent methyl methanesulfonate (MMS) ([Figure 6](#ijms-20-05311-f006){ref-type="fig"}). The results showed that in response to MMS treatment, the expression of *GhHAC1502, GhHAC1503, GhHAG1501* and *GhHAF1501* increased, while the expression of *GhHAC1504* decreased compared with the control, indicating their potential role in DNA damage repair pathways in cotton.

#### Expression of HATs in Response to Phytohormones

Cis-regulatory elements related to phytohormones were found in the promoter of *GhHATs* ([Figure S2, Tables S3 and S4](#app1-ijms-20-05311){ref-type="app"}). We therefore investigated the gene expression of *GhHATs* in response to auxin (NAA) and abscisic acid (ABA) by qRT-PCR. Auxin and ABA treatments differentially regulated the expression of *GhHATs* ([Figure 6](#ijms-20-05311-f006){ref-type="fig"}). The results showed that in response to ABA treatment, the expression of *GhHAC1503, GhHAC1504* and *GhHAG1504* decreased, while the expression of *GhHAG1501* and *GhHAG1503* increased compared with the control. In response to NAA treatment, the expression of *GhHAC1502, GhHAG1501* and *GhHAF1501* increased, while the expression of *GhHAC1504* decreased compared with the control.

2.7. Co-Localization of HATs with QTLs of Fiber Development {#sec2dot7-ijms-20-05311}
-----------------------------------------------------------

To validate the potential function of GhHATs in fiber development, the co-localization of *GhHATs* with reported QTLs/SNPs of fiber development (i.e., fiber length (FL), fiber elongation (FE), fiber micronaire (FM), fiber strength (FS), and fiber uniformity (FU) was analyzed. Nine genes were mapped to eight chromosomes with reported QTLs of FL, FE, FS, FM, and FU ([Figure 7](#ijms-20-05311-f007){ref-type="fig"}). Among nine *GhHATs*, only four genes were co-localized with QTLs of FL, FE, FM, FS, and FU on four chromosomes, i.e., Chr-A05, Chr-A06, Chr-A08 and Chr-A11 of the A subgenome. *GhHAC1502* was located within the qFE-A05-2, *GhHAG1503* gene anchored in qFU-A06-1 and qFE-A06-2, and *GhHAG1504* was mapped in the qFU-A05-2 QTL region, while *GhHAC1504* was 2 Mb from the qFL-A08. Among the D subgenome *HATs*, *GhHAC1502*, *GhHAC1503* and *GhHAG1503* were anchored in the FL-QTL-9, SNP (i20058Gh, i38606Gh), and qFL-D06, respectively, while *GhHAF1501* and *GhHAG1504* were found in surroundings of reported SNP. Interestingly, some genes were co-localized with multiple QTLs related to different fiber development traits. For example, the *GhHAC1502* gene on chromosome A05 and D05 was anchored in FL and FE QTL, and the *GhHAG1503* gene was localized inside the FL, FE and FU QTL. This reveals that these *GhHATs* had pleiotropic effects on fiber development related traits.

3. Discussion {#sec3-ijms-20-05311}
=============

Histone acetylation is a mark of transcriptional activation and has been reported to play an important role in plant development and response to various biotic and abiotic stresses \[[@B12-ijms-20-05311],[@B14-ijms-20-05311],[@B15-ijms-20-05311],[@B16-ijms-20-05311],[@B18-ijms-20-05311]\]. Moreover, levels of histone acetylation are tightly linked with gene expression regulation. HATs carry out histone acetylation and have been classified into different distinct groups in *Arabidopsis* \[[@B7-ijms-20-05311]\], rice \[[@B8-ijms-20-05311]\], *Vitis vinifera* \[[@B9-ijms-20-05311]\], and *Citrus sinensis* \[[@B10-ijms-20-05311]\]. In this study, we revealed that *Gossypium* HATs could also be classified into three major subgroups: CBP, GNAT, and TAF~II~250. This suggests that multiple subgroups of *Gossypium* HATs might play specialized roles in the adaptive evolution of cotton. However, similar to rice \[[@B8-ijms-20-05311]\], the MYST domain containing homologs of HATs in all the three genomes of cotton were absent, while *Arabidopsis* \[[@B7-ijms-20-05311]\] and *Vitis* \[[@B9-ijms-20-05311]\] had MYST domain members. This indicates the specific functional divergence of HATs between the different field crop plants. CBP, GNAT, and MYST HATs are also considered as transcriptional co-activators in addition to their HAT activity. For example, TAZ-type, ZZ-type and PHD-type zinc finger domains of CBP have been reported to play an important role in protein recognition and protein-protein interactions \[[@B36-ijms-20-05311],[@B37-ijms-20-05311]\]. The PHD domain has also been reported to interact with histones and other histone-related proteins \[[@B38-ijms-20-05311]\]. Furthermore, bromodomains are known to bind to acetylated lysine residues \[[@B39-ijms-20-05311],[@B40-ijms-20-05311]\]. Along with functional catalytic domains, all the other conserved domains of *Gossypium* GNAT, CBP and TAF~II~250 were largely similar to their counterparts in monocots and dicots ([Figure 2](#ijms-20-05311-f002){ref-type="fig"}C). These observations suggest that all the *Gossypium* HATs may have similar functions as described in other plant species.

Gene duplication plays a significant role in generating new gene subfamilies in the evolution of genome and genetic systems \[[@B41-ijms-20-05311]\]. Tandem duplication, polyploidy, and segmental duplications primarily contribute to the creation of new gene families \[[@B41-ijms-20-05311]\]. Two segmental duplications were found in *G. hirsutum* (*GhHAC1501A/GhHAC1502A, GhHAG1501A/GhHAG1502A*), *G. raimondii* (*GrHAC1501/GrHAC1502, GrHAC1502/GrHAC1503*), and *G. arboretum* (*GaHAC1501/GaHAC1502, GaHAG1501/GaHAG1502*) ([Figure 3](#ijms-20-05311-f003){ref-type="fig"}). Gene duplications occurred in these genes because the identities of the genes flanking both sides of the paralogous *Gossypium HAT* genes were found to be absolutely conserved and located on duplicated segments on two different chromosomes. Moreover, these duplicated genes were not likely diverged much during the evolution based on their Ka/Ks ratios ([Table S2](#app1-ijms-20-05311){ref-type="app"}), suggesting the functional conservation of duplicated genes. This observation can partially be validated by the overlapping expression patterns of *GhHAG1501* and *GhHAG1502* during the fiber development and in different tissues ([Figure 5](#ijms-20-05311-f005){ref-type="fig"}A,B). However, the expression of *GhHAG1501* did not change in response to Cd, while the expression of *GhHAG1502* increased compared with control ([Figure 6](#ijms-20-05311-f006){ref-type="fig"}), suggesting that these duplicated genes may undergo functional divergence in response to particular stimuli.

*Arabidopsis* HATs play a crucial role in different aspects of plant development \[[@B12-ijms-20-05311],[@B14-ijms-20-05311],[@B15-ijms-20-05311],[@B16-ijms-20-05311],[@B18-ijms-20-05311]\]. *G. hirsutum* has a huge contribution in the textile industry. Therefore, the identification of *HATs* genes in *G. hirsutum* and their role in fiber development will provide the fundamental information for future studies. Our results showed that the expression of *GhHAC1503* (CBP subgroup) was increased during the different stages of fiber development, while the expression of *GhHAG1501/GhHAG1502* (GNAT subgroup) and *GhHAF1501* (TAF~II~250 subgroup) was generally decreased ([Figure 5](#ijms-20-05311-f005){ref-type="fig"}B). This suggests that histone acetylation levels are likely to be dynamic during the course of fiber development. However, further studies are required to investigate the function of histone acetylation in cotton fiber development. Phytohormones, including gibberellic acid (GA) \[[@B42-ijms-20-05311]\], jasmonic acid (JA) \[[@B43-ijms-20-05311]\], abscisic acid (ABA) \[[@B44-ijms-20-05311]\], ethylene \[[@B45-ijms-20-05311]\], and auxin (NAA) \[[@B46-ijms-20-05311]\], are also known to regulate the fiber development. Cis-regulatory elements specific for ethylene (ERE), GA (GARE), ABA (ABRE), and JA (CGTCA) were found in the promoter of four highly expressed *GhHATs* (*GhHAC1503, GhHAG1501, GhHAG1502* and *GhHAF1501*) during the fiber development ([Tables S3 and S4, Figure S3](#app1-ijms-20-05311){ref-type="app"}). This suggests that the expression of these *GhHATs* might be regulated by phytohormones. Consistent with this, the expression of *GhHAG1501/GhHAG1502* slightly increased in response to ABA treatment, while the expression of *GhHAC1503* decreased ([Figure 6](#ijms-20-05311-f006){ref-type="fig"}). The expression of *GhHAF1501* did not change in response to ABA treatment, but the expression was increased in response to NAA treatment ([Figure 6](#ijms-20-05311-f006){ref-type="fig"}). This suggests that these phytohormones directly or indirectly regulate the expression of these *GhHATs*, which may lead to different acetylation levels. The exogenous application of auxin has been shown to cause higher acetylation levels at the promoter of *SKP2B* \[[@B47-ijms-20-05311]\], while the acetylation levels decrease at the promoter of *KRP7* in *Arabidopsis* \[[@B48-ijms-20-05311]\], indicating the crosstalk between phytohormones and histone acetylation in plants. However, further studies are required to investigate the crosstalk between histone acetylation and phytohormones in cotton. Together, our results suggest that histone acetylation might play an important role in fiber development, as well as cotton growth. However, further studies are required to investigate the function of each GhHATs in phytohormone related pathways.

Histone acetylation has been reported to play a key role in DNA damage repair \[[@B49-ijms-20-05311],[@B50-ijms-20-05311],[@B51-ijms-20-05311]\]. Histone acetyltransferase 1 (HAT1) is required for the incorporation of H4K5/H4K12 acetylated H3.3 histones at the sites of double strand breaks (DSB), which then facilitate the recruitment of a key DNA repair factor Rad51 in mammals \[[@B52-ijms-20-05311]\]. Furthermore, DSB-inducing agents have been reported to induce the H4K16ac levels in mammals \[[@B53-ijms-20-05311]\], indicating a key role of histone acetylation in DNA damage and repair pathways. In *Arabidopsis*, the histone acetyltransferases HAM1, HAM2 and HAG3 have been reported to participate in UV-B induced DNA damage \[[@B54-ijms-20-05311],[@B55-ijms-20-05311]\]. Furthermore, a treatment with the histone acetyltransferase inhibitor curcumin increased DNA damage in *Arabidopsis* and maize \[[@B54-ijms-20-05311]\], indicating the conserved role of histone acetylation in DNA damage in plants. Our qRT-PCR results showed that the expression of *GhHAC1502, GhHAC1503, GhHAG1501* and *GhHAF1501* increased, while the expression of *GhHAC1504* decreased compared with the control in response to MMS treatment ([Figure 6](#ijms-20-05311-f006){ref-type="fig"}). In brief, our results suggest the potential implication of GhHATs in DNA damage repair pathways in cotton. However, further studies are required to investigate whether and how GhHATs are involved in DNA damage repair.

Different abiotic stresses, including, cold, drought, salt, and metal stresses, severely affect the cotton growth and yield. Recent studies established the link of histone acetylation and abiotic stresses \[[@B32-ijms-20-05311],[@B56-ijms-20-05311],[@B57-ijms-20-05311],[@B58-ijms-20-05311]\]. For example, salt treatment caused a global increase of H3K9 and H3K4 acetylations \[[@B32-ijms-20-05311]\], while cold treatment decreased the H3K9, H4K5 and H4K4 acetylations compared with the control in maize \[[@B58-ijms-20-05311]\]. Similarly, in response to dehydration treatment, drought responsive genes such as *RD29A*, *RD29B*, *RD20* and *RAP2.4* were found to be differentially acetylated at H3K9, H3K14, H3K23 and H3K27 in *A. thaliana* \[[@B56-ijms-20-05311]\]. We showed that many *GhHATs* were differentially regulated in response to different abiotic stresses, suggesting the dynamic levels of histone acetylation in each abiotic stress adaptation in cotton. Our expression data also suggest the functional diversity and specificity among different subgroups of *Gossypium* HATs (GhHACs, GhHAGs and GhHAFs) in response to stress, as well as hormone treatments ([Figure 6](#ijms-20-05311-f006){ref-type="fig"}). For example, the expression of *GhHAC1502* strongly increased in response to NaCl, while the expression of *GhHAC1503* was strongly induced by Zn. Similarly, in the HAGs subgroup, the expression of *GhHAG1501* was strongly induced by NaCl, while the expression of *GhHAG1502* did not change. These observations suggest that GhHATs likely play an important role in cotton plant adaptation to various abiotic stresses. Furthermore, we also noticed that more than one gene responded to a particular stress. For example, in response to salt stress, the expression of *GhHAC1501, GhHAC1502, GhHAC1503, GhHAG1501, GhHAG1504* and *GhHAF1501* increased, while the expression of *GhHAC1504* and *GhHAG1503* decreased compared to the control. This suggests that different GhHATs may work together in response to particular stimuli and may participate in long-term resistance to different abiotic stresses. However, further molecular and biochemical studies are required to validate GhHATs function and to understand the underlying molecular mechanism.

Four *HAT* genes (*GhHAC1502*, *GhHAC1503*, *GhHAG1503* and *GhHAG1504)* on six chromosomes were anchored in fiber development-related QTL regions. These four genes were identified inside multiple QTLs, suggesting the pleiotropic role of these genes in fiber related traits. These results also support the conclusion that both the A and D genomes of *Gossypium hirsutum* jointly participate in different aspects of cotton fiber trait. Interestingly, *GhHAC1502*, *GhHAC1503*, *GhHAG1503* and *GhHAG1504* also displayed differential expression at different fiber development stages and in response to abiotic stresses, metal stress, MMS and phytohormones. This suggests that these four *HAT* genes are potentially involved in cotton growth and development, fiber-related traits, and plant response to the environment. However, in general, there were few discrepancies between RNA-seq and qRT-PCR data. RNA-seq expression levels were either from the A or D genome of *G. hirsutum* and it would be difficult to properly separate and count the reads from the homologous regions. On the contrary, we regarded *GhHACxA* and *GhHACxD* as one combination, referred to as *GhHACx*, because of the extremely high similarity between the mRNAs of the *GhHACxA-GhHACxD* gene pairs (for example, *GhHAC1502A/GhHAC1502D* and *GhHAG1501A/GhHAG1501D* have 99 and 99.5% similarity, respectively) and their nearly identical transcript sizes. With such a high similarity, we could not distinguish them using the qRT-PCR, and the results of qRT-PCR are theoretically the average of A and D genome's expression. Thus, potential difficulties of read counts on homologous regions in RNA-seq and inability to separate the A and G genome expressions by qRT-PCR may cause the discrepancy between RNA-seq and qRT-PCR results. In brief, further studies are required to comprehensively elucidate the function of HATs in different aspects of cotton plant.

In this study, we highlighted complex and diverse transcriptional regulations of *GhHATs* which significantly broadened our understanding of the underlying epigenetic mechanisms in allotetraploid cotton and unraveled an extra layer of complexity for better allotetraploid cotton adaptation in response to developmental, environmental, and hormonal cues. Furthermore, our results strongly recommend the comprehensive dissection of the biological and cellular function of GhHATS and argue for the potential implication of histone acetyltransferases in cotton molecular breeding in addition to other existing breeding strategies. This will eventually lead to the long-term improvement of stress tolerance in cotton and will allow high fiber yield and quality.

4. Materials and Methods {#sec4-ijms-20-05311}
========================

4.1. Identification of HATs Gene Family {#sec4dot1-ijms-20-05311}
---------------------------------------

The data of three cotton species, *G. raimondii* (JGI, version), *G. arboretum* (BJI, version 1.0), and *G. hirsutum* (NAU, version 1.1) were attained from the COTTONGEN (<http://www.cottongen.org>) \[[@B33-ijms-20-05311],[@B34-ijms-20-05311],[@B35-ijms-20-05311]\]. The HAT protein sequences from *Arabidopsis*, rice, *Brachypodium distachyon*, and *Vitis Vinfera* were downloaded from Phytozome (<http://phytozome.jgi.doe.gov/pz/portal.html>) and then used as queries in BLASTP searches \[[@B59-ijms-20-05311]\] against the *G. raimondii, G. arboretum* and *G. hirsutum* genomes, respectively. Genes with *E*-values \< 1.0 were selected, and redundant sequences were eliminated by following the previously published method \[[@B60-ijms-20-05311]\]. Furthermore, InterProScan (<http://www.ebi.ac.uk/interpro/search/sequence-search>) was used to confirm the presence of the HAT domain \[[@B61-ijms-20-05311]\]. The physicochemical properties were predicted by ExPASy (<http://cn.expasy.org/tools>).

4.2. Analysis of Chromosomal Location and Gene Duplication {#sec4dot2-ijms-20-05311}
----------------------------------------------------------

The loci of *HATs* were deduced from the gff3-files of cotton genome. The localization of *HATs* genes on the chromosomes were visualized using the program Circos \[[@B62-ijms-20-05311]\]. The duplication events of *HATs* and Ka/Ks were calculated using the previously published method \[[@B63-ijms-20-05311]\]. The T = Ks/2λ equation was used to determine the duplication time and deviation of the *HAT* gene pairs, assuming clock-like rates of (λ) 1.5 × 10^−8^ substitutions per synonymous site per year for cotton \[[@B64-ijms-20-05311]\].

4.3. Sequence Alignment and Phylogenetic Analyses {#sec4dot3-ijms-20-05311}
-------------------------------------------------

Multiple sequence alignment was performed for the full-length HAT proteins using Clustal W with standard settings \[[@B65-ijms-20-05311]\]. A neighbor-joining (NJ) phylogenetic tree was constructed using the full length HATs sequences from *G. raimondii, G. arboretum* and *G. hirsutum* by MEGA 6.0 \[[@B66-ijms-20-05311]\], with P-distance and pairwise gap deletion parameters engaged. The bootstrap test was used with 1000 replicates to evaluate the statistical consistency of each node. To confirm the grades from the NJ method, the minimal-evolution method of MEGA 6.0 was utilized with 1000 replicates as well.

4.4. Gene Structure, Protein Motif, and Promoter Cis-Element Analysis {#sec4dot4-ijms-20-05311}
---------------------------------------------------------------------

The exon/intron structures of the *HAT* genes were acquired from bed-file and displayed using the online tool Gene Structure Display Server (<http://gsds.cbi.pku.edu.cn>) \[[@B67-ijms-20-05311]\]. The NJ tree was constructed with MEGA 6.0 as explained above. The deduced HAT protein sequences of three cotton species were submitted to the online Multiple Expectation Maximization for Motif Elicitation (MEME) version 4.11.1 (<http://meme-suite.org/tools/meme>) \[[@B68-ijms-20-05311]\] as described \[[@B60-ijms-20-05311]\]. For cis-element analysis in promoter regions, the 1 kb upstream sequences were analyzed in the PlantCARE database (<http://bioinformatics.psb.ugent.be/webtools/plantcare/html/>) \[[@B69-ijms-20-05311]\].

4.5. Transcriptome Data Analysis and Gene Expression Heatmap {#sec4dot5-ijms-20-05311}
------------------------------------------------------------

The raw data of RNA-seq of *G. hirsutum* were downloaded from the NCBI Sequence Read Archive (SRA: PRJNA248163) to calculate the expression level. TopHat (version: 2.0.13) was used for mapping reads, cufflinks (version: 2.2.1) were used to analyze gene expression levels, and fragments per kilobase million values were used to normalize gene expression levels \[[@B70-ijms-20-05311]\]. Finally, these values of the *GhHAT* candidates were extracted from total expression data and the heatmap was generated by MeV 4.0 (<http://www.tm4.org/>).

4.6. Plant Materials, Stress Treatments, and qRT-PCR {#sec4dot6-ijms-20-05311}
----------------------------------------------------

*G. hirsutum* cultivar 'CRI35' was used for gene expression analysis. All the sampled tissues obtained from cotton plants grown under field condition with standardized cultural practices to determine the expression analysis \[[@B60-ijms-20-05311]\]. For treatments, cotton seeds were surface sterilized and germinated on a moist paper. Young seedlings of same size were selected and exposed to NaCl (200 mM), polyethylene glycol 4000 (PEG4000) (15%), cold (4 °C), methyl methanesulfonate (MMS) (250 ppm), auxin (NAA) (10 µM), abscisic acid (ABA) (10 µM), ZnSO~4~ (1 mM), and CdCl~2~ (1 mM) for 24 h. All treatments were performed in three biological replicates. All samples were frozen quickly in liquid nitrogen and kept at −80 °C. The total RNA was extracted from cotton samples using the RNAprep Pure Plant kit (TIANGEN, Beijing, China). A total of 2 μg of RNA was used as the template, and the first-strand cDNAs were synthesized using the SuperScript III (Invitrogen, Waltham, MA, USA). Quantitative real-time PCR (qRT-PCR) analysis was performed as described previously in \[[@B71-ijms-20-05311]\] using the specific primers for each *GhHAT* gene ([Table S5](#app1-ijms-20-05311){ref-type="app"}). Cotton *UBQ7* (UniProt accession number: AY189972) was used as an internal reference gene for normalization of expression and three biological replicates were performed for each sample. To calculate the relative expression levels, a comparative 2^−∆∆*C*t^ method was used \[[@B72-ijms-20-05311]\]. The heat map for the gene expression profiles was generated with Mev 4.0 (<http://www.tm4.org/>) \[[@B73-ijms-20-05311]\].

4.7. Co-localization of HATs with Fiber Related QTLs {#sec4dot7-ijms-20-05311}
----------------------------------------------------

To identify the localization of QTLs and SNPs for fiber development related traits, QTLs and linked molecular markers were retrieved from the Cotton Gen website (<https://www.cottongen.org>). The sequence of each marker was fetched from Cotton Gen to obtain the physical position information. For this purpose, the sequence of each marker was BLAST against the *G. hirsutum* (AD1) and HAU genome in the CottonGFD database. HAT genes co-localized with QTLs were displayed to show *HAT* gene distribution on chromosomes, along with surrounding loci and QTLs, using mapchart software. Genes identified inside the QTL or ≤500 kb far from SNP were considered as anchored gene in QTL because cotton LD decay was approximately 0.80 Mb \[[@B74-ijms-20-05311]\].

4.8. Statistical Analysis {#sec4dot8-ijms-20-05311}
-------------------------

After performing normal distribution of the data and the homogeneity of variance tests, analysis of variance (ANOVA) was performed, followed by the least significant difference (LSD) test at *p* value ≤ 0.05 for each parameter. Different letters indicate significant difference by LSD test (*p* ≤ 0.05). Statistical analyses were performed using the Statistical Package for Social Sciences (SPSS) software (version 11.5, SPSS Inc., Chicago, IL, USA).

5. Conclusions {#sec5-ijms-20-05311}
==============

In this study, 36 HATs were identified in three genomes of *Gossypium* and clustered into three groups: CPB, GNAT, and TAF~II~250. *Gossypium HATs* are unevenly distributed on the chromosomes, and segmental duplications contributed to the evolution of the *HATs* family. Cis-element analysis discovered several abiotic and biotic stresses and hormonal responsive elements in the promoter region of the *GhHATs*, but each member had peculiar types and numbers. Furthermore, the expression profile of the upland cotton *HAT* gene family exhibited different expression patterns in response to abiotic and hormonal stresses, which disclosed that GhHATs play roles in different aspects of upland cotton abiotic stress tolerance and hormonal signaling. Thus, our study helps to lay the foundation for the functional characterization of the *GhHATs* gene family by overexpression and knockdown/out using RNAi or CRISPR-Cas9 genome editing and provides new insight into the evolution and divergence of *HAT* genes in plants. Furthermore, these results may enhance the understanding of the molecular mechanisms of many agronomic traits of cotton, such as fiber development and other physiological processes.

We would like to thank Yasar Sajjad (Department of Biotechnology, COMSATS University Islamabad, Abbottabad campus, Pakistan) for his valuable suggestions.

###### 

Click here for additional data file.

Supplementary Materials can be found at <https://www.mdpi.com/1422-0067/20/21/5311/s1>.

Conceptualization, S.S.; Formal analysis, M.I., S.S., M.A.F., M.K.N., E.W. and A.B.; Funding acquisition, S.S. and R.R.-C.W.; Methodology, M.I.; Writing---original draft, M.I. and S.S.; Writing---review & editing, S.S., K.B.J. and R.R.-C.W.

This research was supported by the grant PD-IPFP/HRD/HEC/2013/1129 from Higher Education Commission of Pakistan.

The authors declare no conflict of interest.

![Phylogenetic relationships of histone acetyltransferases (HATs) from *Gossypium hirsutum*, *Gossypium raimondii*, *Gossypium arboretum*, *Arabidopsis thaliana*, *Oryza sativa*, *Brachypodium distachyon* and *Vitis vinifera*. The un-rooted phylogenetic tree was constructed using MEGA 6 by the neighbor-joining (NJ) method, and the bootstrap analysis was performed with 1000 replicates. For the phylogenetic tree, amino acid sequences were used and the classification of CPB, GNAT and TAF~II~250 was performed based on the conserved signature domain of each subgroup. The HATs from *G. hirsutum*, *G. raimondii*, *G. arboretum*, *Arabidopsis*, *V. vinifera*, *Brachypodium distachyon*, and *Oryza sativa* were marked with black, green, red, pink, cyan blue, yellow, and blue dots, respectively.](ijms-20-05311-g001){#ijms-20-05311-f001}

![Phylogenetic relationships, conserved motifs, domain organization, and gene structure analysis of HATs in *Gossypium hirsutum*, *Gossypium raimondii* and *Gossypium arboretum*. (**A**) The unrooted phylogenetic tree was constructed using MEGA 6 by the NJ method, and the bootstrap analysis was performed with 1000 replicates. The HAT proteins from *G. hirsutum*, *G. raimondii*, and *G. arboretum* were marked with black, green and red dots, respectively. The branches of each subgroup were indicated in a specific color. (**B**) Motif identification analysis in HAT proteins from *G. hirsutum*, *G. raimondii* and *G. arboretum*. Each motif is shown in unique color, and the regular expression and sequences of the 1--10 motifs are listed in [Figure S1](#app1-ijms-20-05311){ref-type="app"}. (**C**) Domain organization of HAT proteins from *G. hirsutum*, *G. raimondii* and *G. arboretum*. One representative of each subgroup of HAT from *G. hirsutum* is presented. (**D**) The intron/exon structure of *HAT* genes from *G. hirsutum*, *G. raimondii* and *G. arboretum*. The green boxes and grey lines represent exons and introns, respectively.](ijms-20-05311-g002){#ijms-20-05311-f002}

![Chromosomal distribution and gene duplication of HATs in *Gossypium hirsutum* (**A**) and in *Gossypium raimondii* and *Gossypium arboretum* (**B**). The chromosome number was indicated in boxes and represented as Gh1-Gh13A/D (**A**), Ga1-Ga13, and Gr1-Gr13 (**B**) for *G. hirsutum*, *G. raimondii*, and *G. arboretum*, respectively. The orthologous HATs are connected by green, while the segment duplication of HATs is represented by light blue and red colors in different genomes.](ijms-20-05311-g003){#ijms-20-05311-f003}

![Gene expression pattern of *Gossypium hirsutum HATs* by the analysis of RNA-sequencing in different tissues and at different stages of fiber development (**A**) and in response to cold, hot, salt, and polyethylene glycol (PEG) 6000 stresses (**B**). The illumina reads of RNA-seq data were retrieved from the National Center for Biotechnology Information Sequence Read Archive (NCBI SRA) database. The color scale at the bottom of heat map indicates the fragments per kilobase million (FPKM)-normalized log2 transformed counts.](ijms-20-05311-g004){#ijms-20-05311-f004}

![Gene expression validation of *HATs* from *G. hirsutum* by quantitative real-time polymerase chain reaction (qRT-PCR) in different tissues. (**A**) RNA from the root, leaf, stem and flowers was extracted and reverse transcribed. *Ubiquitin 7* (*UBQ7*) was used as an internal control for qRT-PCR. The relative expression is presented in the heat map, and the color scale at the top represents the relative signal intensity. The primer sequences can be found in [Table S5](#app1-ijms-20-05311){ref-type="app"}. (**B**) Gene expression validation of *HATs* from *G. hirsutum* by qRT-PCR at different stages of fiber development. The data were normalized to *UBQ7*. The data presented are the average of three biological replicates. Bar = standard deviation (SD).](ijms-20-05311-g005){#ijms-20-05311-f005}

![Expression pattern of *Gossypium hirsutum HATs* in response to abiotic stresses and hormones. RNA from roots after 24 h of CK (control), cadmium (Cd), zinc (Zn), abscisic acid (ABA), auxin (NAA), DNA damage (MMS), cold, salt (NaCl), and drought (PEG4000) was extracted and reverse transcribed. *UBQ7* was used as an internal control for qRT-PCR. The data presented are the average of three biological replicates. Different letters indicate significant difference by least significant difference (LSD) test (*p* ≤ 0.05). Bar = SD.](ijms-20-05311-g006){#ijms-20-05311-f006}

![Distribution of co-localized *HAT* genes on chromosomes of A and D subgenomes of *G. hirsutum*. The scale represents the physical position of genes and quantitative trait loci (QTL)-linked markers in megabases (Mb). QTLs/Single Nucleotide Polymorphism (SNPs) related to fiber length (FL), fiber elongation (FE), fiber micronaire (FM), fiber strength (FS), and fiber uniformity (FU) are shown. Asterisks indicate that the *HAT* genes co-localized with QTLs/SNPs related to fiber.](ijms-20-05311-g007){#ijms-20-05311-f007}
