Abstract

   Dracaena, a remarkably long-lived and slowly maturing species of plant,
   is world famous for its ability to produce dragon’s blood, a precious
   traditional medicine used by different cultures since ancient times.
   However, there is no detailed and high-quality genome available for
   this species at present; thus, the molecular mechanisms that underlie
   its important traits are largely unknown. These factors seriously limit
   the protection and regeneration of this rare and endangered plant
   resource. Here, we sequenced and assembled the genome of Dracaena
   cochinchinensis at the chromosome level. The D. cochinchinensis genome
   covers 1.21 Gb with a scaffold N50 of 50.06 Mb and encodes 31 619
   predicted protein-coding genes. Analysis showed that D. cochinchinensis
   has undergone two whole-genome duplications and two bursts of long
   terminal repeat insertions. The expansion of two gene classes,
   cis-zeatin O-glucosyltransferase and small auxin upregulated RNA, were
   found to account for its longevity and slow growth. Two transcription
   factors (bHLH and MYB) were found to be core regulators of the
   flavonoid biosynthesis pathway, and reactive oxygen species were
   identified as the specific signaling molecules responsible for the
   injury-induced formation of dragon’s blood. Our study provides
   high-quality genomic information relating to D. cochinchinensis and
   significant insight into the molecular mechanisms responsible for its
   longevity and formation of dragon’s blood. These findings will
   facilitate resource protection and sustainable utilization of Dracaena.

   Key words: genome, Dracaena cochinchinensis, dragon’s blood, longevity,
   flavonoid biosynthesis, ROS
     __________________________________________________________________

   This study reports the chromosome-level genome assembly for Dracaena
   cochinchinensis and demonstrates that the longevity and slow growth of
   this species are associated with the significant expansion of SAUR and
   cZOGT genes. bHLH and MYB transcription factors are the core regulators
   of flavonoid biosynthesis, and ROS act as specific signaling molecules
   during the process of injury-induced dragon’s blood formation.

Introduction

   Dragon’s blood, the red resin produced by tree species in the genus
   Dracaena (Asparagaceae), is a precious traditional Chinese medicine and
   a folk medicine used by many nationalities across the world. It has
   several therapeutic functions and can activate blood circulation,
   dissipate blood stasis, relieve inflammation and pain, and improve
   astringency and hemostasis ([57]Gupta et al., 2008; [58]Fan et al.,
   2014; [59]Sun et al., 2019). Apart from being a famous medicine,
   dragon’s blood is also used as a colorant for artwork in many cultures
   ([60]Zheng et al., 2012). It is called dragon’s blood because the resin
   or sap excreted by some members of the genus Dracaena after external
   injury or microbial invasion is deep red in color and looks like the
   blood of a legendary dragon ([61]Gupta et al., 2008). Therefore,
   Dracaena has become famous all over the world for its red resin and
   belongs to the cultural heritage of humanity.

   Dracaena species are monocotyledonous evergreen plants that are mainly
   distributed in tropical and subtropical regions of Asia and Africa.
   These are remarkably long-lived and slow-growing species and are
   renowned for their longevity ([62]Zheng et al., 2012; [63]Tomlinson and
   Huggett, 2012; [64]Maděra et al., 2018, [65]2020). The age of a famous
   dragon’s blood tree was previously estimated to be several thousands of
   years ([66]Humboldt, 1814; [67]Christ, 1886; [68]Schenck, 1907).
   However, subsequent analysis of several of the largest trees concluded
   that their age was more likely to be several hundred years and that the
   flowering period varies from 10 years to several decades ([69]Pütter,
   1925; [70]Symon, 1974; [71]Mägdefrau, 1975; [72]Zheng et al., 2012;
   [73]Maděra et al., 2018). According to previous research, the juvenile
   stages of D. cinnabari from seedling to first flowering can take
   200–250 years or even more ([74]Maděra et al., 2018), thus
   demonstrating that this species takes an extremely long time to grow
   and mature. Owing to difficulties related to investigating natural
   population aging, several methodologies have been developed for age
   estimation in such trees ([75]Adolt and Pavliš, 2004; [76]Adolt et al.,
   2012; [77]Jura-Morawiec, 2019; [78]Lengálová et al., 2020). These
   methods have all demonstrated that the dragon tree is a long-lived
   species.

   There is much confusion surrounding the origin and identity of dragon’s
   blood trees, and many plants can produce dragon’s blood, including
   Croton, Dracaena, Daemonorops, and Pterocarpus ([79]Gupta et al., 2008;
   [80]Ding et al., 2020). In the wild, Dracaena grows slowly, and only
   trees that are over 30 years of age have been shown to produce dragon’s
   blood ([81]Edwards et al., 1997; [82]Wang et al., 2013; [83]Ding
   et al., 2020). Due to continuous over-exploitation, Dracaena are now
   threatened with depletion. In 1998, many Dracaena species were listed
   on the International Union for Conservation of Nature Red List
   ([84]IUCN Red List, 2017.2). Of these, Dracaena cochinchinensis and
   D. cambodiana were rated as vulnerable in China ([85]Hubálková et al.,
   2015; [86]Kamel et al., 2018). Soon afterward, these two species were
   placed on the List of National Key Protected Wild Plants of China; this
   represents the second highest grade of national protection, and the
   felling of such endangered species is prohibited. Because of the
   endangered status of D. cochinchinensis and the fact that it is now a
   limited genetic resource, it is critical to consider adequate measures
   for its conservation management and to protect against its
   exploitation.

   In China, D. cochinchinensis is the official original species found to
   produce dragon’s blood ([87]Zheng et al., 2009; [88]Luo et al., 2011)
   and was first discovered by Xitao Cai in 1972 in Yunnan province. This
   species is commonly used as a substitute for the traditionally imported
   dragon’s blood, called Long-Xue-Jie (Chinese dragon’s blood). According
   to literature surveys, dragon’s blood has been recorded and used
   medicinally in China for at least 1500 years. The authentic source of
   ancient Chinese Dracaena was thought to have been imported from Africa
   ([89]Cai and Xu, 1979). With increased levels of shipping, the genuine
   form of Dracaena changed from the resin of Dracaena in the agave family
   to the red resin of the rattan fruit in Southeast Asia until the
   discovery of D. cochinchinensis, which is exclusively found in China
   ([90]Cai and Xu, 1979; [91]Zheng et al., 2006; [92]Wang et al., 2011).

   Flavonoids are the main components of dragon’s blood. Flavonoid
   oligomers, composed of one dihydrochalcone unit condensed with one or
   more chalcone, flavane, or homoisoflavane units, constituted the major
   identified components of dragon’s blood from the genus Dracaena
   ([93]Pang et al., 2016). Oligomeric flavonoids account for over 50% of
   the resin produced by Dracaena spp. by weight ([94]Hao et al., 2015).
   Thus far, hundreds of flavonoids have been isolated from the resin wood
   of D. cochinchinensis (for reviews, see [95]Gupta et al., 2008; [96]Fan
   et al., 2014, [97]Sun et al., 2019). In addition, terpenoids and
   steroidal saponins are also effective components of dragon’s blood
   ([98]Fan et al., 2014; [99]Teng et al., 2015; [100]Maděra et al., 2020;
   [101]Vanickova et al., 2020); these components have been reported to be
   cytotoxic and could be used as anticarcinogens ([102]Darias et al.,
   1989; [103]Mimaki et al., 1999; [104]Gonźalez et al., 2003;
   [105]Herńandez et al., 2004; [106]Ding et al., 2020; [107]Maděra
   et al., 2020).

   Although many previous studies have described the chemical constituents
   and pharmacological effects of dragon’s blood (for reviews, see
   [108]Gupta et al., 2008; [109]Fan et al., 2014; [110]Sun et al., 2019;
   [111]Maděra et al., 2020), little is known about the molecular
   mechanisms responsible for its formation. In this study, we report the
   draft genome sequence of D. cochinchinensis at the chromosome level. We
   also conducted transcriptomic and metabolomic analysis. Our findings
   represent an important resource that can enhance our understanding of
   the fundamental biology and molecular basis of D. cochinchinensis, thus
   providing a foundation from which to generate new conservation
   management strategies.

