Gene structure and conserved motifs analysis
In order to better understand the similarity and diversity of different members of the PHD protein in G. hirsutum, we analyzed the phylogenetic tree, exon-intron structure, and conserved motif (Fig. 3). GhPHD genes in the branches with closer homology have similar gene structure and exon-intron number, and their corresponding PHD proteins have highly conserved motif, which is consistent with the conclusion in the phylogenetic tree. Among them, GhPHD49 has the longest genomic sequence and 26 exons, while GhPHD12 has the shortest genomic sequence and only has two exons (Fig. 3B and Table S3). A total of three conserved motifs were identified, named motif 1 to motif 3, respectively (Fig. 3C). Among them, 20 GhPHD proteins all have motif 1, 2, and 3, and 75 GhPHD proteins only have motif 1. Different GhPHD proteins have different number of motifs, but GhPHD proteins with the same clustering relationship have strong motif conservation.
The alignment of 108 GhPHD conserved domain was shown in Fig. S1. The GhPHD motifs have nine conserved sites (indicated with dark blue and pink) and they all have typical Cys4-His-Cys3 model. Conserved histidine (H) was separated by four amino acids from the fourth conserved cysteine (C), and was separated by two amino acids from the latter conserved cysteine C. The third and fourth conserved cysteines (C) before histidine (H) were separated by one or two amino acids, but the interval number between other conserved amino acids was uncertain. However, GhPHD17, GhPHD27, GhPHD71, and GhPHD81 were closely related, but they did not follow the rules about conserved domains mentioned above (Fig. S1 and Fig. 3).
Cis -acting element analysis
Many studies have showed that PHD genes were involved in various stress responses. To elucidate the putative function of GhPHD under different stresses, we identified putative stress-related and phytohormone-related cis-acting elements in the 1500 bp promoter sequence before the start transcription codon (ATG) of 108 GhPHD genes (Fig. 4 and Table S4). Ten kinds of plant hormone-related elements, ABRE (abscisic acid), TGA-element (auxin), AuxRR-core (auxin), ERE (ethylene), TATC-box (gibberellin), GARE-motif (gibberellin), P-box (gibberellin), CARE (gibberellin), CGTCA-motif (MeJA), TCA-element (salicylic acid); and three kinds of stress-responsive regulatory elements, TC-rich repeats (defense and stress responsiveness), MBS (drought), LTR (low-temperature) were identified in the GhPHD promoters. In addition, these GhPHD promoters contain circadian control elements and light-responsive elements (G-box). These results revealed that all GhPHD promoters contain at least one putative stress-responsive element.
In addition to GhPHD52 and GhPHD10, other GhPHD promoters contained elements that respond to single or multiple phytohormones (Fig. 4 and Table S5). There were 85, 67, 64, and 55 GhPHD promoters with ethylene, ABA, GA, and MeJA response elements, respectively. In particular, the 17 GhPHD promoters contain the response elements of four hormones including ethylene, ABA, GA, and MeJA. Further, there were 26 and 37 GhPHD promoters containing auxin and SA response elements, respectively. The promoter region of 27 GhPHD genes contained cis-acting elements in response to four hormones, and the promoter regions of eight GhPHDs contained cis-acting elements in response to five hormones. In particular, GhPHD5, GhPHD47, GhPHD56, and GhPHD65 contained cis-acting elements in response to six hormones, indicating that these genes might be involved in multiple hormone signal pathways. Interestingly, cis-acting element analysis showed that GhPHD may be involved in light response and abiotic stress pathways, because there were 71, 38, and 36 GhPHD promoters with light, drought, and low temperature response elements, respectively. Additionally, some GhPHD genes promoters contain both low temperature and drought response elements. Coupled with these results, we speculated that GhPHDs might play important roles in light response, abiotic stress and hormone signaling pathways.
Tissue-specific expression patterns of GhPHD genes in cotton tissues
To predict the possible functions of GhPHD genes during cotton development, we investigated their expression profiles in 20 different tissues including root, stem, leaf, petal, stamen, pistil, ovule (-3, -1, 0, 1, 3, 5, 10, 20, 25 and 35 days ovule), and fiber (5, 10, 20, 25 DPA) using available transcriptomic data. According to the expression features (Fig. 5), 108 GhPHD genes were mainly clustered into four groups (A-D) based on a hierarchical clustering analysis. Among these, nine GhPHD genes in group A were highly expressed in each tissue and preferentially expressed in ovules, especially GhPHD23 and GhPHD77. The expression of 43 GhPHDs in group B were slightly lower than that of GhPHDs in group A. Groups C and D each contained 28 members and their expression were very low in all the tissues observed. In group D, three GhPHDs (GhPHD6, GhPHD10, and GhPHD63) were not expressed in root, stem, leaf, petal, stamen, and pistil. GhPHD6 was expressed only in 3 DPA ovule, while GhPHD103 was expressed only in root and fiber. GhPHD7 and GhPHD87 were not expressed in fibers, but they are low in other tissues. Different expression patterns of each GhPHD genes revealed that they have different physiological functions in cotton.
Identification of stress-related PHD genes in G. hirsutum
Cis-acting elements analysis presumed that GhPHDs may be involved in the stress responses. To verify this hypothesis, we analyzed the expression profile of GhPHD genes under heat, cold, salt, and drought treatments using the transcriptomic data. Under different stress treatments, each gene responds differently to stress. However, 66 GhPHD genes were highly expressed under stress treatment (Fig. S2). To identify GhPHD genes that play important roles in the stress responses, we verified the expression level of 12 GhPHD genes using qRT-PCR. Under the treatment of heat, cold, salt, and drought (Fig. 6), the relative expression level of GhPHD18 were up-regulated, indicating that GhPHD18 responded positively to all four stress treatments. GhPHD23 was also up-regulated after heat treatment, but was down-regulated after cold, salt, and drought treatment, as the expression level of GhPHD23 decreased dramatically after six hours of drought treatment (Fig. 6). Similarly, under heat and cold treatment, the expression patterns of the six GhPHD genes (GhPHD34, 40, 43, 80, 88, and 107) were consistent with that of GhPHD23, that is, their gene expression levels increased by 2.0-2.5 times under heat treatment, but decreased under cold treatment (Fig. 6).
Under salt treatment, GhPHD5, 24, 34, 40, 43, and 72 were up-regulated, while GhPHD80 and GhPHD88 were down-regulated. Under drought treatment, the relative expression levels of GhPHD5, 18, 80, and 88 increased, while GhPHD23, 40, 43, and 107 decreased (Fig. 6). The relative expression of GhPHD77 increased gradually at different times of heat treatment, but its expression level did not increase significantly under cold, salt, and drought treatments (Fig. 6). These results indicate that different GhPHDs respond differently to different abiotic stresses. These seven GhPHD (23, 34, 40, 43, 77, 80, and 107) showed positive responses to heat stress, while GhPHD23, GhPHD77, and GhPHD80 showed negative responses to cold, salt, and drought.
Identification of the GhPHD genes in response to phytohormone
It is an established fact that phytohormone act as chemical messengers in the physiological processes during the plant’s life cycle, which simultaneously coordinate physiological responses to biotic and abiotic stresses. To further explore the potential functions of GhPHDs in response to stress-related plant hormones, we analyzed the expression pattern of six GhPHDs under GA, MeJA, IAA, SA, and BL treatment by qRT-PCR (Fig. 7). Among 12 GhPHDs mentioned above in response to abiotic stresses (Fig. 6), six GhPHD containing four or more hormone response elements (Fig. 4) in the promoter region were selected and their responses to different hormones were identified.
The GhPHD5 promoter contains responsive elements for six hormones (auxin, ethylene, GA, MeJA, SA, and ABA) (Fig. 4). However, qRT-PCR results showed that the relative expression level of GhPHD5 increased after treatment with MeJA, IAA, and BL. The expression level of GhPHD5 increased at 0.5 hour of SA treatment, while decreased significantly at 1, 3, and 5 hours of SA treatment. However, at different time points of GA treatment, the expression level of GhPHD5 fluctuated, but the difference was not significant. Similarly, after 0.5 hour of GA treatment, the expression level of GhPHD40 increased to 1.2, and then gradually decreased with time. After SA treatment, the expression level of GhPHD40 increased at different time points. However, its expression decreased at 0.5 hour of BL treatment, and then increased to 1.45 at one hour after BL treatment (Fig. 7).
Under the treatment of five hormones at different time points, the relative expression level of GhPHD43 gradually increased, which indicates that GhPHD43 has a strong response to BL and the weakest response to MeJA (Fig. 7). In contrast to GhPHD43, the expression level of GhPHD80 decreased after treatment with five hormones (Fig. 7). The relative expression level of GhPHD80 increased significantly to 2.9 at 0.5 hour after GA treatment, and then rapidly decreased to 0.2 at three hours of GA treatment. Interestingly, the response of GhPHD88 to GA treatment was the same as that of GhPHD80, that is, the relative expression level of GhPHD88 increased to 2.3 at 0.5 hour of GA treatment, and reached a minimum of 0.3 after three hours of GA treatment (Fig. 7). Under SA treatment, the relative expression level of GhPHD88 increased gradually, while under BL treatment, the relative expression level decreased to 0.3 at 0.5 hour, but increased to 1.7 at one hour, and decreased at three and six hours. Under GA and IAA treatments, the relative expression of GhPHD107 increased and reached the maximum at one hour. The relative expression of GhPHD107 decreased significantly at 0.5 hour of BL treatment, reached the maximum at one hour, and then decreased significantly at three hours and six hours after treatment (Fig. 7). The above results indicated that GhPHDs showed diverse expression pattern under different hormone treatment and might play crucial roles in hormone signaling pathways.
Co-expression networks with functional modules for G. hirsutum and G. arboreum
GO ontology analysis of all 108 GhPHDs genes showed that protein binding and zinc ion binding were the most abundant functional terms compared to the regulation of transcription, histone binding, and chromatin binding (Fig. 8B). We then used ccNET software to analyze the co-expression networks of the genes identified above for stress and phytohormone response. After analyzing them one by one, we predicted many co-expressed genes and interacting proteins (Table S6). For example, GhPHD5 depicted a positive co-expression relationship with SLOMO (slow motion), which is a F-box protein required for auxin homeostasis and normal timing of lateral organ initiation at the shoot meristem [39]. GhPHD5 showed positive co-expression with a novel plant-specific DNA ligase, which is involved in seed germination and DNA repair (Table S6).
GhPHD18 may interact with highly hydrophilic proteins that are involved in positively regulating FLC expression. GhPHD18 have positive co-expression with a SHAGGY-related kinase involved in meristem organization, indicating that GhPHD18 play important roles in plant growth and development (Table S6). GhPHD34 was negatively co-expression with ERF subfamily B-1 of the ERF/AP2 transcription factor and participated in the ethylene signal pathway and abiotic stress response (Table S6). Similarly, GhPHD107 exhibited positive co-expression with ARF-GAP and ERF (ethylene response factor), suggesting that GhPHD107 may be involved in auxin and ethylene signaling pathways (Table S6). More interestingly, the physiological function of GhPHD88 is particularly powerful, and the genes predicted to interact with GhPHD88 are leucine-rich repeat protein kinase, late embryogenesis abundant protein (LEA), AP2/B3 transcription factor, R2R3 factor gene, DREB subfamily A-2 of ERF/AP2 transcription factor, cellulose synthase, gibberellin-regulated family protein, ethylene response factor, and other genes (Fig. 8 and Table S6)), which indicates that GhPHD88 may be involved in many physiological processes such as plant growth and development, hormone signal transduction and stress response. These results showed that GhPHDs were involved in the plant growth and development, and phytohormone may mediate the response of the GhPHD to abiotic stress, thereby increasing the plant’s tolerance to adverse environmental stress.