Molecular Defensive Mechanism of Echinacea purpurea (L.) Moench against PAH Contaminations

The understanding of the molecular defensive mechanism of Echinacea purpurea (L.) Moench against polycyclic aromatic hydrocarbon (PAH) contamination plays a key role in the further improvement of phytoremediation efficiency. Here, the responses of E. purpurea to a defined mixture of phenanthrene (PHE) and pyrene (PYR) at different concentrations or a natural mixture from an oilfield site with a history of several decades were studied based on transcriptomics sequencing and widely targeted metabolomics approaches. The results showed that upon 60-day PAH exposure, the growth of E. purpurea in terms of biomass (p < 0.01) and leaf area per plant (p < 0.05) was negatively correlated with total PAH concentration and significantly reduced at high PAH level. The majority of genes were switched on and metabolites were accumulated after exposure to PHE + PYR, but a larger set of genes (3964) or metabolites (208) showed a response to a natural PAH mixture in E. purpurea. The expression of genes involved in the pathways, such as chlorophyll cycle and degradation, circadian rhythm, jasmonic acid signaling, and starch and sucrose metabolism, was remarkably regulated, enhancing the ability of E. purpurea to adapt to PAH exposure. Tightly associated with transcriptional regulation, metabolites mainly including sugars and secondary metabolites, especially those produced via the phenylpropanoid pathway, such as coumarins, flavonoids, and their derivatives, were increased to fortify the adaptation of E. purpurea to PAH contamination. These results suggest that E. purpurea has a positive defense mechanism against PAHs, which opens new avenues for the research of phytoremediation mechanism and improvement of phytoremediation efficiency via a mechanism-based strategy.


Introduction
Polycyclic aromatic hydrocarbons (PAHs) are a family of environmentally persistent, widely distributed organic pollutants produced by the incomplete combustion of carbonbased fuels or originating from the release of petroleum into the environment [1,2]. These substances are readily adsorbed by soil particles and highly resistant to degradation by natural mechanisms, which cause their accumulation in the soil ecosystem and consequent entry into the human food chain [3]. Given their acute toxicity, carcinogenicity, mutagenicity, and teratogenicity, the remediation of PAH contamination is an ongoing endeavor. Among diverse options for remediation of PAHs, phytoremediation, a method that utilizes suitable plants to remove PAHs from contaminated environments, is a cost-effective and environmentally friendly strategy that has considerable advantages, such as low operating costs, non-interference with polluted sites and ecosystems, and the possibility to clean up different pollutions concomitantly [4][5][6].
The success of PAH phytoremediation is largely determined by the selection of suitable plant candidates because phytoremediation efficiency varies remarkably among plant

Growth and PAH Accumulation
The exposure of PHE mixed with PYR at a low level (25 mg kg −1 ) had no remarkable effect on the phenotype, plant height (PH), total leaf area per plant (LA), and shoot dry weight (SDW) of E. purpurea plants (Figure 1a-c). Notably, the exposure of PHE mixed with PYR at a high level (100 mg kg −1 ) resulted in small leaves, leaf senescence, and stunted growth in terms of PH (p < 0.01) and SDW (p < 0.001). Similarly, in PCM (natural mixture of PAH sandwiched with farm soil) treatment, when the plants were exposed to the natural mixture of PAHs derived from aged petroleum-contaminated soils, a significant decrease in growth was observed with respect to the PAH-free control (CK, p < 0.001). These results indicated that E. purpurea plants showed a clear resistance to PAH exposure at low levels but sensitivity to increased concentrations and complex PAH mixture. However, the highest , growth (c), and PAH accumulation (d) of E. purpurea upon PAH exposure. PH, plant height; LA, total leaf area; SDW, shoot dry weight; TPAH, total concentration of polycyclic aromatic hydrocarbon; CK, PAHs-free control; PP25, PP50, PP100, defined PHE and PYR at 25 mg kg −1 , 50 mg kg −1 , and 100 mg kg −1 for each one, respectively; PCM, natural PAH mixture from oilfield site. *, **, and *** represent significance at 0.05, 0.01, and 0.001 test levels, respectively. Different lowercase letters in d indicate significant differences between treatments at p < 0.05 (LSD).