Results

Genome sequencing, assembly, and annotation

   A genome survey using k-mer analysis (k = 17) revealed that the size of
   the D. cochinchinensis genome was approximately 1.257 Gb with a
   heterozygosity of 0.52% and a repeat sequence ratio of 67.8%
   ([112]Supplemental Figure 1; [113]Supplemental Table 1). Based on the
   estimated genome size, using a combination of SMART sequencing
   technology from PacBio, along with 10x Genomics and short-read
   sequencing from Illumina platforms, we obtained a total of 407.06 Gb of
   clean data, thus representing 323.75× coverage ([114]Supplemental
   Table 2). The initial assembly was 1.21 Gb in size, containing 1898
   scaffolds and 3042 contigs, with a scaffold N50 size of 2.07 Mb and a
   contig N50 size of 1.1 Mb ([115]Supplemental Table 3). Assessment of
   the completeness of the genome assembly with CEGMA indicated 96.77%
   coverage of the conserved core eukaryotic genes, and Benchmarking
   Universal Single-Copy Ortholog analyses ([116]Simao et al., 2015)
   showed that the genome was 93.6% complete ([117]Supplemental Table 4).
   These results suggested that the genome assembly was relatively
   complete and of high quality. Next, approximately 137.26 Gb of
   high-quality Hi-C data were used to further improve the genome
   assembly, and 90.26% of the genome was anchored onto 20 chromosomes
   (2n = 40) ([118]Figure 1B; [119]Supplemental Figures 2 and [120]3;
   [121]Supplemental Table 5). The final assembly consisted of 1835
   scaffolds, spanning 1.21 Gb, with a scaffold N50 size of 50.05 Mb
   ([122]Supplemental Table 6).

Figure 1.

   [123]Figure 1
   [124]Open in a new tab

   Morphological and genomic features of D. cochinchinensis.

   (A) Morphology of D. cochinchinensis and dragon’s blood.

   (B) Overview of the D. cochinchinensis draft genome assembly. (a) The
   20 chromosomes of D. cochinchinensis, with a resolution of 5 Mb. (b)
   Gene density, with a sliding window of 500 kb. (c) Percentage of
   repeats, with a sliding window of 500 kb. (d) GC content, with a
   sliding window of 100 kb. (e) Each linking line in the center of the
   circle connects a pair of homologous genes.

   Based on de novo and homology-based predictions and transcriptome data,
   a total of 31 619 protein-coding genes were predicted, and 31 402 genes
   (99.3%) were assigned specific functions ([125]Supplemental Table 7).
   In addition, by comparing with the known noncoding RNA libraries, we
   obtained noncoding RNA information relating to the D. cochinchinensis
   genome ([126]Supplemental Table 8), including 1570 transfer RNAs, 1278
   ribosomal RNAs, 1576 microRNAs, and 722 small nuclear RNAs.

Genomic evolution of D. cochinchinensis

   D. cochinchinensis belongs to the Asparagales clade of monocotyledons
   and is close to Asparagus officinalis, Dendrobium officinale, and
   Apostasia shenzhenica, which have completed genome sequences. We
   clustered the annotated genes into gene families among these species
   along with additional monocot and dicot species, including
   A. officinalis, D. officinale, A. shenzhenica, Elaeis guineensis, Oryza
   sativa, Ananas comosus, Musa acuminata, Arabidopsis thaliana, Populus
   trichocarpa, Aquilaria sinensis, Hevea brasiliensis, Citrus grandis,
   Carica papaya, Eucalyptus grandis, Vitis vinifera, and Amborella
   trichopoda. Family clustering yielded 31 324 gene families (groups),
   including 252 single-copy families. Evolutionary trees were constructed
   based on coding sequences using RAxML software. D. cochinchinensis and
   A. officinalis were gathered in one branch, as were D. officinale and
   A. shenzhenica, thus indicating that Asparaginaceae and Orchidaceae
   underwent independent divergence long ago (138.4 Mya); the time of
   divergence between D. cochinchinensis and A. officinalis was 87.7 Mya
   ([127]Figure 2A).

Figure 2.

   [128]Figure 2
   [129]Open in a new tab

   Comparative genomic analysis.

   (A) Phylogenetic tree of 17 plant species. Divergence time (Mya)
   estimates are indicated by the blue numbers beside the branch nodes.
   The numbers of gene-family contractions and expansions are indicated by
   green and red numbers, respectively.

   (B) Distribution of 4DTv for paralogous genes among D. cochinchinensis,
   A. officinalis, D. officinale, and A. shenzhenica.

   (C) Estimation of LTR activity showed two insertion burst events in the
   D. cochinchinensis genome (Mya).

   (D) In the recent burst of activity, most of the elements responsible
   were autonomous (3359 copies of Ty1-copia and 16 897 copies of
   Ty3-gypsy) and non-autonomous LTRs (5861 copies) with peaks of activity
   <6 Mya.

   To investigate the genome expansion in D. cochinchinensis, we analyzed
   whole-genome duplication (WGD) events. Four-fold degenerate site (4DTv)
   values were estimated based on paralogous gene pairs in collinear
   regions detected in D. cochinchinensis and three other representative
   plant species (A. officinalis, D. officinale, and A. shenzhenica). The
   distribution of 4DTv distances in D. cochinchinensis showed two peaks
   at approximately 0.12 and 0.45 ([130]Figure 2B). The first peak shared
   by D. cochinchinensis and A. officinalis may correspond to the
   Asparagales-β event previously identified in the A. officinalis genome
   ([131]Harkess et al., 2017; [132]Li et al., 2020). The second peak at
   approximately 0.12 indicated that D. cochinchinensis underwent another
   WGD event independently after diverging from A. officinalis
   ([133]Figure 2B and [134]Supplemental Figure 4). To further confirm
   that D. cochinchinensis had undergone two WGD events, we detected
   syntenic blocks across the D. cochinchinensis genome and the
   A. trichopoda genome and performed synteny analysis. We found that
   D. cochinchinensis had experienced the Asparagales-β event and another
   WGD event, as suggested by the 4:1 syntenic relationship between
   D. cochinchinensis and A. trichopoda ([135]Supplemental Figure 5). Our
   analysis thus provides convincing evidence for two rounds of
   whole-genome diploidization events in the D. cochinchinensis genome. In
   total, 29 280 duplicated genes were identified and classified into
   five categories: 6864 WGDs (23.1%), 1506 tandem duplicates (5.1%), 1196
   proximal duplicates (4.1%), 7301 transposed duplicates (24.9%), and
   19 613 dispersed duplicates (66.7%). Next, we investigated the
   duplications of genes involved in flavonoid biosynthesis (i.e., PAL,
   4CL, C4H, CHI, CHS, DFR, F3H, FLS, FAR, PPO, and OMT) and found that
   some genes were duplicated in a different manner. Some genes were
   retained after WGD, thus suggesting that the recent WGD event played an
   important role in the evolution of flavonoid biosynthesis in
   D. cochinchinensis ([136]Supplemental Table 9).

