Epigenetic priming of terpenoid biosynthesis drives cultivar-specific efficient agarwood formation in Aquilaria sinensis

Aquilaria spp. are aromatic tree species in the Thymelaeaceae that constitute the principal botanical source of agarwood, a precious resinous material with longstanding applications in perfumery, ritual incense, and traditional pharmacopeias^[1,2]. Although agarwood is a generic term for resinous products across Asia, the only officially recognized botanical origin in China is Aquilaria sinensis, a species widely distributed throughout Hainan and other southern provinces^[3]. Global demand for agarwood has driven the large-scale plantation of A. sinensis across Southern China. Consequently, improving yield has become an urgent practical goal^[4].

The pathological resin deposition characteristic of agarwood occurs exclusively under stress conditions, such as physical wounding, microbial infection, or other environmental insults^[5]. To meet the rising market demand, artificial induction techniques must be employed to consistently increase yield. Over the past decade, various artificial induction methods have been successively developed and improved, including mechanical wounding, chemical injection, fungal inoculation, and whole-tree drilling systems^[6−9]. Alongside induction technology, cultivar development has produced a diverse collection of A. sinensis lines with varying yield potential. Among these, a category of highly efficient agarwood-inducing cultivars has been widely propagated in recent years^[10−12]. Notably, in recent studies^[13,14], these highly efficient strains are explicitly designated as 'Chi-Nan' or 'Qi-Nan'.

Molecular biology methods are widely applied to authenticate these elite cultivars^[15]. Despite their extremely high efficiency in forming agarwood, population genomics indicates that the genetic differentiation of these cultivars from ordinary A. sinensis (often referred to as Baimu, BM) is exceptionally low^[12]. However, they typically differ substantially in chemical composition from traditional Baimu agarwood. The pharmacopoeia standard specifically mandates the presence of agarotetrol as a diagnostic marker for authentic medicinal agarwood. While this compound is readily generated in traditional Baimu agarwood, it is often absent or present at trace levels in resins produced by typical highly efficient cultivars^[13,14]. Consequently, they offer high yields but have limited medicinal applicability.

To resolve this conflict between yield and medicinal composition, we identified an elite cultivar designated 'Shuxinyou' (SXY) as a unique member among cultivars with high efficiency in forming agarwood. SXY accumulates resin at rates comparable to those of other high-yield cultivars. Notably, it also retains agarotetrol at levels consistent with traditional agarwood chemical profiles. This combination has not been previously reported in this cultivar class. Compared to ordinary Baimu (BM), SXY produces nearly four times the resin yield of BM under identical injury conditions. However, the regulatory basis of SXY's yield superiority remains unknown. The minimal DNA-level differentiation between elite and common cultivars suggests that non-genetic regulatory mechanisms (such as chromatin accessibility) may underlie this phenotypic difference. In other plant systems, constitutive differences in chromatin accessibility have been shown to prime transcription in response to stress, facilitating stronger gene activation upon injury or elicitation^[16−18]. By employing a grafted system to strictly control for environmental variables, we integrated ATAC-seq and RNA-seq analyses to characterize chromatin accessibility and transcriptional responses under standardized injury conditions. We show that SXY maintains a more accessible chromatin landscape together with stronger induction of terpenoid-pathway genes. These results raise the possibility that chromatin accessibility, rather than DNA sequence variation, underlies the yield difference between BM and SXY, with practical implications for cultivar screening.

Grafted A. sinensis plants were used in this study. One-year-old Baimu (BM) seedlings were used as rootstocks. Scions (branches) were collected from an adult Shuxinyou (SXY) mother tree located in Wenru Town, Chengmai County, Hainan Province, China. Grafting was performed by inserting SXY scions onto the BM rootstocks. The grafted plants were maintained under standard nursery conditions for 6 months prior to induction. At the time of mechanical injury, shoots of comparable stem diameter (5–8 mm) and physiological state were explicitly selected for both the BM branch and the SXY branch to ensure morphological consistency. The shoots emerging from the rootstock and the scion were designated as the BM branch and the SXY branch, respectively.

