Transcriptomic and Metabolomic Analyses to Unravel Responses of Populus Simonii × P. Nigra to attack by Asian Gypsy moth Lymantria Dispar (Lepidoptera: Lymantriidae)

Background: Poplar is frequently attacked by herbivorous insects, including the Asian gypsy moth, Lymantria dispar. Here, we combined metabolomic and transcriptomic analysis to identify key genes and metabolites involved into the molecular mechanism of defensive enhancement against L. dispar herbivory on poplar. Results: The 3666 differentially expressed genes (DEGs, 1,799 up-regulated and 1,867 down-regulated) and 1,171 DEGs (695 up-regulated and 476 down-regulated) were specic in L. dispar herbivory and mechanical wounding, respectively. Moreover, the 9,108 and 7,656 ions were detected while 636 and 531 different ions were obtained using positive (pos) mode and negative (neg) mode, respectively. Among these ions, the 33 and 7 different ions were specic in L. dispar herbivory and mechanical wounding in pos mode while 46 and 4 different ions in L. dispar herbivory and mechanical wounding in neg mode, respectively. The 3,666 specic DEGs in L. dispar herbivory group were classied into phenylpropanoid and avonoid secondary metabolism pathways by comprehensive networks between transcriptomes and metabolomes. Conclusions: The current ndings greatly improve our understanding of the induced defensive response in will contribute

abiotic invasions [8]. Asian gypsy moth, Lymantria dispar, is a devastating chewing pest in forestry [9]. Among integrated pest managements, the potential mechanisms of plant-insect interactions and the identi cation of the new strategy to improve plant resistance have been considered as a perspective control strategy.
Nowadays, high-throughput omics techniques, including genomic, transcriptomic, proteomic, metabolomic, phenomic, ionomic, have been extensively employed into plant-insect interactions [10]. To further understanding speci c defense of P. simonii × P. nigra against L. dispar herbviory, we conducted a combined transcriptomic and metabolomics analysis to excavate differences in transcript and metabolic levels of poplar leaves under L. dispar herbivory and mechanical wounding treatments. The resulting data will promote to better understand the underlying molecular defensive mechanisms of poplar and explore more new biological control strategy against forest insect pests.

Transcriptome assembly
To investigate the response of poplar to L. dispar herbivory and mechanical wounding at the transcriptome level, the output of high-quality clean reads was obtained from paired-end BGISEQ-500 Table S1). The clean reads generated from 68.94 to 69.78 Mb among all treatments. Among them, 98.11 ~ 98.36% of bases had a Q ≥ 20, whereas 90.91 ~ 91.86% of bases a Q ≥ 30. The total length, average length, N50 and GC content of unigenes were 109.9 Mb, 1,567 bp, 2,222 bp and 41.68%, respectively. All unigenes was annotated into six functional databases, and 62,326  By the comparative transcriptome analysis, a total of 12,152 DEGs (6,733 up-regulated and 5,419 downregulated genes) were identi ed between the L. dispar herbivory group and the control group, and 6,440 DEGs (3,998 up-regulated and 2,442 down-regulated genes) were found between the mechanical wounding group and the control group. Moreover, the 3,790 up-regulated genes and 4,145 down-regulated genes were identi ed between L. dispar herbivory and mechanical wounding group (Fig. 2). The venn diagram indicated 3,666 DEGs (1,799 up-regulated genes and 1,867 down-regulated genes) and 1,171 (695 up-regulated genes and 476 down-regulated genes) were speci c in L. dispar herbivory and mechanical wounding, respectively (Fig. 3).