LTR insertion and genome size expansion

   Further analysis showed that 60.15% (729.11 Mb) of the
   D. cochinchinensis genome consists of repetitive sequences
   ([137]Supplemental Table 10); this is similar to D. officinale (63.33%)
   ([138]Yan et al., 2015) and E. guineensis (57%) ([139]Singh et al.,
   2013) but higher than A. shenzhenica (42.05%) ([140]Zhang et al.,
   2017a) and A. comosus (43.7%) ([141]D’hont et al., 2012). As in many
   other sequenced plant genomes ([142]Kim et al., 2014; [143]Teh et al.,
   2017; [144]Zhang et al., 2017b), long terminal repeat retrotransposons
   (LTR-RTs) dominated, accounting for 51.06% of the assembly; 41% of the
   genome sequence comprised Gypsy elements (501 175 881 bp), and 5.9%
   comprised Copia elements (72 176 608 bp). By analyzing the insertion
   times of intact LTRs, we found that two insertion bursts had occurred
   in the D. cochinchinensis genome at approximately 0.3–0.4 and 2 million
   years, respectively ([145]Figure 2C); these events were dominated by
   both autonomous and nonautonomous LTRs. Among these, the Gypsy family
   had the largest number of insertions, and there were two obvious
   insertion peaks ([146]Figure 2D), thus demonstrating that there is a
   significant relationship between Gypsy expansion and genome expansion.
   Phylogenetic trees of domains in reverse transcriptase genes from
   complete RTs were constructed for the Copia and Gypsy superfamilies in
   D. cochinchinensis, A. officinalis, and A. shenzhenica. We found that
   the majority of Gypsy elements were species specific in
   D. cochinchinensis, whereas the Copia superfamily exhibited a different
   pattern that consisted of elements from all three species
   ([147]Supplemental Figure 6).

Comparative genomics demonstrated the adaptation of D. cochinchinensis

   Copy numbers of homologous gene families vary greatly among different
   species, caused by differences in rates of gene gain and loss. Analysis
   using Computational Analysis of Gene Family Evolution (CAFÉ) ([148]Bie
   et al., 2006) identified 81 expanded gene families in
   D. cochinchinensis compared with the common ancestor of
   D. cochinchinensis and A. officinalis ([149]Figure 2A). A total of 730
   and 1088 genes in the expanded families were annotated with Gene
   Ontology terms and Kyoto Encyclopedia of Genes and Genomes (KEGG)
   pathways, respectively. Gene Ontology analysis revealed that these
   expanded orthogroups were related to several metabolic processes
   ([150]Supplemental Table 11; [151]Supplemental Figures 7 and
   [152]8) such as “ATP metabolic process,” and “defense response” in the
   biological process category and “binding,” “catalytic
   activity,” “oxidoreductase activity,” and “peroxidase activity” in the
   molecular function category. The enriched KEGG categories were
   “plant-pathogen interaction,” “phenylpropanoid biosynthesis,” “zeatin
   biosynthesis,” “sesquiterpene and triterpenoid biosynthesis,” “flavone
   and flavonol biosynthesis,” and “photosynthesis” ([153]Supplemental
   Table 12; [154]Supplemental Figure 9). We hypothesize that these
   expanded genes are closely related in some way to the functional
   requirements of environmental adaptability in D. cochinchinensis and
   the formation of dragon’s blood.

   We found that 28 genes were under positive selection and showed
   significant enrichment in KEGG terms related to base excision repair,
   nucleotide excision repair, mismatch repair, and DNA replication
   ([155]Supplemental Table 13). These include DNA ligase 1 and DNA
   excision repair protein, which are involved in DNA-repair pathways that
   play crucial roles in maintaining genomic integrity, particularly in
   response to exposure to environmental genotoxicants, including
   ultraviolet light and reactive oxygen species that form in response to
   the harsh rays of the sun ([156]Schuch et al., 2017; [157]Nisa et al.,
   2019). The response of cells to DNA damage restores the DNA structure
   to enable its original function or simply enables the cell to survive
   the damage.

   Next, we identified unique and shared gene families among the four
   Asparagales species.D. cochinchinensis shared 8056 families with
   A. officinalis, D. officinale, and A. shenzhenica, whereas 1367 gene
   families, containing 2954 genes, appeared to be unique to D.
   cochinchinensis ([158]Figure 3A; [159]Supplemental Tables 14 and
   [160]15). Functional annotation revealed that the expansion of two
   classes of genes, small auxin upregulated RNA (SAUR) genes and
   cis-zeatin O-glucosyltransferase (cZOGT) genes, which play a role in
   zeatin biosynthesis, may be closely related to the slow growth and
   longevity of D. cochinchinensis ([161]Figure 3B and 3C;
   [162]Supplemental Tables 16 and [163]17). Interestingly, analysis of
   the D. cochinchinensis genome identified evolutionary gene family
   expansions of the SAUR genes (77 in D. cochinchinensis, 27 in
   A. officinalis, 20 in A. shenzhenica, 33 in D. officinale, 56 in
   O. sativa, and 42 in E. guineensis). SAUR is an auxin-responsive
   protein that acts as the plant’s toolbox for the adaptation of growth
   and development ([164]Ren and Gray, 2015; [165]Stortenbeker and Beme,
   2019). Previous studies in rice confirmed that SAUR inhibits growth by
   negatively regulating auxin synthesis and transport. In addition, SAUR
   has been shown to significantly increase anthocyanin content and ABA
   levels ([166]Kant et al., 2009). Phylogenetic analysis of SAUR proteins
   from D. cochinchinensis and other species identified 11 subgroups
   ([167]Supplemental Figure 10). To investigate the potential mechanisms
   by which SAURs may have expanded, we compared SAUR-containing genomic
   regions in D. cochinchinensis with those of A. officinalis and
   E. guineensis. The two species showed strong synteny to
   D. cochinchinensis ([168]Figure 3D). In addition to the SAURs, cZOGT
   had also undergone unique expansion in D. cochinchinensis (34 in
   D. cochinchinensis, 0 in A. officinalis, 0 A. in shenzhenica, 0 in D.
   officinale, 5 in O. sativa, 2 in E. guineensis, and 0 in A. thaliana)
   and is known to play roles in growth retardation and delay of
   senescence ([169]Rodo et al., 2008; [170]Kudo et al., 2012). Therefore,
   the expansion of SAURs and cZOGTs is predicted to account for the slow
   growth and longevity of D. cochinchinensis, at least in part, and even
   its adaptation to the arid environment of karst landforms. To
   investigate the functions of SAUR and cZOGT genes in
   D. cochinchinensis, we examined their expression levels at different
   time points after injury ([171]Supplemental Figures 11 and [172]12). We
   found that the expression levels of the two genes differed, implying
   that these genes may have other novel functions in D. cochinchinensis.