Physical injury induction was performed using a standard trunk drilling method^[6]. Holes were drilled with a 2.3 mm bit to a depth of 10 mm, with four holes per plant and 1.5 cm spacing between adjacent holes. No liquid or chemical agent was injected. Sampling time points were set at 0 d (CK) and at 15 and 30 d post-injury. For transcriptome analysis, three biological replicates were collected at each time point (CK, 15 d, and 30 d). ATAC-seq sampling was performed at the baseline (0 d, CK) time point. Detailed sample information, including biological replicates, sampling time points, and experimental groups, is provided in Supplementary Table S1.

Phenotyping and quantitative analysis of alcohol-soluble extractive content

Phenotypic characteristics of the induced sites were recorded photographically. Alcohol-soluble extractive content was quantified by ethanol reflux extraction. Approximately 0.5 g of agarwood powder ($ {m}_{0} $) was extracted in 50 mL of 95% ethanol under reflux for 1 h. After filtration, a 25 mL aliquot of filtrate was evaporated to dryness and dried at 105 °C to constant weight.

The alcohol-soluble extractive content (X) was calculated as:

where, $ {{m}}_{1} $ is mass of the empty evaporating dish (g); $ {{m}}_{2} $ is mass of the evaporating dish and the dried extract (g); m₀ is the mass of the sample (g); W is the moisture content of the sample (%). V is total extract volume (50 mL), and V_a is aliquot volume used for drying (25 mL).

For each cultivar and time point, three independent plants were measured (biological replicates, n = 3). For each biological sample, ethanol extraction was performed in technical duplicate and averaged. Statistical significance was determined using a two-sided Student's t-test (p < 0.05).

RNA-seq sequencing and differential expression analysis

Total RNA was extracted using the Trizol method. RNA integrity was verified with an Agilent Bioanalyzer (RIN ≥ 7.0) and OD₂₆₀/₂₈₀ ratios between 1.8 and 2.1. PE150 libraries were sequenced on an Illumina NovaSeq 6000 platform. Raw data were filtered using fastp (v0.23.0)^[19]. Clean reads were aligned to the A. sinensis reference genome^[20,21] using HISAT2 (v2.2.1). Gene-level counts were quantified using featureCounts (v2.0.1).

Differential expression analysis was performed in DESeq2^[22] using raw gene counts. To characterize the overall induction response rather than time-resolved dynamics, samples from 15 and 30 d were pooled into a single 'Treat' group and compared against the CK group within each cultivar. p-values were adjusted using the Benjamini–Hochberg procedure, and significant differentially expressed genes (DEGs) were defined as p-value ≤ 0.05 and |log₂FoldChange| ≥ 1. Raw count matrices and sample-level FPKM values (for visualization) are provided in Supplementary Tables S2 and S3.

ATAC-seq library construction and analysis

Nuclei were isolated from uninduced stem tissues harvested at 0 d (CK), and one ATAC-seq library was generated for each cultivar using a standard Tn5 transposase protocol^[23]. Libraries were purified using AMPure beads and sequenced on an Illumina NovaSeq (PE150). Data were processed by fastp (v0.20.0) for quality control and aligned to the reference genome^[20] using BWA (v0.7.12)^[24]. Mitochondrial and chloroplast reads were removed, and only high-quality unique mapping reads (MAPQ ≥ 13) were retained to prevent false positives. Standard ATAC-seq QC metrics (including fragment-size periodicity, TSS enrichment, and the fraction of reads in peaks) were evaluated and summarized in Supplementary Table S4. Given the single-library design per cultivar, differential peak analysis was strictly interpreted at a systemic pathway-level scale (e.g., broad accessibility vs transcription trends), rather than relying on the isolated statistical significance of any individual peak.

