Mining the candidate transcription factors modulating dendrobine biosynthesis under phosphate deficiency in Dendrobium officinale Kimura & Migo

TL;DR

This study explores how low phosphorus stress affects dendrobine production in Dendrobium officinale and identifies key genes and transcription factors involved.

Contribution

The study identifies candidate transcription factors modulating dendrobine biosynthesis under phosphate deficiency in Dendrobium officinale.

Findings

Transcriptomic analysis revealed 10 phosphorus transport-related transcription factors and 21 TFs associated with dendrobine biosynthesis.

Low phosphorus stress significantly increases dendrobine content in Dendrobium officinale.

DEGs related to phosphate transport and alkaloid biosynthesis were enriched in D. officinale under phosphorus stress.

Abstract

Phosphorus is integral to energy transfer and structural integrity in plants, which plays a significant role in regulating secondary metabolism. Notably, low phosphorus (LP) stress significantly improves dendrobine content in Dendrobium officinale, yet the molecular basis for this induction remains unclear. This study employed transcriptomic analysis to identify the differentially expressed genes (DEGs) related to the dendrobine biosynthesis under LP stress in D. officinale. 1,713, 222, 488, and 174 DEGs were up-regulated among the different phosphorus treatment groups, including the HP (high phosphorus) vs TP (total phosphorus), MP (medium phosphorus) vs TP, LP (low phosphorus) vs TP and NP (no phosphorus) vs. TP, respectively. In contrast, 1,855, 195, 432, and 120 DEGs exhibited a down-regulated expression pattern between each of them, respectively. Gene annotation in public datasets revealed that the DEGs related to phosphate transporter and alkaloid biosynthesis were enriched in D. officinale. By co-expression analysis, 10 phosphorus transport-related transcription factors (TFs) and 21 TFs associated with dendrobine biosynthesis were mined from the D. officinale transcriptome. These above findings provide many candidate TFs related to dendrobines biosynthesis and new insights into dissecting the potential molecular mechanism on regulating dendrobine biosynthesis under LP stress in D. officinale.

Figures

Effects of phosphorus deficiency on the dendrobine accumulation in D. officinale. TP, HP, MP, LP and NP denote total, high, medium, low, and no phosphorus, respectively. Data are presented as mean ± SD (n = 3). Significant differences were determined by Student’s t-test (*P < 0.05; P < 0.01).

Transcriptome sequencing and gene expression analysis. (A) Functional annotation of unigenes; (B) Number of DEGs in reponse to the induction of phosphorus deficiency in four comparative combinations. TP, HP, MP, LP, and NP denote total, high, medium, low and no phosphorus, respectively.

GO and KEGG pathway enrichment analysis of the DEGs in four comparative groups. (A-D) GO enrichment analysis of the DEGs. The visualization presents the significantly enriched GO semantic categories (y-axis) against the magnitude of their enrichment (x-axis), which is expressed as -log10(P-adjust). Therefore, the length of each bar directly reflects the statistical significance of the enrichment for that particular GO term. (E-H) KEGG pathway enrichment analysis of the DEGs is summarized in this visualization. Pathways are arranged along the y-axis. For each pathway, its enrichment significance (rich factor) is indicated on the x-axis, and the scale of the associated gene set is represented by the size of the circle. Log2FC|≥1 and an adjusted P-value<0.05 was defined as the cutoff to select the DEGs. TP, HP, MP, LP, and NP denote total, high, medium, low and no phosphorus, respectively.

Construction of co-expression network to identify the candidate TFs involved in LP stress and heatmap of differentially expressed genes related to phosphorus (P)-transport in response to LP stress. (A) Phosphorus transporter gene expression profiles. (B) Network analysis of transcriptional regulators and phosphorus transporter genes. PHT, phosphate transporter; PHO, phosphate; PHR, phosphate starvation response; PAP, purple acid phosphatase. TP, HP, MP, LP and NP denote total, high, medium, low and no phosphorus, respectively.

Candidate TF mining via co−expression network and expression analysis of DEGs involved in dendrobine biosynthesis. (A) Phosphorus−induced expression profiles of dendrobine synthase genes. (B) Construction of a co−expression network between candidate TFs and dendrobine biosynthetic genes. TP, HP, MP, LP and NP denote total, high, medium, low, and no phosphorus, respectively.