Figure 3.

   [173]Figure 3
   [174]Open in a new tab

   Expansion of the SAUR and ZOGT gene families.

   (A) Venn diagram showing shared and unique orthologs among the four
   Asparagales species D. cochinchinensis, A. officinalis, D. officinale,
   and A. shenzhenica.

   (B) Schematic representation of the D. cochinchinensis chromosomes with
   the positions of the SAUR and ZOGT genes. The genes marked in blue and
   red were generated by WGD and tandem duplication events, respectively.

   (C) Collinearity of SAUR and ZOGT genes in D. cochinchinensis.

   (D) Syntenic analysis of SAUR and ZOGT genes in D. cochinchinensis and
   two other representative plant species (A. officinalis and
   E. guineensis). Paralogous gene pairs from A. officinalis and
   E. guineensis were identified in syntenic blocks using JCVI software
   (v.0.8.4) with a parameter setting of “–cscore = 0.99”. Protein
   sequences were extracted and viewed in TBtools (v.0.98745).

Metabolic profiling and differential metabolite pathways in the formation of
dragon’s blood

   Dragon’s blood does not appear directly after wounding but instead
   forms over time. By observing tissue sections, we noted that the color
   of the resin gradually deepened during the formation of dragon’s blood;
   the resin had filled the vascular bundle and the surrounding
   parenchymal cells after 90 days ([175]Figure 4A). Eventually, the resin
   appeared in the form of tear-like drops or chips that coated the site
   of injury ([176]Figure 1A), protecting the plant from the spread of
   pathogens. Thus, the formation of dragon’s blood is a mechanism by
   which the tree defends itself ([177]Wang, 2007; [178]Jura et al., 2016;
   [179]Maděra et al., 2020).

Figure 4.

   [180]Figure 4
   [181]Open in a new tab

   Histological and metabolic profiling during the formation of dragon’s
   blood.

   (A) Histological images during the formation of dragon’s blood in
   D. cochinchinensis (the yellow and red substances in the image are
   dragon’s blood resin). VB, vascular bundle; PC, parenchymal cell. Scale
   bar: 5 μm.

   (B) Typical total ion chromatograms of D. cochinchinensis (a: positive
   total ion chromatogram [TIC] before dragon’s blood formation; b:
   positive TIC after dragon’s blood formation; c: negative TIC before
   dragon’s blood formation; d: negative TIC after dragon’s blood
   formation).

   (C) Heatmap of identified metabolites from D. cochinchinensis before
   and after the formation of dragon’s blood (BFD, 1–6 represent six
   replicates of the analytical materials before the formation of dragon’s
   blood; AFD, 1–6 represent six replicates of the analytical materials
   after the formation of dragon’s blood).

   To elucidate the changes in metabolite composition during the formation
   of dragon’s blood, we focused on identifying significant changes in
   metabolites by filtering selected peaks. Ultimately, a total of 147
   differential metabolites were identified or tentatively identified
   ([182]Supplemental Tables 18 and [183]19; [184]Figure 4B and 4C). As
   shown in [185]Figure 4C, the levels of most compounds changed
   significantly after injury. The levels of 85 compounds increased,
   including flavonoids, phenols, polyphenols, phenylpropanoids, and
   phenolic acids, whereas the levels of 62 compounds decreased, including
   organic acids, amines, and steroids; normalized peak areas showed the
   same results ([186]Supplemental Figure 13).

   Next, we used transcriptomic data to identify the biological processes
   associated with the formation of dragon’s blood. The upregulated
   differentially expressed genes were mainly enriched in the biological
   processes “metabolic pathway” and “biosynthesis of secondary
   metabolites,” thus showing that different time points correspond to
   different metabolic processes ([187]Figure 5). For example, the first
   metabolic pathways to change (in 2 h) were those related to
   plant-pathogen interaction and plant hormone signal transduction, as
   well as starch and sucrose metabolism. Then, the biosynthesis of amino
   acids and citrate cycle became significantly enriched; this was
   followed by flavonoid biosynthesis at 3 days; this lasted for 17 months
   and represented the most active metabolic pathway. Terpenoid backbone
   biosynthesis became evident at 15 days, sesquiterpenoid and
   triterpenoid biosynthesis at 30 days, and phenylpropanoid biosynthesis
   at 6 months. The differentially expressed metabolic pathways combined
   with the metabolic profiles are shown in [188]Figure 3C. Our results
   indicate that the synthesis of dragon’s blood involves a series of very
   complex signal transduction, gene expression regulation, and metabolic
   processes.

Figure 5.

   [189]Figure 5
   [190]Open in a new tab

   Differentially expressed metabolic pathways during the formation of
   dragon’s blood.

   Transcriptome analysis identified differentially expressed pathways
   during the formation of dragon’s blood. The upregulated pathways are
   shown above the oblique line, and the downregulated pathways are below
   the line. The size of the circles represents the number of genes, and
   the coordinate value represents the mean times of upregulation or
   downregulation.

Genes associated with the biosynthesis of flavonoids in D. cochinchinensis

   Numerous studies have shown that flavonoids are the main components of
   dragon’s blood ([191]Gupta et al., 2008; [192]Fan et al., 2014) and are
   associated with the phenylpropanoid biosynthesis pathway
   ([193]Figure 6). From our genome assembly, we identified 11 gene
   families related to the formation of dragon’s blood ([194]Supplemental
   Table 20), including some common enzymes in the flavonoid synthesis
   pathway: PAL, C4H, 4CL, CHS, CHI, F3H, FLS, and DFR. We also identified
   several downstream biosynthetic genes that produce specific components
   of dragon’s blood. For example, O-methyltransferase forms loureirin A,
   and polyphenol oxidase forms loureirin B, the index component of
   standard dragon’s blood in China. Gene expression analysis showed that
   four of the five CHS genes and two CHI genes were significantly
   upregulated after injury; this effect lasted for 17 months. Two of the
   three F3H genes were not expressed in healthy tissues; the highest
   expression levels occurred 3 and 5 days after injury. Some of the FLS
   and PPO members were also expressed at very low levels in healthy stems
   but were upregulated after injury ([195]Figure 6). We hypothesized that
   these genes may be specifically involved in the formation of dragon’s
   blood and that their expression is significantly induced by injury.