Overview of RNA Sequencing, Functional Annotation, and Identification of DEGs
An overview of RNA-Seq data and a summary of the length distribution of the transcripts and unigenes are presented in Table S1 and S2, which suggest the high quality of RNA-seq data. Totally, 56,161 (61.54%) unigenes were successfully matched to at least one of the databases according to a comprehensive annotation, and 36,541 (94.49%) annotated unigenes matched five species, of which the most significant homology species was Helianthus annuus (L.) (79.22%) ( Figure S1). The total numbers of differentially expressed genes (DEGs) identified across all groups varied from 3176 with PP50/CK, 2606 , growth (c), and PAH accumulation (d) of E. purpurea upon PAH exposure. PH, plant height; LA, total leaf area; SDW, shoot dry weight; TPAH, total concentration of polycyclic aromatic hydrocarbon; CK, PAHs-free control; PP25, PP50, PP100, defined PHE and PYR at 25 mg kg −1 , 50 mg kg −1 , and 100 mg kg −1 for each one, respectively; PCM, natural PAH mixture from oilfield site. *, **, and *** represent significance at 0.05, 0.01, and 0.001 test levels, respectively. Different lowercase letters in d indicate significant differences between treatments at p < 0.05 (LSD).

Overview of RNA Sequencing, Functional Annotation, and Identification of DEGs
An overview of RNA-Seq data and a summary of the length distribution of the transcripts and unigenes are presented in Table S1 and S2, which suggest the high quality of RNA-seq data. Totally, 56,161 (61.54%) unigenes were successfully matched to at least one of the databases according to a comprehensive annotation, and 36,541 (94.49%) annotated unigenes matched five species, of which the most significant homology species was Helianthus annuus (L.) (79.22%) ( Figure S1). The total numbers of differentially expressed genes (DEGs) identified across all groups varied from 3176 with PP50/CK, 2606 with PP100/CK, to 3964 with PCM/CK. The numbers of up-regulated genes (1787 in PP50/CK and 1614 in PP100/CK) were higher than those of the down-regulated ones (1389 in PP50/CK and 992 in PP100/CK), which showed that the majority of genes were switched on after exposure to PHE and PYR in E. purpurea (Figure 2a). Compared with those in the PP50/CK and PP100/CK groups, a larger set of genes in the PCM/CK group showed a response. In addition, the upset plot analysis showed a relatively larger overlap between the three groups with 951 identical DEGs. Meanwhile, 1308 DEGs were shared by at least two but not all groups, and 4317 DEGs were not shared at all (Figure 2b). with PP100/CK, to 3964 with PCM/CK. The numbers of up-regulated genes (1787 in PP50/CK and 1614 in PP100/CK) were higher than those of the down-regulated ones (1389 in PP50/CK and 992 in PP100/CK), which showed that the majority of genes were switched on after exposure to PHE and PYR in E. purpurea (Figure 2a). Compared with those in the PP50/CK and PP100/CK groups, a larger set of genes in the PCM/CK group showed a response. In addition, the upset plot analysis showed a relatively larger overlap between the three groups with 951 identical DEGs. Meanwhile, 1308 DEGs were shared by at least two but not all groups, and 4317 DEGs were not shared at all (Figure 2b).

