Metabolic and Transcriptional Profiling of Fraxinus chinensis var. rhynchophylla Unravels Possible Constitutive Resistance against Agrilus planipennis

: The emerald ash borer (EAB, Agrilus planipennis ), an ash-tree wood-boring beetle, has caused widespread mortality of ash. Asian ash, which coevolved with EAB, is considered more resistant than its North American and European congeners. Although some compounds and proteins related to resistance to EAB have been identified, the underlying ash resistance mechanism to EAB still needs further study. The Asian ash species, Fraxinus chinensis var. rhynchophylla , is highly resistant to EAB. In this study, metabolic and transcriptional profiling of the phloem of this species was investigated, and differentially expressed metabolites and genes were analyzed by comparing them with those of the susceptible F. pennsylvanica . Four hundred and twenty-eight metabolites were detected in both species, and several coumarins and lignans, which were exclusive to F. chinensis var. rhynchophylla , were identified. Compared with susceptible F. pennsylvanica, genes related to phenylpropanoid biosynthesis, ethylene (ET), and jasmonic acid (JA) biosynthesis and signaling in F. chinensis var. rhynchophylla were found to be up-regulated. It was hypothesized that coumarins, lignans, and ET and JA signaling might contribute to greater resistance to EAB in F. chinensis var. rhynchophylla . This study suggests candidate metabolites and genes for biomarker development in future ash-breeding programs.


Introduction
Ash trees are an economically important tree species and are widespread in temperate and sub-tropical regions, from North America to Eurasia. They are important timber species in many countries where wood is widely used for the manufacture of furniture and tools [1]. They are also commonly found as planted trees in urban and suburban landscapes and in forested settings for their aesthetic beauty [2]. The emerald ash borer (EAB), Agrilus planipennis Fairmaire (Coleoptera, Buprestidae), is an ash-tree wood-boring beetle. EAB larvae feed on the phloem and outer xylem of the tree, creating serpentine galleries in these tissues. Serpentine galleries disrupt the flow of water and nutrients throughout the tree and girdle the branches and trunk, eventually leading to canopy decline, followed was more resistant than Manchurian ash [19]. In Eastern Russia, F. chinensis var. rhynchophylla was also identified as an ash species with high resistance to EAB [20,21]. In the genus Fraxinus, a phylogenetic analysis showed that F. chinensis var. rhynchophylla and Manchurian ash belong to different sections [22]. Therefore, this species might contain some species-specific defensive compounds that are different from Manchurian ash. Joint analyses for metabolites and genes in F. chinensis var. rhynchophylla using metabolomics and transcriptomics could reveal potential resistance mechanisms to EAB.
In this study, the metabolomics and transcriptomics of the bark of F. chinensis var. rhynchophylla were analyzed, and differentially expressed metabolites and genes of this species were investigated by comparing them with those of the susceptible species F. pennsylvanica. Some potential defensive compounds and proteins that have been previously reported in Manchurian ash were confirmed and several coumarins and lignans that were exclusive to F. chinensis var. rhynchophylla were identified. Along with the expression levels of genes related to phenylpropanoid biosynthesis and genes involved in hormone biosynthesis and signaling, a possible mechanism underlying resistance to EAB in F. chinensis var. rhynchophylla was discussed.

Metabolic Profiles of F. chinensis var. rhynchophylla and F. pennsylvanica
To obtain a better understanding of differential metabolites between resistant F. chinensis var. rhynchophylla and susceptible F. pennsylvanica, a widely used metabolomics strategy was used for the analysis of bark metabolites. A total of 428 metabolites were identified in the two species. These metabolites included phenolic acids, amino acids and derivatives, saccharides and alcohols, organic acids, nucleotides and derivatives, free fatty acids, coumarins, alkaloids, LPCs (Lysophosphatidyl Cholines), lignans, glycerol esters, flavonoids, flavonols, LPEs (Lyso Phosphatidyl Ethanolamines), monoterpenoids, PCs (Phosphatidyl Cholines), triterpenes, isoflavones, sphingolipids, and tannins ( Figure 1A). A total of 390 metabolites were shared by the two species, while 17 and 21 metabolites were exclusive to F. chinensis var. rhynchophylla and F. pennsylvanica, respectively ( Figure  1B).

