Tissue-cultured seedlings of D. huoshanense were subjected to a sterilization process and then planted on Murashige and Skoog (MS) media. This cultivation occurred in the tissue culture room with long-day conditions (16 h light, 20°C/8 h dark) of West Anhui University’s Anhui Engineering Technology Research Center of Plant Cell Engineering, located in Lu’an City, Anhui Province. The methodology employed in this study was based on previously published articles (Cao et al., 2019). Subsequently, 100 μM MeJA (Aladdin, Shanghai of China), 100 μM SA (Aladdin, Shanghai of China), and 100 μM ABA (Aladdin, Shanghai of China) were added after 0.22 μM microporous filtration. Referring to the method of Cao et al. (2015), 20% PEG6000 was added to the culture medium to simulate drought stress. The leaves were sampled at 1 h, 4 h, and 16 h after treatment. There were three biological replicates per set of samples, and each sample was collected individually for each induction treatment and immediately stored at −80°C to facilitate RNA isolation. In addition, untreated D. huoshanense leaves were used as the control group for comparative purposes.

Using the BLAST tool available on the Dendrobium genome (https://ftp.cngb.org/pub/CNSA/data3) information website, a comparative analysis was conducted based on the KNOX protein sequence in Arabidopsis as the reference sequence (Furumizu et al., 2015), with a BLASTp (E-value = 0.001) and other default parameters. The objective was to identify candidate KNOX genes within the Dendrobium genome. Furthermore, the conserved KNOX1 (PF03790) or KNOX2 (PF03791) domains should be selected for further analysis using the structural properties that were identified using the Pfam website (Cheng et al., 2019). Subsequently, the isoelectric point and molecular weight of the KNOX gene family members of D. huoshanense were analyzed by ProtParam (https://web.expasy.org/protoparam/) (Finn et al., 2013). We predict of subcellular localization of the KNOX protein on the website WoLF PSORT (https://wolfpsort.hgc.jp/).

To elucidate the evolutionary relationship of the DhKNOX gene, we employed MEGA 5.2 software to construct a KNOX gene phylogenetic tree with 1,000 bootstrap values (Tamura et al., 2011). The GenBank accession codes were used to construct the phylogenetic tree (Supplementary Table S1). According to the original annotation file in GFF3 format, the CDS of DhKNOXs was used for gene structure analysis and visualized with TBtools (Chen et al., 2018).

The chromosome location information of the KNOX gene family members was obtained from the Dendrobium genome annotation file. The KNOX gene information from the Dendrobium genome database was utilized to perform chromosome localization using the MapDraw V2.1 software (Lee et al., 2013).

TBtools software was used to obtain the 2,000-bp promoter region upstream of the CDS sequence of each KNOX gene Dendrobium from the genome database and to analyze cis-acting elements in the KNOX gene promoter. The obtained promoter region sequence was submitted to the PlantCARE database to predict cis-acting elements (Cheng et al., 2019). We implement the prediction of protein interactions among members of the DhKNOX gene family through the STRING website (https://cn.string-db.org) (Yang et al., 2022).

The sequence alignment of KNOX proteins was performed using ClustalX software (Cao et al., 2016). Additionally, the conserved motifs of KNOX proteins were identified by the MEME online tools (Bailey et al., 2015).

The total RNA from D. huoshanense leaves was extracted using a plant RNA extraction kit (Tiangen). qRT-PCR was used to detect the expression level of the KNOX gene at different times after treatment. The internal control gene was used as a reference (Fan et al., 2016). The specific primers for each gene are shown in Supplementary Table S2. The qRT-PCR reaction system was 20 µL, and the reaction procedure was 95°C for 3 min, 95°C for 10 s, 60°C for 30 s, and 45 cycles. The melting (65°C~95°C) program is finally run after signal acquisition. Three biological repeats were set for each sample, and each biological duplicate sample was tested three times. Finally, the relative gene expression of KNOX genes was calculated using 2^−ΔΔCt methods (Livak and Schmittgen, 2001).

The query sequence used for this analysis consisted of the KNOX family members from Arabidopsis thaliana, totaling nine sequences. Based on this query, 19 nonredundant KNOX family members were identified from the genome of Dendrobium. These newly identified members were named DhKNOX1 to DhKNOX19. Table 1 presents various characteristics of the DhKNOX genes in Dendrobium. The coding region lengths of the DhKNOX genes range from 369 to 1,872 bp, while the corresponding coding amino acid lengths range from 123 to 624. The molecular weight of the DhKNOX proteins spans from 18.81 to 69.0 kD, and their isoelectric points range from 4.29 to 10.25. Through the examination of predicted subcellular localization, it was ascertained that the nucleus emerges as the predominant site for the subcellular localization of the majority of DhKNOX genes. This observation is in consonance with their established functions as transcription factors.

Multiple sequence alignments revealed that the KNOX protein sequence of Dendrobium has four relatively conserved domains (Figure 1). The KNOX1 and KNOX2 domains are located at the N-terminus of the protein, while ELK and HOX are located at the C-terminus of the protein. The HOX domain demonstrates the highest level of conservation and exhibits the typical structural characteristics found in the TALE homeobox protein superfamily. In addition, the conserved domains of KNOX1, KNOX2, and ELK exhibit notable conservation. For instance, the conserved domain KNOX2 displays a highly conserved E-L-D amino acid sequence, while the conserved domain ELK features a highly conserved E-L-K amino acid sequence.

To elucidate the distribution of individual DhKNOX members within the Dendrobium genome, a distribution map was generated. This map visually depicts the arrangement of the 19 DhKNOX genes across each specific chromosome of Dendrobium, as shown in Figure 2. The 19 DhKNOX genes have been clearly located on chromosomes and are distributed on 11 different chromosomes. The distribution of DhKNOX genes across the 11 chromosomes of Dendrobium is widely scattered. Specifically, chromosomes 16, 17, and 18 harbor two, five, and four DhKNOX genes, respectively, while the remaining eight chromosomes each contain one DhKNOX gene. Furthermore, a notable observation is that the majority of DhKNOX genes are located toward the ends of the chromosomes.

The KNOX protein sequences of well-known model plants such as Arabidopsis and rice, along with the KNOX protein sequence of D. huoshanense were constructed using MEGA 5.2 software (Figure 3). The results showed that these KNOX genes were divided into three categories: class I, class II, and class KNATM. We further analyzed and found that these proteins include six subfamilies: class I includes three subfamilies, class II consists of two subfamilies, and class M contains only the KNATM-like subfamily. In these three large families, most DhKNOX genes are distributed in class II. Among the DhKNOX proteins, both the BP-like subfamily and STM-like subfamily contains two DhKNOX members. While there is no member of the KNAT2/6-like subfamily, the other 13 DhKNOX proteins are classified as belonging to the KNAT7-like subfamily. Previous research reported that class I KNOX genes are related to lignin metabolism (Cheng et al., 2019). In our study, DhKNOX2 and DhKNOX3 belong to the BP-like subfamily and speculated that they may have similar biological functions. Interestingly, Scofield et al. (2014) demonstrated in Arabidopsis that STM-like proteins can regulate the expression level of the AtKNAT1/BP gene. DhKNOX1 and DhKNOX4 belong to the same subfamily, so we speculate that they may have similar biological functions.

To gain deeper insights into the Dendrobium KNOX gene family, we utilized the intraspecific phylogenetic tree (Figure 4A) and protein-conserved domain distribution (Figure 4B), as well as the gene structure of DhKNOX genes (Figure 4C). Subsequently, employing the MEME website (Figure 4B), we identified 10 conserved motifs across the 19 DhKNOX proteins. Significantly, almost all DhKNOX proteins encompass at least one conserved domain, such as motif1, motif3, motif4, and motif10. This observation reinforces the reliability of the DhKNOX gene screening results. However, it is worth mentioning that DhKNOX6 and DhKNOX12 in class M are exceptions, lacking conserved domains. Furthermore, from an evolutionary perspective, the DhKNOX genes within the same subfamily and exhibiting close relationships also demonstrate consistency in gene structure composition and protein domain allocation, which supports their close evolutionary relationships.

The findings indicated that all DhKNOX genes comprised both exons and introns, and the DhKNOX gene family predominantly exhibited a pattern of four to six exons and three to five introns. The number of gene exons in class I was between three and six, while almost all genes in class II contained five or six exons (Figure 4C). Among the DhKNOX family members, DhKNOX3, DhKNOX4, DhKNOX18, and DhKNOX19 contain only one intron. Most DhKNOX members have three introns. Generally, closely related members have similar genetic structures, such as DhKNOX11/DhKNOX14 and DhKNOX15/DhKNOX13. It is noteworthy that except for DhKNOX5, DhKNOX9, and DhKNOX15, class I genes contain a long intron (Figure 4C), which is consistent with the structural characteristics of the class I KNOX gene family (Morimoto et al., 2005).

To analyze the potential expression regulation mechanism of DhKNOX gene family members, we identified the cis-acting elements of each member (Figure 5). For the identified cis-acting elements, the elements related to plant hormones, plant growth, participation in stress responses, and light responses were selected for analysis. First, based on the 2,000-bp promoter sequence, stress-related elements, including LTR, MBS, ARE, and other elements, were identified, and transcriptional regulatory-related elements, including MYB, W-box, and MYB-also, were found. In addition, we discovered and validated a couple of elements related to plant hormones, such as P-box, TGA-element, ABRE, and other elements. It is worth noting that the most common motif is the cis-acting element of the motif related to ABA reactivity, accounting for 36% of the hormone response motifs scanned. Cis-acting elements of the TGACG motif related to MeJA reactivity accounted for 13%. In addition, we found that the TCA element in response to SA, among the 19 DhKNOX gene promoters, appeared 29 times and accounted for 10%. The MBS element was related to drought stress, which accounted for 12%. These findings suggest that the transcription of DhKNOX genes may be affected by these hormones and drought stress.

We conducted protein interaction predictions among the members of the DhKNOX gene family by using the STRING website (Figure 6). The results revealed a significant number of interactions among the family members. We found protein-to-protein interactions among 15 members of the DhKNOX family, including DhKNOX1, DhKNOX2, DhKNOX3, DhKNOX5, and DhKNOX19. There are 17 nodes and 42 sets of interaction relationships in the interaction network. DhKNOX2, DhKNOX6, and DhKNOX12 are located at the center of the entire protein interaction network and can interact with the proteins DhOFP9 and DhMYB6. In summary, the members of the DhKNOX gene family have the ability to form dimers or polymers through protein interactions and to carry out transcriptional regulatory functions.

Various abiotic stresses may affect the growth and development of plants and ultimately affect the regulation of a series of stress-related genes (Grallath et al., 2005). Consequently, it becomes imperative to clarify the regulatory pathway and identify pivotal regulatory factors implicated in hormone responses in D. huoshanense. To better understand the changes in the expression level of the DhKNOX gene under various hormone treatments, we performed qRT-PCR experiments to analyze their expression patterns under MeJA, ABA, and SA treatments.

After treatment with MeJA, we found that the expression of 10 DhKNOX genes was downregulated to different degrees (Figure 7). The results showed that DhKNOX1, DhKNOX2, DhKNOX3, DhKNOX4, DhKNOX5, DhKNOX6, and DhKNOX12 were obviously significantly and rapidly upregulated at all time points. Among these DhKNOX genes, we found that the highest expression levels of DhKNOX1, DhKNOX2, DhKNOX3, DhKNOX4, and DhKNOX5 occurred 16 h after treatment. The expression of one DhKNOX6 and DhKNOX12 gene peaked at 4 h. Notably, it was found that four genes, DhKNOX1, DhKNOX2, DhKNOX3, and DhKNOX4, were all clustered in class I, and DhKNOX6 and DhKNOX12 were clustered in class II. Moreover, eight DhKNOX genes presented reduced expression levels to different degrees. Interestingly, we found that the expression of two DhKNOX genes (DhKNX13 and DhKNOX15) was not obviously changed at any time point, and the two DhKNOX genes lacked the TGACG motif.

ABA induced the transcription of 19 DhKNOX genes. The expression levels of 17 DhKNOX genes increased to varying degrees (Figure 8). Among them, 11 DhKNOX genes were significantly upregulated at the last time point (16 h), such as DhKNOX2 and DhKNOX3, which were upregulated by more than 100- and 500-fold, respectively. Interestingly, we found that the expression of two DhKNOX genes was clustered in class I. However, we found that DhKNOX2 had only two ABREs, and DhKNOX3 lacked the ABREs. In addition, two DhKNOX genes (DhKNOX6 and DhKNOX12) clustered in class II, and their expression levels decreased significantly.

Under SA treatment, in leaves, the expression level of 17 DhKNOX genes increased significantly at three treatment time points (1 h, 4 h, and 16 h) (Figure 9). We further found that the highest expression levels of 12 DhKNOX genes occurred after 16 h of treatment. The expression pattern of most other DhKNOX genes showed a trend of increasing first and then decreasing, and DhKNOX15 and DhKNOX18 had the highest expression levels after 1 h. Four of 17 genes reached their highest expression level after 4 h of SA treatment and 4 h but decreased after 16 h, such as DhKNOX2, DhKNOX13, DhKNOX16, and DhKNOX19. Of these four DhKNOX genes, only DhKNOX2 was clustered in class I and had one TCA element. Notably, DhKNOX13, DhKNOX16, and DhKNOX19 were clustered in class I, and all lacked TCA elements. Only two DhKNOX genes were obviously rapidly and significantly downregulated at all time points.

The expression level of most DhKNOX genes showed a significant increase or decrease under 20% PEG treatment, as shown in Figure 10. We further found that the expression levels of the five genes were clustered in class I, class II, and class III and gradually increased with the prolongation of treatment time, and the expression level of DhKNOX13 increased by 120 times compared to the control group, which contains three MBS elements. In addition, the expression levels of four genes, DhKNOX4, DhKNOX5, DhKNOX6, and DhKNOX18, did not show significant changes, and this gene promoter does not contain MBS elements. Additionally, the expression of DhKNOX7, DhKNOX8, and DhKNOX9 was consistently suppressed under drought treatment. Moreover, three DhKNOX genes (DhKNOX14, DhKNOX15, and DhKNOX16) were clustered in class I, and the expression level of these genes first increased and then decreased, reaching a maximum at 12 h. Furthermore, only two DhKNOX genes were clustered in class I, and the expression of these two genes varied irregularly.