Overview of Metabolomic Data and Identification of DAMs
The good repeatability and instrumental stability of metabolite detection were confirmed by the PCA score plot, which showed that the three quality control samples were overlaid ( Figure 3a). Among identified metabolites, 78, 118, and 208 differentially accumulated metabolites (DAMs) were observed in the three comparison groups, and of these, amino acids, carbohydrates, and their derivatives were the first predominant chemical classes (Figure 3b,c). Furthermore, considering the secondary metabolites detected here, PAH treatments affected metabolites predominantly belonging to the class of phenolic acids, flavonoids, and their derivatives. These results demonstrated that PAH exposure strongly affects the metabolite profiles of E. purpurea, which was associated with PAH exposure levels.

Overview of Metabolomic Data and Identification of DAMs
The good repeatability and instrumental stability of metabolite detection were confirmed by the PCA score plot, which showed that the three quality control samples were overlaid ( Figure 3a). Among identified metabolites, 78, 118, and 208 differentially accumulated metabolites (DAMs) were observed in the three comparison groups, and of these, amino acids, carbohydrates, and their derivatives were the first predominant chemical classes (Figure 3b,c). Furthermore, considering the secondary metabolites detected here, PAH treatments affected metabolites predominantly belonging to the class of phenolic acids, flavonoids, and their derivatives. These results demonstrated that PAH exposure strongly affects the metabolite profiles of E. purpurea, which was associated with PAH exposure levels.

GO Classification
According to GO enrichment analysis, a total of 8 enriched terms in biological process (BP) and 5 enriched terms in molecular function (MF), including 4 sub-categories, namely, "chlorophyll catabolic process", "pigment catabolic process", "porphyrin-containing compound catabolic process", and "tetrapyrrole catabolic process", were shared between PP100/CK and PCM/CK groups ( Figure S2 and Table S3). Totally, 17 genes encoding seven enzymes involved in the chlorophyll metabolic pathway were annotated (Figure 4a). Among these, 1, 3, and 9 genes encoding enzymes related to the chlorophyll cycle, namely, chlorophyllide a oxygenase (CAO), Chl b reductase (non-yellow coloring, NYC1/NOL), and chlorophyllase (CLH), were identified, and most of them were up-regulated under PAH exposure. In addition, the amount of Chla increased, and that of Chlb decreased in PCM, which indicated the enhanced interconversion between Chla and Chlb under complex PAH exposure (Figure 4b). The up-regulated expression of two genes associated with chlorophyll degradation, including Cluster-16663.37827, which encodes pheophorbide a oxidase (PAO) and Cluster-16663.69223, which encodes red Chl catabolite reductase (RCCR), was associated with low Chl level in the PCM leaves ( Figure 4b).

GO Classification
According to GO enrichment analysis, a total of 8 enriched terms in biological process (BP) and 5 enriched terms in molecular function (MF), including 4 sub-categories, namely, "chlorophyll catabolic process", "pigment catabolic process", "porphyrin-containing compound catabolic process", and "tetrapyrrole catabolic process", were shared between PP100/CK and PCM/CK groups ( Figure S2 and Table S3). Totally, 17 genes encoding seven enzymes involved in the chlorophyll metabolic pathway were annotated ( Figure 4a). Among these, 1, 3, and 9 genes encoding enzymes related to the chlorophyll cycle, namely, chlorophyllide a oxygenase (CAO), Chl b reductase (non-yellow coloring, NYC1/NOL), and chlorophyllase (CLH), were identified, and most of them were up-regulated under PAH exposure. In addition, the amount of Chla increased, and that of Chlb decreased in PCM, which indicated the enhanced interconversion between Chla and Chlb under complex PAH exposure ( Figure 4b). The up-regulated expression of two genes associated with chlorophyll degradation, including Cluster-16663.37827, which encodes pheophorbide a oxidase (PAO) and Cluster-16663.69223, which encodes red Chl catabolite reductase (RCCR), was associated with low Chl level in the PCM leaves ( Figure 4b).