Differential Metabolites between F. chinensis var. rhynchophylla and F. pennsylvanica
Multivariate statistical analyses showed that the resistant and susceptible species were clearly divided into two separate groups in PCA analysis ( Figure S1). To identify the different metabolites in the samples of the two groups, supervised PLS-DA (partial least squares discriminant analysis) was performed. As shown in Figure 2A, a clear separation between the two groups was observed. Then, statistically significant differential metabolites between the two species were obtained according to three criteria, namely variable importance in projection (VIP) values from PLS-DA plot >1.0, p-value <0.05 at univariate analysis, and |Log2 fold change |≥1. In total, 183 differential metabolites were determined between F. chinensis var. rhynchophylla and F. pennsylvanica, mainly including simple phenolic acids, coumarins, lignans, flavonoids, and amino acids and derivatives. The top 20 differential metabolites with the largest VIP values are shown in Figure 2B. Fifty-nine differential metabolites were exclusive to or more abundant in F. chinensis var. rhynchophylla. A total of 17 metabolites were detected in the resistant but not in the susceptible species, including 2 coumarins (fraxetin diglucoside and fraxin), 1 lignan (forsythialan B), 1 phenylethanoid (calceolarioside A), 2 amino acid derivatives (L-cystathionine and acetyltryptophan), 1 organic acid (coumalic acid), and 10 phenolic acids (Table  1). For 42 metabolites, the amounts found were significantly greater in F. chinensis var. rhynchophylla than in F. pennsylvanica. Among these metabolites, the amounts of coumarins and lignans, especially fraxidin, esculin hydrate, fraxidin 8-O-glucoside, fraxetin, isofraxetin, syringaresinol-4′-O-β-D-monoglucoside, daphnetin, esculetin, and isofraxidin, were greater in the resistant than in susceptible species (Table 1). The concentration of fraxidin in F. chinensis var. rhynchophylla was 461-fold that in F. pennsylvanica. A total of 21 of metabolites were exclusive to F. pennsylvanica, including eight flavonoids, five phenolic acids, two amino acid derivatives, two lipids, one organic acid, and three other metabolites. For 103 metabolites, the amounts found were greater in F. pennsylvanica than in F. chinensis var. rhynchophylla, and they were mainly amino acids and derivatives, phenolic acids, saccharides and alcohols, flavonoids, and lipids. The pathway enrichment analysis suggested that the most significantly enriched pathway was associated with phenylpropanoid biosynthesis ( Figure S2).

Gene Expression Profiles of F. chinensis var. rhynchophylla and F. pennsylvanica
To investigate the differentially expressed genes (DEGs) between F. chinensis var. rhynchophylla and F. pennsylvanica, barks of the two species were sampled for RNA-seq on the Illumina platform, and their abundance was evaluated by converting read counts to fragments per kilobase per million (FPKM). Compared with gene expression levels in F. pennsylvanica, a total of 8387 DEGs were identified in F. chinensis var. rhynchophylla, and 5907 DEGs were annotated, including 2931 up-regulated and 2976 down-regulated genes. The absolute log2 transformed fold-change values for annotated DEGs ranged from 1 to 17. Gene ontology (GO) analysis was carried out to classify the functions of the DEGs between the resistant F. chinensis var. rhynchophylla and the susceptible F. pennsylvanica. For the biological process categories, the DEGs were mostly enriched in terms of protein ubiquitination and intracellular protein transport. Within the cellular component categories, cytosolic large ribosomal subunit and obsolete intracellular parts were the most enriched terms. In the molecular function categories, DEGs were highly enriched in categories for obsolete coenzyme-binding and monooxygenase activity ( Figure S3). To validate the expression profiles of the RNA-Seq data, six transcripts, including three up-regulated and three down-regulated DEGs in F. chinensis var. rhynchophylla, were selected for analysis using qRT-PCR. The results from the qRT-PCR analysis were in a high level of agreement with those obtained from the RNA-Seq analysis ( Figure S4) (Pearson correlation = 0.91).