RNA-Seq (Supplementary
According to the KEGG annotations, the 3,666 DEGs in L. dispar herbivory group and 1,171 DEGs in mechanical wounding group were involved into ve categories of "metabolism", "genetic information processing", "environmental information processing", "organismal systems" and "cellular processes". In metabolism category, more systemic DEGs in L. dispar herbivory group and mechanical wounding group were classi ed into "global and overview maps" (744 genes and 198 genes), "carbohydrate metabolism" (312 genes and 81 genes) and "energy metabolism" (199 genes and 19 genes), respectively. In genetic information processing, the DEGs in L. dispar herbivory group and mechanical wounding group were group into subcategory "translation" (169 genes and 41 genes) and "folding, sorting and degradation" (124 genes and 46 genes), respectively. In "environmental information processing" and "organismal systems", the DEGs in L. dispar herbivory group and mechanical wounding group were involved into "signal transduction" (175 genes and 57 genes) and "environmental adaptation" (105 genes and 45 genes), respectively (Fig. 4). These results indicated that herbivory triggered a large scale of transcription recon guration, and the induced molecular mechanism of the herbivory response was different from those of the mechanical wounding responses.
Non-targeted metabolic pro les of poplar leaves associated with L. dispar herbivory To overview the metabolic shifts over the process in response to herbivory and wounding, non-targeted metabolic pro ling was performed using liquid chromatography-mass spectrometry (LC-MS). Nine samples (three control, three L. dispar herbivory and three mechanical wounding) were analyzed after 24 h treatment. The differential metabolites by fold change ≥ 1.2 or ≤ 0.8333 and p-value ≤ 0.05 were obtained in pos mode and neg mode, respectively ( and mechanical wounding conditions compared with the control condition by the PC2 (Fig. 6). Under herbivore condition, the PC2 could separate the difference between L. dispar herbivory and control groups, which explains 28.51% of the total metabolome variation among samples in neg mode, and 26.91% of the total metabolome variation in pos mode. The dominated components for the PC2 were pantothenol, 4'-cinnamoylmussatioside and others in neg mode, and vinblastine sulfate, coenzyme B and others in pos mode, correspondingly, ubiquinone-6 (37.62-fold) and cimicifugoside (53.89-fold) were obviously up-accumulated in neg mode and pos mode, respectively ( Fig. 6A and Fig. 6B).
Under mechanical wounding condition, the PC2 could separate the difference between mechanical wounding and control groups, which explains 19.85% of the total metabolome variation of among samples in neg mode, and 20.54% of the total metabolome variation in pos mode. The dominated components for the PC2 were phylloquinone, eriojaposide B and others in neg mode, and 1D-1-Guanidino-1-deoxy-3-dehydro-scyllo-inositol, strictosamide and others in pos mode, correspondingly, arbutin (59.91fold) and strictosamide; (2.17-fold) were obviously up-accumulated in neg mode and pos mode, respectively ( Fig. 6C and Fig. 6D).

Comprehensive Networks Of Transcriptomes And Metabolomes
To explore the speci c transcriptional characterization of the herbivory responses, 3,666 speci c DEGs in L. dispar herbivory group were classi ed into secondary metabolism pathways according to KEGG pathway analysis, which showed a large amount of DEGs involved into the phenylpropanoid and avonoid biosynthesis ( Fig. 7 and Supplementary Table S2). Thirty eight and thirty DEGs were involved into phenylpropanoid biosynthesis pathway and avonoid biosynthesis pathway, respectively.
In the phenylpropanoid biosynthesis pathway, a phenylalanine ammonia-lyase (PAL) gene was upregulated in the L. dispar herbivory group and down-regulated in the mechanical wounding group compared with control group. Of six 4-coumarate-CoA ligase (4CL) genes, all 4CL genes were signi cantly induced in the L. dispar herbivory group compared with the mechanical wounding group. Of ten cinnamoyl-CoA reductase (CCR) genes, expression levels of six CCR genes in the L. dispar herbivory group were 2.31 ~ 12.89-fold of those in the mechanical wounding group. Compared with the mechanical wounding group, three of four cinnamyl-alcohol dehydrogenase (CAD) genes were overexpressed in the L. dispar herbivory group. Of fourteen peroxidase (POD) genes, six POD genes were up-regulated in the L. dispar herbivory group and in the mechanical wounding group while eight POD genes were suppressed compared with control group. A shikimate O-hydroxycinnamoyltransferase (HCT) gene was downregulated in the L. dispar herbivory group and the mechanical wounding group, a beta-glucosidase (BGL) gene was up-regulated in the L. dispar herbivory group compared with the mechanical wounding group.
In the avonoid biosynthesis pathway, ve chalcone synthase (CHS) genes and two chalcone isomerase (CHI) genes were signi cantly regulated in the L. dispar herbivory group with 10.53 ~ 23.41-fold and 6.70 7.54-fold of those in mechanical wounding group, respectively. Compared with the control group, ve avanone 3-hydroxylase (F3H) genes were overexpressed in the L. dispar herbivory group and two of them were suppressed in the mechanical wounding group. Moreover, the expression levels of F3H in L. dispar herbivory group were 2.82 ~ 10.88-fold of those in mechanical wounding group. A avonoid 3'hydroxylase (F3'H) gene was induced in the L. dispar herbivory group with 6.07-fold of that in the mechanical wounding group. Except for one dihydro avonol 4-reductase (DFR) gene, four DFR genes were up-regulated in the L. dispar herbivory group with 8.24 ~ 15.41-fold of those in the mechanical wounding group. Of three avonol synthase (FLS) genes, one FLS gene was overexpressed in the L. dispar herbivory group and was suppressed in the mechanical wounding group. Two anthocyanidin synthase (ANS) genes were signi cantly induced in the L. dispar herbivory group compared with those in the mechanical wounding group. Five of six leucoanthocyanidin reductase (LAR) genes were overexpressed in the L. dispar herbivory group compared with those in the mechanical wounding group. An anthocyanidin reductase (ANR) gene was down-regulated in the L. dispar herbivory group and in the mechanical wounding group (Fig. 7).
DiscussionPlants can reorganize their transcriptome, proteome and metabolome to deal with abiotic and biotic stress [11]. Plant resistance traits are extremely plastic against insect herbivores [12]. However, the protein and metabolite constituents of host plants are determined by the plant species, herbivore species, infested time, and detection techniques [13]. The C. suppressalis larvae herbivory results in 4729 genes and 151 metabolites differently regulated into rice plants. The defense-related phytohormones, transcript factors, and secondary metabolism mediated by shikimate and terpenoid were activated while the growth-related counterparts were suppressed after C. suppressalis herbivory [14]. Latexes protect the Ficus carica from microbes and herbivores attack by toxic proteins and metabolites produced in immature fruit, young petioles and ligni ed trunks [15]. In present study, the speci c DEGs and differential metabolites were con rmed according to transcriptomes and metabolomes of P. simonii × P. nigra in response to attack by L. dispar larvae. Gene expression analyses revealed that differently DEGs were involved into L. dispar herbivory.
Variable gene expressions resulted in recon guration of primary and/or secondary metabolism after insect infestation [16,17]. After pre-infested by the O. furnacalis, DEGs in maize plants are associated with resistance to herbivores by affecting some metabolic pathways such as benzoxazinoids, phytohormones, volitale organic compounds (VOC) pathways [18]. Our results showed that a large amount of genes related to the biosynthesis of phenylpropanoid and avonoid were activated in the herbivory group and some secondary substances obviously induced such as anethole, methylchavicol, chavicol, p-coumaryl alcohol and leucodelphinidin. Interestingly, twenty metabolites, such as lolicine A, cimicifugoside and vinblastine sulfate were specially associated with L. dispar infestation (Supplementary Table S3). In the phenylpropanoid biosynthesis pathway, one PAL, six 4CL, six CCR, six CCR, three CAD, six POD, and one BGL genes in poplar were up-regulated after L. dispar herbivory. For example, tobacco (Nicotiana tabacum) overexpressing Arabidopsis transcription factor, AtMYB12, enhanced expressions of key genes involved in the phenylpropanoid pathway, resulting in higher accumulation of avonols [19]. PAL, CAD and 4CL are key enzyme and/or rate-limiting enzyme in lignin biosynthesis, which play an important role in plant responses to biotic and abiotic stresses [21,22,23]. As a scavenger of reactive oxygen species (ROS), POD is usually generated in early resistance responses of plants to a variety of pathogens through blocking oxidant-mediated programmed cell death [24,25]. The 4CL3 in Isatis indigotica hairy roots elicited by methyl jasmonate (MeJA) was con rmed to esterizate caffeic acid associated with lignan production using the transcriptional and metabolic analyses [26]. Oviposition by Spodoptera exigua increased Solanum dulcamara resistance to the subsequent feeding larvae using transcriptome and metabolome analyses. The genes involved in phenylpropanoid metabolism signi cantly overexpressed in previously oviposited S. dulcamara plants by lepidopteran herbivore S. exigua [12]. Moreover, pro les of volatile emission did vary in the indirect defenses of Populus nigra to leaf herbivory by lepidopteran species (Lymantria dispar, Laothoe populi and Amata mogadorensis) and beetles species (Phratora vulgatissima and Chrysomela populi) [27]. The transcriptome and metabolome analyses showed the ROS and calcium in early signaling, and salicylic acid (SA) and jasmonic acid (JA) mediated the defence regulation of Arabidopsis thaliana attacked by Brevicoryne brassicae. Secondary metabolites changes revealed defense against aphids, such as accumulation of 4-methoxyindol-3-ylmethyl glucosinolate [28]. Flavonoids synthesized by the phenylpropanoid pathway participate in multiple physiological and biochemical processes in plants. In the avonoid biosynthesis pathway, our results showed ve CHS, two CHI, ve F3H, a F3'H, four DFR, a FLS, two ANS, ve LAR, and an ANR genes in poplar were signi cantly induced in response to L. dispar herbivory. However, insect can change biochemical and physiology to combat with plant secondary metabolites. Our group reported effects of six characteristic poplar secondary metabolites (caffeic acid, salicin, rutin, quercetin, avone, and catechol) and three mixtures ( avone + salicin, salicin + caffeic acid + catechol, avone + catechol + caffeic acid) on performance and GST and P450 activities of L. dispar [29].
In conclusion, the transcriptome and metabolome analyses were conducted to generate a series of data set concerning P. simonii × P. nigra defense attacked by L. dispar. Speci c genes and metabolic substances in poplar were induced by lepidopteran L. dispar. The phenylpropanoids and avonoid biosynthesis pathway involved into the defense responses to L. dispar. The comprehensive networks of transcripts and metabolites provide novel insights into poplar defense mechanisms in this study.

Conclusions
In this study, the speci c mechanism of poplar was investigated in response to L. dispar herbivory based on transcriptomic and metabolomic analyses. The defense responses of poplar are partially similar by L. dispar herbivory and mechanical wounding treatments. The 3666 speci c DEGs were caused by L. dispar herbivory, and these speci c DEGs are mainly enriched in phenylpropanoid biosynthesis pathway (ko00940) and avonoid biosynthesis pathway (ko00941) by the analysis of KEGG metabolic pathway. The ubiquinone-6 and cimicifugoside showed signi cantly different expression in metabolome caused by L. dispar herbivory. In summary, the speci c different genes and metabolites are involved in the poplar defense caused by L. dispar herbivory, and phenylpropanoid biosynthesis (ko00940) and avonoid biosynthesis (ko00941) may be key defense pathways in poplar activated by L. dispar herbviory damage stress.

Plants and insects
The P. simonii × P. nigra was kindly provided by Prof. Jing Jiang from State Key Laboratory of Tree Genetics and Breeding, Northeast Forestry University, Harbin, Heilongjiang, China. Brie y, P. simonii × P. nigra was clonally propagated in tissue culture, transferred to soil and grown in a mixture of peat and sand (2:1 v/v) in a greenhouse with day temperatures of 24°C and night temperature of 12°C conditions. No chemical sprays were used into the poplar seedlings during the all experiments.
Lymantria dispar eggs were initially obtained from the Research Institute of Forest Ecology, Environment, and Protection, CAF, and stored at 4°C before hatching. Larvae were reared on arti cial diets at (25±1)°C under 14/10 h light/dark conditions, as described previously [30]. Healthy and hatching on the rst day of 4th instar L. dispar larvae were randomly selected for subsequent experimental processing.

Herbivory and mechanical wounding treatments
Thirty healthy 4th instar L. dispar larvae were conducted into P. simonii × P. nigra seedlings without pests and diseases with the similar diameter (mean: 18.31mm±2.15mm) and height (mean: 110.27cm±3.79cm). After the poplar plants were placed into the greenhouse for 48 hours, the L. dispar larvae pretreated with starvation for 6 h were uniformly transferred to the poplar plants, and the L. dispar larvae on the poplar were evenly re-distributed for four times per every day. To prevent larvae from escaping, the poplars were caged with a nylon net (Fig. 1). Corresponding to the treatment time and damage leaves of the L. dispar herbivory group, the poplar leaves was cut as mechanical wounding to simulate the L. dispar herbivory corresponding four times per day. After 24 h of treatment, all poplar leaves were used liquid nitrogen for quick freezing and stored in a -80°C refrigerator. The L. dispar herbivory group was abbreviated as GH, the mechanical wounding group was abbreviated as M, and the control group was abbreviated as CK. All treatments are performed in three biological replicates. One replicate was mixture of three technical replicate at equal ratios.
RNA extraction, library construction, and sequencing The CTAB method was used to extract total RNA according to the manufacturer's instructions [31]. DNase I was used to clean out DNA contaminant in RNA samples. The mRNA molecules in total RNA were isolated using oligo (dT)-attached magnetic beads and were broken into small pieces using fragmentation reagent. First-strand cDNA was generated by random hexamer-primed reverse transcription, followed by synthesis of second-strand cDNA. The synthesized cDNA was subjected to endrepair and following 3 -adenylated. And the ends of these 3 -adenylated cDNA fragments ligated adaptors. PCR products were puri ed with the solid phase reverse immobilization (SPRI) beadsand then were dissolved in EB solution. The double stranded PCR products were heating denatured and circularized by the splint oligo sequence. The single-strand circle DNA (ssCir DNA) was formatted as the nal library. The library was con rmed by Agilent Technologies 2100 bioanalyzer, and was ampli ed with phi29 to form DNA nanoball (DNB) which have more than 300 copies of one molecular. The DNBs were loaded into the patterned nanoarray and single end reads of 50 bp were generated in the way of sequenced by synthesis. The library was used for sequencing through the sequencing BGISEQ-500 platform (BGI, Shenzhen, China).

De novo assembly and annotation
Prior to assembly, the reads containing adaptor, unknown base N content > 5%, and low-quality bases (>20% of the bases with a quality score ≥10) were removed by SOAPnuke (https://github.com/BGIexlab/SOAPnuke, the BGI, Shenzhen, China). After ltering, the remaining reads were called "clean reads" and stored in FASTQ (A format that stores biological sequences and corresponding quality assessments) format. De novo assembly of clean reads was performed using Trinity (v2.0.6, min_kmer_cov:3). Then, we used the TGI Clustering Tool (TGICL) to cluster transcripts to unigenes [32].
To acquire comprehensive information on gene functions, assembled unigenes were annotated using four databases, including the NCBI non-redundant (Nr) protein sequences, eukaryotic ortholog groups (KOG), SwissProt, and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases using BLASTx and BLASTn with an E-value <10 -5 . The gene ontology (GO) annotation was performed by Blast2GO with an E-value <10 -5 based on the protein annotation results of the Nr database [33]. InterProScan5 was used for interPro functional annotation [34]. The expression level of unigenes was calculated using FPKM (fragments per kilobase pairs per million mapped reads) values, which were indicated by running RSEM (RNA-Seq by Expectation Maximization) [35].

Metabolomics sample preparation
Chemical extraction was carried out on nine samples (three groups × three biological treatments) for metabolic analyses. Each sample aliquots (25 mg) were treated with 800 μL extraction liquid (V methanol : V H2O = 1:1) and mixed well. The samples were ground using a tissue lyser (55 Hz, 5 min) into tubes where two magnetic beads were put into each tube, then were precipitated at -20°C for 2h. Samples were centrifuged at 25000 g for 20 min at 4°C, and supernatants were used for metabolomics analysis by ultraperformance liquid chromatography coupled to Q-TOF mass spectrometry (UPLC/Q-TOF MS).
The high-resolution tandem mass spectrometer was used to detect the metabolites eluted from the column. The Q-TOF was conducted in both positive and negative ion modes. The capillary and sampling cone voltages were set at 3.0 kV and 40.0 V in the positive (POS) ion mode, respectively. The capillary and sampling cone voltages were set at 2.0 kV and 40.0 V in the negative (NEG) ion mode, respectively. The mass spectrometry data were obtained in Centroid MSE mode. The TOF mass ranged from 50 to 1200 Da with 0.2 s of scan time. All precursors were fragmented using 20-40 eV for the MS/MS detection with 0.2 s of scan time. During the acquisition, the LE signal was acquired every 3 s to calibrate the mass accuracy. In addition, to evaluate the stability of the LC-MS during the whole acquisition process, a quality control (QC) sample was taken from every 10 samples.

Metabolomic data analysis
Data pretreatment of raw UPLC-MS data including peak picking, peak grouping, retention time (RT) correction, second peak grouping, and annotation of isotopes and adducts was performed using Progenesis QI software (Waters, Manchester, U.K) and the freely available R statistical language (metaX). LC-MS raw data les were initially converted into net CDF format, which is further processed by the XCMS and CAMERA toolbox. RT and the m/z data pairs were used as identi ers for each ion.
To obtain consistent variables, the resulting matrix was further reduced by removing peaks with more than 80% missing values (ion intensity = 0) and those with isotope ions from each group. To ensure the high quality of the data within the analysis process, an approach based on the periodic analysis of a standard biological QC sample and the true samples was used as a quality assurance strategy in metabolic pro ling. Each retained peak was normalized to the QC sample by using robust Loess signal correction. The relative standard deviation (RSD) value of metabolites in the QC samples was set a threshold of 30% to assess the repeatability of the metabolic data sets. The logarithmic transformation and pareto scaling were applied to the datasets before statistical analysis.
Principal components analysis (PCA) was used to investigate and visualize the changes of metabolite in an unsupervised manner. Then, t test and fold change (FC) analysis (nonparametric univariate method [Mann-Whitney-Wilcoxon test]) were performed to measure the importance of metabolic characteristics from the preprocessed UPLC-MS datasets. The p-values generated by statistical test were adjusted to obtain q-values by using false discovery rate (FDR) correction. The screening requirement for differential metabolites was FC ≥ 1.2 or ≤0.8333.

Metabolite annotation and identi cation
The progenesis QI software was used for annotations and identi cation of metabolites. Brie y, the exact molecular masses and their fragmentation traits were annotated using public databases, Human Metabolome Database (HMDB; http://www.hmdb.ca/) and KEGG (www.genome.jp/kegg/) databases. Metabolites searched from databases were fragmented to obtain theoretical fragment spectrum. Then, the similarity scores were determined by comparing ion secondary spectrum from deconvolution of observed and theoretical spectra. The speci c metabolites were screened according to the scores and difference between observed and theoretical mass. The isotopic distribution measurement was further used to identify the molecular formulae of matched metabolites. All plant materials used in our study were permitted and obtained from the State Key Laboratory of Tree Genetics and Breeding, Northeast Forestry University. Experimental research on plant, including collection of plant materials, complied with institutional, national, or international guidelines.

Consent for publication
Not applicable.

Availability of data and materials
The relevant datasets in this study are included within this article and the supplementary les.

Competing interests
The authors declare that they have no competing interests.