Transcriptome and Metabolome Analysis of the Synthesis Pathways of Allelochemicals in Eupatorium adenophorum.

Eupatorium adenophorum (Crofton weed) is an invasive weed in more than 30 countries. It inhibits the growth of surrounding plants by releasing allelochemicals during its invasion. However, the synthetic pathways and molecular mechanisms of its allelochemicals have been rarely reported. In this study, the related genes and pathways of allelochemicals in E. adenophorum were analyzed. Transcriptome analysis showed that differentially expressed genes (DEGs) were mainly enriched in the phenylpropanoid biosynthetic pathway and flavonoid biosynthetic pathway. Thirty-three DEGs involved in the synthesis of allelochemicals were identified, and 30 DEGs showed significant differences in blades and stems. Six allelochemicals were identified from blades and stems by ultraperformance liquid chromatography-tandem mass spectrometry (UPLC-MS/MS). Correlation analysis of genes and metabolites showed a strong correlation between the five genes and allelochemicals. In addition, this study supplemented the biosynthetic pathway of Eupatorium adenophorum B (HHO). It was found that acyclic sesquiterpene synthase (NES), δ-cadinene synthase (TPS), and cytochrome P450 (P450) were involved in the synthesis of HHO. These findings provide a dynamic spectrum consisting of allelochemical metabolism and a coexpression network of allelochemical synthesis genes in E. adenophorum.

Introduction

Eupatorium adenophorum (Crofton weed) is a perennial shrub of the genus Eupatorium in the dicotyledonous family, which belongs to the family Asteraceae.[1] It is native to Mexico and Costa Rica, has been introduced to Europe, Australia, and Asia as an ornamental plant, and is gradually expanding to the wild.[2] The distribution of E. adenophorum ranges from Spain at 37°N to South Africa and Australia at 35°N. The invasion of this plant has caused great damage to the economy and biodiversity, seriously affecting the ecological balance of these areas. In the Convention on the Control of Exotic Pests and International Biological Control, E. adenophorum has been listed as one of the four pernicious weeds. The invasion of E. adenophorum causes a decrease in the Simpson diversity and Shannon–Wiener diversity indices of species under different habitat conditions.[3] It is estimated that the invasion of E. adenophorum causes annual losses of 99 billion CNY (1.55 billion USD) in the grazing industry in China and 263 billion CNY (4.13 billion USD) in the service function of grasslands and forests.[4]

Allelopathy is an important weapon for the successful invasion and rapid spread of invasive alien plants. As a worldwide malignant weed, E. adenophorum also has strong allelopathy. It releases allelochemicals into the external environment to affect the growth and development of surrounding plants, thereby gaining an advantage over the competition and allowing its populations to grow and expand rapidly.[5−8] For example, aqueous extracts of E. adenophorum strongly inhibit seed germination and the growth of ryegrass (Lolium perenne L.), wheat (Triticum aestivum L.), maize (Zea mays L.), and dry rice (Oryza sativa L.).[9−12] The blades and stems have been shown to have significant allelopathy, with different organs having different allelopathy potentials. It has also been shown that sesquiterpenes (HHO), phenolic acids (cinnamic acid, ferulic acid), coumarins (coumarin), and flavonoids (catechins, epicatechin) are potential allelochemicals in E. adenophorum.[13,14] Among them, HHO is defined as the main allelochemical. The malondialdehyde content and peroxidase activity of plants are significantly increased at certain concentrations of HHO. In addition, the concentrations of some endogenous hormones, such as abscisic acid (ABA), indole-3-acetic acid (IAA), and zeatin (ZR), are also affected by HHO.[15,16] However, the biosynthetic pathways of these allelochemicals in E. adenophorum have not been elucidated. Given the severe damage caused by E. adenophorum invasion, investigating the synthesis pathways of the allelochemicals in E. adenophorum has become an urgent task for comprehensive management. This study aims to explore the biosynthetic pathways of allelochemicals and their associated regulatory networks in E. adenophorum through a combined transcriptome and metabolome analysis. The results can contribute to the further understanding of the biosynthesis process of allelochemicals in E. adenophorum.

Results and Discussion

Statistics of Transcriptome Sequencing Data

Twelve cDNA libraries (blades 1–3, petioles 1–3, stems 1–3, and roots 1–3) were sequenced using the Illumina HiSeq platform.[17] After removing the low-quality reads, 83.19 Gb of clean data was obtained. The clean data of each sample reached more than 6 Gb, and the Q30 base was above 86%, indicating that the sequencing quality was acceptable (Table ).

Table 1

Data Statistics of Filtered E. adenophorum Transcriptome Sequencing Samples

sample	raw reads	clean reads	clean base (G)	error rate (%)	Q20 (%)	Q30 (%)	GC content (%)
blade-1	49 948 118	44 938 498	6.74	0.04	93.49	86.67	44.20
blade-2	47 103 618	43 584 774	6.54	0.05	93.32	86.31	44.31
blade-3	62 254 956	56 443 184	8.47	0.04	93.52	86.70	44.43
petiole-1	52 923 698	48 428 776	7.26	0.05	93.34	86.30	44.05
petiole-2	53 819 518	49 142 258	7.37	0.04	93.55	86.63	44.05
petiole-3	48 054 422	41 227 364	6.18	0.04	94.25	87.80	43.96
root-1	47 703 320	43 811 214	6.57	0.04	93.50	86.62	43.57
root-2	46 501 406	42 586 210	6.39	0.04	94.14	87.51	43.30
root-3	52 838 988	48 958 758	7.34	0.04	93.93	87.00	43.22
stem-1	55 142 368	50 324 762	7.55	0.04	93.60	86.70	43.74
stem-2	46 604 292	42 154 430	6.32	0.05	93.41	86.46	43.86
stem-3	47 547 798	43 093 014	6.46	0.04	94.09	87.43	43.49