Experimental validation of the candidate TFs co-regulating LP stress and dendrobine biosynthesis. (A) Venn diagram of TFs from two co-expression networks. (B) qRT−PCR validation of the seven candidate TFs in (D) officinale, compared with transcriptome−derived expression levels. The TP group was taken as the reference. Data are presented as mean ± SD (n = 3). Significant differences were determined by Student’s t-test (*P < 0.05; P < 0.01). TP, HP, MP, LP and NP denote total, high, medium, low, and no phosphorus, respectively.

Tables

Table 1. Summary of the D. officinale transcriptome.

Sample	Raw reads	Clean reads	Clean bases	Q30 (%)	GC content (%)
TP_1	49,300,714	48,913,774	7,306,208,623	96.74	45.35
TP_2	41,807,246	41,474,796	6,202,588,927	96.62	45.2
TP_3	46,143,358	45,731,912	6,824,993,444	96.67	45.16
HP_1	47,627,626	47,251,262	7,047,670,065	96.71	44.93
HP_2	45,157,276	44,782,028	6,691,124,939	96.5	45
HP_3	41,758,584	41,427,312	6,189,244,159	96.52	45.04
MP_1	72,159,922	71,501,996	10,663,062,508	96.63	45.19
MP_2	74,014,092	73,419,776	10,969,199,898	96.62	45.16
MP_3	47,542,948	47,016,608	6,999,118,657	96.62	46.59
LP_1	43,111,822	42,733,908	6,380,666,647	96.67	45.12
LP_2	42,019,830	41,663,600	6,217,752,257	96.64	45.24
LP_3	43,423,554	43040112	6,427,049,986	96.66	45.24
NP_1	47,766,948	47,381,206	7,079,937,680	96.69	45.4
NP_2	44,193,734	43,807,358	6,565,192,759	96.69	44.97
NP_3	45,707,294	45,354,450	6,804,259,651	96.55	45.07
Average	48,782,330	48,366,673	7,224,538,013	96.63	45.24

Full text

Introduction

The medicinal plant Dendrobium officinale Kimura & Migo (commonly known as “Tiepi Shihu” in China) is a perennial epiphytic herb of the Orchidaceae family and the Dendrobium genus. It is a rare and valuable medicinal plant in traditional Chinese medicine (Cakova et al., 2017; Wang, 2021; Li et al., 2025). D. officinale contains various bioactive compounds, with polysaccharides and dendrobines as the major criteria for quality assessment in Chinese pharmacopoeia (Cakova et al., 2017; Wang, 2021; Li et al., 2024). Among them, the dendrobines present in D. officinale exhibit a variety of therapeutic effects against cancer, cardiovascular diseases, and gastrointestinal disorders. Additionally, they possess analgesic and antipyretic properties (Inubushi et al., 1963; Li et al., 2017; Liu et al., 2020b; Zhang et al., 2022; Huang et al., 2023; Okoro et al., 2023; Yang et al., 2023; Guo et al., 2024; Huang et al., 2024). Many studies have validated that environmental factors significantly influence dendrobine levels in D. officinale. Notably, phosphorus has been identified as a confirmed regulator of secondary metabolite biosynthesis in medicinal plants (Zuo et al., 2020; Chen et al., 2025).

Dendrobine is a picrotoxane-type sesquiterpene alkaloid whose biosynthesis follows a typical sesquiterpene pathway (Gong et al., 2021; Li et al., 2017). Its building blocks, isopentenyl diphosphate (IPP) and dimethylallyl diphosphate (DMAPP), are derived from both the mevalonate (MVA) pathway (starting from acetyl−CoA) and the methylerythritol phosphate (MEP) pathway (starting from glyceraldehyde−3−phosphate and pyruvate) (Shi et al., 2016; Huang et al., 2023). IPP and DMAPP are then further converted by the enzymes of geranyl diphosphate synthase (GPPS) and farnesyl diphosphate synthase (FPPS) to generate farnesyl pyrophosphate (FPP). In the downstream pathway of dendrobine biosynthesis, FPP is used as the substrate and undergoes two catalytic reactions by multiple terpene synthases (TPSs) to form the molecular skeleton of copacamphane (Gong et al., 2021). This skeleton is then subjected to multiple catalytic reactions by cytochrome P450 (CYP450) enzymes (including CYP71D55, CYP94C1 and other unknown CYP450 members) to produce picrotoxane-lactone, which is the basic skeleton molecule of dendrobine (Li et al., 2017). This basic skeleton molecule subsequently undergoes amination, methylation, and cyclization-decarboxylation reactions to ultimately yield dendrobine.

