To identify PHD gene family members from wheat, whole genome data for T. aestivum (IWGSC RefSeq_v1.1) were obtained from the Ensembl plant database (http://plants.ensembl.org/info/website/ftp/index.html), and the PHD-finger domain (PF00628) was downloaded from the PFAM database (https://pfam.xfam.org/). The PHD protein sequences from A. thaliana (70) and O. sativa (59) ( Supplementary Table S1 ) (Sun et al., 2017) were used as query sequences to search against the wheat protein dataset using the BLASTP program, and the threshold was set as E-value < 1e-5. The NCBI-Batch CD-Search (Marchler-Bauer et al., 2017) (https://www.ncbi.nlm.nih.gov/Structure/bwrpsb/bwrpsb.cgi), PFAM database, and SMART database (http://smart.embl.de/) were used to further confirm the candidate PHD-finger genes of T. aestivum. There were other spliced transcripts in the candidate genes of these species, and we selected the first splice variant as a representative for subsequent analysis.

The protein sequences of TaPHDs were computed using the ExPASy server (Artimo et al., 2012) to obtain the theoretical isoelectric point (pI), molecular weight (MW), instability index (II), aliphatic index (AI), and grand average hydrophobicity (GRAVY). Plant-mPLoc (Chou and Shen, 2010) (http://www.csbio.sjtu.edu.cn/cgi-bin/PlantmPLoc.cgi) and BUSCA (Savojardo et al., 2018) (Bologna Unified Subcellular Component Annotator, http://busca.biocomp.unibo.it) were used to predict the subcellular localization of the TaPHD proteins.

The PHD-finger protein sequences of T. aestivum, A. thaliana, and O. sativa were used for phylogenetic analysis. Jalview 2.11 software (http://www.jalview.org/) with the MUSCLE method with default parameters was utilized to conduct multiple sequence alignment. Evolutionary analysis involved 342 amino acid sequences (all wheat PHD genes, and most rice and Arabidopsis PHD genes). These analyses were conducted in MEGA X (Kumar et al., 2018) using the neighbor-joining method (Saitou and Nei, 1987). The percentage of replicate trees in which the associated taxa clustered together in the bootstrap test (1000 replicates) is shown next to the branches. The evolutionary distances were computed using the Poisson correction method and were expressed as the number of amino acid substitutions per site. The iTOL website (http://itol.embl.de/) was used to visualize the phylogenetic tree.

MCScanX software (Wang et al., 2012) was used to detect collinear regions between TaPHD genes as well as collinear blocks of TaPHDs with three monocotyledons (H. vulgareto, Z. mays, and O. sativa) and three dicotyledons (A. thaliana, B. rapa, and G. raimondii). Whole genome data for H. vulgareto, Z. mays, O. sativa, A. thaliana, B. rapa, and G. raimondii were obtained from the Ensembl plant database (http://plants.ensembl.org/info/website/ftp/index.html). All TaPHD genes were mapped to their respective loci in the wheat genome in a circular diagram using shinyCircos (Yu et al., 2018). Gene duplication events of TaPHDs and synteny relationships between the aforementioned species were visualized using TBtools (v1.082) (Chen et al., 2020). The Ka/Ks values (non-synonymous substitution rate/synonymous substitution rate) were calculated after identification of duplicated genes, using the method of Nei and Gojobori as implemented in KaKs_calculator (Zhang et al., 2006) based on the coding sequence alignments. Subsequently, the divergence time of collinear gene pairs was calculated using the duplication events formula T = Ks/(2λ × 10^-6) in millions of years (Mya), with λ = 6.5 × 10^-9 (Wang et al., 2015b).

GO annotation of TaPHD proteins was available from the KOBAS database (http://kobas.cbi.pku.edu.cn/kobas3) (Xie et al., 2011). The full-length amino acid sequences of TaPHD proteins were uploaded to the original program, followed by drawing and annotation. GO annotations were performed for three types of analyses: biological processes, molecular functions, and cellular composition. The GO annotation results were visualized using the online tool OmicStudio (https://www.omicstudio.cn/tool) (Ye et al., 2018). All the predicted TaPHD proteins were submitted to the STRING database (https://string-db.org/cgi/input.pl). The minimum required interaction score was set to a high confidence (0.700). The maximum number of interactors was no more than 10 on the first shell.

Transcriptional data for TaPHDs were obtained from the wheat expression website (http://www.wheat-expression.com/download) (Borrill et al., 2016; Ramírez-González et al., 2018) and were used to explore the potential biological functions of TaPHD genes in growth and development, abiotic and biotic stress, and other conditions. Systematic clustering analysis was performed based on the log2 of transcripts per million (TPM) values for the 244 TaPHD genes. R was used to display the expression patterns in a heat map, and OmicStudio (https://www.omicstudio.cn/tool) was used to display the histogram, volcano plot, and Venn diagram.

In this study, the seeds of the hexaploid common wheat variety “Zhengmai 7698” were surface-sterilized with 2% hydrogen peroxide, rinsed thoroughly with distilled water, and germinated with water saturation at 25°C for two days in Petri dishes on three layers of filter paper. The young seedlings were transformed and grown in 1/2 Hoagland’s culture solution under a 14 h light (25°C)/10 h dark (20°C) photoperiod. When the wheat grew to two leaves and one heart, the plants were subsequently treated with 16% polyethylene glycol 6000. For cold stress, wheat seedlings were exposed to 4°C for 12 h. For heat stress, wheat seedlings were exposed to 40°C for 12 h. New leaves of the three seedlings were collected as biological replicates, and each treatment had three replicates.

Total RNA was extracted using RNAiso Reagent (TaKaRa, Beijing, China) and Cdna was synthesized using the RT Master Mix Perfect RealTime kit (TaKaRa, Beijing, China). Quantitative real-time PCR was performed using the CFX Touch™ Real-Time PCR Detection System (Bio-Rad Laboratories, Hercules, CA, USA) and the SG Fast Qpcr Master Mix (Sangon Biotech, Shanghai, China). Relative expression levels were determined using the 2^(-ΔΔCt) method (Livak and Schmittgen, 2001), and β-actin was used as the internal control to normalize the expression levels of TaPHD genes. Specific primers used for qRT-PCR are listed in Supplementary Table S2 .

Full-length open reading frames of TaPHD11, TaPHD19, and TaPHD133 were obtained from “Zhengmai 7698” Cdna ( Supplementary Table S2 ). The Coding sequence (CDS)of TaPHD11, TaPHD19, and TaPHD133 were cloned into the pJIT16318 vector at the BamHI site using specific primers ( Supplementary Table S2 ). The pJIT16318 vector contained a CaMV 35S promoter and C-terminal GFP. Transient expression assays were conducted as described by Cui et al. (2019). Approximately 4 × 10⁴ mesophyll protoplasts were isolated from 12-day-old wheat seedlings. The transfected protoplasts were incubated at 23°C for 12 h. The GFP fluorescence in the transformed protoplasts was imaged using a confocal laser-scanning microscope (LSM 700; Zeiss).

In this study, 244 T. aestivum genes were designated PHD genes with two query methods; HMM and BLASTP were used for identification, and three websites, NCBI-Batch CD-Search, PFAM database, and SMART database, were used for confirmation ( Supplementary Table S3 ). These PHD genes were renamed TaPHD1 to TaPHD244, based on the order of their chromosomal locations and physical positions.

To further determine the characteristics of TaPHD genes, the ExPASy Server online tool was used to analyze the protein characteristics ( Supplementary Table S3 ). The shortest protein contained 216 amino acids (TaPHD158, TaPHD175) and the longest protein contained 2853 amino acids (TaPHD204); the molecular weight was between 24567.82 Da (TaPHD158) and 310347.53 Da (TaPHD204). The protein instability index showed that all PHD genes were unstable proteins. The isoelectric point of TaPHD genes varied markedly from 4.42 (TaPHD36) to 9.65 (TaPHD78), and the aliphatic index varied significantly from 48.13 (TaPHD26/39/51) to 97.51 (TaPHD42). The GRAVY of TaPHD proteins in wheat varied from 0.016 (TaPHD160) to -1.285 (TaPHD23), indicating that they were all hydrophilic proteins, except for TaPHD160 ( Supplementary Table S3 ). We used two methods (Plant-mPLoc and BUSCA) to predict the subcellular localization of the TaPHD proteins. The results showed that a few TaPHDs may be localized in the chloroplast, mitochondrion, or cytoplasm, and most members were predicted to be located in the nucleus ( Supplementary Table S3 ).

Multiple sequence alignments of PHD domains were performed ( Figure 1 ). Approximately 60 amino acids (aa) comprised a PHD domain containing basic Cys4-His-Cys3 sequence motifs in each TaPHD.

To evaluate the evolutionary relationships of PHD genes in T. aestivum, O. sativa, and A. thaliana, a neighbor-joining phylogenetic tree was constructed using full-length PHD proteins ( Figure 2 and Supplementary Table S1 ). Phylogenetic analysis showed that PHD family proteins can be divided into four clades (clades 1 to 4). TaPHD members were found in all clades. Clade 1 was the largest, with 95 TaPHD members, and clade 4 was the smallest, with only 38 members. The results showed that there were many small branches under each clade, and almost every small branch had corresponding genes of rice and Arabidopsis. This indicates that the TaPHD gene is not an evolutionary characteristic of monocotyledonous and dicotyledonous plants, and that the PHD gene family was formed before the differentiation of these two types of plants.

Protein domains are often functional carriers. According to phylogenetic and domain analyses (NCBI-Batch CD-Search, PFAM, and SMART database), 30 dominant types were identified in all wheat PHD proteins ( Table 1 ). The results showed that among all wheat PHD proteins, 43 contained a typical PHD domain. The next most common, the jas-PHD and alifn-PHD domains, had 28 and 25 members, respectively; the PHD-Oberon_cc domain and the PHD-RING domains had 11 members, and the remaining domain types had less than ten members. The results showed that wheat PHD proteins contained a canonical PHD domain or double PHD domains. Owing to their different domains, differentiation in function was achieved.

To better understand why PHD-finger genes are abundant in the wheat genome, we analyzed the homoeologous groups in detail ( Table 2 ). A total of 35.8% of wheat genes were present in homoeologous groups of three, also termed triads (A:B:D = 1:1:1) (Consortium et al., 2018). In contrast, 84.8% of the PHD-finger genes identified were present in triads ( Table 2 ). Also, the percentage of PHD-finger genes with homoeolog-specific duplications was lower for PHD-finger genes than for all wheat genes (1.6% vs 5.7%; Table 2 ). Loss of one homoeolog, on the other hand, was less pronounced in PHD-finger genes (6.6% vs 13.2%; Table 2 ). Only four PHD-finger genes were orphans/singletons. Thus, the high homoeolog retention rate could partly explain the high number of wheat PHD-finger genes.

Based on the reference GFF3 files, the physical positions of PHD genes on the corresponding chromosomes are shown in Figure 3 . The identified TaPHDs could be mapped on every chromosome and evenly across the three sub-genomes. The map shows that chromosomes 5B and 5D harbor the largest number of TaPHD genes (18), whereas chromosome 1D contains the least (6).

Gene duplication is an indispensable mechanism by which organisms create new genes with similar or different functions (Song et al., 2019). Therefore, we analyzed the duplication events that occurred in the TaPHD gene family. A total of 230 PHD gene pairs from wheat were identified as duplicated ( Figure 4 and Supplementary Table S4 ). These similar PHD gene pairs had the same domain type and appeared in the same branch of the phylogenetic tree. Tandem and segment duplications are critical for the evolution of gene families to adapt to different environmental conditions. Interestingly, all the TaPHD gene pairs were associated with segmental duplication events. This suggests that this was the main route for expanding PHD genes in wheat and the many homologous genes on different wheat chromosomes suggest the high conservation of the family. To further infer the evolutionary origin and homology of the wheat PHD family, we constructed a collinear chart comparing six species with wheat, including three monocotyledons (H. vulgareto, Z. mays, and O. sativa) and three dicotyledons (A. thaliana, B. rapa, and G. raimondii) ( Figure 5 and Supplementary Table S4 ). We identified pairwise homologues of the TaPHD genes and detected 119, 186, 168, 7, 2, and 6 pairs of homologous genes from H. vulgareto, Z. mays, O. sativa, A. thaliana, B. rapa, and G. raimondii, respectively ( Figure 5 and Supplementary Table S4 ). This implies that TaPHD genes share a strong evolutionary relationship with ZmPHDs, HvPHDs, and OsPHDs. Furthermore, these results indicated that the PHD gene family was differentiated between monocotyledonous and dicotyledonous plants. This also indicated that TaPHD genes had a strong evolutionary relationship with ZmPHDs, HvPHDs, and OsPHDs. The average differentiation time was: barley (12.78 Mya) < rice (22.09 Mya) < maize (60.87 Mya).

Ka/Ks, the non-synonymous substitution ratio, determines the selection pressure for duplicated genes. According to the results ( Supplementary Table S4 ), only a very few TaPHD gene pairs had Ka/Ks ratios >1, suggesting that the evolution of TaPHD genes was accompanied by strong purifying selection. The Ka/Ks ratios between wheat and three monocotyledonous plants were calculated based on the collinear gene pairs. Except for very few genes, the values of the other collinear gene pairs were all below 1, which confirmed that the evolution of the wheat PHD gene family underwent strong purifying selection. However, the Ka/Ks ratios of the collinear gene pairs between wheat and the three dicots could not be calculated properly. This is because most synonymous mutation sites have synonymous mutations; that is, the degree of sequence divergence and evolutionary distance is too large. Some TaPHD genes have formed at least five homologous gene pairs, such as TaPHD9, which may have played key roles in the evolution of the PHD gene family ( Figure 5 and Supplementary Table S4 ).

We performed GO annotation analysis of the 244 TaPHD proteins, revealing that they may participate in a range of cellular components, molecular functions, and biological processes ( Figure 6 and Supplementary Table S5 ). The 244 TaPHD proteins were assigned a total of 105 GO terms. In biological processes, the three most highly enriched categories were related to the regulation of DNA-templated transcription, heat acclimation, and chromatin organization. Developmental growth and jasmonic acid-mediated systemic resistance were also particularly enriched. In the cellular component category, the most highly enriched categories were related to the nucleus, and 85% of the TaPHDs could participate in this process, whereas less than 10% of TaPHDs were involved in plasmodesma. Regarding molecular functions, the 65 most enriched TaPHDs were involved in histone binding, 28 TaPHDs were involved in chromatin binding, and 81 TaPHDs were related to protein binding.

To understand protein-protein interactions between TaPHDs and other proteins in wheat, we constructed a protein-protein interaction network ( Figure 7 and Supplementary Table S6 ). A total of 89 TaPHD proteins and 548 interacting protein branches were identified. According to the strength of the interaction, we divided the 89 proteins into four interaction regions, which are represented by different colors, as shown in Figure 7 . Some TaPHDs, such as TaPHD15, TaPHD145, and TaPHD162, could interact with up to 28 proteins, suggesting that these TaPHD proteins play a significant role in the regulation of protein networks. Notably, we found that these proteins had a PHD domain or a PHD-SWIB-Plus3-GYF domain. Therefore, we believe that such domains are likely to play an important role in the PHD family.

RNA-sequencing is a powerful tool for exploring certain gene transcription patterns using high-throughput sequencing methods (Wang et al., 2009). Systematic clustering analysis was performed based on the log2 of TPM values for 244 TaPHD genes ( Figure 8A and Supplementary Table S7 ). The data showed that TaPHD gene expression showed great differences with the change in the growth period. In general, the expression of TaPHDs can be divided into three categories: the first group contains members that are widely expressed in many tissues under multiple developmental stage conditions; the second group contains those that are highly induced only at specific growth and development stages; and the last group includes members that do not appear to be expressed during growth and development. For example, TaPHD100, TaPHD108, and TaPHD122 had high expression during most growth and developmental processes, except in the endosperm. There were also some genes (TaPHD222 and TaPHD232) that had higher expression only in shoots and roots. Furthermore, some genes, such as TaPHD68, TaPHD78, and TaPHD86, were not expressed, which implies that these genes may have functional redundancy.

To further study the expression differences of this family in different stages and organs of wheat, we counted the number of high, medium, and low expression genes in each period and organ ( Figure 8B ). The data showed that the number of highly expressed genes was the largest in the stigma and ovary, reaching as high as 60, followed by a spike in the boot period, reaching 41. The lowest number of highly expressed genes (none) was found in the flag leaf blade at night in the flag leaf stage. Our results suggest that some TaPHDs may play important roles in many biological processes during wheat growth, especially during anthesis.

The differential expression of TaPHDs under different conditions is shown in Figures 9A–F and Supplementary Table S8 . During biological stress, we found that inoculation with Fusarium, powdery mildew, pathogen associated molecular patterns (PAMP), crown rot, Septoria, or stripe rust caused few changes in the expression of TaPHD genes. This suggests that TaPHD family members may not be associated with disease resistance.

Under abiotic stress, there are many TaPHD genes whose expression changes are more obvious under high-temperature, drought, and cold conditions ( Figures 9G–K and Supplementary Table S8 ). For example, after high-temperature treatment, the expression levels of many TaPHD genes (TaPHD26, TaPHD75, TaPHD100, TaPHD115, TaPHD117, and TaPHD167) were significantly altered compared to those in the experimental control group. In the drought starvation treatment, TaPHD11, TaPHD19, TaPHD99, TaPHD141, TaPHD153, and TaPHD171 expression levels changed significantly. However, in the phosphorus starvation treatment, there were few changes in the expression of TaPHD genes. To further understand whether there is an intersection between the differential genes of the PHD family under drought, high-temperature, and low-temperature treatments, we drew a Venn diagram of DEGs in TaPHD genes during the four different transcriptomes ( Figure 10 , Supplementary Table S9 ). The data showed that TaPHD215 and TaPHD223 were significantly altered in every treatment. TaPHD30, TaPHD96, TaPHD180, TaPHD174, and TaPHD239 gene expression varied greatly between the two drought and heat treatments. In addition, in cold and heat stress environments, the expression levels of five genes (TaPHD109, TaPHD118, TaPHD120, TaPHD167, and TaPHD178) were significantly changed.

To elucidate the possible regulatory mechanisms of TaPHD genes under cold, drought, and heat conditions, we performed qRT-PCR analysis of 20 genes ( Figure 11 ). The results showed that all 20 TaPHDs responded to different stress conditions and had different manifestations. Under low temperature stress induced by 4°C, the expression of five TaPHDs was significantly upregulated at different time points, and the expression of six TaPHDs was significantly downregulated at different time points compared with the control. In contrast, under 40°C-induced high-temperature stress, the expression of 12 TaPHDs was significantly upregulated at different time points compared with the control. The expression of five TaPHDs was inhibited at different time points. This indicated that compared with low temperature stress, high temperature stress could induce more changes in the expression of TaPHDs and could upregulate the expression more. In wheat under 16% PEG stress, the expression of ten TaPHDs was significantly upregulated at different time points. The expression of seven TaPHDs was inhibited at different time points. Among them, TaPHD72 was most significantly inhibited, and it was downregulated four-fold at 6 and 12 h after treatment. The expression levels of TaPHD69 and TaPHD135 significantly increased after the three treatments. However, the expression levels of TaPHD23 and TaPHD141 significantly decreased after the three treatments. In addition, TaPHD99 was strongly upregulated or downregulated by high temperature, low temperature, and PEG, and we speculated that this might be a key regulator of abiotic induction. In conclusion, we verified the effect of PHD-finger gene expression on the effect of three abiotic stresses in wheat using qRT-PCR. These results indicate that PHD-finger genes play an important role in coping with abiotic stress in wheat.

Previous studies have shown that most PHD finger proteins are localized in the nucleus, and only a few are localized in the membranes or other organelles (Gozani et al., 2003; Wu et al., 2016; Sun et al., 2017). For example, ZmPHD14 and ZmPHD19 are localized to the nucleus (Wang et al., 2015a). Also, GmPHD1 to GmPHD6 target the nucleus, and their nuclear localization requires the PHD domain (Wei et al., 2009). To better understand the functions of TaPHDs, we used Plant-mPLoc and BUSCA to predict their subcellular localization. The results showed that more than 90% of the TaPHD proteins were localized in the nucleus ( Table S1 ). In Arabidopsis thaliana, the PHD genes AL5 and AL6 play a very important role in improving the resistance of plants to abiotic stress. Therefore, we selected TaPHD11 and TaPHD19, which are highly homologous to AtALs, for subcellular localization of wheat protoplasts. As shown in Figure 12 , this suggests that, in wheat, the proteins TaPHD11 and TaPHD19 not only function in the nucleus but also in the membrane. In addition, research has shown that PHD finger ING2 is a phosphoinositide binding module and a nuclear PtdInsP receptor and suggests that PHD-phosphoinositide interactions directly regulate nuclear responses to DNA damage (Gozani et al., 2003). However, we studied the protein TaPHD133, which is highly homologous to ING1, and found that it is localized not only in the nucleus but also in the membrane. In summary, the subcellular localization of PHD proteins in wheat differs from that in other species.