Through the Trinity assembly program,[18] short-read sequences were assembled into 380 569 transcripts with an average length of 819 bp, and the length of N50 was 1210 bp. After further filtering of the low-expressed transcripts, 331 049 unigenes with an average length of 902 bp were obtained, and the length of N50 was 1267 bp. Unigenes with lengths of between 300 and 400 bp accounted for the largest percentage (46 847, 14.15%). The percentages of unigenes with lengths of between 200 and 300, 400 and 1000, and 1000 and 2000 bp were 13.11% (43 391), 42.23% (39 813), and 21.89% (72 453), respectively. In addition, 28 545 (8.62%) unigenes were greater than 2000 bp in length (Figure A).

Figure 1

Transcriptome sequencing analysis of E. adenophorum. (A) Length distribution of transcripts and unigenes in E. adenophorum. (B) Distribution of sequence alignment results in the NR database. (C) Functional annotation of unigenes based on the GO database. (D) Functional annotation of unigenes based on the KOG database. (E) Statistics of DEGs in four organs of E. adenophorum.

Annotation and Classification of Transcriptome Sequencing Results

To speculate on the function of unigenes of E. adenophorum, unigene sequences were compared with the KEGG, NR, Swiss-Prot, GO, COG, KOG, and Trembl databases using BLAST software.[19−26] Unigenes were translated into amino acid sequences and then compared with the Pfam database using HMMER software to obtain annotation information on unigenes.[27,28] The results showed that 216 287 unigenes (65.33%) were annotated in at least 1 database. The highest annotation rate was obtained in the NR database, which assigned 214 036 (64.65%) unigenes. In other databases, 150 874 (45.57%) unigenes were annotated in KEGG, 140 483 (42.44%) were annotated in SwissProt, 206 164 (62.28%) were annotated in Trembl, 120 729 (36.47%) were annotated in KOG, 177 740 (53.69%) were annotated in GO, and 149 261 (45.09%) were annotated in Pfam.

Compared to the data in the NR database, the gene sequences of E. adenophorum were most similar to those of a sunflower (Helianthus annuus, 60.47%), followed by quercus variabilis (Quercus suber, 10.33%), lettuce (Lactuca sativa, 7.28%), artichoke (Cynara cardunculus var. scolymus, 6.97%), japonica rice (O. japonica group, 2.27%), and rice (O. indica group, 0.77%) (Figure B).

The GO annotation indicated that 177 740 unigenes were categorized into 60 functional terms. In the group of cellular components, the cell (119 889), cell part (119 598), and organelle (90 361) were mainly involved. In the molecular function group, binding (106 387), catalytic activity (92 678), and transporter activity (12 156) were mainly involved. In the biological process group, cellular processes (104 109), metabolic processes (12 156), and the response to stimulus (46 249) were mainly involved (Figure C).

The unigenes obtained by sequencing were compared with the KEGG database, and 150 874 annotated unigenes assigned to 145 biological pathways were obtained, including transcriptional regulation, signal transduction, translational regulation, substance metabolism, and secondary metabolite biosynthesis. Among these pathways, d-arginine and d-ornithine metabolic pathways (KO00472) were analyzed, and they involved the largest number of unigenes. Several pathways related to substance metabolism were identified, including carbon metabolism, amino acid metabolism, lipid metabolism, and secondary metabolism. Because most allelochemicals in E. adenophorum are secondary metabolites, metabolic pathways related to secondary metabolite synthesis were further analyzed, including the phenylpropanoid biosynthesis pathway (KO00940), the flavonoid biosynthesis pathway (KO00941), and the sesquiterpene and triterpenoid biosynthesis pathway (KO00909). A total of 17 744 genes were involved in secondary metabolism.

On the basis of the KOG databases, 120 729 unigenes of E. adenophorum were categorized into 25 functional groups. The largest proportion of the grouping was “general function prediction only”, followed by “signal transduction mechanisms”, “post-translational modifications”, “protein turnover”, and “chaperones”. In addition, only 57 unigenes were categorized as having “cell motility” (Figure D).

In this study, transcriptome sequencing and the preliminary analysis of data from blades, petioles, roots, and stems of the invasive plant E. adenophorum were assembled de novo into 331 049 unigenes. Compared to other species, including Euphorbia fischeriana (18 180), Paeonia suffruticosa (72 997), and Picrorhiza kurrooa (74 336), the number of unigenes and possible protein sequences obtained in this study was higher.[29−31] Possible reasons are (a) the genes with specific coding sequences in E. adenophorum may differ significantly from currently available protein sequences and cannot be effectively annotated; (b) E. adenophorum lacks a reference for genomic information, and no genome of a single species has been resolved even in the Asteraceae. Therefore, when short sequence assembly is performed without reference genome guidance, splicing errors are prone to occur in short contigs joined into long contigs. This phenomenon results in sequences belonging to the same gene not being spliced together, causing an increase in the number of unigenes. The transcriptomes of four organs of E. adenophorum were studied, and their complete transcriptome sequencing data were obtained, providing a basis for further molecular biology and genomics studies of E. adenophorum.

Differential Expression Analysis of E. adenophorum Transcripts in Different Organs

On the basis of FPKM values, changes in the transcript levels of unigenes in each group were analyzed to identify DEGs between groups (FDR < 0.01 and |log2FC| ≥ 2). There were 20 900–49 974 DEGs, 7477–24 821 up-regulated genes, and 13 423–27 277 down-regulated genes in each group. In addition, there were 21 007 and 18 555; 22 390 and 27 277; 24 821 and 25 153; 12 135 and 20 891; 7477 and 13 423; and 11 597 and 151 124 up-regulated and down-regulated differential genes in blade vs petiole, blade vs root, blade vs stem, petiole vs root, petiole vs stem, and stem vs root, respectively. Compared to other groups, blade vs stem had the highest number of up-regulated and down-regulated genes (24 821 and 25 153). The group of petiole vs stem was the lowest (7477 and 13 423), and blade vs stem was better enriched for DEGs (Figure E).

