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₂PO₄ 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.

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.

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).

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.

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₂ fold change (|log₂FC|)≥1 and a false discovery rate (FDR) < 0.05 were identified as significant DEGs (Love et al., 2014).

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.

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₂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.

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.

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).

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.

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).

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).

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).

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.

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.