Figure 6.

   [196]Figure 6
   [197]Open in a new tab

   Flavonoidsbiosynthesis in D. cochinchinensis.

   The flavonoids biosynthesis pathway and the expression profiles (reads
   per kilobase per million reads mapped) of the genes involved in
   dragon's blood biosynthesis. PAL, phenylalanine ammonia-lyase; C4H,
   cinnamate 4-hydroxylase; 4CL, 4-coumarate CoA ligase; CHS, chalcone
   synthase; CHI, chalcone isomerase; FNS, flavone synthase; F3H,
   flavanone 3-hydroxylase; DFR, dihydroflavonol 4-reductase; FLS,
   flavonol synthase; LAR, leucoanthocyanidin reductase; PPO, polyphenol
   oxidase; OMT: O-methyltransferase.

   Furthermore, based on the co-expression model, two types of
   transcription factors (bHLH and MYB) expressed in wounded stems were
   predicted to be core regulators of flavonoid synthesis; these
   transcription factors are likely to target PAL, CHS, and FLS and
   participate in the regulation of the flavonoid biosynthesis pathway
   ([198]Supplemental Figure 14; [199]Supplemental Table 21). Analyzing
   the 2000 bp upstream of the gene coding region, we found many
   cis-acting elements that were specifically recognized by bHLHs and MYBs
   ([200]Supplemental Table 22), further demonstrating that these
   enzyme-encoding genes are likely to be regulated by the two
   transcription factor types. Moreover, an expression heatmap of bHLH and
   MYB genes ([201]Supplemental Figures 15 and [202]16) showed that they
   were significantly induced just 2 h after injury. We hypothesize that
   bHLHs and MYBs not only regulate the key enzyme-encoding genes PAL,
   CHS, and FLS but also regulate the entire flavonoid synthesis pathway.

Specific regulation of ROS signaling during the formation of dragon’s blood

   Analysis of genes involved in different signaling pathways showed that
   genes related to reactive oxygen species (ROS) production and
   scavenging were substantially upregulated ([203]Supplemental
   Figure 17). The oxidative burst is a common early event in plant
   defense responses, during which ROS rapidly accumulate following
   injury. The generation of ROS depends on NADPH oxidase and apoplastic
   peroxidases ([204]Bolwell and Wojtaszek, 1997). Peroxisomes are the
   major sites of intracellular H[2]O[2] accumulation as a consequence of
   their oxidative metabolism ([205]del Río et al., 2006). From the
   transcriptome data, we found that several enzymes, such as respiratory
   burst oxidase homolog, amine oxidases, NAD(P)H oxidases, and
   peroxidase, contributed to ROS generation and were significantly
   upregulated 2 h after injury ([206]Supplemental Figure 17). Thereafter,
   ROS scavenging enzymes, such as superoxide dismutase, ascorbate
   peroxidase, glutathione peroxidase, and catalase, were also upregulated
   ([207]Supplemental Figure 17). Obviously, the degree of oxidative
   stress is determined by the level of ROS, and the balance between ROS
   and antioxidants is essential if a plant is to maintain a balanced
   redox state.

   H[2]O[2] is the most stable type of ROS within plants. To preliminarily
   validate the H[2]O[2] response during the formation of dragon’s blood,
   we detected the changes in H[2]O[2] content within stems of
   D. cochinchinensis over time after mechanical damage. The concentration
   of H[2]O[2] increased just 1 h after mechanical damage and continued to
   increase to maximal levels at 6 h. After 8 h, the concentration of
   H[2]O[2] began to decrease slightly until the minimum value was reached
   at 5 days ([208]Supplemental Figure 18). By contrast, other signals,
   such as ethylene, nitric oxide, and brassinosteroids, showed no
   specific changes in response to injury ([209]Supplemental
   Figures 19–21). We infer that ROS act as the first sensor and
   represent the core signal regulator during the formation of dragon’s
   blood ([210]Figure 7). This is consistent with the conclusion that ROS
   generation represents a point at which various signaling pathways
   converge ([211]Camejo et al., 2016; [212]Sewelam et al., 2016).

Figure 7.

   [213]Figure 7
   [214]Open in a new tab

   A proposed model for ROS as a specific regulatory signal in the
   adaptation of D. cochinchinensis.

   An external wound leads to a burst of ROS. Excessive levels of ROS
   cause DNA damage. Most of the positively selected genes in the genome
   were enriched in DNA-repair pathways such as base excision repair,
   nucleotide excision repair, and mismatch repair; these processes play
   crucial roles in maintaining genomic integrity. ROS, as signaling
   molecules, also activate transcription factors, thus activating the
   flavonoid synthesis pathway to synthesize a large amount of flavonoids.
   Flavonoids are the main components of dragon’s blood and act as a
   defensive substance to resist further damage. In turn, flavonoids,
   together with enzymes in the ROS scavenging system, scavenge the excess
   ROS that have accumulated in vivo to maintain homeostasis.