Peaks were called using MACS2 (v2.1.0) with parameters: -q 0.05 --call-summits --nomodel --shift -100 --extsize 200 --keep-dup all^[25]. To support integrative comparison with the RNA-seq dataset, cultivar-biased peak sets were defined using a binary overlap strategy in BEDTools^[26]. Peaks overlapping (> 1 bp) only with SXY peaks were defined as SXY-up, and peaks overlapping only with BM peaks were defined as BM-up. Given the lack of biological replication in the ATAC-seq dataset, these categories were used for descriptive comparison rather than as statistically validated differential peaks.

Peak annotation and gene integration

Peaks were annotated to genomic features and associated genes using ChIPseeker^[27]. Promoter regions were defined as TSS ± 3 kb. For integration with RNA-seq, cultivar-biased promoter peaks (SXY-up or BM-up) were linked to genes based on promoter overlap. We defined cultivar-specific peak-linked genes as those associated with peaks unique to one cultivar in the baseline ATAC-seq comparison. These accessibility categories were then related to transcriptional changes observed after injury. Integrated analysis related accessibility state (SXY-up vs BM-up) to transcriptional changes (RNA log₂FC). A quadrant analysis classified genes based on the concordance between accessibility state and expression change.

Pathway and network analysis

Key genes in the MVA/MEP and downstream terpenoid pathways were curated from genome annotation^[28]. Their expression profiles were integrated with ATAC signals to summarize pathway-level patterns. Where relevant, metabolite-level interpretations are discussed in light of prior metabolomic studies in Aquilaria^[29].

Functional enrichment analyses were performed on DEGs to summarize pathway-level trends. Over-representation analysis was performed against the expressed-gene background, and p-values were adjusted using the Benjamini–Hochberg procedure. In addition, gene set enrichment analysis (GSEA) was conducted using preranked gene lists derived from RNA-seq differential expression results to assess coordinated shifts in terpenoid-related gene sets. Volcano plots were generated directly from DESeq2 outputs. Hierarchical clustering and heatmaps were generated in R using sample-level FPKM values (Supplementary Table S3) for visualization.

For regulatory network construction, motif enrichment analysis was performed on SXY-up and BM-up peak sets using HOMER (v4.9.1)^[30]. TFs were prioritized using three criteria. First, TFs were significantly upregulated in SXY. Second, the corresponding TF family motif was enriched in SXY-up peaks. Third, the TF locus showed high promoter accessibility. Putative TF-target edges were defined when an enriched motif was present in a cultivar-biased promoter peak linked to the target gene. Edges were retained only for TFs that met the prioritization criteria. The resulting network is reported as a predicted regulatory model.

Statistics and visualization

Statistical analyses were performed in R. Alcohol-soluble extractive content differences were assessed using two-sided Student's t-tests (p < 0.05). RNA-seq differential expression was performed using DESeq2 with Benjamini–Hochberg adjustment for multiple testing (p-value). ATAC-seq peak calling used MACS2 with q-value thresholding, and motif enrichment significance was reported from HOMER outputs. Plots were generated using standard R visualization workflows. ATAC-seq cultivar comparisons (peak counts, genomic feature distributions, and peak-linked gene sets) were performed as descriptive analyses without statistical modeling of biological replication, consistent with the single-library design described in ATAC-seq library construction and analysis.

We evaluated BM and SXY responses to identical drilling injury utilizing a grafted-plant design (Fig. 1a). At 30 d post-injury, SXY developed distinct dark resinous zones, whereas BM exhibited only lighter discoloration (Fig. 1b). Consistent with this visual disparity, the accumulation of alcohol-soluble extractives in SXY was nearly four times higher than in BM (13.2% vs 3.6%, p < 0.05) (Fig. 1c). The SXY extracts retained the diagnostic chemical profile of authentic traditional agarwood, including the agarotetrol peak (Supplementary Fig. S1).

Figure 1. 