KEGG Pathway Analysis
The KEGG enrichment analysis showed that 40, 58, and 176 DEGs were significantly enriched in 3, 3, and 9 KEGG pathways, respectively, in the PP50/CK, PP100/CK, and PCM/CK groups (Table S4). The pathway of "circadian rhythm of plant", which was co-enriched between three groups, and the pathways of "plant hormone signal transduction" and "starch and sucrose metabolism", which showed the top two largest numbers of DEGs, were selected to perform the subsequent analysis ( Figure 5). We identified

KEGG Pathway Analysis
The KEGG enrichment analysis showed that 40, 58, and 176 DEGs were significantly enriched in 3, 3, and 9 KEGG pathways, respectively, in the PP50/CK, PP100/CK, and PCM/CK groups (Table S4). The pathway of "circadian rhythm of plant", which was coenriched between three groups, and the pathways of "plant hormone signal transduction" and "starch and sucrose metabolism", which showed the top two largest numbers of DEGs, were selected to perform the subsequent analysis ( Figure 5). We identified and annotated nine key components (Figure 6a), including circadian clock-associated 1 (CCA1), timing of cab expression 1 (TOC1), late elongated hypocotyl (LHY), pseudo-response regulator 5/7 (PRR5/7), early flowering 3 (ELE3), lux arrhythmo (LUX), phytochrome interacting factor 4 (PIF4), and gigantea (GI), which are involved in the "circadian rhythm of plant". The genes encoding CCA1, LHY, and PIF4 were down-regulated, and other genes showed up-regulated expression patterns in response to PAH exposure (Figure 6a,b). The genes in jasmonic acid (JA) and salicylic acid signaling pathways, which are in crosstalk with the circadian clock, were also identified, and most of them in JA were significantly downregulated in three groups, including coronatine insensitive 1 (COI1)-encoding gene, JA ZIM Domain proteins (JAZ)-encoding gene and transcription factor (MYC2). Meanwhile, in SA, genes encoding non-expressor of pathogenesis-related gene 1 (NPR1) were upregulated. More genes were coregulated in PCM than PP50 and PP100. We also identified three metabolites, namely, JA, JA-isoleucine conjugate (JA-Ile), and methyl JA (MeJA), which showed an increased trend from 1.2-to 2.2-fold (Log2FC), especially under PHE and PYR exposure at 100 mg kg −1 , which means an enhancement of JA biosynthesis in E. purpurea leaves after PAH exposure (Figure 6c). In concert with this, most genes in the JA synthesis pathway were up-regulated in response to PYR and PHE exposure. Together with JA biosynthetic genes, genes in JA-catabolic pathways, namely, JA-Ile oxidation and JA-Ile deconjugation, were coregulated differently to maintain JA homeostasis in response to PAH exposure.  Furthermore, together with the pathway of "starch and sucrose metabolism", which was enriched in PP100 and PCM groups, two pathways, namely, "galactose metabolism" and "fructose and mannose metabolism", were enriched in the PCM/CK ( Figure S3). In response to PAH exposure, a majority of genes (27 out of 32 in PP50, 39 out of 48 in PP100, and 37 out of 57 in PCM) encoding critical enzymes involved in the three pathways were up-regulated; these enzymes included alpha/beta-glucosidase (malZ) and β-fructofuranosidase (INV) ( Figure S3). Consistent with the gene expression, 7, 7, and 12 DAMs classified as carbohydrates and derivatives, including galactose 1-phosphate, D-glucose 1-phosphate, galactinol, melibiose, sucrose, cellobiose, trehalose, trehalose 6phosphate, maltose, isomaltose, D-fructose 6-phosphate, and D-mannose 6-phosphate, were mapped to three of three pathways in PP50, PP100, and PCM, respectively. In addition, the abundance of these metabolites increased more under PP100 than PP50. In contrast to more metabolites accumulated in PP50 and PP100 in a concentration-dependent manner, more metabolites were depleted in PCM.

Association Analysis of Metabolome and Transcriptome
Association analysis of KEGG based on metabolome and transcriptome indicated that "starch and sucrose metabolism" was mainly enriched in all three treatments ( Figure  S4a-c), which was consistent with the results obtained in Figure 4. To explore the links