Discussion

   Dragon’s blood is a precious traditional medicine with multiple
   therapeutic applications ([215]Masaoud et al., 1995; [216]Gupta et al.,
   2008; [217]Yang and Yao, 2011; [218]Xin et al., 2013; [219]Fan et al.,
   2014; [220]Ding et al., 2020; [221]Maděra et al., 2020).
   D. cochinchinensis is one of the most important genuine resources of
   Chinese dragon’s blood. Here, we constructed a high-quality reference
   genome for this valuable species at the chromosome level by combining
   PacBio sequencing and Hi-C data. The genome sequence of Dracaena
   reported herein will facilitate our understanding of the fundamental
   biology of this species and provide a foundation from which to generate
   strategies for protection and regeneration.

   Dracaena is world famous for its longevity ([222]Gupta et al., 2008;
   [223]Tomlinson and Huggett, 2012; [224]Zheng et al., 2012; [225]Maděra
   et al., 2018). However, the molecular mechanisms that underlie its
   unique biological traits remain unclear. Our analysis revealed that
   D. cochinchinensis has experienced two genome duplication events, one
   shared with A. officinalis and another that occurred independently
   after their divergence ([226]Figure 2B). In addition, two LTR insertion
   bursts occurred recently ([227]Figure 2C). Comparative genomic analysis
   revealed that expansions of SAUR and cZOGT genes are likely to be the
   main factors responsible for the longevity and slow growth of
   D. cochinchinensis ([228]Figure 3). In addition, genes in the DNA
   repair pathway were found to be positively selected, thus contributing
   to genome maintenance and environmental adaptability ([229]Supplemental
   Table 13).

   Dragon’s blood, the red resin of Dracaena, is a defensive substance
   produced by Dracaena when it is damaged or infected by pathogenic
   microorganisms ([230]Wang, 2007; [231]Jura et al., 2016; [232]Ding
   et al., 2020; [233]Maděra et al., 2020). Flavonoids are known to
   exhibit pharmacological effects and represent the main component of
   dragon’s blood; furthermore, flavonoids are highly beneficial to
   Dracaena during growth and development. For example, flavonoids are
   associated with plant defense against pathogens and microbes and the
   absorption of free radicals and ultraviolet light ([234]Neill and
   Gould, 2003; [235]Crozier et al., 2006; [236]Agati et al., 2009,
   [237]2011; [238]Agati and Tattini, 2010; [239]Landi et al., 2015). In
   the present study, we aimed to identify differentially expressed
   pathways and differential metabolites during the formation of dragon’s
   blood. We also attempted to identify genes responsible for regulating
   the biosynthesis of flavonoids in D. cochinchinensis. We found that
   several secondary metabolite-related pathways, including those related
   to flavonoids, the terpenoid backbone, and phenylpropanoid
   biosynthesis, were upregulated 0–6 months after injury ([240]Figure 5).
   Correspondingly, metabolomic analysis showed that the levels of
   flavonoid compounds increased significantly, along with those of
   phenylpropanoids, phenols, and phenolic acids. By contrast, the levels
   of coumarins, organic acids, and amine compounds were significantly
   reduced ([241]Figure 4C; [242]Supplemental Figure 13; [243]Supplemental
   Table 18). Upregulated genes related to the biosynthesis of secondary
   metabolites included PAL, CHS, and methyltransferase; these genes are
   all involved in the flavonoid biosynthesis pathway. The dynamic
   expression profiles of genes in the flavonoid biosynthesis pathway
   during dragon’s blood formation revealed that most of them were
   induced, including CHS, CHI, F3H, FLS, and PPO; these were only
   expressed at very low levels in healthy stems but were rapidly and
   significantly induced after injury. These changes are considered to be
   specifically involved in the formation of dragon’s blood
   ([244]Figure 6). This pattern of gene expression is essentially
   consistent with our results relating to differential metabolites and
   metabolic pathways ([245]Figures 4C and [246]5). Furthermore, bHLH and
   MYB transcription factors were predicted to be core regulators of the
   flavonoid biosynthesis pathway ([247]Supplemental Figures 14–16); these
   results are consistent with previous studies showing that MYB, bHLH,
   and WD40 transcription factors participate in flavonoid metabolism by
   regulating many structural genes ([248]Hichri et al., 2011;
   [249]Schaart et al., 2013). The corresponding relationship between
   these transcription factors and their targets was also predicted, thus
   providing important clues for studying the function of these genes and
   their regulatory mechanisms. Further studies are now required to
   comprehensively elucidate the process by which the red resin is formed
   in D. cochinchinensis and thereby enable controlled and reasonable
   extraction of dragon’s blood for commercial purposes.

   To reveal the regulatory network of dragon’s blood formation induced by
   injury, we systematically analyzed a range of signaling pathways and
   considered ROS to be the most crucial signal stimuli. Genes associated
   with ROS generation and scavenging were all upregulated after injury
   ([250]Supplemental Figure 17). Moreover, the concentration of H[2]O[2]
   showed regular dynamic changes at different time points after injury
   ([251]Supplemental Figure 18), and genes related to the main site of
   ROS production were also expanded ([252]Supplemental Tables 12 and
   [253]23). D. cochinchinensis grows in limestone areas that are exposed
   to strong sunlight and high UV-B; these stimuli can cause the
   degradation of CP43 and CP47 and reduce photosynthetic capacity, thus
   leading to ROS production, DNA damage, and protein denaturation
   ([254]Takshak and Agrawal, 2019). KEGG pathway analysis indicated that
   many of the positively selected genes were related to these deleterious
   processes; they participated in pathways such as base excision repair,
   nucleotide excision repair, mismatch repair, and DNA replication
   ([255]Supplemental Table 13), implying that D. cochinchinensis has
   adapted to the environment during evolution to maintain genomic
   integrity. In addition, gene families involved in phenylpropanoid
   biosynthesis and the photosynthesis pathway were all expanded
   ([256]Supplemental Tables 12 and [257]23). An increasing body of
   evidence supports the fact that flavonoids have multiple
   photoprotective effects ([258]Winkel-Shirley, 2002; [259]Agati and
   Tattini, 2010) and can function as ROS scavengers ([260]Yamasaki
   et al., 1997; [261]Agati et al., 2009) as well as moderators of auxin
   transport ([262]Peer and Murphy, 2007). We believe there are some
   potential associations in D. cochinchinensis relating to ROS signaling,
   flavonoid biosynthesis, and environmental adaptability ([263]Figure 7).
   This system is the result of adaptation to the environment during the
   process of evolution. We hypothesize that the longevity of
   D. cochinchinensis depends on its strong self-healing ability and
   environmental adaptability.

   In summary, this study reports a genome assembly for D. cochinchinensis
   at the chromosome level, the first genome for Dracaena Vand. ex L. It
   provides a valuable genetic resource and creates significant scope for
   studying Dracaena. We revealed the genetic basis of the longevity and
   slow growth of Dracaena at the genomic level. We also provided a
   preliminary analysis of the molecular mechanisms associated with
   injury-induced flavonoid biosynthesis to form dragon’s blood by
   performing transcriptomic and metabolomic analysis. These findings are
   important, as they improve our understanding of Dracaena adaptation,
   resource protection, and utilization. Further research is likely to
   have significant benefits for human medicine and health.

Methods

Plant materials, library construction, and genome sequencing

   All plant materials used in this study were 10-year-old adult
   D. cochinchinensis trees growing in the germplasm bank of Dracaena at
   the Yunnan Branch of the Institute of Medicinal Plant Development,
   Chinese Academy of Medical Sciences (Jing-hong City; N:22.0058°,
   E:100.7885°); the plants were formally identified by Professor Hai-tao
   Li.

   Total genomic DNA was extracted from fresh young leaves using the
   DNeasy Plant Mini Kit (QIAGEN Bio-Tec, Hilden, Germany). The DNA was
   sheared by a Covaris M220 focused ultrasonicator (Covaris, Woburn, MA,
   USA) to create fragment sizes of 350 bp. Four types of libraries were
   constructed: (1) a short-read sequencing library with 350-bp insertions
   sequenced with the Illumina HiSeq X Ten instrument platform; (2)
   linked-read sequencing libraries constructed on the 10x Genomics
   GemCode platform and sequenced with the Illumina platform; (3) a SMRT
   cell library with an insert size of 20 kb for PacBio sequencing; and
   (4) a chromatin interaction mapping (Hi-C) library with 350-bp
   insertions for Illumina sequencing.

Genome assembly and assessment of assembly quality

   De novo assembly of PacBio single-molecule long reads from
   SMRT sequencing was performed using FALCON
   ([264]https://github.com/PacificBiosciences/FALCON/) ([265]Chin et al.,
   2013). BWA-MEM ([266]Li, 2013) was used to align the 10x Genomics data
   to the assembly using default settings. Scaffolding was performed by
   fragScaff (v.1.1) with barcoded sequencing reads. These contigs were
   used to form super scaffolds. BWA (v.0.7.8) ([267]Li and Durbin, 2009)
   software was then used to align the clean Hi-C data to the preceding
   assembly. According to the linkage information and restriction enzyme
   sites, the string graph formulation was used to construct a scaffold
   graph with LACHESIS ([268]Burton et al., 2013).