Phosphorus is a key nutrient element that regulates both primary and secondary metabolism in plants. Notably, phosphorus has been confirmed as a regulatory factor for the biosynthesis of secondary metabolites in medicinal plants. Low-phosphorus (LP) stress can promote the biosynthesis of phenolic acids in Salvia miltiorrhiza (Hao et al., 2020) and enhance the accumulation of tanshinones in its hairy roots (Zheng et al., 2023), while it inhibits alkaloid biosynthesis in Anisodus tanguticus (Zhang et al., 2025). Similarly, in D. officinale, applying LP stress within a comparable concentration range (typically spanning from deficiency to sufficiency, e.g., 0 to 2.5 mM) is known to increase alkaloid content (Liu et al., 2021). However, the transcriptional regulatory network, especially the key transcription factors, remains largely unexplored in D. officinale.

To dissect such complex regulatory mechanisms, identifying key transcription factors (TFs) is essential. In medicinal plants, transcription factors play a fundamental role in regulating secondary metabolism by controlling the transcription of biosynthetic genes. In Camptotheca acuminata, OpWRKY1 directly downregulates the expression of cytochrome P450 reductase (CYP) gene in hairy roots, thereby inhibiting the camptothecin (CPT) biosynthesis (Xu et al., 2020). OpERF2, when suppressed via RNA interference (RNAi), leads to a reduced expression pattern of those genes in MEP and secologanin-strictosidine pathways, indicating that OpERF2 promotes CPT biosynthesis (Udomsom et al., 2016). In Catharanthus roseus, CrWRKY1 positively regulates the biosynthesis of bisindole alkaloids by activating the TDC gene through binding to the W-box element in its promoter, while simultaneously suppressing the expression of several terpenoid indole alkaloids (TIA) biosynthesis activators, such as ORCA2, ORCA3, and CrMYC2 (Suttipanta et al., 2011). In Anisodus acutangulus, AaWRKY11 activates the expression of hyoscyamine 6β-hydroxylase (H6H1) gene, leading to tropane alkaloids accumulation in A. acutangulus (Zhou et al., 2024). These cases establish TFs as central molecular switches connecting environmental cues to metabolic outputs. Nevertheless, the specific TFs that mediate the response to LP stress and concurrently regulate dendrobine biosynthesis in D. officinale have not to be systematically identified and characterized.

To address this, the present study employed a phosphorus gradient (0–2.5 mM KH_2_PO_4_) centered around the standard MS medium concentration. An integrated approach of transcriptome sequencing and gene co−expression network analysis was implemented. First, transcription factors (TFs) responsive to low−phosphorus (LP) stress were identified. Subsequently, candidate TFs potentially coregulating phosphorus−signaling adaptation and dendrobine biosynthesis were mined. Finally, the expression patterns of key candidates were validated via qRT−PCR detection. This work provides a transcriptional regulatory framework for LP−enhanced dendrobine accumulation and offers genetic targets for breeding phosphorus−efficient D. officinale varieties.

Materials and methods

Experimental design

For RNA-sequencing and dendrobine quantification, ten-month-old sterile D. officinale seedlings were used. Seedlings were grouped (ten per replicate) and treated with a gradient of inorganic phosphorus (Pi) concentrations. This was achieved by supplementing phosphorus-free Murashige and Skoog (MS) medium with KH_2_PO_4_ to the following final concentrations: 1.25 mM (total phosphorus, TP), 2.5 mM (high phosphorus, HP), 0.625 mM (medium phosphorus, MP), 0.0625 mM (low phosphorus, LP), and 0 mM (no phosphorus, NP). The TP treatment group served as the mock group. All plants were cultivated in a greenhouse for 40 days under controlled conditions: 25 °C, 60% relative humidity, and a 12-h light (200 μmol·m^-^²·s^-^¹)/12-h dark photoperiod. After the treatment period, the seedlings were harvested for subsequent transcriptome analysis and dendrobine content measurement (Ren et al., 2022; Liu et al., 2022). All collected samples were snap-frozen and stored at -80 °C, with three biological replicates used for all experiments.