Association Analysis of Metabolome and Transcriptome
Association analysis of KEGG based on metabolome and transcriptome indicated that "starch and sucrose metabolism" was mainly enriched in all three treatments ( Figure S4a-c), which was consistent with the results obtained in Figure 4. To explore the links among DEGs, DAMs, and growth properties, the Mantel test was used ( Figure S4d). The results indicated that both TPAH (Mantel's r = 0.546, 0.610, p < 0.05) and SDW (Mantel's r = 0.616, 0.710, p < 0.05) showed a stronger correlation with DEGs and DAMs, respectively. Furthermore, Chla was significantly correlated with DEGs and DAMs; LA was correlated with DAMs. These results support that metabolite accumulations and transcription regulations contribute to the growth response of E. purpurea to PAH exposure.

Responses of Growth and Chlorophyll Metabolism in E. purpurea Leaves to PAH Exposure
The advantages of phytoremediation have been mentioned previously; however, the most considerable disadvantages are the toxicity of PAH for the plants over a certain concentration threshold, which can limit the phytoremediation efficiency [31]. Tolerance of plants to contaminants is an essential criterion for successful phytoremediation [32,33]. Our results showed that E. purpurea exhibited a tolerant phenotype under PAH treatments (as a mixture of PHE and PYR) within the tested dosages, especially at lower dosages of 50 mg kg −1 and 100 mg kg −1 (Figure 1). Conversely, E. purpurea was less resistant to the natural PAH mixture in PCM treatment and exhibited reduced growth and decreased physiological parameters (Figures 1 and 4). Supporting this result, PCC analysis among leaf parameters also showed that TPAH in the leaf was negatively correlated with SDW, LA, and Chla ( Figure S5). In this work, members of genes involved in the chlorophyll cycle, such as CAO, NYC1/NOL, and CLH, displayed similar expression patterns and were significantly up-regulated and then promoted the cycle of chlorophyll in E. purpurea leaves when PAH-exposed ( Figure 4). PaO is a key enzyme for chlorophyll breakdown, and it catalyzes open the porphyrin ring of pheophorbide to form red chlorophyll metabolite; then it is converted into a primary fluorescent chlorophyll metabolite by another key enzyme, RCCR [34]. Here, the high expression of genes encoding PaO and RCCR showed similar trends to those observed in Chl content analysis ( Figure 4). Hence, the low chlorophyll content in PCM could be explained partly by the promotion of the chlorophyll cycle and degradation under exposure to the natural PAH mixture.

Effects on Circadian Rhythm and Hormone Signal Transduction in E. purpurea Leaves
In A. thaliana, the core oscillator of the circadian clock mainly consists of CCA1, LHY, TOC1, PRR9/7/5, ELF3, ELF4, LUX, and GI, which comprises three interlocking negative feedback loops (central/morning/evening loop) [27,28,35,36]. Given that the clock transcriptional network in A. thaliana is partly conserved among angiosperms [37], the expression of associated genes in E. purpurea was characterized ( Figure 6). After PAH exposure, circadian genes, including CCA1, LHY, and PIF, were downregulated, whereas other genes, such as ELF3, LUX, GI, PPR5, and PPR7, showed up-regulated expression patterns. To our best knowledge, this work provides the first observation of circadian regulation in plants under PAH exposure, and the results may help in elucidating the function of circadian clock regulation in the defense responses of plants to PAH contamination. As another component of the circadian clock molecular system, circadian output pathways, which include a series of downstream events, can contribute to the ability of plants to adapt to environmental alterations [38,39]. Emerging evidence shows that as key components of the core oscillator, CCA1, PRR7, and ELF3 are associated with plant growth and leaf senescence [36,40]. Therefore, we concluded that during PAH phytoremediation, the reduced plant growth and enhanced leaf senescence in E. purpurea leaves are modulated by the circadian clock (Figures 1, 4 and 6). In addition, the circadian clock also functions as a master regulator in phytohormone biosynthesis and signaling [40]. In this concern, JA biosynthesis, homeostasis, and signaling are regulated by the circadian clock, whereas, in a feedback manner, JA also participates in the regulation of circadian speed [41,42]. Furthermore, JA can induce secondary compounds biosynthesis and participate in the plants' adaptation to environmental stresses [29,[43][44][45]. Here, the expressions of genes associated with the pathways of JA and the circadian network were the most affected, and the biosynthesis of hormones, such as JA, JA-Ile, and Me-JA, was remarkably enhanced (Figures 4 and 5), which contributed to the ability of E. purpurea to adapt to PAH stress.

