Leaf materials of the borneol-type C. burmannii were collected from a six-year-old plant (Figure 1A) in the Tianluhu Forestry Park, Guangzhou, China, which permitted sample collection. Leaves at four developmental stages, including 7 days (CBS1), 14 days (CBS2), 21 days (CBS3), and 28 days (CBS4) (Figure 1B), were collected. For each developmental stage, three replicates were prepared, resulting in a total of 12 samples. The samples were immediately stored in a -80°C freezer and used for leaf transverse sectioning, GC-MS analysis, ONT sequencing, and qRT-PCR analysis.

The C. burmannii leaf samples were soaked in a formaldehyde-acetic acid-alcohol mixture (FAA) for 48 h and subsequently embedded in paraffin wax. Then, 8 μm thick transverse sections were obtained from these paraffin blocks using a rotary microtome (Leica, RM2265, Germany). The freshly cut sections were dried on a slide drier (KEDEE, KD-P, US) at 40°C and stained with the safranin O-fast green (Getzy et al., 1982). The prepared sections were observed using a light microscope (Nikon TE2000-U, Nikon, Japan).

The leaf samples were ground into a powder with liquid nitrogen and dissolved in 4 ml of hexyl hydride, ethanol, and ethyl acetate mixture (2:1:1). The organic solution containing the tissue powder was treated with ultrasonic waves at 20 kHz for 20 min and incubated in a water bath at 56°C for 30 min. After centrifugation, the supernatant was collected, mixed with cyclohexanone (Sigma-Aldrich, Shanghai, China), and analyzed on a gas chromatograph-mass spectrometer (GC-MS) with the SH-Rxi-5Sil.MS quartz capillary column (30 m × 0.25 mm × 0.25 µm). To assess relative content of d-borneol with GC-MS, the initial column temperature was set at 70°C, then raised to 160°C at a rate of 2°C/min, maintained for 2 min, then to 220°C at a rate of 10°C/min, and maintained for 5 min; it took 51 min for the warming program to detect one sample. The ionization temperature was set at 200°C, the interface temperature at 250°C, and the scanning range (m/z) at 50–500. Besides, to measure absolute content of d-borneol, 200 mg of standard d-borneol (purity ≥99%, Sigma-Aldrich, Shanghai, China) were measured and transferred into a 10 ml volumetric flask, added with cyclohexanone to the constant volume, yielding 20 mg/ml of d-borneol standard solution. 0.005, 0.025, 0.1, 0.5, and 2.5 ml of the above solution was transferred into five 10 ml volumetric flasks, correspondingly. Cyclohexanone was added to each constant volume, yielding 10, 50, 200, 1,000, and 5,000 mg/L of d-borneol standard solution. 1 ml of each sample was transferred into the injection bottle for GC-MS testing. Under the same GC-MS conditions shown above, the absolute content of borneol in the leave samples were calculated according to the linear equation.

Total RNA was extracted from the leaf samples using a RNeasy Plant Mini Kit (Qiagen, Valencia, CA, United States, product No. 74903). After removal of the relic using RNase-free DNase (Qiagen), the quality of the RNA sample was assessed using 1% agarose gel electrophoresis, a NanoDrop spectrophotometer (ThermoFisher Scientific, Wilmington, DE, United States), and an Agilent 2,100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, United States). The RNA was then reverse-transcribed to cDNA according to the protocol by Oxford Nanopore Technologies (Jain et al., 2016). The sequencing libraries were created using a preparation kit (Ultra-Long, Grandomics, Wuhan, China) by fragmenting the cDNA and adding PCR adapters to both ends. After PCR amplification (14 cycles) with LongAmp Taq (NEB, New England Biolabs LTD., Beijing, China), the ONT adaptors were ligated onto the PCR products using T4 DNA ligase (NEB). The ligated PCR products were sequenced on a MinION Mk1B sequencer (Oxford Nanopore, Oxford, United Kingdom). All sequencing data have been deposited in the NCBI Sequence Read Archive (SRA) under the data have been released from the SRA database since November 23, 2022.