Genome annotation

   Genome repeat sequences were annotated de novo by RepeatMasker
   ([269]http://www.repeatmasker.org). Protein-coding genes were predicted
   based on the repeat-masked genome, homologous conserved proteins, and
   assembled transcripts. microRNA and small nuclear RNA genes were
   predicted by INFERNAL software using the Rfam database (detailed
   methods related to genome annotation are given in [270]Supplemental
   Method 5).

Gene family construction and phylogenetic analysis

   Gene families were generated by OrthoMCL (v.12.0.9) ([271]Li et al.,
   2003). First, nucleotide and protein data from 17 species, including D.
   cochinchinensis, A. officinalis ([272]Harkess et al., 2017), D.
   officinale ([273]Yan et al., 2015), A. shenzhenica, E. guineensis, O.
   sativa (phytozomev10), A. comosus ([274]Ming et al., 2015), M.
   acuminata ([275]D’hont et al., 2012), A. thaliana (phytozomev10), P.
   trichocarpa (release-32), A. sinensis, H. brasiliensis, C. grandis, C.
   papaya (Phytozome.ASGPBv0.4), E. grandis, V. vinifera (phytozomev10),
   and A. trichopoda (release-32) were downloaded from Ensembl (Release
   70), NCBI, and [276]http://marinegenomics.oist.jp/gallery/. Before
   running an “all against all” blastp (E value ≤ 1E−07) program, the
   longest transcript was selected from among the alternative splicing
   transcripts of each gene, and genes predicted to encode ≤50 amino acids
   were removed. To identify homologous gene pairs, more than 30% coverage
   of the aligned regions in both homologous genes was required. Finally,
   the alignments were clustered into gene families using OrthoMCL
   (v.2.0.9) ([277]Li et al., 2003) with a 1.5 inflation index.

   Protein sequences from 252 single-copy gene families were used for
   phylogenetic tree construction. MUSCLE (v.3.7) ([278]Edgar, 2004) was
   used to generate multiple sequence alignment for protein sequences in
   each single-copy family with default parameters. RAxML (v.8.0.19)
   ([279]Stamatakis, 2006) was used to construct the phylogenetic tree.
   Divergence times were estimated using MCMCTree in PAML (v.4.9)
   ([280]Yang, 2007). Five calibration points were selected and applied:
   V. vinifera and O. sativa (130–200 Mya) ([281]Iorizzo et al., 2016);
   O. sativa and A. shenzhenica (137–156 Mya) ([282]Zhang et al., 2017a);
   V. vinifera and C. papaya (114–134 Mya) (Ming et al., 2008);
   P. trichocarpa and A. thaliana (107–109 Mya) ([283]Tuskan et al.,
   2006); and H. brasiliensis and P. trichocarpa (65–86 Mya) ([284]Tuskan
   et al., 2006; [285]Tang et al., 2016).

Expansion and contraction of gene families

   Gene gain and loss in different gene families were determined using
   CAFÉ (v.2.2). The enriched KEGG pathways in the expanded gene families
   of D. cochinchinensis were summarized and visualized using KOBAS
   software.

WGD

   To investigate the genomic evolution of D. cochinchinensis, we
   identified collinearity within the D. cochinchinensis genome and
   between the D. cochinchinensis genome and the genomes of
   A. officinalis, D. officinale, and A. shenzhenica. We searched for
   putative paralogs and orthologs within and between genomes with MCScanX
   based on BLASTP alignment (with an E value ≤1E−5) for each genome pair.
   4DTv values were calculated based on each gene pair in a syntenic
   block, and values of all gene pairs were plotted to identify putative
   WGD events in D. cochinchinensis. We also identified different modes of
   gene duplication as WGDs, tandem duplicates, proximal duplicates,
   transposed duplicates, or dispersed duplicates using DupGen_finder
   ([286]Qiao et al., 2019) with default parameters.

Dynamics of LTR-RTs

   We investigated the dynamics of LTR-RTs in the genomes of
   D. cochinchinensis, A. officinalis, and A. shenzhenica. De novo
   detection of LTR-RTs was performed using LTRharvest (v.1.5.8) (-motif
   tgca -motifmis 1) and LTRfinder (-minlenltr 100 -maxlenltr 5000
   -mindistltr 1000 -maxdistltr 20000) ([287]Ellinghaus et al., 2008).
   Primer-binding site motifs were annotated using LTRdigest (v.1.5.8)
   ([288]Steinbiss et al., 2009), and only elements containing
   primer-binding site motifs were retained for further domain annotation.
   The domains of RTs were identified by searching against HMM profiles
   collated by the Gypsy Database ([289]http://gydb.org/), and the
   internal features of LTR-RTs were annotated with LTRdigest (v.1.5.8)
   ([290]Steinbiss et al., 2009). Intact LTR-RTs that contained domains
   (e.g., the gag domain, protease domain, reverse transcriptase domain,
   and integrase domain) were clustered, and those belonging to families
   with fewer than five members were discarded. The nucleotide divergence
   rate (λ) was calculated from the MUSCLE alignment of the 5′ and 3′ LTR
   sequences of LTR-RTs with the EMBOSS program distmat. The insertion
   date (T) was computed for each LTR-RT (T = K/2r, K: genetic distance,
   K = −0.75ln(1 − 4λ/3)), which was set to a substitution rate (r) of
   1.3E−08^8 substitutions per site per year. The translated amino acid
   sequences were then searched for Ty1/Copia (PF07727) and Ty3/Gypsy
   ([291]PF000078) domains using HMMER (v.3.1b2) ([292]Johnson et al.,
   2010) (E value ≤1E−5). The amino acid sequences of the two
   superfamilies were then aligned using MAFFT (v.7.310) ([293]Katoh and
   Standley, 2013), and phylogenetic trees for Copia and Gypsy LTR-RTs
   were constructed using FastTree
   ([294]http://www.microbesonline.org/fasttree/).

Detection of key candidate functional genes

   We used hmmsearch (HMMER v.3.1b2) ([295]Johnson et al., 2010) and
   BLASTP (E value ≤ 1E−05) (BLAST+ 2.71) ([296]Altschul et al., 1997;
   [297]Johnson et al., 2008) to identify genes involved in the
   biosynthesis of flavonoids based on amino acid sequences from O. sativa
   and genes involved in ROS signaling pathways based on amino acid
   sequences from A. thaliana. The candidates were further curated using
   the NCBI nonredundant protein database, and incorrect proteins were
   removed manually.

RNA sequencing and analysis

   Five tissues (root, stem, leaves, flower, and fruit) were harvested
   from D. cochinchinensis, as well as 33 cut and injury-treated stems at
   different time points (0, 2, 6, 12, and 24 h, 2, 3, 4, 5, 10, 15, and
   30 days, and 6 and 17 months). Three biological replicates of each
   tissue were collected. Raw reads from different samples were filtered
   to obtain high-quality clean reads following standard protocols. For
   gene expression analysis, clean reads were mapped to the
   D. cochinchinensis genome using TopHat (v.2.0.8) ([298]Kim et al.,
   2013). Then, the read counts for each gene were normalized for
   downstream analysis. A gene was defined as expressed if its value of
   reads per kilobase per million reads mapped was ≥1 in at least one of
   the transcriptomes.

   Differential expression analysis between a treatment sample and the
   control sample was performed using the DESeq R package (v.1.18.0). The
   resulting p values were adjusted using Benjamini and Hochberg’s
   approach to control the false discovery rate. Genes with an adjusted p
   <0.05 identified by DESeq were classified as being differentially
   expressed. KEGG pathway enrichment analysis of the differentially
   expressed genes was performed using KOBAS software ([299]Mao et al.,
   2005).

Metabolite profiling based on ultra-high-performance liquid chromatography
(LC) Q-Exactive Orbitrap mass spectrometry (MS)UHPLC-Q-Orbitrap-MS

   Six healthy ∼10-year-old adult trees of D. cochinchinensis were
   selected, and a 6-cm^3 (length: 3 cm; width: 2 cm; depth: 1 cm) wound
   was artificially cut using a small knife in the smooth xylem stem area
   of each tree. Then, 5% (v/v) benzoic acid was sprayed on the wound to
   prevent microbial effects (5 ml each day), and 3 months later, the
   xylem stem with red resin was collected from each tree. We also
   collected a sample of healthy xylem stem surrounding each wound. The
   freeze-dried samples were dissolved in 50% methanol-water (v/v)
   solutions at a concentration of 5 g·ml^−1. The solutions were filtered
   and centrifuged at 13 200 × g for 10 min, and 10 μl of the resulting
   supernatant was directly injected into the LC-MS system.

   Untargeted metabolomics was performed using a U3000 UHPLC equipped with
   a Q-Exactive Q-Orbitrap MS through an HESI source (Thermo Fisher
   Scientific, Waltham, MA, USA). A Waters Acquity UPLC HSS T[3] column
   (2.1 × 100 mm, 1.8 μm) was used for chromatographic separation with a
   column temperature of 35°C. The flow rate was set to 0.4 ml·min^−1, and
   the autosampler temperature was maintained at 4°C.

   The mobile A and B phases were 0.1% formic acid-water (v/v) and
   acetonitrile, respectively. The following elution gradient was used:
   10% B in 0–3 min; 10%–22% B in 3–5 min; 22%–27% B in 5–7 min; 27%–45%
   B in 7–25 min; 45%–57% B in 25–35 min; 57%–100% B in 35–36 min; 100%
   B in 36–37 min; 100%–10% B in 37–37.5 min; and 10% B in 37.5–40 min.

   The HESI source parameters were as follows: the scan mode was
   full-ms/dd-ms mode, and the scanning range was m/z 150–1500 at a
   resolution of 70 000; the spray voltage was 3.0 kV in the positive mode
   and 2.8 kV in the negative mode; the capillary temperature and probe
   heater temperature were 320°C and 350°C, respectively; the sheath gas
   (N[2]) was 35 arbitrary units, and the auxiliary gas was 10 arbitrary
   units (99.9% pure N[2]). The S-lens level value was 50 V. The
   normalized collision energy was 20/40/60 V. Data acquisition was
   performed using Xcalibur 4.1 software (Thermo Fisher Scientific).

   The metabolites were identified or tentatively identified as follows:
   first, we used a full MS/dd-MS[2] (TopN) scan method that included a
   precursor ions list involving all of the secondary metabolites ever
   isolated from Dracaena, and the “If idle-pick others” function was used
   to sensitively characterize the metabolites within one injection
   analysis. Literature on the phytochemistry of Dracaena retrieved from
   CNKI ([300]www.cnki.net), PubChem
   ([301]http://pubchem.ncbi.nlm.nih.gov), Web of Science
   ([302]www.webofknowledge.com), ChemSpider ([303]www.chemspider.com),
   and Chemicalbook ([304]www.chemicalbook.com) was summarized (up to
   January 2022), and all the isolated compounds were included in the
   precursor ions list. Second, for compounds not annotated by the above
   method, we used Compound Discoverer 3.1 (Thermo Fisher Scientific), a
   software tool for untargeted metabolomics, to achieve unknown
   identification using the following parameters: polarity,
   negative/positive; maximum shift, 0.2 min; mass tolerance, 5 ppm;
   minimum peak intensity, 100 000. The compounds were annotated with
   different types of databases, including m/z Cloud, ChemSpider, and the
   CD internal database. After data preprocessing, a peak table was
   generated containing m/z, retention time, and peak areas of each
   compound in all samples. The identification levels for all compounds
   were in accordance with the Metabolomics Standards Initiative
   ([305]Blaženović et al., 2018).

Data processing of metabolomic assays

   Data processing of the metabolomic assays was performed using SIEVE 2.2
   (Thermo Fisher Scientific) for background subtraction and component
   extractions. The obtained peak list was further processed by principal
   component analysis and Student’s t-tests to compare changes in compound
   composition before and after wounding stress. Principal component
   analysis was performed with the SIMCA-P program v.14.1 (Umetrics, Umea,
   Sweden). The paired-sample Student’s t-tests were performed using
   Office Excel 2010 (Microsoft, Redmond, WA, USA), and the obtained p
   values were corrected using a false discovery rate (FDR) calculated in
   the R programming language. The null hypothesis was that the content of
   a compound showed no significant difference before and after the
   formation of dragon’s blood (FDR < 0.01). SIMCA-P 14.1 was also used to
   transform data for orthogonal partial least squares discriminant
   analysis. To ensure the reliability of the data analysis, variable
   importance in projection values, which are the weighted values of
   variation obtained from orthogonal partial least squares discriminant
   analysis, were used as another measure to screen for differential
   compounds. Metabolites with a variable importance in projection >1.0
   and FDR <0.01 were chosen for further analysis of differential
   compounds before and after wounding stress. Multiple Experiment Viewer
   software ([306]http://mev.tm4.org/) was used for the generation of
   heatmaps. Prism 6 (GraphPad, San Diego, CA, USA) was used to produce
   box plots and line graphs.

Funding

   This work was supported by the program of the CAMS Initiative for
   Innovative Medicine (2021-1-I2M-032), the National Natural Science
   Foundation of China (82173925 and 81573525), and the National Key
   Research and Development Program of China (2018YFC1706401).

Author contributions

   Y.X. and J.W. conceived and designed the project. Y.X. organized the
   manuscript and supervised and participated in all the related analyses.
   K.Z. managed the bioinformatics analyses. Z.Z. conducted metabolome
   analysis. Y.L. carried out the H[2]O[2] determination experiment. Y.X.,
   Z.Z., and Y.L. contributed to sample preparation for genome and RNA
   sequencing. F.L., P.S., S.G., Q.W., C.Y., J.J., C.L., M.S., Z.G., C.S.,
   H.L., Y.J., and X.G. provided constructive suggestions. All authors
   contributed to the writing of the paper.

Acknowledgments