Determination of dendrobine content

Total dendrobine content was determined according to a previous method (Wang et al., 2016). Each 0.5 g sample of dried stem was finely powdered and poured into a distillation flask. The sample was moistened with 10% ammonia solution and kept for 30 minutes. Subsequently, 10 mL of chloroform was added into sample, and the distillation flask was initially weighed. Next, the dendrobine was extracted for 2 h with a condensation reflux device at 65 °C. After cooling, an appropriate volume of chloroform was added into the crude extract to restore the total weight to its initial value. The crude extracts were collected and filtered. Finally, 2 mL of filtered extract was mixed with 8 mL of chloroform to prepare the final tested solution, which was used to measure the absorbance value at a wavelength of 620 nm using an ultraviolet-visible spectrophotometry device. The total dendrobine content in each D. officinale sample was calculated based on the standard curve.

RNA extraction and Illumina sequencing

D. officinale seedlings subjected to five Pi treatments were collected. Following extraction of total RNA with TRIzol^®^ Reagent and quality confirmation (Nanodrop 2000 and gel electrophoresis), mRNA was purified using Oligo (dT) magnetic beads. Subsequently, a library was prepared from 1 µg of mRNA. Paired-end sequencing was conducted on the Illumina NovaSeq X Plus platform by Shanghai Majorbio Bio-pharm Technology Co., Ltd (Zhang et al., 2024).

De novo assembly and functional annotation

To ensure data quality, raw sequencing reads were processed with fastp (Version 1.0.1) with default parameters. to remove adapter sequences and filter low-quality reads (Chen et al., 2018). The resulting clean reads were aligned to the reference genome (ASM160598v2) in orientation mode using HISAT2 software (Version 2.2.1) to generate mapped reads for subsequent transcript assembly and expression level calculation (Kim et al., 2015; Zhang et al., 2016). Functional annotation of the assembled transcripts was then performed by scanning six public databases (EggNOG, Swiss-Prot, GO, Pfam, KEGG, and NR) with Diamond and HMMER, using an E-value cutoff of 1×10^-5^ (Cheng et al., 2025a). The highest-scoring annotation for each gene was retained to compile the final annotation list.

Identification of DEGs

Differentially expressed genes (DEGs) were identified between the mock (TP) group and each treatment group (HP, MP, LP, and NP). Expression levels were quantified as FPKM (Fragments Per Kilobase of transcript per Million mapped reads) using RSEM software (Version 1.3.3) (Pertea et al., 2015). Differential expression analysis was performed with DESeq2. Genes with an absolute log_2_ fold change (|log_2_FC|)≥1 and a false discovery rate (FDR) < 0.05 were identified as significant DEGs (Love et al., 2014).

Mining the candidate TFs by weighted gene co-expression network analysis

To identify candidate transcription factors involved in modulating dendrobine biosynthesis in response to low-phosphorus (LP) stress, expression correlation analysis was performed. Spearman’s rank correlation was used to calculate pairwise associations between genes, with an absolute correlation coefficient threshold of |r| ≥ 0.9. Correlations were considered significant if the adjusted P-value (padj) was < 0.05 after Benjamini–Hochberg correction. The expression profile of candidate TFs was presented in heatmaps generated using the Cytoscape software.

Quantitative real-time quantitative PCR

For experimental verification of transcriptome data, qRT−PCR assays were performed. First−strand cDNA was generated from 100 ng total RNA. The Applied Biosystems 7500 system and Taq Pro Universal SYBR qPCR Master Mix (Vazyme, China) were employed, with primers listed in Supplementary Table S1. The 10−µL reaction volume comprised 1 µL cDNA, 0.5 µL of each primer, 5 µL 2× Master Mix, and 3 µL ddH_2_O. Using the 2^−ΔΔCT^ method and DoActin as the reference gene, relative expression of candidate gene was derived. All data were obtained from 3 independent biological replicates.

Results

Low Pi improves dendrobine biosynthesis in D. officinale