Effects on Sugar Metabolism and Secondary Metabolism in E. purpurea Leaves
In plants experiencing different environmental stress conditions, saccharide accumulation in the cytoplasm is an important mechanism for maintaining osmotic balance [46,47]. Moreover, by providing more energy resources to the plant, carbohydrates can improve the synthesis of particular anti-stress compounds, such as secondary metabolites, for survival under adverse conditions [48]. An earlier study confirmed that galactose metabolism is the most remarkable metabolic pathway disrupted by BaP and PYR exposure in Zea mays (L.) leaves [25]. Equally, we expanded on those findings and confirmed that most DEGs, i.e., malZ and INV, in galactose metabolism, were significantly up-regulated in response to PAH exposure, either as a defined or natural mixture ( Figure S3a). Sucrose, which is the basic stored form of saccharides, can function both in plant growth and signal transduction as well as stress mitigation [47]. In the present study ( Figure S3), the plants subjected to PAH exposure presented a significant accumulation of sucrose, together with galactinol and melibiose, which provided the E. purpurea plant with high energy levels and protected it against PAH stress. This was confirmed by association analysis of KEGG pathways ( Figure S4), which also indicated that the "starch and sucrose metabolism" pathway was mainly enriched based on metabolomic and transcriptomic data.
The Mantel test revealed that TPAH was strongly correlated to DEGs and DAMs, which indicated that PAH concentration is tightly associated with gene regulation and metabolite accumulation in E. purpurea leaves and may explain the metabolic responses caused by PAH exposure ( Figure S4). In this concern, the accumulation of several phenolic metabolites, such as phenolic acids, flavonoids, and coumarins, which were produced via the phenylpropanoid pathway, was observed (Figures 3 and 7). Furthermore, the regulation of phenylpropanoid pathway genes (COMT, CCoAOMT, BGLU, and CHS), which were associated with the metabolism of phenolic acids and synthesis of coumarins and flavonoids, was found ( Figure 7). Phenolic secondary metabolites have different functions in plant growth and stress response [49][50][51]. Spinedi et al. (2021) [6] and   [21] demonstrated that the up-regulation of genes in flavonoid biosynthesis and the accumulation of flavonoid compounds could contribute to the antioxidant activity against anthracene toxicity in Marchantia Polymorpha (L.) or is the key to enhancing tolerance to PHE stress in A. thaliana, respectively. Among numerous structural modulations, glycosylation can alter the homeostasis of polyphenols and glycosidic phenols with antioxidant activity and function in plant defense response [52]. In this study, glycosidic flavonoids, such as syringetin 3/5/7-O-hexoside, accumulated higher contents in E. purpurea leaves of PP50 and PP100 compared with the CK (Figure 7). Although our understanding of the roles of secondary metabolites and their derivatives in PAH phytoremediation remains limited, the accumulation of phenolic secondary metabolites and the regulation of phenylpropanoid and flavonoid biosynthesis pathways partially function in the biochemical adaptability of E. purpurea after PAH exposure. However, the ecological aspect of invasive exotic plant species [53] needs to be considered in future research in order to achieve sustainable phytoremediation of PAHs using E. purpurea.