The KEGG-based enrichment analysis showed that all DEGs were enriched in 143 metabolic pathways, with metabolic pathways (46.41%) and secondary metabolic pathways (26.38%) being in the top 2 of each group (Supporting Information Figure 1). Flavonoid biosynthesis and phenylpropanoid biosynthesis were significantly enriched in the top 20 KEGG pathways of each group. Sesquiterpene biosynthesis and triterpene biosynthesis were enriched in petiole vs root and stem vs root. GO enrichment was divided into three major categories: biological processes, cellular components, and molecular functions (Supporting Information Figures 2 and 3). In biological processes, the enrichment of DEGs was higher for the metabolic process of blade vs stem (18 628), indicating that some important metabolic activities differed among the six groups. Many DEGs were enriched in the metabolic pathways known to be associated with allelopathy, flavonoid biosynthetic pathways were enriched in blade vs stem and stem vs root, and phenylpropanoid catabolic processes were enriched in petiole vs root and petiole vs stem.

There are 39 562, 49 667, and 49 974 DEGs in blade vs petiole, blade vs root, and blade vs stem, respectively, with 24 747 overlapping DEGs (Figure A). The DEGs in petiole vs root, petiole vs stem, and stem vs root are 33 026, 20 900, and 26 721, respectively (Figure B), with 4813 overlapping DEGs. The most genes were shared among the three groups of blade vs petiole, blade vs root, and blade vs stem, and the least genes were shared among the two groups of blade vs root and petiole vs stem, with a total of 852 DEGs in the 6 groups. In addition to common DEGs in four different organs, each organ has its own specific genes (e.g., the blade vs stem grouping has 4805 unique DEGs), indicating that many differential genes were detected in the grouping of different sampled organs with gene specificity and organ specificity (Figure C).

Figure 2

Distribution of DEGs in E. adenophorum. (A) DEGs in blade vs petiole, blade vs root, and blade vs stem. (B) DEGs in petiole vs root, petiole vs stem, and stem vs root. (C) DEGs in blade vs petiole, blade vs root, blade vs stem, petiole vs root, petiole vs stem, and stem vs root. (Connected black dots represent common genes within groups.)

Analysis of DEGs in the Biosynthetic Pathway of Allelochemicals

Cinnamic acid, coumarin, and ferulic acid are phenylpropanoid metabolites, catechin and epicatechin are flavonoids metabolites, and HHO is a sesquiterpene metabolite.[32−34] Five DEGs were identified in the phenylpropanoid synthesis pathway (p < 0.05): two genes (PAL1 and PAL2) annotated as phenylalanine amino lyase (PAL), six β-glucosidase (bglx) genes (bglx1, bglx2, bglx3, bglx4, bglx5, and bglx6), two cinnamate-4-hydroxylase (C4H) genes (C4H1 and C4H2), two 4-coumaroyl-CoA ligase (4CL) genes (4CL1 and 4CL2), and three caffeic acid O-methyltransferase (COMT) genes (COMT1, COMT2, and COMT3). Most genes were highly expressed in roots and blades. The expression patterns of bglx2, bglx3, 4CL2, COMT1, COMT2, and C4H2 genes tend to be consistent in the four organs, and the expressions were statistically different (p < 0.01) in both blades and roots, with relatively high expressions in blades. The expression patterns of PAL1, PAL2, bglx4, bglx5, and C4H1 genes tend to be consistent in the four organs, showing statistically significant differences in the expression in roots and blades (p < 0.01), with relatively high expression in roots and a high expression of bglx1 in petioles (Figure A). The expressions of the PAL1 gene are 3.86, 3.33, and 2.57 times higher in roots, stems, and petioles than in blades, respectively. In addition, C4H2 gene expression in blades is 4.33, 3.96, and 1.78 times higher than in roots, stems, and petioles, respectively.

Figure 3

Allelochemical synthesis pathways of E. adenophorum. (A) Biosynthesis pathways of cinnamic acid, coumarin, ferulic acid, catechin, and epicatechin. (B) Speculation pathways of HHO biosynthesis. Red and green represent high and low gene expressions, respectively.

In the flavonoid pathway, a total of six DEGs were identified in the synthesis of allelochemicals (p < 0.05): a gene annotated as chalcone synthase (CHS), two chalcone isomerase (CHI) genes (CHI1 and CHI2), a naringenin 3-dioxygenase (F3H) gene, a flavonoid 3′-monooxygenase (CYP75B) gene, and a flavanone 4-reductase (DFR) gene. The expression pattern of DEGs of the flavonoid pathway tends to be consistent in the four organs, and the expression was statistically different in stems and blades (p < 0.01), with a relatively high expression in stems. Compared to the DFR expression in the blades, the expressions were 8.88, 7.78, and 2.08 times higher in the stems, petioles, and roots, respectively.