To examine how low Pi affects dendrobine biosynthesis in D. officinale, the total content of dendrobine in each D. officinale plantlet treated with five different Pi concentrations was detected by spectrophotometry, respectively. Under LP treatment, the dendrobine content elevates to 0.07% of the total dry weight, representing a 1.85-fold increase relative to the Mock control (TP) (Figure 1). It implies that low Pi can improve the dendrobine biosynthesis in D. officinale.

*Effects of phosphorus deficiency on the dendrobine accumulation in D. officinale. TP, HP, MP, LP and NP denote total, high, medium, low, and no phosphorus, respectively. Data are presented as mean ± SD (n = 3). Significant differences were determined by Student’s t-test (*P < 0.05; *P < 0.01).

Transcriptome sequencing and de novo assembly of D. officinale transcripts

Through transcriptome sequencing, a total of 731 million raw reads were obtained from the 15 samples, with the average clean data more than 6.19 Gb (Supplementary Table S2). The GC content of each sample ranged from 44.93% to 46.59%, and the Q30 of clean reads was greater than 96.5% (Table 1). The high-quality reads from every sample were aligned against the D. officinale reference genome (ASM160598v2), yielding alignment rates between 86.75% and 89.87% (Supplementary Table S3). Evaluation of sequencing coverage and transcript length distribution indicates that the quality of assembled transcriptome data are sufficient for further data mining (Supplementary Figure S1).

Transcriptome sequencing and gene expression analysis

By homologous sequence similarity alignment, a total of 25,286 coding genes were successfully annotated in the EggNOG, Swiss-Prot, GO, Pfam, KEGG, and NR databases, accounting for 96.56% of the total 26186 splicing unigenes. Among the successfully annotated genes, 20,959 unigenes (84.89%) were annotated in the EggNOG database, 18,294 unigenes (74.07%) in the Swiss-Prot database, 18,829 unigenes (78.72%) in the GO database, 18,289 unigenes (74.02%) in the Pfam database, 9,417 unigenes (37.74%) in the KEGG database, and 22,712 unigenes (96.56%) in the NR database (Figure 2A; Supplementary Table S2). Among them, 7,744 unigenes were matched to the known sequences in multiple databases.

Transcriptome sequencing and gene expression analysis. (A) Functional annotation of unigenes; (B) Number of DEGs in reponse to the induction of phosphorus deficiency in four comparative combinations. TP, HP, MP, LP, and NP denote total, high, medium, low and no phosphorus, respectively.

Through differentially expressed gene (DEG) analysis, the DEGs in response to phosphorus deficiency in D. officinale were identified (Figure 2B; Supplementary Table S4). In total, 4,049 DEGs were identified. Among these, 1,713, 222, 488 and 174 DEGs were up-regulated in the HP vs. TP, MP vs. TP, LP vs. TP and NP vs. TP comparisons, respectively.

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis were introduced to predict the potential function of DEGs. GO enrichment analysis was performed to annotate the DEGs in the LP vs TP group. In the cellular component (CC) category, extracellular region was the most enriched subcategories. In the biological process (BP) category, tyrosyl-tRNA aminoacylation peaked at all the function classifications. In the molecular function (MF) category, copper ion binding, heme binding, small molecule binding and oxidoreductase activity were the top four enrichments (Figure 3; Supplementary Table S5-8).

GO and KEGG pathway enrichment analysis of the DEGs in four comparative groups. (A-D) GO enrichment analysis of the DEGs. The visualization presents the significantly enriched GO semantic categories (y-axis) against the magnitude of their enrichment (x-axis), which is expressed as -log10(P-adjust). Therefore, the length of each bar directly reflects the statistical significance of the enrichment for that particular GO term. (E-H) KEGG pathway enrichment analysis of the DEGs is summarized in this visualization. Pathways are arranged along the y-axis. For each pathway, its enrichment significance (rich factor) is indicated on the x-axis, and the scale of the associated gene set is represented by the size of the circle. Log2FC|≥1 and an adjusted P-value<0.05 was defined as the cutoff to select the DEGs. TP, HP, MP, LP, and NP denote total, high, medium, low and no phosphorus, respectively.