Soil Preparation, Experimental Design, and Plant Cultivation
Two soil types were used in this study. One (farm soil) was clay loam soil with no history of contamination and contained only background levels of PAH or no detectable PAH. This soil was collected from an experimental farm near Shenyang City, Liaoning Province, China, at a depth of approximately 20 cm. The other (oilfield soil) was aged PAH-contaminated soil obtained from Daqing Oilfield, which has been operating for more than 30 years in Daqing City, Heilongjiang Province, China. After collection, the soils were air-dried and homogenized by sieving to remove stones and residual plant materials. Table S5 summarizes the basic characteristics of the soil samples. Five treatments based on a component-based approach and a natural-mixture approach were set up as follows: farm soil without PAHs (CK); farm soils spiked with PHE and PYR mixtures at 25 (PP25), 50 (PP50), and 100 mg kg −1 (PP100) for each one; oilfield soil containing natural mixture of PAH sandwiched with farm soil (PCM). The selected PAH concentration was the range described in previous studies [25,54]. In the PP25, PP50, and PP100 treatments, 10% total quantity of farm soils was spiked with highly-pure PHE and PYR in acetone. Afterward, the soils were left under the fume hood to evaporate any traces of acetone. Then, the spiked soils were thoroughly mixed with uncontaminated farm soils to obtain initial total PAH concentrations of 50, 100, and 200 mg kg −1 . The sum of seven PAHs concentrations in the oilfield soil was about 150 mg kg −1 , and chrysene (CHR, 65.3%), PHE (20.6%), and fluoranthene (5.7%) were the major species (analysis method is shown in Section 2.3). The aged oilfield soil was sandwiched with uncontaminated farm soil in accordance with the method of Dai et al. (2020) [55] for the PCM treatment. The contaminated soil in the PCM treatment was diluted to a PAH concentration of about 50 mg kg −1 through the addition of uncontaminated farm soil.
The experimental soil was placed in 1.5 L pots with an underpan to avoid chemical and soil loss when watering. Then, 2-month-old E. purpurea seedlings with similar appearances, plant height (PH), and leaf quantity were selected and transferred to these pots for the phytoremediation experiment stretched over a period of 2 months (60 d). After transplanting, the pots with twelve replicates of each treatment were cultivated in a growth chamber (a 16/8 h light/dark cycle and temperatures of 18 • C-27 • C.) following the method of Heidari et al. (2018) [11]. Each pot was irrigated regularly to avoid excessive drainage. After 2 months of phytoremediation, E. purpurea plants were collected. Then, the plants were divided into two parallel groups for growth measurement (twelve biological replications), PAH analysis, metabolome analysis, and transcriptome sequencing (three biological replications). After collection, the leaf samples were immediately frozen and then stored at −80 • C until subsequent metabolite identification, transcriptome sequencing, and quantitative real-time polymerase chain reaction (qRT-PCR) validation.

Growth Measurement
The plant height (PH), total leaf area per plant (LA), Chlorophyll a (Chla), Chlb, and total Chl contents in leaves were measured in accordance with the method of Sun et al. (2018) [56]. Afterward, the shoot dry weight (SDW) was determined via oven-drying at 80 • C until constant weight.