Detection of the agarotetrol peak in SXY resin extracts via HPLC. (a) Representative grafted A. sinensis plants used for induction. (b) Representative induced stems at 30 d post-injury, showing stronger resinous discoloration in the SXY branch than in the BM branch. (c) Alcohol-soluble extractive content at 0 and 30 d post-injury (mean ± SD, n = 3; two-sided Student's t-test, p < 0.05). (d) ATAC-seq mapping and unique mapping rates (0 d). (e) Total number of called peaks in BM and SXY (0 d). (f), (g) Genomic feature distribution of ATAC-seq peaks in (f) BM, and (g) SXY. (h) Metagene plot of ATAC-seq signal across gene bodies with transcription start site (TSS) and transcription end site (TES) landmarks (± 3 kb flanks).

We profiled chromatin accessibility in baseline stem tissues using ATAC-seq (Fig. 1d; Supplementary Table S4). Peak calling identified substantially more accessible regions in SXY (71,680) than in BM (51,489) (Fig. 1e−h). At baseline, the SXY genome exhibited a more broadly open conformation.

Cultivar-specific chromatin signatures favour transcriptional potential

SXY's accessible chromatin was also more heavily enriched at promoters (31.4% vs 26.8% in BM; Fig. 2a, b). A count-based Fisher's exact test yielded p = 1.48e−68, but because this comparison reflects peak-count asymmetry rather than replicate-aware inference, it is interpreted here as descriptive support only. Visualizing the signal around transcription start sites (TSS) likewise indicated a broader, more intense accessibility signature in SXY (Fig. 2c).

Figure 2. 

Landscape of accessible chromatin regions (ACRs). (a), (b) Distribution of peaks relative to genomic features in (a) BM, and (b) SXY. (c) Heatmap of ATAC-seq signal intensity within ± 3 kb around the TSS of all genes. (d) Venn diagram showing cultivar-specific (SXY-only, BM-only) and shared ACRs. (e) Distributions of peak length and peak intensity for BM and SXY.

When we filtered for cultivar-specific peaks (defined by binary peak overlap rather than statistical differential analysis (see ATAC-seq library construction and analysis), we found 5,355 genes with SXY-specific accessibility, compared to just 1,523 in BM (Fig. 2d). SXY thus maintains a much larger set of accessible chromatin regions (ACRs) (Fig. 2e).

Cultivar-specific transcriptional reprogramming

RNA-seq data indicated that the two cultivars mount highly divergent transcriptional programs following injury. BM exhibited a broad response comprising 2,653 differentially expressed genes (DEGs). In contrast, SXY's response was narrower and more targeted, involving 1,779 DEGs (Fig. 3a, b). Hierarchical clustering of these DEGs cleanly separated the two cultivars (Fig. 3c; Supplementary Tables S5, S6).

Figure 3. 

Gene set enrichment analysis (GSEA) of terpenoid-related gene sets. (a), (b) Volcano plots of differentially expressed genes (DEGs) in (a) BM, and (b) SXY (Treat vs CK; |log₂FC| ≥ 1, p-value ≤ 0.05). (c) Hierarchical clustering heatmap of all identified DEGs across cultivars and treatments.

Functional enrichment analysis further showed that wound-responsive genes in SXY were mainly associated with sesquiterpenoid biosynthesis and secondary metabolism, whereas BM showed enrichment for RNA processing and more general stress-related terms (Supplementary Fig. S2). This divergence was further supported by gene set enrichment analysis (GSEA), which showed positive enrichment of gene sets involved in terpenoid biosynthesis in SXY but not in BM (Supplementary Fig. S3).

Integration of chromatin and transcriptional landscapes

To relate baseline chromatin state to the injury-induced transcriptional response, we integrated the cumulative expression changes (pooled 15 d and 30 d RNA-seq) with the chromatin state captured at 0 d (prior to induction). To see if this accessibility predicts transcription, we classified the cultivar-specific peaks as SXY-up (27,353) or BM-up (10,980) (Fig. 4a) and plotted them against RNA-seq fold changes (Fig. 4b).

Figure 4. 

Integration of chromatin accessibility and gene expression. (a) Classification of cultivar-biased peaks (SXY-up vs BM-up). (b) Quadrant analysis integrating accessibility direction (SXY vs BM) and expression change (Treat vs CK). (c) Violin plot showing the distribution of expression bias (log₂(Treat/CK) in Treat samples) for genes associated with cultivar-specific or shared peaks.

The 'concordant upregulation' quadrant (genes with both open chromatin and increased expression after injury) was heavily skewed toward SXY. Genes linked to SXY-specific peaks generally demonstrated a corresponding positive shift in expression (Fig. 4c; Supplementary Tables S7, S8).

Coordinated activation of the terpenoid pathway

When we mapped these integrated data onto the pathway for terpenoids, several nodes that limit the rate of reaction exhibited clear activation specifically in the SXY cultivar (Fig. 5a). Expression profiles for representative genes in the terpenoid pathway at each time point are shown in Supplementary Fig. S4. For example, DXS (the MEP pathway entry enzyme) and IDI (the core isomerase) were strongly upregulated in SXY but not in BM.

Figure 5. 

Multi-omics evidence for the terpenoid biosynthetic pathway. (a) Simplified pathway map showing cultivar-biased expression changes at key nodes (e.g., DXS, HMGS, TPSa1). (b) Genome browser view of the TPSa1 locus (corresponding to AsTPS1 in the text; evm.model.Scaffold11.100), showing higher promoter accessibility in SXY alongside stronger expression after injury.

The AsTPS1 locus provides a particularly clear example of this coupling. BM and SXY exhibited opposing transcriptional differentiation: AsTPS1 was downregulated in BM (log₂FC = −1.13, p-value = 0.025), whereas SXY showed a positive log₂FC (1.08), although this change did not reach the DESeq2 significance threshold in Table S5 (p-value = 0.0837). Looking at the locus itself, we found a prominent ATAC-seq peak at the SXY promoter that was largely absent in BM (Fig. 5b). Here, localized promoter accessibility was associated with a more positive transcriptional response in SXY.

A candidate bHLH/bZIP regulatory module

We searched the SXY-specific peaks for enriched transcription factor binding sites and found strong signals for bHLH, bZIP, and MYB families (Fig. 6a). To prioritize candidates, we focused on TFs that were upregulated in SXY, associated with enriched motifs, and linked to accessible promoters (Fig. 6b).

Figure 6. 

Candidate transcription factors and a putative regulatory framework. (a) Motif enrichment analysis for SXY-up and BM-up peaks. (b) Integrated evidence matrix for candidate TFs, summarizing differential expression, promoter accessibility, and motif enrichment. (c) Putative regulatory network linking candidate TFs to terpenoid biosynthesis genes.

Two factors, BHLH137 and HYH (a HY5 homolog), emerged as candidates. Both were highly expressed in SXY, and enriched bHLH/bZIP motifs were present in accessible promoter regions associated with DXS, HMGS, and AsTPS1. Based on these integrated observations, we linked the prioritized TFs to candidate terpenoid pathway targets to construct a putative regulatory network (Fig. 6c; Supplementary Tables S9, S10). Each predicted edge is a testable hypothesis.

Under standardized induction conditions, SXY produced nearly four times the alcohol-soluble extractives of BM (13.2% vs 3.6%). SXY is therefore a useful model for understanding efficient agarwood formation. Given the minimal population differentiation between the two lineages and the rigorous control of environmental plasticity through our paired grafting design, this discrepancy in yield is consistent with underlying regulatory differences between the two cultivars.

Chromatin accessibility differences are associated with the yield disparity

Our multi-omics comparison suggests that chromatin accessibility differences are associated with the yield disparity between SXY and BM. Even at the uninduced baseline, SXY displayed a more globally open chromatin landscape than BM. In SXY, accessibility signals were also more frequently associated with promoter regions. In plant stress and defense responses, chromatin accessibility is closely linked to the ability of transcription factors to engage cis-regulatory elements^{[16,17,31,32]}. Broader chromatin accessibility in SXY may thus facilitate a stronger metabolic response after injury. However, because our ATAC-seq data used only one biological replicate for each cultivar, these findings should be seen as evidence for a candidate regulatory state linked to the accessibility of chromatin, rather than final proof of a stable epigenetic mechanism. We also cannot exclude the possibility that transcriptional differences drive chromatin remodeling rather than the reverse; future time-resolved ATAC-seq experiments would be needed to establish causality.

Chromatin differences are associated with terpenoid biosynthetic pathways

These chromatin differences are closely associated with key metabolic branches centered on terpenoid biosynthesis^[29,33]. In particular, SXY-biased accessible regions coincided with stronger transcriptional activation of multiple rate-limiting enzyme genes, including DXS, HMGS, and AsTPS1. The AsTPS1 locus provides a representative example, in which higher promoter accessibility in SXY was accompanied by a more positive transcriptional response after injury. Motif enrichment and network inference further suggested that these accessible promoters may be linked to transcription factor families such as bHLH and bZIP, including candidates such as BHLH137 and HYH^[34]. Baseline chromatin accessibility and transcription factor availability may thus jointly shape pathway activation in SXY.

Implications for understanding agarwood synthesis and renewing germplasm Resources

These results help explain why agarwood formation varies by cultivar. Furthermore, they provide a new basis for evaluating future germplasm. Rather than focusing only on acute transcriptional responses or single-gene expression^[4,35−37], our results suggest that pre-existing differences in chromatin accessibility may influence the magnitude of wound-induced metabolic activation. From an applied perspective, this raises the possibility that accessibility states at key regulatory loci, such as the AsTPS1 promoter, could be explored as candidate markers in future screening or breeding efforts. However, such applications remain prospective and will require validation across additional elite lines and broader germplasm resources.

Limitations and future directions

This study is limited by the single-replicate nature of ATAC-seq. This limits variance estimation and constrains inference on differential accessibility. In addition, conclusions are derived from one elite cultivar (SXY) and one common cultivar (BM). Broader sampling across independent SXY lineages and additional cultivars will be necessary to assess generality. Furthermore, pooling 15 d and 30 d RNA-seq samples to define a compound 'Treat' group captures broad, persistent transcriptional responses but sacrifices the resolution of distinct temporal dynamics.

Several lines of evidence remain to be established to test causality and refine the permissive-landscape interpretation. Although our ATAC-seq data capture the uninduced baseline, time-resolved accessibility measurements after injury would help determine whether the observed constitutive differences are further amplified during induction or remain static. Differences in tissue composition at the wound site may contribute to apparent accessibility and expression differences. Wound healing dynamics and local microenvironmental factors may also contribute. These factors should be evaluated with expanded sampling and controlled tissue collection. Locus-focused validation of prioritized promoters/enhancers (e.g., ATAC-qPCR/FAIRE-qPCR) and targeted metabolite quantification would strengthen links from regulatory DNA to pathway output. The inferred TF-target relationships remain predictive in nature. Direct tests such as ChIP-seq/CUT&Tag, transient reporter assays, and perturbation of prioritized TFs will be required to validate binding and regulatory function.

Chromatin accessibility is unlikely to be the only epigenetic layer involved. DNA methylation at promoter regions is one plausible component^[38]. Integrating methylome assays with functional TF tests will help clarify how epigenetic modifications and TF coordination jointly shape the high-yielding agarwood phenotype.

We compared BM and SXY under identical drilling induction and found that SXY showed higher resin accumulation together with a more open, promoter-enriched chromatin landscape. Integrating ATAC-seq and RNA-seq further identified a set of terpenoid-pathway genes, including AsTPS1, whose accessibility and transcriptional patterns were more closely aligned in SXY. Motif enrichment and network analysis also highlighted candidate regulators from the bHLH and bZIP families, including BHLH137 and HYH. The results point to a model in which chromatin accessibility-associated regulatory differences contribute to stronger terpenoid-pathway activation in SXY.