By KEGG enrichment analysis, cutin, suberin, and wax biosynthesis, phenylpropanoid biosynthesis, β-alanine metabolism, pentose and glucuronate interconversions represented the top four pathways in the HP vs. TP comparative group (Supplementary Table S9). In the NP vs. TP comparative group, the DEGs were preferentially enriched in the pathways of cutin, suberin, and wax biosynthesis, tyrosine metabolism, fructose and mannose metabolism, and isoquinoline alkaloid biosynthesis (Supplementary Table S10). In the MP vs. TP group, pathways were predominantly enriched in brassinosteroid biosynthesis, glutathione metabolism, photosynthesis, and sesquiterpenoid and triterpenoid biosynthesis (Supplementary Table S11). The LP vs. TP groups exhibited significant pathway enrichment in peroxisome, cutin, suberin and wax biosynthesis, and fatty acid degradation. Notably, the biosynthesis of tropane, piperidine, and pyridine alkaloids (terpenoid alkaloids) exhibited higher rich factors compared to other groups. This finding indicates that under phosphorus deficiency, enzyme-coding genes associated with terpenoid alkaloid synthesis become substantially enriched (Supplementary Table S12).

Mining of the candidate transcription factors in response to low phosphorus stress

To identify candidate transcription factors (TFs) responsive to low-phosphorus (LP) stress, a gene co-expression network was constructed. As the core drivers of phosphorus uptake and allocation, key phosphorus transport and signaling genes serve as ideal entry points for discovering their upstream regulators. Based on functional annotation,six such “bait” genes were selected, each representing a pivotal component of the conserved plant phosphate response network: DoPHT2 (a phosphate transporter central to Pi uptake, LOC110101372), DoPHO1;3 (LOC110106399) and DoPHO1;4 (mediators of phosphate allocation from roots to shoots, LOC110098888), DoPHR1 (LOC110110494) and DoPHR2 (master transcription factors of phosphate starvation response, LOC110113993), and DoPAP (a purple acid phosphatase involved in phosphorus deficiency adaptation, LOC110111096) (Figure 4A; Supplementary Table S13). These six bait genes, together with 194 TFs identified from the DEGs, were then subjected to co-expression correlation analysis. A co-expression network containing 15 nodes and 25 links was constructed with the Pearson correlation coefficients >0.90 as the cutoff. In total, 10 candidate TFs were identified that exhibited strong correlations with five of the above six P-transport genes (DoPHO1;3, DoPHO1;4, DoPHR1, DoPHR2, and DoPAP). Among the 10 correlated TFs, DoNAC68 (gene-LOC110092932) and DoMYC3 (gene-LOC110099101) exhibited the highest correlations with the above five P-transport genes, followed by DoERF5 (LOC110093822), DoERF8 (LOC110095855), DobHLH57 (LOC110111619), DoERF6 (LOC110106428), DobZIP25 (LOC110116565), DoBEL1 (LOC110111122), DoHOX14 (LOC110107591) and DoZF-HD4 (LOC110109979), which were named based on their homologs in the genome of Arabidopsis thaliana (Figure 4B; Supplementary Table S14).

Construction of co-expression network to identify the candidate TFs involved in LP stress and heatmap of differentially expressed genes related to phosphorus (P)-transport in response to LP stress. (A) Phosphorus transporter gene expression profiles. (B) Network analysis of transcriptional regulators and phosphorus transporter genes. PHT, phosphate transporter; PHO, phosphate; PHR, phosphate starvation response; PAP, purple acid phosphatase. TP, HP, MP, LP and NP denote total, high, medium, low and no phosphorus, respectively.

Mining of the candidate transcription factors involved in dendrobine biosynthesis

To identify candidate TFs involved in regulating dendrobine biosynthesis, a co-expression network containing 26 nodes and 63 links was constructed using 194 selected TFs and 5 known dendrobine biosynthetic genes including DoDMAPP (LOC110095726), DoFPPS (LOC110096432), DoMVK (LOC110113415), DoDXR (LOC110096522) and DoCYP71D55 (LOC110114012). Pearson correlation coefficients >0.9 and an adjusted P-value <0.05 were used as cutoffs (Figure 5A; Supplementary Table S15). In this co-expression network, DoRAP2.4 (LOC110110379) and DoGRF10 (LOC110104491) exhibited the highest Pearson correlation coefficients reaching to 0.983, followed by DoFAR1 (LOC110105203), DoGRF6 (LOC110107122), DoNAC68 (LOC110092932), DoERF8 (LOC110095855), and DoBEL1 (LOC110111122) with the Pearson correlation coefficients greater than 0.95 (Figure 5B; Supplementary Table S16). These TFs were thought to be the potential regulators of dendrobine biosynthesis in D. officinale.

