The D. nobile plants used in this study were cultivated at the Anhui University of Chinese Medicine (Hefei, Anhui, China). D. nobile samples are stored in the Herbarium of Anhui University of Chinese Medicine (code: 20220601). The plant was identified as D. nobile by Professor Yang Qingshan of Anhui University of Traditional Chinese Medicine. In this study, the root, stem, leaf, and flower of D. nobile were tested for gene expression tissue specificity. Each tissue was collected from three different plants, and the collected samples were immediately stored at –80 °C. For the stress treatment, the plants were exposed to 100 µM abscisic acid (ABA), 100 µM methyl jasmonate (MeJA), and low-temperature (4 °C). After treatment, the leaves were collected at different time points (0, 2, 6, 12, and 24 h) and frozen in liquid nitrogen.

The genomes of D. nobile were downloaded from the National Center for Biotechnology Information (NCBI) (https://www.ncbi.nlm.nih.gov/genome/?term=dendrobium+nobile) to identify the candidate R2R3-MYB genes. The hidden Markov model (HMM) profile of the MYB domain (PF00249) from the Pfam database was used to search the protein sequences [17, 18]. The putative DnMYB sequences were checked by the National Center for Biotechnology Information (NCBI). Redundant sequences were removed using CD-HIT online software (http://www.bioinformatics.org/cd-hit/). The remaining sequences were submitted to HMMscan online analysis (https://www.ebi.ac.uk/Tools/hmmer/search/hmmscan) to determine the presence and quantity of MYB domain, and screen for R2R3-MYB sequences were screened for structural features. A total of 126 R2R3-MYB TF sequences were retrieved from the Arabidopsis Information Resource (TAIR) database (https://www.arabidopsis.org/).

The physicochemical properties of the protein sequences were analyzed using the ProtParam tool of ExPASy (https://web.expasy.org/protparam/). TMHMM software (v 2.0) (DTU Health Tech, Lyngby, Denmark) was used to determine transmembrane domain numbers. The subcellular localization of the MYB proteins was predicted using the PSORT tool.

MEGA X software (https://www.megasoftware.net/) was used for sequence alignment, and a phylogenetic tree based on D. nobile and Arabidopsis R2R3-MYB proteins was constructed using the neighbor-joining method. The 126 R2R3-MYB protein sequences of Arabidopsis were used as references. Bootstrap analysis with 1000 repetitions was performed to compute the reliability of the phylogenetic tree.

For gene structure analysis, the D. nobile genome sequences corresponding to the DnMYB genes were downloaded from NCBI. The Gene Structure Display Server (v.2.04) (http://gsds.gao-lab.org/) was employed to analyze the exons-introns structures of the obtained full sequences. The MEME suite (Multiple Expectation Maximization for Motif Elicitation, version 5.4.1) (https://meme-suite.org/meme/tools/meme) was used for conserved motif analysis [19, 20].

The physical positions of the DnMYB genes on the chromosome were mapped using MG2G online software (http://mg2c.iask.in/mg2c_v2.1/) [21]. The synteny of DnMYB genes between D. nobile and its homologs in Dendrobium chrysotoxum and Arabidopsis thaliana was analyzed and visualized using TBtools software (https://github.com/CJ-Chen/TBtools/releases) [22].

The sequences of the upstream 2 kb regulatory regions (from the translation start site) of the DnMYB genes were obtained from the D. nobile genome database, and were analyzed using PlantCARE software (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/).

RT-qPCR analysis was performed using a SuperReal PreMix Plus SyBr Green PCR kit (Qiagen, Shanghai, China) on a Cobas z480 Real-Time PCR System as previously described by Wang et al. [23]. Primers were designed with Primer 5.0 (Supplementary Table 1). The reactions mixtures contained 2 µL of diluted cDNA, 0.6 µL of each primer, 10 µL of 2 $\times$ SuperReal PreMix Plus, and 6.8 µL of RNase-free double-distilled water. The conditions for all the RT-qPCRs were performed as follows: denaturation at 95 °C for 15 min, followed by 45 cycles of 95 °C for 10 s, 58 °C for 20 s, and 72 °C for 30 s. Each sample contained three biological replicates. The 18S housekeeping gene was used as a reference and the relative expression level of each gene was calculated using the 2${}^{-\mathrm{\Delta }\mathrm{\Delta }\text{Ct}}$ method [24].

All experiments were performed with three replicates, and each replicate was repeated three times. The significance of differences in the expression levels of selected 8 DnMYB genes was assessed using the student’s t-test. p-value 0.05 was considered significant.

To identify D. nobile MYB proteins, the genomic sequence of the plant was downloaded from NCBI. After a searching performed through the Biolinux system using PF00249 as a seed sequence to perform HMMER alignment, the redundant sequences were removed by the CD-HIT online analysis tool, yielding a total of 276 candidate MYB genes. The number of MYB domains in the sequence of each candidates was assessed using HMMscan online analysis software, and 125 sequences were screened out that matched the characteristics of R2R3-MYB genes. The R2R3-MYB genes of D. nobile were named DnMYB1 to DnMYB125 (Supplementary Table 2). These DnMYBs were classified into four distinct subfamilies: R1-MYB (123 proteins, 44.57%), R2R3-MYB (125 proteins, 45.29%), 3R-MYB (27 protein, 9.78%), and 4R-MYB (1 protein, 0.36%). The observation that the R2R3-MYB family encompassed the greatest number of DnMYB genes is consistent with that previously reported, namely, that the R2R3-MYB subfamily is the largest subfamily of TFs in plants [25]. The physicochemical properties of the125 DnMYB genes, including their chromosomal locus and the length, molecular weight (MW), isoelectric point (pI), number of transmembrane regions, and subcellular localization of MYB proteins, are listed in Supplementary Table 2. The lengths of the DnMYB proteins ranged from 88 (DnMYB50) to 2424 amino acids (DnMYB64), with an average number of 348 amino acids. The MW and pI of predicted DnMYB proteins ranged from 10.32 to 270.01 kDa and 4.88 to 10.04. None of DnMYBs were predicted to be transmembrane proteins. Subcellular localization results show that all predicted DnMYB proteins are localized in the nucleus, which is consistent with their functions as TFs [26].

In order to explore the evolutionary relationship among the DnMYB proteins, a phylogenetic tree was constructed based on 125 and 126 MYB proteins from D. nobile and A. thalian, respectively (Fig. 1). The results showed most of the DnMYB could be classified into different subclasses, which was supported by high bootstrap values. The phylogenetic tree of the 125 DnMYB were divided into 26 subgroups with the number of members ranging from 1 to 26. The DnMYB proteins has the most members in the S21 subfamily, including 14 DnMYBs. This result is consistent with the classification results of D. officinale and P. aphrodite. The S14 and S18 subfamily contained 7 DnMYBs respectively, followed by 5 DnMYBs were both divided into S4 and S17 subfamilies. The S13, S22, and S25 subfamilies all contained 4 DnMYBs, while the S1 and S11 subfamilies encompassed 3 DnMYBs each. Four subfamilies, including S2, S6, S10, and S20 all harbored 2 DnMYBs. Each of the remaining subfamilies (S5, S9, S19, S23 and S24) contained only one DnMYBs. Interestingly, no DnMYB proteins were classified into the S12 subfamily. Arabidopsis R2R3-MYB TFs are classified into 25 subclasses. Members of subclasses S1, S2, S20, and S22 are involved in regulating plant responses to abiotic stresses [27], while R2R3-MYB TFs within the S4, S5, and S6 subfamilies are involved in plant anthocyanin synthesis [27]. The fact that some AtMYB proteins with similar functions are clustered in the same clade may indicate that DnMYB proteins within a clade may also have similar functions [28]. However, the nature of these functions requires further investigation.

Fig. 1.

Phylogenetic tree of R2R3-MYBs in D. nobile and A. thaliana. The figure was generated by using MEGA X software coupled with a neighbor-joining method and a bootstrap of 1000 replicate.

Exon/intron analysis showed that 125 DnMYB genes had introns in their coding regions. The 125 DnMYB genes displayed a high degree of diversity in the number and relative positions of introns and exons. It is well-established that genetic structural diversity can explain the evolution of polygenic families [28]. The gene structure map of these 125 DnMYB genes was constructed, and the results showed that the clustering pattern of 125 DnMYB genes was not clearly consistent with the exon/intron structure (Fig. 2) (Supplementary Table 3). Interestingly, most genes in the same subclass showed similar exon/intron structure, which is in line with that reported for Liriodendron [28]. The number of exons ranged from 1 to 13, and most of these genes contained either two (27 genes) or three exons (67 genes) based to the 125 DnMYB gene structure.

Fig. 2.

Phylogenetic relationships, conserved motifs, and gene structures of R2R3-MYBs family member. (a) The phylogenetic tree was constructed based on the full-length sequences of 125 DnMYB proteins using the neighbor-joining method in MEGA X, with 1000 bootstrap replicates. (b) Conserved motifs of DnMYB proteins. Different motifs are represented by various colored boxes. (c) Exon/intron structures of DnMYB genes. Exon(s), and intron(s) are represented by black boxes and black lines, respectively.

There are many conserved sequences in MYB TFs and these conserved sequences may combine with some parts of DNA and play a role in expression regulation [29]. To comprehensively investigate the potentially conserved motifs in 125 DnMYB proteins of D. nobile, we systematically evaluated their evolutionary relationship and motif composition using MEME online analysis software. A total of 10 motifs were identified in the 125 proteins. Motif 3, located in N-terminal domain, was present in almost all the DnMYBs. Compared with the N-terminus, the C-terminus had greater motif and sequence variability. These differences suggest that different subgroups may have different functions [30] (Fig. 2). Moreover, most of the 125 DnMYBs contained motifs 1, 2, 3, 4, and 6, and most members in the same subgroup harbored more than one identical motif. Motifs 3, 6, and 2 constitute the complete R2 domain, which has three highly conserved tryptophan (W) residues. Motif 1 comprises the complete R3 domain, in which the first tryptophan (W) residue is often replaced by phenylalanine (F), while the remaining two tryptophan (W) residues, which are separated by 19 amino acids, are highly conserved. These characteristics are consistent with those of the MYB domains identified in soybean [30], Ginkgo biloba [31], and Arabidopsis [32], among other species. Interestingly, the results also showed that most members of the same subclass have similar motif composition.

To understand the genomic distribution of the R2R3-MYB genes, we employed MG2G online software to identify the chromosomal location of each of the 125 DnMYB genes. D. nobile has 19 chromosomes (Chr1 to Chr19; Supplementary Table 4) and DnMYBs were present on all of them, but displayed uneven the arrangement and density (Fig. 3). Chromosome 9 had the highest percentage of DnMYBs (12%) and chromosome 4 the lowest. There was a high density of DnMYBs located on chromosomes 6, 11, 12, 13, and 17.

Fig. 3.

Distribution of 125 DnMYB genes on chromosomes of D. nobile.

A collinearity analysis was undertaken to detect DnMYBs gene duplication events. A total of 46 pairs of gene duplications were identified among the DnMYBs (Fig. 4). The DnMYB genes were homologous to genes in other plants, and syntenic conservation was found among D. chrysotoxum (174 orthologous gene pairs dispersed on all chromosomes). 28 collinear gene pairs between D. nobile and A. thaliana (Fig. 5). The number of orthologous events of D. nobile-D. chrysotoxum was far greater than that of D. nobile-Arabidopsis, and the closer evolutionary distance between D. nobile and D. chrysotoxum was confirmed.

Fig. 4.

Synteny analysis of interchromosomal relationships of DnMYB genes. Red lines indicate duplicated gene pair.

Fig. 5.

Synteny analysis of DnMYBs genes in the genomes between D. nobile and D. chrysotoxum or A. thaliana. (a) D. nobile and D. chrysotoxum. (b) D. nobile and A. thaliana. The gray lines show collinear blocks. The red and blue lines indicate the syntenic gene pairs between D. nobile and D. chrysotoxum or A. thaliana, respectively.

Cis-acting elements in promoters are bound by TFs and play an important role in the regulation of gene expression [33]. Here, we used the PlantCARE online prediction tool to identify cis-acting elements located within the first 2,000 bp upstream of the transcription start site of 125 DnMYB genes (Supplementary Table 5). Besides the basic elements, the predicted cis-acting elements could be divided into three groups, including hormone-responsive elements (ABA, salicylic acid, gibberellin, MeJA, auxin), abiotic stress elements (low temperature-responsive elements and drought-responsive elements), and metabolism-related elements (regulation of flavonoid biosynthesis) (Supplementary Fig. 1). The MeJA-responsive element was the most abundant element with 350 instances. It is well known that regulates a wide variety of physiological processes such as plant growth and development, especially the response of plants to biotic and abiotic stresses [34]. DnMYB23 contained 24 MeJA-responsive elements in its 2-kb upstream regulatory region. The hormone-related elements identified included 210 ABA-responsive elements, 87 gibberellin-responsive elements, 64 auxin-responsive elements, and 58 salicylic acid-responsive elements. Among abiotic stress elements, there were 60 low temperature responsive elements, 3 of which were in the promoter region of DnMYB100. A total of 80 drought responsive elements were also predicted, with DnMYB34, DnMYB37, DnMYB58, DnMYB6, DnMYB74, and DnMYB83 harboring 3. Finally, 14 response elements associated with the regulation flavonoid biosynthesis were also identified in some of the 2-kb sequences analyzed. These results suggested that the expression of DnMYB genes is controlled by complex regulatory networks.

DnMYB genes in different subgroups may have different functions in the growth and development of D. nobile. To determine the spatial expression pattern of MYB genes in D. nobile, we measured the expression levels of 8 genes from eight different subtypes in four tissues (roots, stems, leaves, and flower) using RT-qPCR. DnMYB8, DnMYB14, DnMYB26, DnMYB27, DnMYB81, DnMYB90, DnMYB108, and DnMYB116 belong to the S22, S2, S21, S1, S7, S4, S19, and S20 subfamilies, respectively. As shown in Fig. 6, 8 genes were expressed in all the tissues. The results indicated that DnMYB8, DnMYB14, DnMYB26 and DnMYB81 had similar expression patterns and were predominantly expressed in root. Meanwhile, the expression of DnMYB27was higher in flowers than in the other three tissues. Similarly, DnMYB108 showed a tissue-specific expression pattern, showing 30-fold higher expression levels in flowers than in other tissues. DnMYB90 was primarily expressed in the leaf and stem. DnMYB116, classified in the S20 subfamily, exhibited relatively high expression levels in the leaf and/or stem. The results suggested that these genes may play a role in the growth and development of D. nobile.

Fig. 6.

Relative expression levels of 8 selected DnMYB genes in different tissues. Error bars indicate the standard error from three independent replicates. Different lowercase letters (a, b, c, and d) are significantly different (p-value 0.05).

Based on the reported functions of R2R3-MYB in other plants, we selected 8 DnMYBs belonging to eight subfamilies to study the roles of R2R3-MYB in the regulation of secondary metabolism and the responses to abiotic stress for D. nobile. We measured the expression levels of these genes by RT-qPCR after treatment with ABA, MeJA and low-temperature (4 °C). Six of the 8 DnMYB genes-DnMYB8, DnMYB14, DnMYB26, DnMYB27, DnMYB90, and DnMYB116-showed increased expression levels at different times under the low temperature, ABA, and MeJA treatment (Fig. 7). While DnMYB81 was the most responsive to low temperature, and showed reduced expression patterns under ABA- and MeJA- mediated stress. DnMYB108 showed downregulated expression patterns under both low-temperature and MeJA stress. With ABA treatment, the expression of DnMYB108 reached its lowest level at 12 h and peaked after 24 h.

Fig. 7.

Expression of 8 DnMYB genes at 0, 2, 6, 12, and 24 h after various stresses treatment. LT means Low temperature, ABA means abscisic acid, MeJA means methyl jasmonate. Error bars indicate the standard error from three independent replicates. Different lowercase letters (a, b, c, d and e) are significantly different (p-value 0.05).

R2R3-MYB TFs have been identified in many plants, including Arabidopsis [32], maize [35], and soybean [30], etc. To further explain the evolutionary relationship of DnMYB genes, we performed a phylogenetic analysis based on 125 DnMYB proteins and 126 AtMYB proteins. Within the phylogenetic tree, these proteins were divided into 26 subgroups. Among them, S21 subgroup had the highest number of DnMYB proteins. In addition, we also found that some DnMYB genes had not been classified, indicating that their protein products may have special functions. We further analyzed the conserved motifs of the DnMYB proteins based on evolutionary relationships and identified a total of 10 conserved motifs. Most DnMYB proteins shared 3 motifs (motifs 1, 5, and 8; Fig. 2), which suggested that these motifs may have species specific functions. Our results indicated that although motif composition varied according to subgroup, the similarity in motif composition within the same subgroup was very high. The motif composition of closely related members in phylogenetic trees is consistent, and there have also been reports of tobacco [36]. Gene structure usually reflects the evolution of gene family [37]. In our gene structure analysis, we found that most genes in the same subgroup had either two or three exons, which was consistent with that reported for D. officinale and P. aphrodite [15]. These results indicate that the DnMYB protein is relatively conserved during evolution.

Low temperature, drought, and other stresses can seriously affect the growth and variety of D. nobile. Secondary metabolites such as flavonoids and anthocyanins contribute to plants adaptations to hostile environments. Meanwhile, flavonoid biosynthesis can be induced under a wide range of abiotic stresses. Accordingly, stress resistance of D. nobile can be enhanced via the regulation of genes related to secondary metabolite production. MYB genes contribute to the development of trichomes and root hairs, and play a crucial role in regulating the synthesis of plant secondary metabolites, such as flavonoids and anthocyanins, that help plants resist stress [38]. To identify DnMYB genes that may be involved in regulating growth, development, and secondary metabolic synthesis of D. nobile, we selected 8 genes, each from a different subgroup, and measured their expression levels in different tissues and under different stresses. Gene functions are closely related to tissue-specific expression. Flavonoids, the major pigments in plants, have attracted substantial research attention owing to their important biological properties. It has been reported that flavonoid biosynthesis in plants can be induced under a wide range of abiotic stresses [38]. For instance, in wheat, the key enzymes in the flavonoid biosynthesis pathway were found to be upregulated by drought stress [39]. AtMYB4 and AtMYB7, which belong to the S4 subfamily, are related to the flavonol synthesis pathway [40]. Similarly, the AtMYB11, AtMYB12 and AtMYB111 genes within the S7 subfamily are involved in the regulation of flavonol biosynthesis in Arabidopsis [41]. Genes with high homology in the same branch of a phylogenetic tree generally have high sequence similarity and may also have similar functions. The DnMYB90 gene, which was also categorized into the S4 subfamily, is a homolog of AtMYB4 and AtMYB7. RT-qPCR analysis showed that although DnMYB90 was expressed in all the tissues of D. nobile, the highest expression level was detected in the leaf. DnMYB90 was also expressed at very high levels under phytohormonal and abiotic stresses (ABA, MeJA, low-temperature stress), but it was most responsive to low temperature. DnMYB81, belonging to the S7 subfamily, exhibited close homology with AtMYB11, AtMYB12, and AtMYB111. DnMYB81, belonging to the S7 subfamily, was highly homologous to AtMYB11, AtMYB12, and AtMYB111. The results of the RT-qPCR analysis showed that DnMYB81 was expressed in all the tissues, but was most highly expressed in the root. DnMYB81 was also responsive to phytohormonal and abiotic stresses. This gene was most responsive to low temperature and showed downregulated expression patterns under ABA and MeJA stress. Future, studies should focus on whether DnMYB81 and DnMYB90 can regulate flavonoid synthesis.

Several R2R3-MYB genes were reported to be involved in regulating responses to biotic and abiotic stresses. For example, the overexpression of AtMYB94 and AtMYB96, the S1 subfamily genes, enhances drought tolerance in plants [42]. AtMYB14, a member of the S2 subfamily, regulates cold tolerance in Arabidopsis [43]. Under conditions of phosphate starvation, the S2 subfamily AtMYB20, which plays a role in plant responses to salt and drought stress, directly regulates the expression of miR399f [44]. Meanwhile, the Overexpression of the S21 subfamily member AtMYB52 confers ABA hypersensitivity and drought tolerance [45]. AtMYB44, belonging to the S22 subfamily, responds to salt, MeJA and drought [46]. To identify DnMYB genes related to hormonal and abiotic stresses, we detected the expression levels of DnMYB genes from the S1, S2, S20, S21, and S22 subgroups following low temperature, ABA, and MeJA treatments. We found that DnMYB27 in the S1 subfamily, DnMYB14 in the S2 subfamily, DnMYB116 in the S20 subfamily, DnMYB26 in the S21 subfamily, and DnMYB8 in the S22 subfamily, which we found were homologous to stress-response genes, were expressed in in roots, stems, leaves, and flowers. To further identify candidate genes involved in hormone signaling pathways or abiotic stress responses, we measured their expression levels under low temperature, ABA and MeJA stress. The results showed that the five genes showed increased expression levels at different times under different stress, which was consistent with the presence of the relevant cis-acting elements in their promoters (Figs. 2,7). Interestingly, the expression level of DnMYB108, classified in the S19 subfamily, gradually decreased after treatment with cold and MeJA stress. These results may suggest that members of this subfamily play an important role in the response of D. nobile to abiotic stress. Some of the DnMYB genes screened in this study may play a role in responses to abiotic stresses such as low temperature, ABA and MeJA stress; however, the underlying mechanism remains unclear and requires future clarification.