Raw sequencing data were quality filtered to remove the adapter sequences, reads with quality scores ≤7, and reads with length ≤500 bp using MinKNOW (version 2.2) (Oxford, United Kingdom). Reads similar to ribosomal RNA were also deleted by searching against the Silva rRNA database (https://www.arb-silva.de). We then identified the full-length reads and clustered them after mapping to the C. burmannii reference genome (Hou et al. unpublished) using minimap2 (version 2.16) with the setting (-t 6 -ax splice--split-prefix -uf--secondary = no) (Li 2018). Each cluster of full-length reads was collapsed to obtain a consensus sequence compared with the C. burmannii reference genome using pinfish (https://github.com/nanoporetech/pinfish, version 0.1.0) in the minimap2 software. The full-length reads with a minimum coverage of 85% and a minimum identity of 90% were saved using the cDNA_Cupcake package (version 5.80), and the remaining redundant reads were deleted. The full-length reads with different 5′ ends were not considered redundant.

Gene expression levels were quantified using the non-redundant full-length reads by calculating the CPM (counts per million) values as follows: CPM = reads mapped to transcripts/total reads aligned to one sample × 1,000,000. The differential expression analysis was performed using the DESeq2 R package (version 1.10.1) to identify the genes differentially expressed between two developmental stages (Love et al., 2014). Genes with fold change ≥2 and false discovery rate (FDR) < 0.01 were defined as the differentially expressed genes (DEGs). Adjusted p-values with the Benjamini–Hochberg approach were used to assess the DEGs accurately. Further, the gene ontology (GO, http://www.geneontology.org) annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG, http://www.genome.jp/kegg/) enrichment analysis were performed for the DEGs. Genes enrichment were performed using the Fisher’s exact test.

After GC-MS analysis, a trait matrix of based on relative content was prepared for the most 10 abundant monoterpenoids and sesquiterpenoids (Supplementary Table S1). Then, the correlation of co-expressed genes and these metabolites was analyzed in a weighted correlation network analysis (WGCNA) implemented in the R package WGCNA (version 1.42) (Langfelder and Horvath 2008). A hierarchal clustering tree based on the co-expressed genes was constructed using the Dynamic Tree Cut R package (Langfelder et al., 2008). Genes significantly correlated with the abundance of the ten terpenoids were grouped into different modules and further annotated based on the KEGG database.

The coding sequences (CDS) of the non-redundant full-length reads were identified with TransDecoder based on log-likelihood score and open reading frames (ORFs) length. After searching against the CPC (Coding Potential Calculator) (Kong et al., 2007), CNCI (Coding-Non-Coding Index) (Sun et al., 2013), Pfam, and CPAT (Coding Potential Assessment Tool) (Wang et al., 2013) datasets, the remaining non-protein-coding reads with two exons and at least 200 nt long were identified as lncRNAs. These lncRNAs were further divided into four groups: lincRNA, antisense lncRNA, intronic lncRNA, and sense lncRNA. We further predicted the target genes regulated by the identified lncRNAs using LncTar (version 1.0) (Li et al., 2015). The lncRNA expression was quantified as CPM values, and the differentially expressed lncRNAs were identified with the setting (-p 2 -d -0.1 -s F). The target genes regulated by lncRNAs through cis- or trans-action were predicted separately using lncTar (Kung et al., 2013; Yang et al., 2014).

We selected nine DEGs (the TPS genes) for transcript expression validation using qRT–PCR. Primers were designed using the Primer Premier 5 software (Lalitha 2000). The information related to qRT-PCR is shown in Supplementary Table S2. The PCR was carried out using the following program: an initial denaturation at 95°C for 30 s, followed by 40 cycles of 10 s at 95°C, 30 s at 60°C, and 5 s at 95°C, and a final extension at 60°C for 5 s and 90°C for 5 s; the melting curves were drawn for every 0.5°C. The Actin gene was used as the endogenous control, and the relative expression levels were determined using the ΔΔCt method (Livak and Schmittgen 2001). Three replicates were maintained per sample, and the gene expression levels were represented as mean ± standard deviation. Differences of expression levels between CBS1 and the other three stages of C. burmannii leaves were analyzed by one-way analysis of variance (ANOVA), followed by an independent sample t-test with the R software. p values <0.05 was thought to be statistically significant.

In order to predict functions of TPS genes, the sequences of the TPS genes identified in this study were aligned using MUSCLE (version v3.8.31) (Edgar 2004). The sequence alignment comprises the sequences of seven TPS genes reported earlier (CbTPS1-S7 in Ma et al., 2021) and those mentioned by Yang et al., 2020 (see Supplementary_data18, DOI: 10.7717/peerj.9311/supp-18). Besides, the sequences of the IDS genes were also aligned, as mentioned by Yang et al., 2021 (See Supplementary Table S3). We performed multiple alignments for the IDS and TPS genes of C. burmannii using MUSCLE (version 3.8.31) (Edgar 2004). Then, RAxML-HPC2 (version 8.2.12) (Stamatakis 2014) implemented in CIPRES Gateway (version 3.3) (www.phylo.org) was used to construct two maximum likelihood (ML) trees, one for IDS genes and the other for TPS genes. Before the phylogenetic analysis, the GTR + CAT model was identified as the best-fit model for amino acid using jModeltest2 (version 2.1.6) (Darriba et al., 2012). A rapid bootstrapping was performed to search for trees with the highest score, and statistical support for the ML trees for IDI and TPS genes was derived from 1,000 pseudo-replicates of simulated bootstraps.

The cross sections showed that C. burmannii leaves are mainly composed of epidermis, palisade tissues, spongy tissues, vascular bundles, and hypodermis (Figures 1C,D). Glands were found at the lateral sides of vascular bundles in both the middle and lateral veins. In addition, dramatic changes were observed in the anatomical structures of the main veins (Figure 1C) and lateral veins (Figure 1D) across four developmental stages of C. burmannii leaves. At the CBS1 and CBS2 stages, protoxylem and metaphloem were present in both the main and lateral veins. At CBS3 and CBS4, the xylem and phloem substituted protoxylem and metaphloem. Besides, protuberances, composed of parenchymatous cells, were found above the vascular bundles of the middle and lateral veins from CBS2 to CBS4 (Figures 1C,D). In addition, a few eleocytes, which produce essential oil, emerged in the peripheral areas of the vascular bundles in the middle and lateral veins at CBS1 and CBS2. Subsequently, numerous eleocytes were found at CBS3 and CBS4. Besides, some idioblasts emerged at CBS2 and more found at CBS4 in the middle vein (Figure 1C), while many were merely found in the lateral veins at CBS4 (Figure 1D).

We further analyzed the secondary metabolites in the leaves of C. burmannii using GC-MS. The relative content of borneol significantly increased (p-value < 0.01) from CBS1 to CBS2 (Figure 2A). The relative content of d-borneol increased from 0.2 at CBS1 to 53.34 at CBS4, while the absolute content increased from 100.64 mg/g to 4762.67 mg/g. Other terpenoids in the C. burmannii leaves also increased from CBS1 and CBS4 (Figure 2B). We found the top 10 terpenoids from CBS2 to CBS4 were eucalyptol, d-limonene, alpha-pinene, caryophyllene, bornyl acetate, beta-myrcene, camphene, germacrene B, beta-pinene, and alpha-phellandrene.

Sequencing on the ONT platform generated 2,341,602 to 3,941,015 raw reads per sample (Supplementary Table S3), with a mean length ranging from 988 to 1,228 bp and N50 length ranging from 1,092 to 1,329 bp. After removal of adapters and low-quality reads, we obtained 2,218,026 to 3,725,516 clean reads, comprising 1,935,584 to 3,316,638 full-length reads (Supplementary Table S4). Of this, 1,920,760 to 3,284,474 reads were mapped to the C. burmannii reference genome, with mapping ratios ranging from 98.72% to 99.35% (Supplementary Table S5). After collapsing each cluster, the data generated 23,587 to 37,176 consensus sequences, with N50 length ranging from 1,257 to 1,497 (Supplementary Table S6). Subsequently, after deleting the redundant reads, 11,615 to 20,320 non-redundant full-length transcripts were obtained, with N50 lengths ranging from 1,415 to 1700 bp and mean lengths ranging from 1,283 to 1,463 bp (Supplementary Table S6).

Further analysis identified 881 (CBS2 vs. CBS3) to 6,051 (CBS1 vs. CBS4) differentially expressed genes (Supplementary Table S8). Among the 6,051 DEGs between CBS1 and CBS4, 2,749 were upregulated, and 3,302 were downregulated (Figure 3A). Venn diagram showed that 1,504 DEGs were commonly shared among the three comparison groups (CBS1 and CBS4, 6,051 DEGs; CBS1 and CBS3, 5,505 DEGs; CBS2 and CBS4, 2,687 DEGs; Figure 3B). We then annotated these 1,504 DEGs based on GO and KEGG databases. GO annotation showed that 21 DEGs were involved in the terpene metabolic process (GO ID: 0,042,214), 24 in terpene synthase activity (GO:0010333), six in the monoterpene metabolic process (GO:0043692), and six in monoterpene biosynthetic process (GO:0043692) (Figure 3C). Meanwhile, KEGG annotation showed that 28 DEGs were involved in the terpenoid backbone biosynthesis (KEGG orthology number: ko00900), 14 in sesquiterpenoid and triterpenoid biosynthesis (ko00909), six in diterpenoid biosynthesis (ko00904), and four in monoterpenoid biosynthesis (ko00902).

We then performed a weighted correlation network analysis (WGCNA) for correlations between the contents of top 11 terpenoids and co-expressed genes involved in the biosynthesis of various terpenoids. A phylogeny of co-expressed genes indicates 12 modules (Figure 4A), of which four modules (represented by the colors plum, salmon, dark olive green, and pale turquoise) accounted for most co-expressed genes involved in the synthesis of the top 11 terpenoids (Figure 4B). The co-expressed genes belonging to the salmon module showed positive correlations with the top 11 terpenoids (p < 0.01), while those of the plum module showed negative correlations (p < 0.05; Figure 4B). Co-expressed genes in the dark olive-green module and pale turquoise module exhibited significantly positive correlations with a few exceptions (Figure 4C). KEGG annotation of the co-expressed genes in the four modules revealed that 18 genes (salmon-6, pale turquoise-1, and plum-11) were involved in terpenoid backbone biosynthesis, 3 (salmon-2 and pale turquoise-1) in monoterpenoid biosynthesis, and 23 (salmon-11, dark olive green-1, pale turquoise-6, and plum-5) in sesquiterpenoid and triterpenoid biosynthesis.

Further, 15,201 ORFs were identified from the non-redundant full-length reads, obtained by TransDecoder, of which the majority were shorter than 1,000 bp (Figure 5A). Data filtering based on the four databases (CNCI, CPC, Pfam, and CPAT) identified 922 lncRNA with at least two exons and longer than 200 nt (Figure 5B, Supplementary Table S9). These lncRNAs were further classified into four groups, including 708 lincRNAs (74.9%), 49 antisense lncRNAs (4.3%), 28 intronic lncRNAs (2.1%), and 137 sense lncRNAs (18.7%) (Figure 5C). In addition, a total of 22 differentially expressed lncRNAs were commonly shared among the abundant DEGs of the three comparison groups (CBS1 and CBS4, 141 lncRNAs; CBS1 and CBS3, 1,365 lncRNAs; CBS2 and CBS4, 35 lncRNAs) (Figure 5D). Further analysis using LncTar identified 3,337 cis-regulated target genes/full-length reads for the lncRNAs and 536 trans-regulated genes. These cis-regulated target genes included five involved in monoterpenoid synthesis regulated by 12 differentially expressed lncRNAs (Figure 5E).

We further identified 44 DEGs (22—upstream and 22—downstream) associated with the molecular pathways regulating the synthesis of the top 10 monoterpenoids and sesquiterpenoids (Figure 6). The sequences of these DEGs are shown in Supplementary Table S10. Of the 44 DEGs detected, in the MVA and MEP pathways, our analysis identified more DEGs of the MVA pathway (12 genes) than the MEP pathway (ten genes). Moreover, the expression profiles of DEGs dramatically differed between the two pathways. The DEGs associated with the MVA pathway were highly expressed at CBS1 except for acetoacetyl-coenzyme A thiolase enzyme (AACT, Cbur04G014790 and Cbur04G014840) while the DEGs in the MEP pathway showed high expression at CBS3 except for 1-deoxy-d-xylulose 5-phosphate synthase (DXS, Cbur02G033440, Cbur05G01050, and Cbur04G023450). Besides, two IDI genes, Cbur07G020680 and Cbur07G020710, probably involved in inter-conventing between DMAPP and IPP, were highly expressed at CBS4. In the downstream, two DEGs were annotated as FPPS genes and four as GPPS genes, of which Cbur10G020080 (a FPPS gene) and Cbur10G002890 (a GPPS gene) were highly expressed at CB1, while the remaining four genes were highly expressed at CBS3. We found nine genes annotated as sesqui-TPS, which were highly expressed at CBS3 and CBS4, except for Cbur09G007320. Meanwhile, 14 genes were annotated as mono-TPS genes, of which 10 genes (Cbur03G002680, Cbur0G021000, Cbur03G002160, Cbur03G002310, Cbur09G004970, Cbur0G006370, Cbur09G005050, Cbur09G005010, Cbur09G005210 and Cbur09G005220) showed high expression at CBS3, while Cbur03G002300 and Cbur09G007320 were highly expressed at CBS4 and CBS1, respectively. The high expression of the mono-terpene synthase genes at CBS3 are most likely to explain the increasing content of monoterpenoids and sesquiterpenoids from CBS2 to CBS4. We further validated the expression levels of nine TPS genes at four developmental stages using qRT-PCR are consistent with the gene expression generated by the ONT sequencing (Figure 7).

In order to characterize TPS and IDI genes identified in the present study, two phylogenetic trees were established using the maximum likelihood method. In the TPS tree, we found the 14 TPS genes were resolved into seven clades (Figure 8); seven were assigned to the TPS-a subfamily, five to the TPS-b subfamily, and two to the TPS-g subfamily. In the IDS tree, the six IDI genes were resolved into four clades (Figure 9); two each were nested within the two clades representing the IDS-a and IDS-d subfamilies. Meanwhile, Cbur10G002890 alone was assigned to the IDS-c clade, while Cbur12G014500 was assigned to the IDS-d clade.