Candidate TF mining via co−expression network and expression analysis of DEGs involved in dendrobine biosynthesis. (A) Phosphorus−induced expression profiles of dendrobine synthase genes. (B) Construction of a co−expression network between candidate TFs and dendrobine biosynthetic genes. TP, HP, MP, LP and NP denote total, high, medium, low, and no phosphorus, respectively.

Experimental validation of differential expressed TFs by qRT-PCR analysis

To identify potential co-regulators involved in both dendrobine biosynthesis and LP stress response, this study conducted an integrative analysis combining 10 candidate TFs implicated in the regulation of dendrobine biosynthesis with 21 potential candidate TFs responsive to LP stress (Figure 6A). This analysis revealed seven overlapping candidate TFs: DoNAC68 (LOC110092932), DoERF5 (LOC110093822), DobHLH57 (LOC110111619), DobZIP25 (LOC110116565), DoBEL1 (LOC110111122), DoERF8 (LOC110095855) and DoERF6 (LOC110106428). The expression patterns of all seven candidate TFs were subsequently validated by qRT−PCR under a gradient of phosphorus concentrations. The results confirmed that their expression trends were highly consistent with the transcriptomic data, with correlation coefficients ranging from 0.92 to 0.99 (Figure 6B; Supplementary Figure S2). These results demonstrate the reliability of the transcriptomic data in mining the candidate TFs involved in modulating dendrobine biosynthesis under phosphate deficiency in D. officinale.

*Experimental validation of the candidate TFs co-regulating LP stress and dendrobine biosynthesis. (A) Venn diagram of TFs from two co-expression networks. (B) qRT−PCR validation of the seven candidate TFs in (D) officinale, compared with transcriptome−derived expression levels. The TP group was taken as the reference. Data are presented as mean ± SD (n = 3). Significant differences were determined by Student’s t-test (*P < 0.05; *P < 0.01). TP, HP, MP, LP and NP denote total, high, medium, low, and no phosphorus, respectively.

Discussion

Phosphorus acts as an essential nutrient in regulating plant growth, development, and metabolism (Gu et al., 2016). In wild environment, D. officinale prefers to grow as epiphytes on tree bark and rocks. This species characteristic enables the D. officinale to adapt to the poor nutrient environment (Xiang et al., 2016; Hou et al., 2017; Ren et al., 2020; He et al., 2024). These findings prompted an investigation into whether LP stress genuinely enhances the quality of D. officinale by promoting dendrobine accumulation in its stems (Figure 1). In Anisodus tanguticus, phosphorus deficiency enhances arbuscular mycorrhiza-mediated coordination of carbon-nitrogen-phosphorus metabolism, thereby promoting alkaloid accumulation (Zhang et al., 2025). LP stress has been validated to improve the tanshinone biosynthesis in S. miltiorrhiza hairy roots (Zheng et al., 2023). In summary, the judicious application of phosphorus during the production of traditional Chinese medicinal materials is critical, as it not only ensures product quality but also minimizes phosphate fertilizer usage, thereby alleviating environmental discharge pressures.

