D. Huoshanense materials were collected at the Plant Cell Engineering Center of West Anhui University (Luan, China). The culture conditions were 25°C ± 2°C, with a 12 h day and 12 h dark cycle (Zhang et al., 2023). MS medium (without hormone addition) was used for the tissue culture of seedlings. When the seedlings grow to a height of 5–8 cm in subculture, hormone treatments with different concentrations are added. A 100 mol/L ABA solution was supplied to the medium, and the expression levels of Dof genes were measured on days 0, 2, and 4. The samples were treated with a 50 mol/L MeJA solution and collected at 0, 2, and 8 days. The expression level of the Dof genes was measured on days 0, 2, and 4 after adding the 0.1 mg/L IAA solution and the 100 mol/L SA to the medium. After the application of the 0.5 mmol/L GA solution to the plants, samples were collected on days 0, 3, and 6 to measure the expression of Dof genes. During the low-temperature treatment, the samples were first stored at a temperature of 4°C in a refrigerator. To measure the expression of Dof genes, samples were obtained on days 0, 1, and 4.

The genome sequence and annotation of D. huoshanense were obtained from the National Center for Biotechnology Information database (accession: PRJNA597621). The Hidden Markov Model (HMM) of the Dof domain (PF02701) is downloaded from the InterPro database. The Arabidopsis Dof homologs were downloaded from the Arabidopsis Information Resource (https://www.arabidopsis.org/), and the Blastp method was used for sequence alignment to obtain non-redundant Dof genes. Then, Pfam (http://pfam-legacy.xfam.org/), InterPro (http://www.ebi.ac.uk/interpro/), and SMART (https://smart.embl.de/) were applied for the verification of the Dof candidates. Using genome and gff annotations, chromosomal location analysis of 22 DhDof genes was completed, and then TBtools software was used to perform visual analysis of Dof genes (Chen et al., 2020). By using the Expasy web service (http://expasy.org/), it was possible to predict the molecular weight, isoelectric points, instability index, aliphatic index, and hydropathicity. The Plant-mPLoc plug-in of Cell-PLoc 2.0 software (http://www.csbio.sjtu.edu.cn/bioinf/Cell-PLoc-2/) was used for the subcellular localization prediction.

To determine the composition and classification of DhDof, a phylogenetic tree of the two species was constructed using the Arabidopsis nomenclature. MEGA (v.6.0.6) software was used to construct the evolutionary tree using the neighbor-joining method. First, ClustalW software was used to perform protein sequence alignment with the Poisson model and pairwise deletion mode. The bootstrap test was set at 1,000 replicates. IQ-TREE (v.1.6.12) software was used to construct the evolutionary tree using the maximum likelihood method. First, IQ-TREE was used to calculate the optimal model of these Dof sequence alignment files. “VT + F + R5” was chosen as the best-fit model under the Bayesian information criterion. Then, using this model and the “iqtree.exe -s./bidui.fas -m VT + F + R5 -bb 1000 -alrt 1000 -nt AUTO” command, we constructed an unrooted consistency evolutionary tree. The two evolutionary trees passed the identity (>40%) and SH-aLRT support (>80%)/ultrafast bootstrap support (>95%) thresholds to be manually classified.

The conserved domain database (https://www.ncbi.nlm.nih.gov/Structure/bwrpsb/bwrpsb.cgi) was used for domain searches of Dof proteins. The MEME-suite web service (http://meme.nbcr.net/meme/cgi-bin/meme.cgi) was applied for motif searching and the result was visualized using TBtools software (https://github.com/CJ-Chen/TBtools). Motifs identified by MEME were further retrieved in the InterPro database (http://www.ebi.ac.uk/interpro/). The gene structure was performed using Gene Structure Display Server (http://gsds.cbi.pku.edu.cn/) and visualized by TBtools.

ClustalW was used to align DhDof protein sequences first, and then the alignment files were imported into GeneDoc software (http://psc.edu/biomed/genedoc) for visualization of protein sequences. Black background, gray background, and white background represent different matching values of amino acids. WebLogo 3 software (https://weblogo.threeplusone.com/create.cgi) was used to visualize the sequence logo of the Dof domain.

The MCScanX software (https://github.com/wyp1125/MCScanX) was applied to examine collinearity between D. huoshanense and the other four species. The genome sequences and annotations of D. nobile, D. chrysotoxum, O. sativa and A. thaliana are downloaded from NCBI (accession: PRJNA725550), NCBI (accession: PRJNA664445), Ensembl Plants (genome assembly: IRGSP-1.0), and TAIR (genome assembly: TAIR10). The microsynteny circus of the gene duplication event was visualized by TBtools software. Then, each Dof sequence of D. huoshanense, D. nobile, D. chrysotoxum, A. thaliana, and O. sativa was aligned against itself using BLASTp with an E-value threshold of e⁻¹⁰ to obtain the syntenic blocks between each two species.

The 2 kilobase (kb) regions flanking the promoters were obtained from the genomic sequences. The cis-acting elements in Dof promoters were obtained from the PlantCARE database (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/). TBtools software was utilized to visualize the basic elements associated with growth, phytohormones, and environmental responses.

The previous transcriptome datasets (https://bigd.big.ac.cn/gsa/browse/CRA006607) were used to compute and normalize the expression of Dof genes as fragments per kilobase per million mapped fragments (FPKM). The sample reads were aligned to the D. huoshanense reference genome using the HISAT software, which could be available at https://daehwankimlab.github.io/hisat2/. The FPKM values for mRNA expression analysis were computed using StringTie (https://github.com/gpertea/stringtie). Salmon was employed to quantify expression levels from RNA-seq data (https://combine-lab.github.io/salmon/). The heatmap construction utilizes the normalized FPKM values (log₂ (FPKM)) of the unigenes.

RNA was isolated and extracted using the RNA extraction reagent (QIAGEN, Germany) according to the manufacturer’s protocol. The genomic DNA (gDNA) was cleaved using the RNAse-Free DNase Set from QIAGEN, a company based in Germany. The RNA samples were assessed for purity using the NanoDrop 2000c (Thermo Fisher Scientific, United States) and gel electrophoresis. The reverse transcription PCR was performed using the PrimeScriptTM II 1st Strand cDNA Synthesis Kit from TaKaRa, a company based in China. The SYBR Premix Ex TaqTM II, manufactured by TaKaRa in China, was employed for quantitative reverse transcription polymerase chain reaction (qRT-PCR) on a 7,500 series real-time fluorescence quantitative cycler manufactured by Bio-RAD in the United States of America. The primers for the qRT-PCR test were prepared using the Primer Premier 5.0 software. The ACTIN gene served as the reference gene (Song et al., 2021). The primers utilized for qRT-PCR analysis are listed in Supplementary Table S2. Each experiment was repeated three times in triplicate, and a total of three biological replicates were undertaken. The gene expression levels were determined using the 2^-△△CT method.

A total of 22 Dof genes were screened from D. huoshanense genome, and their conserved domains were further compared and verified through Pfam, InterPro, and SMART databases. Using the genome and annotation files, the chromosomal locations of these Dof genes were determined. The original Dofs were renamed based on their order and location on different chromosomes. The results showed that 22 Dof genes were dispersed on 12 chromosomes (Figure 1). The Dof genes anchored on chromosome 8 genetically formed a cluster at a physical distance. However, there was only one Dof gene on chromosomes 11, 12, 13, 15, 16, and 19. To clarify the evolutionary relationships and classification of these DhDofs, the phylogenetic tree was initially built using the neighbor-joining method (Figure 2A). According to the branch value and the default classification, DhDof proteins were categorized into 5 groups and 11 subgroups, namely, A, B1, B2, B3, C1, C2.1, C2.2, C3, D1, D2, and E. Among them, DhDof4, DhDof5, DhDof17 and AtDof1.4 were on the B3 subgroup, which was different from the previous result. There were no DhDofs that belonged to the C3 subgroup. Group E is a new branch identified for the first time. The four genes DhDof6, DhDof7, DhDof8, and DhDof9 are independent from branches A and D2. In addition, the ML phylogenetic tree of the two species was also constructed (Figure 2B). The only difference between two phylogenetic trees was that AtDof1.7 and AtDof3.1 were classified in the D2 subgroup. The properties and subcellular localization of Dof genes were further predicted (Table 1). Subcellular localization prediction showed that all the DhDofs were localized in the nucleus.

To clarify the gene structure and composition of different types of Dof genes, the number and composition of motifs of DhDofs were first analyzed. Motif 1 was a conserved zinc-finger motif common to all Dof genes (Figure 3A). For the same subgroup, their conserved motif composition was the same, like DhDof10 and DhDof11, which both contained motif 1, motif 3 and motif 9. DhDof14 and DhDof15 included motif 1, motif 2, motif 8 and motif 10. The conserved domain database search results showed that all DhDofs contained the zf-Dof superfamily domains, while the N-ternimus of DhDof4, DhDof5 and DhDof21 also contained an InsA superfamily domain (Figure 3B). Gene structure analysis revealed that the number and composition of introns and exons in different DhDof genes varied greatly. Some DhDofs (such as DhDof2, DhDof20, DhDof10 and DhDof1) also contained 3’ UTR regions (Figure 3C). There was a long intron insertion between the two exons of DhDof18. The zinc finger structure in Dof domain is required to recognize AAAG/CTTT elements. A sequence comparison of the Dof proteins revealed that DhDofs contained a conserved CPRC(X)21CK(X)1C motif, which formed part of the Dof domain (Figure 4).

Comparative genomics analysis provided clues for the evolutionary relationships among species of D. spp. To clarify the genetic relationship of Dof genes between D. huoshanense and its closely related species, a collinearity analysis of Dof genes was employed. Microsynteny analysis showed that out of all the DhDof genes, only one pair of genes (Dhu000025304 and Dhu000009040) had collinearity, and dispersed and proximal duplication contributed to the expansion of DhDof genes. This implied that the formation of DhDof homologs was not due to whole genome duplication (Figure 5A). Synteny analysis from interspecies reveals that D. huoshanense and D. nobile exhibited the highest level of collinearity, with a total of 21 gene pairs, followed by D. huoshanense and D. chrysotoxum. However, there were only 7 and 6 syntenic blocks compared with O. sativa and A. thaliana, respectively. This was in accordance with the species’ relationship (Figure 5B). Artificial selection pressure provides the driving force for genome evolution and domestication. Through calculating Ka and Ks values, the Ka/Ks ratios of all DhDof genes were shown less than 1, which meant that DhDof genes were subject to purifying selection (Supplementary Table S1).

Based on the transcriptome sequencing and annotation results, the expression profile of Dofs was examined under MeJA treatment. Almost half of the DhDof genes were not expressed, and only a few genes were highly expressed after treatment (Figure 6A). DhDof20 and DhDof21 were more sensitive to MeJA stimulation and showed high expression throughout the entire period. This suggests that they are the core JA-responsive genes and participate in the downstream regulation of JA signaling. The expression of DhDof2 continued to increase with treatment time. The expression level reached its maximum at Time 7, showing a delayed JA-induced effect. On the contrary, the expression of DhDof17 was the highest at the beginning and decreased significantly after MeJA treatment, which suggested that DhDof17 was a negative effector gene of JA signaling. The expressions of DhDof6, DhDof8, and DhDof13 were highest at Time 3 with MeJA treatment and then gradually decreased, indicating that these DhDofs were early-responsive genes in JA signaling. To find out which Dof genes are co-expressed, a correlation analysis was performed based on their expression levels. DhDof6 showed a significantly positive correlation with DhDof8 and DhDof13. DhDof8 and DhDof13 also showed a significantly positive correlation, followed by DhDof15. DhDof6, DhDof8, and DhDof13 were co-expressed during MeJA treatment (Figure 6B). In addition, the expression of DhDof21 and DhDof17 showed a significant negative correlation, which was consistent with the previous result.

To determine which kind of environmental stresses these Dof genes could respond to, the composition of cis-acting elements was analyzed. The DhDof gene contained numerous types of regulatory elements. DhDofs from the same branch or subgroup contained similar cis-acting elements but differed in the distribution of seldom elements. For instance, DhDof18 contained a wound-induced element at the 5′end, but not in DhDof19. DhDof19 contained a GA-responsive element at the 5′end, but not in DhDof18. A total of 37 functional classifications were present here, which were mainly involved in growth, hormone induction, and abiotic stress response (Figure 7A). The most distributed element was related to light responsiveness, which was present in all DhDof genes. Among them, DhDof2, DhDof14, and DhDof21 were Dof genes that contained more functional categories. The second most common component was the MYB-binding element. DhDof1, DhDof12, and DhDof16 contained a large number of these elements. The third common component is the MYC-binding element. A total of 76 different cis-acting elements were discovered, with Myb, MYC, and Box4 elements being the most prominent (Figure 7B). The regulatory elements related to plant growth and development mainly included GCN4_motif, circadian, A-box, CAT-box, CCGTCC motif, re2f-1, HD-Zip 1 and RY-element. Cis-acting elements involved in hormone responses were deferentially distributed in the DhDof genes. DhDof2 and DhDof17 were the key genes involved in ABA responsiveness. DhDof21 was the key gene involved in GA responsiveness. DhDof2 and DhDof13 were the key genes involved in MeJA responsiveness. DhDof10 and DhDof11 were the key genes involved in SA responsiveness. DhDof22 was the key gene involved in ethylene responsiveness. The cis-acting elements involved in abiotic stress mainly included ARE, GC-motif, TC-rich repeats, MBS, DRE1, DRE core, LTR, STRE, WUN-motif, and WRE3, which were relevant to drought, low temperature, defense stresses, anaerobic induction, mechanical injury, etc. (Figure 7C).

The expression pattern of some of Dofs was confirmed by qRT-PCR method. The results showed that on the second and fourth days of ABA treatment, the expression of DhDof2 dramatically increased by approximately 92-fold and 49-fold compared with the control, respectively. DhDof2 was the main ABA-responsive gene and could be strongly induced. The expression levels of DhDof1 and DhDof12 exhibited an initial increase followed by a later drop under ABA treatment. The expression of DhDof4 decreased significantly after ABA treatment (Figure 8A). The expression level of DhDof13 decreased significantly at the 2nd hour after GA treatment and remained at a low level at 8 h. The expressions of DhDof21 and DhDof16 were also inhibited by GA. However, the expression of DhDof11 and DhDof17 decreased sharply at the 2nd hour of GA treatment but basically returned to the initial level at the 8th hour (Figure 8B). IAA treatment could significantly inhibit the expression of all DhDof candidates, and this inhibitory effect persisted until day 4. TGA-element (AACGAC) and AuxRR-core element (GGTCCAT) in the promoter region of DhDof were the negative auxin-responsive elements (Figure 8C).

Dof TFs are ubiquitously found across the full spectrum of the plant kingdom, encompassing both lower plants, such as unicellular green algae, and higher plants, like angiosperms and gymnosperms (Yu et al., 2020; Song et al., 2024). D. huoshanense did not have as many Dof genes as other angiosperms, which meant it probably had not experienced large-scale genome duplication (Melgar and Zelada, 2021). DhDofs were unevenly distributed on 12 pseudochromosomes and not equivalent to chromosome size (Figure 1). We infer this result from the unequal gene duplications of chromosomal segments. Microsynteny analysis revealed that only one pair of DhDof genes had collinearity, which proved that other duplication may be relevant (Figure 5). In addition to the Dof gene, some other transcription factors, like bHLH and WRKY, have been shown shared syntenic blocks among D. huoshanense, D. chrysotoxum, and D. nobile. This provided evidence for their close genetic relationship (He et al., 2023; Pan et al., 2023). D. huoshanense and D. nobile had the most collinear gene pairs, followed by D. huoshanense and D. chrysotoxum (Figure 5). This result was consistent with the other monocots, such as Z. mays, Brachypodium distachyon, Saccharum officinarum, etc. (Zou and Sun, 2023).

Phylogenetic analysis helps us understand the evolution and genetic relationship between Dof genes and Dof in other species. Based on the composition of domains and motifs, Dof genes are usually divided into four groups: A, B, C, and D. The B, C, and D groups can be further divided into several subgroups. There are also some studies that classify Dof genes into groups I-VII, which are obtained through domains and motifs. For example, 30 Dof proteins of O. sativa were divided into four groups and further subdivided into seven subgroups by constructing an unrooted phylogenetic tree with Arabidopsis, sorghum, and maize (Khan et al., 2021). Dof gene family was divided into six subfamilies (Group A-F) in sorghum, of which Group B contained only AtDof4.2, AtDof4.3, AtDof4.4, and AtDof4.5. 96 Dof genes from the wheat genome were divided into five subfamilies (Groups A-E), but only four AtDofs were categorized in Group A. Our results indicated that these four AtDofs were classified into the C3 subgroup, and no DhDof was assigned to this subgroup (Figure 2). Based on the phylogenetic tree, DhDof genes were divided into 11 subgroups, among which B3 and E subgroups are new branches identified (Figure 2). In addition to the E group, the D1 subgroup, also known as the CDF subfamily, contains the most Dof genes, including DhDof2, DhDof20, DhDof12, and DhDof16. AtDof1.4 was classified into Group IV or B2 subgroups in previous studies (Yanagisawa, 2002; Lijavetzky et al., 2003; Zhang et al., 2023). However, the phylogenetic trees in our study indicated that AtDof1.4 did not group with AtDof1.7, AtDof3.1, AtDof3.4, and AtDof5.8 in the D2 subgroup but as an independent subgroup. Three DhDof genes also belong to this subgroup (Figure 2). Remarkably, the E group, to which DhDof6/7/8/9 belonged, was a new branch that originated from the paralogs in D1. All the DhDof proteins possessed the highly conserved zinc-finger Dof domain (Figures 3, 4). Gene structure analysis revealed notable variations among distinct subgroups, whereas similar structures were observed in the common subgroup, as demonstrated in the motif analysis (Figure 3). Some subgroups contained their own unique motifs. Studies indicated that the introns of IbDof genes were relatively small, and most genes had only one intron or even no intron. Such intron-free genes may play a role in the accelerated stress response.

The involvement of Dof genes in response to biotic stress has been well documented. However, the regulatory function of Dof on abiotic stress responses was only reported in a limited number of plants (Zhang et al., 2018). RcDOFs expressed at two distinct levels varied in response to ABA (Waschburger et al., 2022). CmDOFs had a role in the response to ABA and SA, which contributed to differential expression patterns. Exogenous ABA specifically dramatically increased the expression of CmDOF12 and CmDOF20. OsDOF15 mediated the growth of the main root in rice when exposed to high salinity by releasing ethylene (Qin et al., 2019). The expressions of OsDof1 and OsDof19 increased at low temperatures in the cold-tolerant cultivar. Overexpression of OsDof1 resulted in a greater seed setting rate (Liu et al., 2021). An investigation and analysis of DhDofs were conducted in several abiotic stresses. Expression profiles showed the spatiotemporal profile of DhDofs after MeJA induction. DhDof20 and DhDof21 were identified as the main responsive genes of JA signaling, and DhDof17 was considered to be an early-responsive gene (Figure 6). DhDofs showed distinct expression patterns under different abiotic stresses. Under a series of abiotic stresses, the expression patterns of the DhDofs were determined at different phases. DhDof2 was strongly induced by ABA, while DhDof7 was significantly induced by low temperatures. The expression of most DhDof genes showed an inhibitory state after treatments. Most genes were inhibited after IAA and SA treatment, which were probably related to the cis-acting elements in their promoter (Figures 7, 8). Our results revealed how Dof genes work and how Dendrobium plants interact and coordinate with the external environment.