DEGs Involved in Hormone Biosynthesis
Plant hormones play important roles in regulating plant defense against herbivore insects [24]. In this study, 17 DEGs were found to be related to the biosynthesis of common stress-related phytohormones, including JA (jasmonic acid), SA (salicylic acid), ET (ethylene), and ABA (abscisic acid). Five DEGs that were annotated as enzymes catalyzing the production of JA, i.e., LOX, AOC, and AOS, were expressed higher in F. chinensis var. rhynchophylla than in F. pennsylvanica. Seven DEGs were annotated as ethylene-forming enzymes (ACO and ACO-like protein) and were up-regulated in F. chinensis var. rhynchophylla. One DEG (ICS) and four DEGs (two NCEDs, ZEP, and β-CH), which control the biosynthesis of SA and ABA, respectively, were down-regulated in F. chinensis var. rhynchophylla. Two, three, and two DEGs were annotated as DELLA proteins GAI1-like, EIN3 proteins, and genes related to SA signaling (NPR1 and EDS), respectively, showing greater levels in the resistant species than in the susceptible one. Eight WRKYs were differentially expressed between the two species, five of which were higher in F. chinensis var. rhynchophylla than in F. pennsylvanica (Figure 4).

Discussion
F. chinensis var. rhynchophylla is a native Asian ash species and is widely planted in China. Several studies have identified that it is more resistant to EAB than Manchurian ash, which is generally used in the research on ash resistance mechanisms to EAB. Previous studies on F. chinensis var. rhynchophylla have mainly focused on the description of resistance to EAB. Although one study analyzed potentially defensive compounds in the bark, it only carried out a simple analysis for total tannin, total polyphenols, and total flavonoids [19]. The chemical composition and transcriptomic profile related to EAB resistance in F. chinensis var. rhynchophylla still remain unknown. In this study, a high throughput analysis of phloem tissues from this resistant species and susceptible green ash was carried out using comparative metabolomics and transcriptomics approaches for the purpose of identifying potential constitutive resistance-related metabolites and genes against EAB.
The metabolic analysis showed that phenolics, including phenolic acids, lignans and coumarins, and flavonoids were the largest group in the metabolite profiling of F. chinensis var. rhynchophylla, followed by amino acids and derivatives, saccharides and alcohols, organic acids, nucleotides and derivatives, free fatty acids, and coumarins. Compared with F. pennsylvanica, 59 metabolites were exclusive to or present in greater amounts in F. chinensis var. rhynchophylla, most of which were hydroxycoumarins, lignans, and phenolic acids. Nine hydroxycoumarins were reported, for the first time, in the comparative chemical analysis between EAB-resistant and -susceptible species. The metabolite fraxetin diglucoside, which was found to be exclusive to F. chinensis var. rhynchophylla in the present study, might be especially worthy of further study for the development of biomarkers for EAB-resistant species. A hydroxycoumarin, fraxin, was detected in the resistant F. chinensis var. rhynchophylla but not in susceptible F. pennsylvanica. This result was consistent with previous studies that found high levels of this metabolite in Manchurian ash, but its detection in green and white ash failed [11][12][13], indicating that this metabolite may be one of determinants of resistance to EAB in Asian ash species. Four hydroxycoumarins, fraxidin, fraxetin, esculetin, and esculin, which were considered exclusive to Manchurian ash in previous studies [12][13][14], were detected in the present study and were shown to be present in substantially higher levels in F. chinensis var. rhynchophylla than in F. pennsylvanica. Hydroxycoumarins have been suggested to play a role in plant resistance. The defensive role of some hydroxycoumarins has been confirmed by some experiments to consist in the inhibition of insect development or in anti-microbial activities [25,26]. Although the role of differentially accumulated hydroxycoumarins needs to be further confirmed, the present results broaden the knowledge of the chemical composition associated with ash resistance to EAB.
The lignan, forsythialan B, were found to be exclusive to F. chinensis var. rhynchophylla in this study. Another lignan, pinoresinol diglucoside, was detected in greater amounts in F. chinensis var. rhynchophylla than in F. pennsylvanica, in accordance with several previous studies on resistant Manchurian ash and susceptible North American ash [12,13]. Lignans have various effects on insects by exerting antifeedant activity or by inhibiting larval growth and molting, and these effects have been well documented for some lignans [27][28][29][30]. The present results confirm previous research and indicate a potential role of lignans in F. chinensis var. rhynchophylla resistance against EAB. In addition to hydroxycoumarins and lignans, several other phenylethanoids, such as calceolarioside A, that were found in Manchurian ash but were absent in susceptible green and white ashes [12,14,31] were detected only in the resistant species, further corroborating the contribution of phenylethanoids to Asian ash species resistance to EAB.
Pinoresinol, verbascoside, and oleuropein, which were considered Manchurian ashspecific by previous research [12,13], exhibited similar concentrations in the resistant and susceptible species in the present study. This suggests that these compounds might not contribute to the resistance to EAB in F. chinensis var. rhynchophylla, unless they take part in synergistic interactions with other metabolites or processes unique to the resistant species [14]. Several metabolites, such as tyrosol, syringing, luteolin and its glucosides, and apigenin and its glucosides, were present in low levels or absent in F. chinensis var. rhynchophylla, which is consistent with previous studies [13,14]. Although some of these metabolites, such as syringin, have shown a certain degree of antifeedant effects against insects, they may not confer resistance to EAB, given their high amounts in several susceptible species. It is possible that some of these metabolites, especially those found exclusively in susceptible ashes, act as feeding-stimulants or growth-promoting traits for EAB larvae.
Although many phenylethanoids were not assigned to the KEGG pathway due to the limited molecular interactions and reactions, the pathway of phenylethanoid biosynthesis was still significantly enriched in the pathway enrichment analysis of differential metabolites ( Figure S2). Transcriptomics analysis showed that 19 DEGs might be related to this pathway. The genes coding for 4CL, CAD, CAGT, CCoAOMT, and C4H, which catalyze the hydroxylation of trans-cinnamic acid to p-coumaric acid, the sequential enzymatic reactions in monolignol biosynthesis, the synthesis of coniferin from coniferyl alcohol, the methylation of caffeoyl-CoA into feruloyl-CoA and the hydroxylation of cinnamic acid into coumaric acid, respectively [32], were up-regulated in F. chinensis var. rhynchophylla. The enhanced expression of these genes might provide sufficient substrates for the production of various phenylpropanoids and their derivatives and might lead to the accumulation of phenylpropanoids in F. chinensis var. rhynchophylla.
Three and four DEGs were annotated as β-glucosidases and peroxidases, respectively, and showed higher levels in F. chinensis var. rhynchophylla. β-Glucosidases catalyze the hydrolysis of the β-glucosidic bond between a carbohydrate moiety and the basic coumarin core, resulting in bioactive coumarin aglycone forms, and are involved in the production of hydroxycoumarins [33,34]. Peroxidases catalyze the last step of lignin biosynthesis, i.e., oxidization of monolignols to produce p-hydroxyphenyl (H), guaiacyl (G), and syringyl (S) units [33]. Constitutively elevated levels of these enzymes may be responsible for the greater amounts of hydroxycoumarins and lignin and may contribute to the higher level of resistance in F. chinensis var. rhynchophylla.
Plants are sessile, so immediate recognition of the biotic stimuli in the early stages of infection is the key to determining plant resistance [35]. Phytohormones and their hormonal signaling networks connect perception and early signaling to broad transcriptional reorganization and defense induction, thereby playing a vital role in plant resistance to insects [24]. Although hormone signaling pathways are generally thought to play an important role in induced resistance, constitutive expression levels of genes related to plant hormone synthesis and signal transduction reflect the ability to sense and react to attacking insects [17,36]. In the present study, 17 DEGs were found to be involved in the biosynthesis of JA, ET, SA, and ABA. All DEGs that participate in JA and ET biosynthesis were up-regulated in F. chinensis var. rhynchophylla compared with F. pennsylvanic. Additionally, one DEG (FRAX09_000207320), which might be involved in the transport of the JA precursor 2-oxophytodienoic acid (OPDA) from the chloroplast [37], and EIN3, which mediates JA and ET signaling synergy [24], showed greater levels, suggesting their constitutive enhancement of the JA and ET pathways and the synergy between the two pathways in F. chinensis var. rhynchophylla.
Isochorismate synthase, which controls SA biosynthesis [38], and NPR1, which is involved in SA signaling [39], were expressed at lower levels, indicating antagonistic interactions between SA and JA pathways in this resistant species. Interestingly, the expression of four DEGs that were identified as enzymes controlling ABA biosynthesis was lower in F. chinensis var. rhynchophylla, as opposed to that of enzymes controlling JA biosynthesis. Generally, ABA has an antagonistic effect on JA in plant disease response but acts synergistically with JA in response to insects [40][41][42][43]. However, since ABA is also a vital signal for responses to desiccation, an effect accompanying insect attack in many cases, the role of this stress hormone in resistance to insects still remains to be determined [24]. Therefore, ABA and its interaction with JA in F. chinensis var. rhynchophylla resistance to EAB will be further studied.
Constitutive resistance appears to be particularly important in long-living trees due to the lack of short vegetation period or generation time typical of short-lived plants; thus, they frequently suffer from herbivore-caused damage [8,11,44]. Using metabolomics and transcriptomics technologies, qualitative and quantitative differences between F. chinensis var. rhynchophylla, which is highly resistant to EAB, and F. pennsylvanica, which is highly susceptible, were investigated. Several coumarins and lignans that were exclusive to F. chinensis var. rhynchophylla were detected, and nine hydroxycoumarins were reported, for the first time, in the comparative chemical analysis of EAB-resistant and -susceptible species. Combined with the higher levels of enzymes involved in phenylpropanoid biosynthesis, it seems that phenylpropanoids, especially hydroxycoumarins, play an important role in F. chinensis var. rhynchophylla constitutive resistance to EAB. Genes related to ET and JA biosynthesis and signaling were up-regulated in F. chinensis var. rhynchophylla when compared with F. pennsylvanica. These data suggest that the ET and JA signaling pathways contribute to F. chinensis var. rhynchophylla constitutive resistance, although they are generally thought to participate in plant-induced resistance. The present results are helpful for a fundamental understanding of ash defense mechanisms against EAB and provide some candidate metabolites and genes for biomarker development in ash-breeding programs.