Among the sesquiterpene synthesis pathways, the synthesis pathway of HHO was speculated on the basis of the relevant metabolic pathways in the KEGG database and related references.[35] HHO is a cadinene-type sesquiterpene. In the sesquiterpene synthesis pathway, farnesyl pyrophosphate (FPP) is catalyzed by NES1 to produce nerolidol and by TPS1 to produce δ-cadinene. After a series of oxidations and hydroxylations, HHO is finally obtained.[36] The enzymes involved in this process remain unclear. Eight DEGs were identified during the synthesis of allelochemicals in the sesquiterpene pathway (p < 0.05): one gene annotated as acetyl-coenzyme A acyltransferase (ACAT), one hydroxymethylglutaryl-coenzyme A synthase (HMGS) gene, three methylglutaryl-coenzyme A reductase (hydroxymethylglutaryl-CoA reductase (HMGCR) genes (HMGCR1, HMGCR2, and HMGCR3), a phosphomevalonate kinase (mvaK2) gene, a diphosphomevalonate decarboxylase (diphosphomevalonate decarboxylase (MVD) gene, a 2Z,6Z-farnesyl diphosphate synthase (ZFPS) gene, a NES1 gene, and a TPS1 gene. The expression patterns of ACAT, HMGS, HMGCR2, mvaK2, MVD, and ZFPS tend to be consistent in the four organs and were statistically different (p < 0.01) in both the stems and blades, with relatively high expression in the stems (Figure B). HMGCR1 and HMGCR3 were highly expressed in the root, NES1 and TPS1 were highly expressed in the blades, and P4501 and P4502 were highly expressed in the petioles. As a key step in the synthetic HHO biosynthetic pathway, TPS1 was expressed 9.20, 5.71, and 5.64 times higher in blades than in roots, stems, and petioles, respectively.

Transcripts from the four organs were explored by a two-by-two comparison of transcriptome data. The results of the comparative analysis of DEGs show that E. adenophorum is enriched in specifically expressed genes involved in the phenylpropanoid synthesis pathway, flavonoid synthesis pathway, and sesquiterpene synthesis pathway. Most ZFPS upstream genes in the sesquiterpene synthesis pathway follow a trend of high expression in stems and low expression in blades. However, three upstream genes of HHO synthesis show a relatively high expression in blades, such as NES, TPS, and P4501. The upstream genes of ZFPS may be involved in the biosynthesis of other terpenoids. In contrast, the expression of MVA pathway genes is relatively high in the roots of the perennial plant Cyanotis arachnoidea. In ginsenoside-producing plants, the core genes of the MVA pathway show a relatively high expression in the roots.[37,38] These two findings indicate that the general biology of these plants differs from that of E. adenophorum.

Verification of the Expression Levels of Selected Genes by qRT-PCR

A total of 20 E. adenophorum DEGs were selected for qRT-PCR analysis. Specifically, there are nine genes of the phenylpropanoid synthesis pathway (PAL1, PAL2, bglx1, bglx3, bglx6, 4CL1, C4H1, C4H2, and COMT2), four genes of the flavonoid synthesis pathway (CHS, CHI1, F3H, and CYP75B), and seven genes of the sesquiterpene synthesis pathway (ACAT, HMGS, HMGCR1, NES, TPS, P4501, and P4502). Pearson correlation coefficients were calculated by SPSS. The results show that the relative expression pattern of DEGs is similar to that of transcriptome sequencing data, proving that the transcriptome data is reliable (Supporting Information Table 1). It can also be seen that NES1 and TPS1 are highly expressed in the blades (Figure ).

Figure 4

Changes in gene relative expressions in E. adenophorum. The gene expression levels of PAL1, PAL2, bglx1, bglx3, bglx6, 4CL1, C4H1, C4H2, COMT2, CHS, CHI1, F3H, CYP75 B, ACAT, HMGS, HMGCR1, NES, TPS, P4501, and P4502 in four organs were analyzed by qRT-PCR. Different lowercase letters represent significant differences in organs (p < 0.05), and different uppercase letters represent highly significant differences in organs (p < 0.01). (Error strips are used to describe the deviation among three biological repeats in the same organ.)

Identification of Metabolites

The most differential genes were enriched in the blade vs stem group of E. adenophorum. To further understand the molecular mechanisms of allelochemicals in E. adenophorum, a metabolomic analysis of blades and stems was performed with the UPLC-MS/MS system, with 667 metabolites identified. Moreover, PCA showed that blades, stems, and mixed samples were significantly separated, and the differences among organs were significantly greater than the differences between the mixed samples and each organ (Figure A). In addition, biological replicates were projected to be spatially close to each other, indicating a good correlation between replicates. For variables with low correlation, sensitive OPLS-DA maximizes the distinction between the blade and stem differences (Figure B).

Figure 5

Metabolomics analysis of E. adenophorum. (A) PCA of blade vs stem metabolism. (B) OPLS-DA of blade vs stem metabolism. (C) Proportion of identified metabolites in different metabolic pathways. (D) Proportion of identified metabolites in secondary metabolic pathways. (E) DAMs classification map of the E. adenophorum blade vs stem metabolic group. (F) First 20 KEGG pathways enriched by DEGs in blade vs stem for E. adenophorum. (G) First 20 KEGG pathways for DMA enrichment in E. adenophorum blade vs stem.

Using the KEGG database and the Plant Metabolic Pathway Databases (https://plantcyc.org/),[39] metabolites were classified into 12 categories according to metabolic pathways: carbohydrate metabolism (68), lipid metabolism (43), energy metabolism (13), nucleotide metabolism (28), amino acid metabolism (138), metabolism of terpenoids and polyketides (7), biosynthesis of other secondary metabolites (120), membrane transport (36), signal transduction (5), translation (17), metabolism of cofactors and vitamins (30), and metabolites that cannot be classified (162) (Figure C). The top three metabolic pathways for secondary metabolite enrichment are flavonoid biosynthesis (Ko00941) (22), phenylpropanoid biosynthesis (Ko00940) (20), and flavonoid and flavonol biosynthesis (Ko00944) (19) (Figure D). The metabolites were classified into various metabolic pathways: catechin and epicatechin were involved in the flavonoid biosynthetic pathway, and cinnamic acid, coumarin, and ferulic acid were involved in the phenylpropanoid biosynthetic pathway.

Differential Metabolite Analysis

In the comparative analysis of blade and stem metabolomes of E. adenophorum, there are 394 differential metabolites. Compared to the blade, 209 and 185 metabolites are up-regulated and down-regulated in the stem, respectively. More metabolites are up-regulated than down-regulated, indicating that most metabolites are efficiently accumulated in the stems (Figure E).

The categories with the largest number of differential metabolites are flavonoids, lipids, and phenolic acids. There are 1.73 times more up-regulated metabolites of flavonoids than down-regulated metabolites, and steroids show almost no difference between blades and stems. The phenylpropanoid pathway and flavonoid pathway are significantly enriched in the transcriptome and metabolome of blades and stems (Figure F,G). The numbers of DAMs identified in the metabolism of phenylpropanoids, flavonoids, and sesquiterpenes are 13, 17, and 5, respectively (Figure A–C). The accumulated contents of the allelochemical cinnamic acid, coumarin, ferulic acid (all involved in phenylpropanoid biosynthesis), catechin, and epicatechin (all involved in flavonoid biosynthesis) and the accumulated contents among blades and stems are 3709.73, 23 204; 7104.97, 41  826.33; 73 553.33, 164 183.33; 169 806.67, 6 851 166.67; and 363 393.33, 12 347,000, respectively (Supporting Information Figure 4). The accumulation in the stems is 2.64, 2.56, 1.16, 5.33, and 5.09 times higher than in the blades, and the accumulation of ferulic acid shows the smallest difference. The intermediate products of p-coumaric acid and methyl eugenol of the phenylpropanoid pathway have the largest differences in relative content accumulation between the blade and stem groups of E. adenophorum (13.50 and 12.25 times higher, respectively). The accumulation of these two metabolites is higher in blades. Additionally, five intermediate metabolites of the phenylpropanoid pathway are accumulated at a higher level in the stem (such as coniferyl alcohol), and eight other metabolites are accumulated at a higher level in the blade (such as 2-hydroxycinnamic acid). Six intermediate metabolites of the flavonoid pathway are better accumulated in stems (e.g., luteolin), and 11 metabolites are better accumulated in blades (e.g., naringenin). The accumulation of the sesquiterpene pathway HHO in blades and stems is 30 140 666.67 and 315 583.33, respectively. It increases 6.58 times more in blades than in stems and has a higher accumulation in blades, indicating that the blade is the main effective organ for allelopathy. The accumulation of the other four sesquiterpenes is also higher in blades.

Figure 6

Analysis of DAMs and DEGs in E. adenophorum. (A–C) DAMs identified in the phenylpropanoid metabolic pathway, DAMs identified in the flavonoid metabolic pathway, and DAMs identified in the sesquiterpene metabolic pathway, respectively. (Red and blue circles indicate that the contents of metabolites in the stem metabolism group are up-regulated and down-regulated, respectively, compared to the content in the blade metabolism group (VIP > 1).) The circle size indicates the difference multiple (FC). (D) Protein interaction of differential genes in the metabolic pathways of allelochemicals. (E) Correlation network analysis among six differential sensing substances. Purple circles represent DAMs involved in phenylpropanoid metabolism, yellow circles represent DMAs involved in flavonoid metabolism, and green circles represent DMAs involved in sesquiterpene metabolism. The round and triangular borders represent the up- and down-regulation of metabolites, respectively. Red and blue lines represent positive (PPC > 0) and negative (PPC < 0) correlations of metabolites with genes. (F–H) Correlation network diagrams of phenylpropanoid metabolic pathway DAMs and DEGs, flavonoid metabolic pathway DAMs and DEGs, and sesquiterpene metabolic pathway DAMs and DEGs. Circles and pentagons represent genes and metabolites, respectively. Pink circles represent DEGs involved in phenylpropanoid metabolism, yellow circles represent DEGs involved in brass metabolism, and brown circles represent DEGs involved in sesquiterpene metabolism. Red and black borders indicate the up- and down-regulation of genes and metabolites, respectively. Red and black lines represent metabolites and genes that are positively correlated (PPC > 0) and negatively correlated (PPC < 0), respectively.

Protein Interactions of Differential Genes and the Regulatory Networks between Differential Metabolites

A total of 21 proteins of DEGs were matched in the String database,[40] with an average interaction score of 0.83. Two interaction networks were identified. PAL1, PAL2, COMT2, 4CL1, CYP75B, CHI2, F3H, CHS, CHI1, DFR, HMGS, ZFPS, NES, HMGCR1, MVD, ACAT, mvak2, and C4H1 proteins interacted with each other; bglx2, bglx4, and bglx6 proteins interacted with each other. The interaction score of NES and ZFPS proteins is 0.661, and there is a protein interaction relationship (Figure D). Cinnamic acid, coumarin, ferulic acid catechin, and epicatechin are all positively correlated with each other and are negatively correlated with HHO (Figure E).

Correlation Analysis of Differential Genes and Differential Metabolites

To explore the association between gene expressions and metabolite accumulation patterns in the blade and stem groups of E. adenophorum, PCC for the correlation between DEGs and DAMs involved in the metabolic pathways of phenylpropanoids, flavonoids, and sesquiterpenes were calculated, and a network diagram of correlations between genes and metabolites was constructed (Figure F–H). The 30 identified DEGs involved in the synthesis of allelochemicals are significantly correlated with chemosensitive DAMs (p < 0.05, PCC > 0.7). When cinnamic acid, coumarin, ferulic acid, catechin, and epicatechin are positively correlated with DEGs in the phenylpropanoid, flavonoid, and sesquiterpene metabolic pathways, respectively, gene expression and metabolite accumulation patterns will be up-regulated for PAL1, PAL2, bglx4, bglx5, bglx6, C4H1, COMT3, CHS, CHI1, CHI2, F3H, CYP75B, ACAT, HMGS, HMGCR1, HMGCR2, HMGCR3, mvak2, MVD, and ZFPS. All were highly accumulated in stems. The expression of genes and metabolite accumulation patterns are inconsistent when a negative correlation exists. Genes (such as bglx2, 4CL1, 4CL2, C4H2, COMT2, NES, TPS, and P4502) have a higher accumulation in blades, and metabolites have a higher accumulation in stems. The genes with a high correlation with cinnamic acid are bglx4 (correlation coefficient 0.9906), HMGCR1 (0.9850), CYP75B (0.9841), HMGCR3 (0.9839), and ZFPS (0.9824). The genes that have a high association with coumarins are bglx4 (correlation coefficient 0.9982), CYP75B (0.9964), ZFPS (0.9960), HMGCR2 (0.9948), and HMGCR1 (0.9947). The genes with a high correlation with ferulic acid are PAL1 (correlation coefficient 0.9839), ZFPS (0.9818), bglx4 (0.9809), HMGCR1 (0.9804), and F3H (0.9803). The genes with a high correlation with catechins are ZFPS (correlation coefficient 0.9993), PAL1 (0.9993), CYP75B (0.9992), HMGCR2 (0.9987), and CHS (0.9986). The genes with a high correlation with epicatechin are CHI2 (correlations efficient 0.9993), HMGCR2 (0.9992), HMGS (0.9986), ACAT (0.9986), and CHS (0.9978). When HHO is positively correlated with DEGs in the metabolic pathways of phenylpropanoids, flavonoids, and sesquiterpenes, the expression of genes and the accumulation pattern of metabolites are both down-regulated and have higher accumulation in blades. When HHO is negatively correlated with DEGs, the expression of genes and the accumulation pattern of metabolites are inconsistent. Genes positively associated with HHO are negatively correlated with the above five metabolites, and genes negatively associated with HHO are positively correlated with the above five metabolites. There is a strong correlation between HHO and TPS1 (0.9929), NES1 (0.8912), and P4502 (0.8787) genes. The genes with a strong correlation with HHO are ranked in terms of the correlation coefficient as HMGCR2 (0.9999), HMGS (0.9995), CYP75B (0.9994), CHS (0.9992), and CHI2 (0.9991). Genes with a strong correlation with at least three allelochemicals are bglx4, CHS, HMGCR1, HMGCR2, CYP75B, and ZFPS, whose expressions are up-regulated 1.17, 2.00, 2.00, 2.46, 2.17, and 3.14 times in the blade and stem groups, respectively, and a higher accumulation is found in the stem. DFR is the most up-regulated gene, with an 8.88 times higher relative expression in stems than in blades, and is strongly correlated with all six allelochemicals (PCC > 0.9).

Currently, comparative analysis of the expression trends of six allelochemicals in the blades and stems of E. adenophorum has been poorly reported in the literature. Cinnamic acid, coumarin, ferulic acid, catechin, and epicatechin accumulate at relatively high levels in stems. The main effect of the allelochemicals’ HHO is high accumulated content in blades, which is consistent with the previously reported result that the blade is the main affected organ in allelopathy.[41] In this article, a complete synthetic pathway of allelochemicals was reported for the first time. HHO belongs to the cadinene-type sesquiterpene and is generated through the terpenoid pathway. The HHO prediction steps are speculated on according to the literature. The predicted pathway is that FPP is catalyzed by NES1 to produce (S,E)-nerolidol and then catalyzed by TPS1 to produce δ-cadinene, which in turn is catalyzed by a series of P450 enzymes to produce HHO.[42] Genes in the predicted steps (NES1, TPS1, P4501, and P502) were also identified in the biosynthetic pathway of KEGG sesquiterpenes and triterpenes. Among them, NES1, TPS1, and P4501 have a positive correlation with HHO and all have a high accumulation in blades. Real-time fluorescence validation results are consistent with transcriptomic data. NES1 and terpenoid pathway ZFPS genes have protein interactions. The above results support the predicted HHO synthesis pathway.

It is found that the accumulation and gene expression levels of metabolites in the blades and stems of E. adenophorum are not completely consistent, suggesting that genes involved in substance synthesis are regulated by a complex network. In the gene and metabolite correlation analysis, strong correlations were found among bglx4, CHS, HMGCR1, HMGCR2, CYP75B, ZFPS, and allelochemicals. It is assumed that these genes play a crucial role in allelopathy. It has been shown that the expression levels of two genes, CHS and HMGCR, differ in the blades of E. adenophorum mays treated with different doses of HHO. This result indicates that the expression levels of CHS and HMGCR are correlated with the expression levels of HHO, which is consistent with the result of this study.[43]

Under natural conditions, the secretions of the aboveground part of E. adenophorum dropped onto the surrounding environment through the leaching of rainwater, dew, and fog. When accumulated to a certain threshold, they would have allelopathy on adjacent plants, inhibit the seed germination of adjacent plants, and affect the growth of the radicle and hypocotyl of seedlings.[5] In this study, we conducted a joint analysis of the molecular mechanism of allelochemical synthesis in transcriptome and metabolomics, which will provide a unique opportunity for us to obtain candidate genes related to allelochemical synthesis so as to ultimately reveal the molecular mechanism of E. adenophorum invasion. Moreover, in application, the reasonable prevention and control of E. adenophorum invasion can be realized through the detection of genes, and the optimal management and control strategies can be formulated to reduce the economic losses caused by E. adenophorum and the damage to the ecosystem.

Materials and Methods

Plant Materials and Extraction of Total RNA

The test material was purple-stemmed zelenia distributed in Baise City, Guangxi Zhuang Autonomous Region, China. The three replicate samples were all grown in a sugar cane field for 2 years and have a height of about 40 cm. Total RNA was extracted from the blades, petioles, stems, and roots of E. adenophorum mays using an RNAprep Pure Plant Total RNA Extraction Kit (Tiangen Biochemical Technology Co., Ltd.), and the mRNA was purified. Three biological replicates were set for each sample.

Transcriptome Sequencing

The cDNA library was constructed by the reverse transcription of mRNA from each sample of E. adenophorum and sequenced with a Illumina HiSeq sequencing platform. After removing splice sequences and low-quality reads from the sequencing data and performing data filtering, clean reads were obtained, and Q20, Q30, and GC contents in the clean data were calculated. Clean reads were assembled by Trinity assembly software to obtain high-quality unigenes.

Functional Annotation of Unigenes and Screening of DEGs

The unigenes obtained by transcriptome sequencing were compared with KEGG, NR, Swiss-Prot, GO, COG/KOG, and Trembl databases by BLAST software to obtain the annotated information on unigenes. DESeq2 was used to analyze DEGs between groups. For the detection of DEGs, changes ≥2- or ≤1/2-fold and a false discovery rate (FDR) < 0.01 were set as the screening criteria. Go enrichment analysis and KEGG pathway enrichment analysis of DEGs were performed using GOseq and KOBAS with p ≤ 0.05, respectively.

qRT-PCR Primer Design and Amplification Reaction

The mRNAs from the blade, petiole, root, and stem of E. adenophorum were reversely transcribed into cDNA using a reverse transcription kit (FastKing RT Kit, Tiangen Biochemical Technology Co., Ltd.). Twenty genes were selected for real-time quantitative PCR (RTQ-PCR) using specific primers designed by primer premier 5.0 software (Table ). Specifically, RTQ-PCR was performed on an ABI Prism 7900-HT Sequence Detection System. The qRT-PCR reaction volumes were 10 μL, including 0.5 μL of cDNA, 0.3 μL of forward primer, 0.3 μL of reverse primer, 1 μL of 5×ROX, 5 μL of 2×Talent, and 2.9 μL of RNase-free H2O. The qPCR condition was set as 3 min at 95 °C, followed by 40 cycles at 95 °C for 5 s and at 60 °C for 5 s. Each sample was set up with three biological repetitions. The relative expression was determined after normalization against β-actin as an internal reference and was calculated using the 2 method.

Table 2

List of Primers for qRT-PCR of Functional Genes in E. adenophorum

gene name	amplified gene	sequence name	sequence (5′ −3′)	expected amplification length (bp)
phenylalanine amino lyase	PAL1	cluster-8876.107461	CCCCTCCGTGGAACCATTACC	193
			CCAGTGCCAGCCCTTCTTTAG
	PAL2	cluster-8876.107361	GATGAGGTGAAGAAGATGGTGG	201
			ATCCGTCCCTTTATTCATACTC
β-glucosidase	bglx1	cluster-8876.247895	CAAGCTTATCAGCCATGGAATA	223
			GAGACCCACATCATAACCACCTA
	bglx3	cluster-8876.131581	GCATTAGGTGGTTATGATGCG	219
			AGGTTCTAACCAATAGGCAAAG
	bglx6	cluster-8876.38078	GAAGATGAATATGGAGGATGGC	225
			ACCCGCATCATAACCACCTAA
4-coumarate-CoA ligase	4CL1	cluster-8876.57421	GTTGATTTGCGTGTTACCGCTG	214
			GTCGTATTTATCCACCACTTCTTC
cinnammate-4-hydroxylase	C4H1	cluster-8876.122929	TAAAGAGAAGAGGTTGAAGCTG	215
			ATTCGATAGACCATAGGGTTG
	C4H2	cluster-8876.138119	CAATCGAAACAACTCTATGGTCG	198
			GTGGGATAGCCATACGAAGACG
caffeic acid O-methyltransferase	COMT2	cluster-8876.146519	CCACATGTTATTGAAGATGCCA	222
			CGAGTCGGGTGCCTCAGGAAG
chalcone synthase	CHS	cluster-8876.138533	GCCTTCGGTCAAACGCTTCAT	202
			CCCGTCACCAAACAAAGCCTG
chalcone isomerase	CHI1	cluster-8876.130777	CGATTAGCAGCCGATGACAAG	208
			TCTCCACCACATTCCCGTTCT
naringenin 3-dioxygenase	F3H	cluster-8876.145359	AGCAATGGGCGGTCCAAGAAC	197
			AGATCGGTACTCATCTTCTTC
flavonoid 3′-monooxygenase	CYP75B	cluster-8876.138817	ACGTTGATCGGACTCAAGGAC	180
			TTGGGCTTGCTTTAGTAGACG
acetyl-CoA C-acetyltransferase	ACAT	cluster-8876.144261	ATTTGCTGTTGTGGCTCTTGC	209
			CAAATGCAGATGCACCTCCTC
hydroxymethylglutary l-CoA synthase	HMGS	cluster-8876.148013	ACTTGTATCCAGCAGGACACCC	207
			TACCGTTTCACTACCAACTTCC
hydroxymethylglutary l-CoA reductase	HMGCR1	cluster-8876.100846	TTGTGCGGACAAGAAACCTACC	200
			TAGCATGTGCATTAAAGCCTCC
acyclic sesquiterpene synthase	NES	cluster-8876.265980	GAAACTTTGCGTCAACCATTAC	203
			TATCTCCCTTGAACTTGCCATC
δ-cadinene synthase	TPS	cluster-8876.177344	GTGGTGGAAAGACTTGGGTGC	201
			ACGCCTTCACAGCGGTGGTAC
cytochrome P450	P4501	cluster-8876.151347	AGTCGTAATTGAATGGCTGATGC	231
			TGGTAGCCTCCAACCTCACAG
	P4502	cluster-8876.161609	ATATCTCCCGTTTGGTTCAGGG	215
			AATCGGTGTTGTAAGGATGTGG

Preparation and Extraction of Samples for Metabolomic Analysis

Metabolome analysis was performed by Wuhan Maiwei Biotechnology Co., Ltd.

The samples of E. adenophorum blades and stems were freeze-dried with a vacuum freeze-dryer (Scientz-100F). The freeze-dried samples were crushed using a mixer mill (MM 400, Retsch) with a zirconia bead at 30 Hz for 1.5 min. Then 100 mg of lyophilized powder was dissolved into 1.2 mL of 70% methanol solution. The mixture was swirled for 30 s every 30 min six times and placed in a refrigerator at 4 °C overnight. Centrifugation was performed at 12 000 rpm for 10 min, and then the extracts were filtered (SCAA-104, 0.22 μm pore size; ANPEL, Shanghai, China, http://www.anpel.com.cn/) before UPLC-MS/MS analysis.

UPLC Conditions

The sample extracts were analyzed using a UPLC–ESI–MS/MS system (UPLC, Shimadzu Nexera X2, http://www.shimadzu.com.cn/; MS, Applied Biosystems 4500 Q TRAP, www.appliedbiosystems.com.cn/). The analytical conditions were as follows. UPLC: column, Agilent SB-C18 (1.8 μm, 2.1 mm × 100 mm). The mobile phase consisted of solvent A (pure water with 0.1% formic acid) and solvent B (acetonitrile with 0.1% formic acid). Sample measurements were performed with a gradient program with the starting conditions of 95% A, 5% B. Within 9 min, a linear gradient of 5% A, 95% B was programmed, and the composition of 5% A, 95% B was held for 1 min. Subsequently, the composition was adjusted to 95% A, 5.0% B within 1.1 min and was held for 2.9 min. The flow velocity was set at 0.35 mL/min. The column oven was set to 40 °C. The injection volume was 4 μL. The effluent was alternatively connected to an ESI triple-quadrupole linear ion trap (QTRAP)-MS.

ESI-Q TRAP-MS/MS

LIT and triple quadrupole (QQQ) scans were acquired on a triple-quadrupole linear ion trap mass spectrometer (Q TRAP; AB4500 Q TRAP UPLC/MS/MS system) equipped with an ESI turbo ion-spray interface. The spectrometer was operated in positive and negative ion modes and controlled with Analyst 1.6.3 software (AB Sciex). The ESI source operation parameters were as follows: ion source, turbo spray; source temperature, 550 °C; ion spray voltage (IS), 5500 V (positive ion mode)/–4500 V (negative ion mode); ion source gas I (GSI), gas II(GSII), and curtain gas (CUR) set at 50, 60, and 25.0 psi, respectively; and a high collision-activated dissociation (CAD). Instrument tuning and mass calibration were performed with 10 and 100 μmol/L polypropylene glycol solutions in QQQ and LIT modes, respectively. QQQ scans were acquired as MRM experiments with the collision gas (nitrogen) set to medium. Optimized DP and CE were performed for individual MRM transitions. A specific set of MRM transitions were monitored for each period according to the metabolites eluted within this period.

Identification and Analysis of Metabolites

On the basis of the self-built database MWDB (Metware database) of Wuhan MetWare Biotechnology Co., Ltd.[44] and substance characterization based on secondary spectral information, a quantitative analysis of metabolites was performed with triple-quadrupole mass spectrometry in the multiple reaction monitoring (MRM) mode. Analyst 1.3 software was used to process the mass spectrometry data,[45] perform the integration and calibration of peaks, and export the peak area integration data for storage. Principal component analysis (PCA) and orthogonal partial least-squares discriminant analysis (OPLS-DA) were used for metabolites. The variable importance in a project (VIP) of the OPLS-DA model was obtained. Significantly regulated metabolites in groups were determined by VIP ≥ 1 and absolute Log2FC(fold change) ⩾ 1. VIP values were extracted from OPLS-DA results, and score plots and permutation plots were generated using R package MetaboAnalystR. The data was log-transformed (log 2) and mean centered before OPLS-DA. To avoid overfitting, a permutation test (200 permutations) was performed.

Metabolite Annotation

Metabolites identified by KEGG annotation and enrichment analysis were annotated using the KEGG Compound Database (http://www.kegg.jp/kegg/compound/). Annotated metabolites were then mapped to the KEGG Pathway Database (http://www.kegg.jp/kegg/pathway.html). Pathways with significantly regulated metabolites were then fed into MSEA (metabolite sets enrichment analysis), and their significance was determined by p values from a hypergeometric test.

Construction and Analysis of Protein Interaction Networks

The differential gene protein interaction network was analyzed using the String (https://cn.string-db.org/) protein interaction website. Then Cytoscape (version 3.8.0) was used to visualize the network and calculate the network-related topological properties by Network Analyze.[46]

Statistical Analysis

One-way ANOVA (Duncan’s test) was used to compare the differences among the four organ samples. p < 0.05 indicated a significant difference. Pearson correlation coefficient (PCC) was applied to analyze the correlation among a pair of metabolites, genes, and metabolites. PCC > 0.7, p< 0.05 and PCC > 0.7, p < 0.01 were considered to indicate statistically significant differences between the correlation of a pair of metabolites, metabolites, and genes, respectively. All analyses were performed using the SPSS 25.0 software.[47]