Co-expression network analysis is thought to be a reliable strategy for identifying candidate TFs in plant species (Liu et al., 2022; Han et al., 2024; Zhang et al., 2024). For example, integration of transcriptome data and co-expression network analysis identified the core transcription factor MabHLH355 in banana, which regulates reactive oxygen species (ROS) scavenging under cold stress ( (Lin et al., 2024). Under low−nitrogen stress, six transcription factors implicated in the biosynthesis of tanshinones and salvianolic acids were identified in S. miltiorrhiza (Cheng et al., 2025b). By gene co-expression network analysis, 4 important TFs including PpZAT6, PpPMZ, PpMYB-like, and PpONAC077 in association with the terpenoid metabolites biosynthesis were identified in Poa pratensis (Poirier et al., 2022). In this study, based on transcriptomic sequencing data in combination with co-expression correlation analyses, 21 TFs related to dendrobine biosynthesis and 10 TFs involved in phosphorus absorption and transport were identified, respectively. By integrating the above two co-expression networks, seven TFs are identified: DoNAC68 (LOC110092932), DoERF5 (LOC110093822), DobHLH57 (LOC110111619), DobZIP25 (LOC110116565), DoBEL1 (LOC110111122), DoERF8 (LOC110095855) and DoERF6 (LOC110106428). As the homologous gene of DoERF8 (LOC110095855), AtERF8 functions as a transcriptional repressor that modulates light-responsive anthocyanin biosynthesis in A. thaliana (Koyama and Sato, 2018; Isidra-Arellano et al., 2021; Ren et al., 2022). DoGRF6 (LOC110107122) exhibits a high sequence similarity with ClNAC68 in Citrullus lanatus, which has been verified to regulate sugar accumulation during maturation stage in watermelon (Pertea et al., 2015). DoBEL1 (LOC110111122) is a homolog of MaBEL1 in banana, and it was verified to activate the expression of genes involved in starch degradation thereby promoting fruit ripening (Song et al., 2023). Although the functions of these genes still require experimental validation in D. officinale, they provide critical clues for elucidating the underlying biological functions of the seven candidate transcription factors identified in this study.

Phosphate transporters (PHTs), acting as the high-affinity phosphate transporters, localize to the plasma membrane and they determine phosphorus uptake and redistribution via plant root (Paz-Ares et al., 2022). Among them, PHT1;1 and PHT1;4 exhibit a root-specific upregulation expression pattern in response to LP stress (Tao et al., 2024). PHO1 members were generally thought to mediate phosphorus transport from roots to shoots in plant species (Isidra-Arellano et al., 2021). Therefore, to identify candidate TFs responsive to LP stress, DoPHO1;3, DoPHO1;4, DoPHR1, DoPHR2, and DoPAP were employed as ‘baits’ to construct a co-expression network to prey the candidate TFs among the 194 differentially expressed TFs. Ultimately, 10 TFs including DoNAC68, DoMYC3, DoERF5, DoERF8, DobHLH57, DoERF6, DobZIP25, DoBEL1, DoHOX14, and DoZF-HD4 were strongly correlated (r > 0.9) with the above 6 phosphorus transport-related genes (Figure 4B; Supplementary Table S14). In apple, MdMYB306 interacts with MdMYB17 and MdbHLH33 to inhibit anthocyanin synthesis (Wang et al., 2022). The CUL3-BPM E3 ubiquitin ligase AtPUB10 has been reported to regulate jasmon acid signaling by modulating the stability of AtMYC2, AtMYC3, and AtMYC4 proteins in A.thaliana (Chico et al., 2020). In Catharanthus roseus, CrERF5, an AP2/ERF transcription factor, modulates the production of bisindole alkaloids through activating the expression of tryptophan decarboxylase (TDC) gene (Teng et al., 2025). AtbHLH57 has been validated to interact with AtODR1 (REVERSAL OF RDO5 1) to regulate ABA synthesis in A.thaliana (Liu et al., 2020a). The ABA pathway integrates environmental signals through a signal transduction network to regulate metabolic processes within plants, enabling them to adapt to adverse conditions. Studies on S. miltiorrhiza reveal that plants adapt to environmental changes by improving anthocyanin synthesis under LP stress (Tao et al., 2024). Therefore, elucidating how these 10 TFs participate in phosphate uptake and secondary metabolite synthesis in D. officinale, with reference to previous reports, is worth exploring.

Conclusions

The present study has employed transcriptome sequencing data of D. officinale in combination with gene co-expression network analysis to identify candidate TFs involved in regulation of dendrobine biosynthesis under LP stress. In total, 7 candidate TFs are mined. qRT-PCR analysis showed that their expression levels are drastically induced under LP stress compared to the Mock (TP), and this induction trend is consistent with the transcriptome sequencing data. In general, this study offers many valuable candidate TFs for elucidating the underlying molecular mechanism on modulating dendrobine biosynthesis under LP stress. These findings will promote scientific production with the condition of low phosphorus consumption in D. officinale.