Plant Materials
The EAB-resistant F. chinensis var. rhynchophylla and the susceptible green ash (F. pennsylvanica) were grown in Laoniujuan, which is located in Jingyuetan National Forest Park (43°52′ N, 125°21′ E) of Changchun City in Jilin Province. The resistance phenotypes of these two species have been confirmed in several previous studies [3,[18][19][20][21]. The ash trees had mean stem diameters of 3.12 ± 0.14 (SE) and 2.92 ± 0.12 cm at 15 cm above the soil when sampling. The trees sampled were healthy and did not have D-shaped exit holes or vertical splits on the trunk, which are indicators of EAB infestation [45]. Four individual trees from each species were selected for metabolomic analysis, and three of them were also used for transcriptomic analysis. In July 2019, phloem and outer bark samples were taken from the main stems at 1.5 m above the soil. A 10 cm long section of bark was cut from each tree and placed immediately in liquid nitrogen; then, it was stored at −80 °C until extraction and analysis.

Metabolite Extraction and UPLC Conditions
Eight bark samples, including four from F. chinensis var. rhynchophylla and four from F. pennsylvanica, were freeze-dried by vacuum freeze-dryer (Scientz-100F) and ground using a mixer mill (MM 400, Retsch). A total of 100 mg of powder was dissolved in 1.2 mL of 70% methanol solution. The solution was vortexed for 40 s every 30 min for 3 h and then stored at 4 °C overnight. The mixture was centrifuged at 12,000 rpm for 10 min. The supernatant was collected and filtered through a 0.22 micron membrane filter (SCAA-104; ANPEL, Shanghai, China). The extracts were analyzed using a UPLC-ESI-MS/MS system (UPLC, SHIMADZU Nexera X2, www.shimadzu.com.cn/ accessed on 15 January 2020; MS, Applied Biosystems 4500 Q TRAP, www.appliedbiosystems. com.cn/ 15 January 2020). The UPLC system was equipped with a Waters ACQUITY UPLC HSS T3 C18 column (1.8 µ m, 2.1 mm × 100 mm). The mobile phase included pure water with 0.04% acetic acid (solvent A) and acetonitrile with 0.04% acetic acid (solvent B). Separations for UPLC analysis were performed in gradient elution from 5% B. The proportion of phase B increased linearly to 95% within 10.00 min, stayed steady at 95% for 1 min, then decreased to 5% within 0.10 min and was held at 5% for 14 min. The injection volume was 4 μL, and the flow rate was set at 0.35 mL/min at a column temperature of 40 °C .

Qualitative and Quantitative Analyses of Metabolites
Qualitative analyses of the metabolites were performed by referencing primary and secondary mass spectrometry data against the self-built database MWDB V2.0 (Metware BiotechnologyCo., Ltd. Wuhan, China) and other available mass spectrometry databases such as MassBank (http://www.massbank.jp accessed on 14 February 2020), KNAPSAcK (http://kanaya.naist.jp/KNApSAcK accessed on 16 February 2020), HMDB (http://www.hmdb.ca accessed on 17 February 2020), and METLIN (http://metlin.scripps.edu/index.php 17 February 2020). Quantitative analyses were carried out by the multiple reaction monitoring mode (MRM) using triple quadrupole mass spectrometry according to the method described previously [46]. The relative contents of metabolites in each sample were represented with chromatographic peak area integrals.

Multivariate Data Analysis and Statistical Analysis
Raw data signals were processed by Analyst 1.6.3 (AB Sciex, Framingham, MA, USA) and normalized using log-transform. Multivariate data analysis, such as principal component analysis (PCA) and partial least squares regression discriminant analysis (PLS-DA), were carried out by R (http://www.r-project.org/ accessed on 25 February 2020). In the loadings plot of PLS-DA model, metabolites with variable importance for the project (VIP) >1.0 were regarded as potential differential metabolites between F. chinensis var. rhynchophylla and F. pennsylvanica. Then, nonparametric tests (Mann-Whitney U tests) were performed on abundance of the above metabolites. Using samples of F. pennsylvanica as the control group, the metabolites meeting the following three criteria were selected as differential metabolites between the two species, namely VIP >1.0, Mann-Whitney p-value < 0.05, and |Log2 fold change|≥1. Differential metabolites were annotated and classified using the Kyoto Encyclopaedia of Genes and Genomes (KEGG) database (http://www.kegg.jp/kegg/pathway.html accessed on 26 February 2020). Significantly enriched pathways in which the metabolites in a module were involved were compared with the background and defined by both a hypergeometric test and a threshold of p-value < 0.05.

RNA Sequencing and Differential Expression Analysis
Six individual trees, including three of resistant F. chinensis var. rhynchophylla (R1, R2 and R3) and three of susceptible F. pennsylvanica, were used for transcriptome analysis. The total RNA of each sample was extracted using the RNeasy kit (Qiagen) as indicated by the manufacturer. The yield and purity of RNA were measured using a nanodrop spectrophotometer. Quality and concentration of the isolated RNA were measured using the NanoDrop ND-1000 spectrophotometer, Qubit 2.0 fluorometer, and Agilent 2100 Bioanalyzer. Unstranded cDNA libraries were constructed using the NEBNext ® Ultra™ RNA Library Prep Kit for Illumina (NEB, Beverly, MA, USA) according to the manufacturer's protocol, and six cDNA libraries with an average insert size of 300 bp were obtained. Then, each constructed cDNA library was subjected to Illumina HiSeq™ 4000 (San Diego, CA, USA) sequencing. A total of 50,753,180, 52,618,848, 57,110,948, 60,186,388, 52,947,722, and 46,568,142 raw reads were obtained from the R1, R2, R3, S1, S2, and S3 libraries, respectively.
The adapter, primer sequences, and reads with poly-N were removed from the raw reads, and low-quality data (>50% bases with Q-value ≤ 20) were also filtered by calculating the Q20 and Q30 scores and sequence duplication level. Finally, a total of 49810040, 51728718, 56054150, 59228236, 52148240, and 45870910 clean reads were obtained in the R1, R2, R3, S1, S2, and S3 libraries, respectively, with Q30 values ≥ 93.56%. The genome Version 0.2 assembly and Version 0.2 annotation of F. pennsylvanica (accession PE_00248) [37] were downloaded from http://www.ashgenome.org (accessed on 5 June 2021) and used as reference. The alignment of the paired-end reads to the reference was performed using the HISAT2 software [47]. Transcript assembly, GTF document mergence, and transcript abundance estimation were performed using StringTie (Version 2.0) [48]. The read count information of each transcript was extracted by Python script (prepDE.py). Differential expression analyses were performed on the count data using the R package DESeq2. Gene expression was quantified in FPKM units. P-values were adjusted for multiple testing using the Benjamini-Hochberg linear step-up procedure and controlling the false discovery rate (FDR). Compared with gene expression levels in F. pennsylvanica, differentially expressed genes (DEGs) in F. chinensis var. rhynchophylla were screened based on the cutoff-adjusted p-values < 0.01 and |LFCShrink| ≥ 1. GO and KEGG pathway enrichment analyses were conducted to analyze the biological functions of DEGs, and clus-terProfiler package in R was used to conduct the functional enrichment analysis. Significantly enriched GO terms and pathways were defined by a hypergeometric test and a threshold of FDR < 0.05.

Quantitative Real-Time PCR (qRT-PCR)
To validate the accuracy of the RNA-seq data, qRT-PCR analyses were performed on a Roche LightCycler 480 (Roche Applied Science, Penzberg, Upper bavaria, Germany) using the SYBR Premix Ex Taq™ II Kit (TaKaRa, Dalian, China) following the manufacturer's instructions. Three up-regulated and three down-regulated DEGs were randomly selected for expression verification. The primer sequences used in qRT-PCR were listed in Table S1. PCR reactions were prepared in 20 µ L volumes containing 10 μL of SYBR Premix Ex TaqTM (TaKaRa, China), 2 μL of cDNA template, 0.8 µ L 10 µ M of forward primer, 0.8 µ L 10 µ M of reverse primer, and 6.4 μL of ddH2O. The PCR conditions were set according to the manufacturer's instructions (pre-incubation of 5 min at 95 °C , amplification of 40 cycles of 95 °C for 10 s, 60 °C for 10 s, and 72 °C for 10 s). The specificity of amplicons was verified by melting curve analysis (60 to 95 °C ) after 40 PCR cycles. All experiments were performed in three biological replicates and four technical replicates. The actin genes were used as internal controls.

Supplementary Materials:
The following are available online at www.mdpi.com/article/10.3390/f12101373/s1, Figure S1: Principal component analysis (PCA) of metabolites detected in F. chinensis var. rhynchophylla and F. pennsylvanica, Figure S2: KEGG enrichment analysis of differential metabolites between F. chinensis var. rhynchophylla and F. pennsylvanica, Figure S3: GO enrichment analysis of DEGs between F. chinensis var. rhynchophylla and F. pennsylvanica, Figure S4: The qRT-PCR validation of transcriptome data, Table S1: Primers used for qRT-PCR.
Author Contributions: L.Q., L.W., and X.W. initiated and designed the research; L.Q., R.W., J.L., T.Z., and Y.C. performed the experiments; L.W. analyzed the data; L.Q. and L.W. wrote the manuscript. All authors have read and agreed to the published version of the manuscript.