PAH Extraction and Analysis
The PAH concentration in leaf samples was determined in accordance with the method described by Wang et al. (2012) and  with modifications [32,57]. Briefly, 5 g each of the leaf sample and quartz sand (blank sample) were placed in a centrifuge tube with 10 mL dichloromethane and 20 µL 200 mg mL −1 deuterated triphenyl solution as the recovery standard. Subsequently, the tube was ultrasonicated thrice for 1 h each at 35 • C. The extracts were combined, concentrated with gentle nitrogen blowing to 1 mL, and then passed through a preconditioned magnesium silicate SPE column (1000 mg 6 mL −1 ) with 11 mL 1:1 (v/v) elution of dichloromethane and n-hexane. After filtration, the extracts were re-concentrated with gentle nitrogen blowing, transferred to a sample vial containing 50 µL 20 µg mL −1 intermediate standard solutions, and diluted with nhexane to 1 mL. The internal standard was the mixture of deuterated PAHs, including Acenaphthylene d10, PHE d10, CHR d12, and perylene d12 (J&K Scientific Inc.; Beijing, China). A gas chromatograph mass spectrometer (HP6890-5975B, Agilent Co.; Santa Clara, CA, USA) equipped with an HP-5MS column was used for the quantification of target PAH concentrations in the extracts.

Transcriptomic Sequencing and qRT-PCR
RNA isolation was carried out following a method described in our previous work [58], and RNA-seq was carried out using an Illumina HiSeq 6000 system (Illumina, San Diego, CA, USA) at Novogene Biotech (Beijing, China). The sequenced reads (NCBI SRA accession number: PRJNA905488) were filtered by the removal of low-quality reads for subsequent analysis. Transcriptome assembly sequences were annotated in seven databases, including NR, NT, KO, Swiss-Prot, Pfam, KOG, and GO. Differentially expressed genes (DEGs) were identified at a cutoff of |Log2FC)| > 1 (false discovery rate and adjusted p < 0.05). Ten DEGs with different putative functions were selected to confirm the RNA-Seq results by performing qRT-PCR following our previous work [58]. The employed primers are summarized in Table S6, and the correlation between RNA-seq and qRT-PCR data is presented in Figure S6.

Metabolomic Analysis
The widely targeted metabolomic analysis was performed using the high-performance liquid chromatography-tandem mass spectrometry with an ExionLC™ AD system and an Xselect HSS T3 column (150 × 2.1 mm 2 , 2.5 µm) coupled to a QTRAP ® 6500+ mass spectrometer (SCIEX) at Novogene Biotech (Beijing, China) following a method described in our previous work [57]. Compound Discoverer software was used to search the online libraries of mzCloud, mzVault, and ChemSpider to identify the metabolites. After that, the differentially accumulated metabolites (DAMs) were screened considering variable importance in the projection > 1, p < 0.05, and |Log2FC)| > 1 and then subjected to KEGG enrichment analysis.

Data and statistical analysis
One-way analysis of variance was carried out to test the effect of PAH treatments on growth parameters using SPSS 25.0 (IBM Co.; Armonk, NY, USA). A volcano plot was performed to filter DEGs based on the −Log10p values of genes. Principal component analysis (PCA) was performed to visualize the differences among samples using the Vegan package in R (v4.0.3, https://www.r-project.org/ (accessed on 9 March 2023)). Heatmaps were generated to identify the expression profiles of DEGs or alterations in DAMs under treatments using the pheatmap package in R. To integrate transcriptomic and metabolomic data, the Pearson correlation coefficients (PCCs) analysis among growth parameters and partial Mantel test between growth parameters, DEGs and DAMs were performed using the Vegan package in R (9999 permutations), respectively.

Conclusions
Integrating transcriptomic data with metabolomic data, we concluded that PAH exposure induced the regulation of multiple pathways and processes, including chlorophyll cycle and degradation, circadian rhythm, JA biosynthesis and signaling, and starch and sucrose metabolism. The accumulation of metabolites, especially secondary metabolites produced via the phenylpropanoid pathway (coumarins and flavonoids), is also responsible for the adaptation of E. purpurea to PAH contamination. The identified specific gene and metabolite targets are associated with primary processes set up to sense PAH as well as with molecular and metabolic mechanisms to tolerate PAH exposure in E. purpurea. Collectively, this novel information opens a new avenue for the mechanism research of PAH phytoremediation, which will be valuable in the improvement of phytoremediation efficiency with E. purpurea via a mechanism-based strategy.