Integrated Transcriptomic and Metabolomic Analysis of Lymantria dispar multiple nucleopolyhedrovirus (LdMNPV) infected Lymantria dispar in Tree-top Disease

The pioneer record described the phenomenon of nucleopolyhedroviruses (NPVs) manipulating host behavior was named as Tree Top Disease or “Wipfelkrankheit” by Hofmann in 1891. Following publics advised multipartite effect during infection progressing that NPVs adjusting manipulated mode to actualize crossing barriers, immune escape and hormone regulation. In this study, the molecular mechanism of Lymantria dispar Multiple embedded Nucleopolyhedrovirus(LdMNPV)-infected was investigated by a comparison of transcriptomic and metabolomic analysis at four selected time points of LdMNPV-infected and healthy Lymantria dispar (Ld) larvae. Among data collecting, the notably up-regulated and down-regulated genes and metabolomics were used for further analysis. Integrative analysis revealed the content changes of caffeine, glutamate and gamma-Aminobutyric acid involved in LdMNPV inducing behavior change, and the alteration of myo-inositol, phosphatidyl-1D-myo-inositol and relevant genes involved in larval extended instar duration and apoptosis avoiding. versus Ld6 the enriched gene sets are numerous in quantity, means larger variations in these two groups. Between NB3 and CK3, the enriched gene sets involved material metabolism and surroundings signal receptor, like phototransduction-y and neuroactive ligand-receptor interaction, reecting the larval physiological change of growing process. In the combination as NB3 versus Ld3 and NB3 versus Ld6 the associating gene sets, most of which present an up-regulated trend, mainly included cellular nutrition utilization, energy metabolism and material metabolism, which is shown as complicated network painted in violet background. In this two comparison, the up-regulated gene sets were peroxisome, phototransduction-y and insect hormone biosynthesis, the down-regulated gene sets were lysosome, apoptosis, proteasome, ubiquitin-mediated proteolysis, endocytosis regulation of autophagy, dorso ventral axis formation, protein export and ABC transporters. Compared with the group NB3 versus Ld3, the trend of group NB3 versus Ld6 enriched gene sets integrally appears down-regulation and the scanty up-regulated gene sets focused on energy and material metabolism. As the infection progresses further, it is revealed that in Ld3 versus Ld6 the alteration mainly consists of material metabolism and protein export, especially in fatty acid metabolism. Other gene sets were peroxisome, protein export and neuroactive ligand-receptor interaction. And the down-regulation emphasized the apoptosis. Between CK3 and Ld3, only the Citrate cycle (TCA cycle) and tyrosine metabolism showed the up-regulation. The down-regulated gene sets covered maintenance and development signal and genetic information processing. The decreasing extensively pathways of host genetic information processing accords with the viral manipulation in early replication. biosynthesis pathway, ABC transporters pathway, FoxO signaling pathway, phosphatidylinositol signaling system, neuroactive ligand-receptor interaction pathway, mTOR signaling pathway, protein digestion and absorption pathway, bile secretion pathway, parkinson disease pathway and alcoholism pathway. Among them, the ABC transporters pathway only keeps invariant in NB3 versus CK3 in NEG mode, which showed severe change. The aminoacyl-tRNA biosynthesis pathway was enriched in NB3 versus Ld3, NB3 versus Ld6 and CK3 versus Ld3 in both POS and NEG mode. The protein digestion and absorption pathway was enriched in NB3 versus Ld3 and NB3 versus Ld6 in both POS and NEG mode. The bile secretion pathway shows slight variation in NB3 versus CK3 and NB3 versus Ld6 in POS mode, and in NB3 versus Ld6 in NEG mode. The mTOR signaling pathway, parkinson disease and alcoholism were only changed in NB3 versus CK3,and the FoxO signaling pathway was only changed in Ld3 versus Ld6. It is remarkable that the changes of neuroactive ligand-receptor interaction pathway and phosphatidylinositol signaling system pathway, which showed incompatibility in POS and NEG mode. The change of neuroactive ligand-receptor interaction pathway appears at NB3 versus CK3, NB3 versus Ld3, Ld3 versus Ld6 and CK3 versus Ld3 in POS mode, and at NB3 versus Ld6 and CK3 versus Ld3 in NEG mode. And the change of phosphatidylinositol signaling system pathway occurs at CK3 versus Ld3 in POS mode, and at NB3 versus Ld3, Ld3 versus Ld6 and CK3 versus Ld3 in NEG mode. By comparing the different metabolites, it is reveal that in neuroactive ligand-receptor interaction pathway the incompatibility causing by detection of dopamine, acetylcholine, adenosine and trace amine were different in POS and NEG mode, and in phosphatidylinositol signaling system pathway Based on the interaction of KEGG pathways and dynamic change between each groups, it is lter out that the change of neurotransmitter, cyclohexanehexol and secondary metabolism may associate with host-parasites interaction. The neurotransmitter like dopamine and acetylcholine, increase in NB3 versus CK3 but remain stable in other compares. The caffeine metabolism pathway was enriched in NB3 versus Ld3 and CK3 versus Ld3. In CK3 versus Ld3, the caffeine was increase, and xanthine dehydrogenase/oxidase (XDH) gene was decrease, which relevant to purine metabolism by oxidizing xanthine. In C.elegans and Drosophila melanogaster cyclohexanehexol was reported to extend lifespan 23, 24 . The inositol phosphate metabolism pathway and phosphatidylinositol signaling system pathway were both enriched in CK3 versus Ld3 and Ld3 versus Ld6, which two pathways link by cyclohexanehexol. In CK3 versus Ld3, the myo-inositol (MI) was increase and PI was decrease. The inositol polyphosphate 5-phosphatase (INPP5A) and phosphatidylinositol-4-phosphate 3-kinase (PI3K) gene were upregulated and inositol 1,4,5-triphosphate receptor (IP3R) gene was downregulated. And in Ld3 versus Ld6, the PI was increase and MI was decrease. Meanwhile, the phosphatidylinositol-3,4,5-trisphosphate 3-phosphatase (PTEN/MMAC1), multiple inositol-polyphosphate phosphatase/2,3-bisphosphoglycerate 3-phosphatase (MINPP1) and ip3r gene were up-regulated, and the pi3k gene was down-regulated. Other metabolisms like pyruvate, N-acetylglucosamine (GluNAc) and NAD + also relevant to life prolonging. Among them, only GluNAc was increase, which was enriched in Ld3 versus Ld6. In Bombyx mori, after BmNPV infection the trehalose, malate, fumarate, citrate, succinate, α-ketoglutarate, phenylalanine and tyrosine showed accumulate in the resistant strain YeA, but in susceptible strain YeB the gluconate, glutamate and aspartate were upregulation 25 . The results showed the phenylalanine and tyrosine were kept increasing in CK3 versus Ld3 and Ld3 versus Ld6. In Arginine biosynthesis pathway the different metabolites were glutamine, glutamate and citruline, which show changing consistency. In NB3 versus CK3 this three all decreases, with glutamine synthetase (GS) gene, glutamate dehydrogenase (GDH) gene, acylamide amidohydrolase (ASL) gene and amido acid deacylase gene were downregulated, and NADPH-diaphorase gene was upregulated. In CK3 versus Ld3 and Ld3 versus Ld6, the glutamine, glutamate and citruline keep accumulate, but genes expression were altered. In CK3 versus Ld3, gs gene, gdh gene, asl gene and arginosuccinate synthetase (ASS1) gene were upregulated, amido acid deacylase gene was downregulated. In Ld3 versus Ld6, gs gene, and NAD(P) + oxidoreductase gene were upregulated, NADPH-diaphorase gene was downregulated. Moreover, in Ld3 versus Ld6 glutamate as the principal excitatory amino acid neurotransmitters, accumulate indicate targeted utilization by glutamine synthetase (GLUL) gene upregulated, which product was glutamine showed increase. And the glutamate decarboxylase gene was downregulated keeping the content of it product, gamma-Aminobutyric acid (GABA) working as a chief inhibitory neurotransmitter in reducing neuronal excitability, change after In CK3 versus Ld3, the alpha-D-glucoside glucohydrolase gene and alpha-D-galactoside galactohydrolase gene and Than in Ld3 versus Ld6, the expression level of alpha-D-glucoside glucohydrolase gene and alpha-D-galactoside galactohydrolase gene remained stable, but the galactinol, ranose, sucrose, mannose and MI present lack. rapid fashion to meet the high energy costs and avoid cell stress resistant and apoptosis of virus replication, and the increasing of lipids by altering the fatty acid synthesis, is needed to afford the additional membrane material for the viral particle envelopment, and immune function relevant pathways were not signicantly affected. Moreover, maintaining infected cell stability is important to virion production and assembly at terminal of infection and the extended instar duration of host was benetted to virus production. In summary, we explored the complicated interactions between LdMNPV and 3th-instar L. dispar larvae at the metabolome and transcriptome levels, conducted the preliminary analysis, uncovered that the content patterns of caffeine, glutamate and GABA involved in LdMNPV inducing behavior change, and the MI, PI and relevant genes involving in larval extended instar duration and apoptosis avoiding. It is anticipated that further research in relevent pathways will reveal in the complex reactions in conict between L. dispar and LdMNPV.


Introduction
In various parasites, host manipulating in protein transcription and expression, cellular metabolism, signal transduction and even behavior is substantive phenomenon, which is bene ted to increase fecundity and transmission rates [1][2][3] . The baculoviruses provide an example of parasite-induced host manipulation, by which infection resulting in a series of obvious behavior change in the hosts.
In retrospect, the intriguing phenomenon was rst recorded in 1891 and was named "Tree Top Disease" or "Wipfelkrankheit" by Hofmann 4 .
Till the late 1940s, it was con rmed the diseases causing by virus 5 . Baculoviruses are classi ed as a family of large rod-shaped circular double-stand DNA viruses, which speci cally infect arthropods and have a unique replication 6 . Relying on the adaption for insect growth pattern and life cycle, the baculoviuses acquired a biphasic replication cycles, which including two individual virion phenotypes, the occlusion derived viruse (ODV) and the budded viruse (BV) 7,8 . The two viral forms, being genetically identical but morphologically different, preformed a synergistic action, giving expression to functions in both cell-to-cell (BV) and insect-to-insect (ODV) transmission, and also bene ted for increasing e ciency of infection and virus titer 5 . On the other hand, the unique replication model causing the interaction mechanism between virus and host was time-ordered, dynamic and complicated, which hardly to be revealed.
Baculoviruses lead typically abnormal behaviors on scheduling, which including increasing food intake, ELA, extended instar duration and climbing to elevated location before death 5 . As nonliving entity, virus needs to manipulate host cell program upon successfully infecting to survival and repilication. The effects of parasitifer infection were shown as behavior change in macroscopic aspect, which leading by manipulation of signaling transduction and alteration of cellular metabolism in microscopic aspect. The Enhanced Locomotory Activity (ELA) is typical early sign of baculoviruses infection 8 . In 3rd-instar Lymantria dispar larvae, the ELA was documented as appearing at average of 3.75 days post infection (dpi) and the phenomenon dying out at average of 4.75 dpi 9,10 . With systemic infection further expanding, ELA was the typical description of usurping caterpillar behavior, which normally occurred at healthy larval last instar and was administered by codeterminantion of larval physiological indexs like body size, hormone, phototaxis and geotaxis 11 . It is reported that in Bombyx mori the hormone metabolism manipulation is important to virus leading ELA 12 . The virus-induced climbing behavior also was refereed to such physiological alterations as geotaxis, phototaxis, vision, olfaction and rhythm 1 . Strinkly, the virus-induced climbing behavior happens in terminal, when the virus vastly replicating at all tissues, but the sick larvae still keep strong climbing ability. Goulson et al was hypothesized that the climbing behavior might be driven by sick larvae to get behavioral fever to respond to infection 4 .
Investigations have suggested that the protein tyrosine phosphatase (ptp) gene is required for inducing the climbing behaviors in Group I Alphabaculovirus, including bombyx mori Nucleopolyhedrovirus(BmNPV) and AcMNPV 9,13 . In 2011, Hoover et al. reported that the ecdysteroid uridine 5'-diphosphate (UDP)-glucosyltransferase (egt) gene of LdMNPV was necessary to trigger the host behavior changes 14 .
Similarly, the egt gene of Spodoptera exigua Multiple embedded Nucleopolyhedrovirus(SeMNPV) has been con rmed to have the same effect in hyperactivity and virus-induced climbing behavior 5,15 . By analyzing the transcriptome of silkworm larvae brains, it was reported that BmNPV-induced hyperactive behaviors were related to circadian rhythms, synaptic transmission and the serotonin receptor signal pathway 16 . Baculovirus were described to be equipped with photolyase-like genes, which have circadian clock regulation function and the UV-induced DNA repair activity 17 . Previous study was revealed that in LdMNPV-infected hyperactived larvae the photosensitivity was remarkably increased and the rhythmicity relevant genes, especially per and tim, was violently impacted 18 . The phenomenona occurring in host may operate by more complicated interaction as gene expression change, hormone regulation and metabolism alteration, which created by competition of control between parasite and host 19, 20 . However, the mechanism by which NPV leading behavior changes is unclear.
In the current study, based on the biphasic life pattern of virus and experience of larvae behavior alteration, four different time points were selected for the integrated analysis of transcriptomic and metabolomic comparison to be carried out. In order to reveal interactions between LdMNPV and Lymantria dispar and the interaction changes in chosen time points.

Results
Transcriptome sequencing and data analysis Before the analysis of transcriptome and metabolome, we drew the survival curve in LD90 (Lethal Dose, 90%) LdMNPV-infected and the result shown in Supplemental Figure S1. The curve shown two turning points in 72 and 144 hours post infenction (hpi). It is reported that the primary infection usually proceeded between 24 and 48 hpi, and the system infection occurred at later than 72 hpi 21,22 . Refer to infected larvae behavior change, LdMNPV biphasic infection mode and in ection points of survival curve are considered together, the 72 and 144 hpi (3 and 6 dpi) were chosen for analysis. Totally, we choose four time points as the newborn 3th-instar larvae (NB3), 72 hours 3th-instar larvae post mock-infected (CK3), 72 hours 3th-instar larvae post infected (Ld3) and 144 hours 3th-instar larvae post infected Four groups (NB3, CK3, Ld3 and Ld6) from twelve samples of L. dispar larvae transcriptomic sequencing respectively generated 63,602,453, 61,981,656, 60,543,314 and 60,633,947 raw reads, and was produced with an average insert size of 298.98bp. In all libraries more than 93.37% of the sequences had a quality score greater than Q30 (Supplemental Table. S1), which indicates that all libraries obtained high-quality raw reads. After qualities were ltered and de novo assembled, the clean reads were totally assembled into 56,923 unigenes by Trinity software, which represents 73.78 Gb of the transcription data. The average length of the unigenes was 1083.78bp, and the N50 was 1877bp (Supplemental Table S2; Supplemental Fig. S2A).
The functional identity of all unigenes was investigated. Homologous annotation for all unigenes was based on information within the GO, Pfam, Swissporot, KOG, eggNOG, Nr, and KEGG databases. Consequently, 26,562 unigenes were annotated at least in one database (Supplemental Table S3). The box-plot gure of all samples is exhibited in Supplemental Figure S2B, of which the FPKM (Fragments Per Kilobase of transcript per Million mapped reads) was used to describe the degree of discretization as evaluation criteria of the gene expression level. The Pearson correlation between the groups is an important parameter, which re ects the experimental reliability and indicates whether the sample selection is reasonable or not. The Pearson correlation and the correlation plot between two samples are exhibited in Supplemental Figure S2C.
By using the DESeq2 softwave, parameters were set as FDK (false discovery rate) < 0.01 and FC (fold change) ≥ 2; the unigenes from each of the two groups were compared to acquire the different expression genes (DEGs). The results list in Supplemental Table S4. The Cluster analysis was shown in Fig. 1A. It is indicated that the CK3, Ld3 and Ld6 were grouped together and separated with NB3, and the CK3 and Ld3 were similary than others groups also reveal the difference in infection proceeding. The MA-plot (Supplemental Fig. S3A-E) clearly reveals the overall distribution of gene expression levels and fold changes. And the results of MA-plot indicated that the unchanged genes distribution in Ld3 versus Ld6 was exhibited in fusiform, which was different with other groups. Based on the gene expression in the different samples, to identify the annotation of DEGs function, the number of annotated DEGs in different databases is listed as in Supplemental Table S5.
In GO annotations, these unigenes were divided into 53 categories belonging to biological processes, cellular components and molecular functions. The DEGs were speci ed to GO categories to determine functional classi cations. Related to biological process, most DEGs were involved in metabolic processes, cellular processes, single-organism processes, biological regulation, localization and response to stimulus. In the cellular component, the DEGs were involved in cells, cell parts, membranes, membrane parts and extracellular regions. In molecular functions, most GO terms concentrated upon binding and catalytic activity. The results of GO database annotated DEGs between groups were shown in Supplement Figure S4.
In KEGG annotations, a total of 10,157 unigenes were mapped onto 303 different KEGG pathways found in human diseases, genetic information processing, metabolism, organismal systems, environmental information processing and cellular processes in six categories.
The stacking bar chart of DEGs was shown in Fig. 1B. In general, in NB3 versus CK3 and Ld3 versus Ld6, the number of up regulation DEGs is more than down, and the amount of DEGs is smallest in Ld3 versus Ld6. It is also indicated that the immune system was not obviously activate in infection larvae, but Signal transduction, Xenobiotics biodegradation and metabolism, Carbohydrate metabolism and Lipid metabolism were signi cant enriched. Additionally, Replication and repair and Folding, sorting and degradation were drastic change with infection process. The Supplemental Table S6 list out the top 3 of drastic changed, up-regulated and down-regulated pathways, which were annotated between groups. And the bubble diagrams of pathway enrichment are shown in Fig. 1C-1G. According to that, we summarized that between two un-infection groups, NB3 versus CK3, the DEGs mainly enriched in small bioactive molecule metabolism like vitamin and cofactor, the extracellular matrix (ECM)-receptor interaction and the defense to external environment. The two compares, NB3 with Ld3 and NB3 with Ld6, mostly give expression to carbon and amino acid metabolism. Between Ld3 and Ld6, the DEGs enrichment pathways tend to focus on cell supporting in material. Especially the carbon metabolism pathway and biosynthesis of amino acids pathway both were drastic down regulated. But the glycosphingolipid biosynthesis pathway was increased. In CK3 versus Ld3, the enriched pathways involved ribosome cytochrome P450 and energy metabolism, which implicates alteration provide the condition of viral genome replication and satisfy the energy supply.
Furthermore, the DEGs were compared with the COG database(Cluster of Orthologous Groups of proteins)and resulted in 26 different functional categories (Supplemental Fig. S5). The DEGs also were annotated to the eggNOG database in different functional categories (Supplemental Fig. S6). Each of them provides panoramic account DEGs between groups from different angles.
Comparing with the DEGs annotations, the GSEA emphasizes the gene, which plays an important biological effect but only have a subtle change of expression, by enrichment the gene sets. The GESA result, which visualization by Cytoscape, was shown in Figure (Fig. 2A-2E). The topological structure describes the correlation characteristic of gene set, which combines the Jaccard similarity coe cient standard by distance and the overlap coe cient standard by green line. To separate different patterns, we pained the nodes background as violet represents metabolism, aqua represents genetic information processing, and pink represents environmental information processing. Generally, it is shown that in combination as NB3 versus Ld3 and NB3 versus Ld6 the enriched gene sets are numerous in quantity, means larger variations in these two groups. Between NB3 and CK3, the enriched gene sets involved material metabolism and surroundings signal receptor, like phototransduction-y and neuroactive ligand-receptor interaction, re ecting the larval physiological change of growing process. In the combination as NB3 versus Ld3 and NB3 versus Ld6 the associating gene sets, most of which present an up-regulated trend, mainly included cellular nutrition utilization, energy metabolism and material metabolism, which is shown as complicated network painted in violet background. In this two comparison, the up-regulated gene sets were peroxisome, phototransduction-y and insect hormone biosynthesis, the down-regulated gene sets were lysosome, apoptosis, proteasome, ubiquitin-mediated proteolysis, endocytosis regulation of autophagy, dorso ventral axis formation, protein export and ABC transporters. Compared with the group NB3 versus Ld3, the trend of group NB3 versus Ld6 enriched gene sets integrally appears down-regulation and the scanty up-regulated gene sets focused on energy and material metabolism. As the infection progresses further, it is revealed that in Ld3 versus Ld6 the alteration mainly consists of material metabolism and protein export, especially in fatty acid metabolism. Other gene sets were peroxisome, protein export and neuroactive ligand-receptor interaction. And the down-regulation emphasized the apoptosis. Between CK3 and Ld3, only the Citrate cycle (TCA cycle) and tyrosine metabolism showed the up-regulation. The down-regulated gene sets covered maintenance and development signal and genetic information processing. The decreasing extensively pathways of host genetic information processing accords with the viral manipulation in early replication.
The transcriptomic analysis shown that after being infected, the metabolism, signal transductions and detoxication relevant metabolism of L. dispar larvae were visibly affected, but immune function relevant pathways, like the JAK/STAT signaling pathway, NF/kB signaling pathway and Toll-like Receptors pathway, were not signi cantly affected. Further more, with infection progressing, the apoptosis was continuous down-regulation.

Metabolite extraction and analysis
Other than protein-interaction, metabolism is constituted by chemical reactions, which is described metabolites as nodes and functional proteins as edges. Metabolomic analysis of twenty-four samples from four groups (NB3, CK3, Ld3 and Ld6) of L. dispar larvae was performed via UHPLC-QTOF-MS. As a result, 1,930 and 1,831 metabolins were respectively obtained using POS and NEG mode detection.
Based on the PCA Analysis (Supplement Fig. S7) and Cluster Analysis ( Fig. 3A-B), the quali ed metabolite data was subjected to differential analysis. After the metabolites were detected by the qualitative analysis, the quantities were compared between groups to show the extent the metabolites changed. The standard P-value < 0.05 of Student's T tests and VIP-values > 2 of OPLS-DA were used to screen the different metabolites. The amount of different metabolites between groups was shown in Table 1. To shown the variability we clustered the different metabolins, the result was shown in Supplemental Figure S8. And the Wayne diagram exhibits in all groups the relationship of the change in metabolites ( Fig. 3C-D). Taken together, the results re ected that in all comparatives the metabolite changes in NB3 versus CK3 and Ld3 versus Ld6 were easement, but in CK3 versus Ld3 was erce.
To identify the pathways of metabolites involved, a total of 1,359 metabolites in the POS mode and 1,341 metabolites in the NEG mode were mapped to KEGG metabolic pathways; in POS and NEG mode, the stacking bar charts of different metabolites in KEGG enrichment were shown in Fig. 3E-F; the details are listed in Table 1; and the bubble diagram of different metabolites in KEGG enrichment were shown in Fig. 3G-P. The stacking bar charts were shown that LdMNPV leading a largescale alteration in metabolism and organismal systems of infected larvae. The results revealed that the different metabolites annotated into KEGG database were concentrated in some key pathways. Generally, in POS mode, metabolic pathway, ABC-transporters and biosynthesis of amino acid were indicated to undergo a violent change in all comparisons. An exception to the above is that protein digestion and absorption emerged in NB3 versus CK3, NB3 versus Ld3,Ld3 versus Ld6, and CK3 versus Ld3. Meanwhile in NB3 versus Ld6, it is shown as glycerophospholipid metabolism. In NEG mode, metabolic pathway, biosynthesis of amino acid, 2-Oxocarboxylic acid metabolism and ABC-transporters were indicated to a violent change in all comparisons. An exception to above, the central carbon metabolism in cancer emerged in NB3 versus Ld3, NB3 versus Ld6 and CK3 versus Ld3; the protein digestion and absorption emerged in NB3 versus CK3, Ld3 versus Ld6 and CK3 versus Ld3; moreover in CK3 versus Ld3, it is shown as aminoacyl-tRNA biosynthesis and biosynthesis of unsaturated fatty acids.
The metabolomics analysis has shown that the LdMNPV induced a largescale alteration in cellular metabolism of infected larvae. In normal physiology status, from NB3 to CK3, the different metabolites, which were enriched into metabolism pathways evenly distribute in amino acid, lipids, and other small molecular metabolistes. After be infected, from CK3 to Ld3, the pathways including amino sugar and nucleotide sugar metabolism, biosynthesis of unsaturated fatty acid biotin metabolism, D-glutamine and D-glutamate metabolism, nitrogen metabolism and galactose metabolism were enriched dramatic change. From Ld3 to Ld6, arginine biosynthesis, galactose metabolism, purine metabolism, amino sugar and nucleotide sugar metabolism, fatty acid biosynthesis and galactose metabolism pathways were enriched further change. Though different viruses appear to activate different metabolic altering pattern, the reprogramming of host cell is determinant to viral growth. In this context, it was pertinent to conclude that the alteration in fatty acid synthesis, carbon source utilization and biosynthesis of amino acid was changed to guarantee optimal microenvironment for LdMNPV's replication and spread. Other wise, after infection the GABAergic synapse, cholinergic synapse, doparninergic synapse, glutarnatergic synapse, neuroactive ligand-receptor interaction and phosphatidylinositol signaling system also constant with infection processing, indicating the viral infection may affect on host's nervous and signal transduction.
Integrative analysis of metabolome and transcriptome To explore further dynamic variation of the infection-caused behavior change in L. dispar, an integrated and detailed analysis was performed on the transcript and metabolic levels. Before screening data, the Pearson correlation coe cient (PCC) and corresponding Pvalue were used to ltrate the data simultaneously, with values set to PCC > 0.80 and PCCP < 0.05.
Totally, the integrative analysis of metabolome and transcriptome screens 52 pathways in the POS mode and 56 pathways in the NEG mode. The amount of different expression genes and metabolites in KEGG pathway was list on Table. 2. The Wayne diagram exhibits in all groups the relationship of the change in covering KEGG pathways ( Fig. 4A-B). The stacking bar charts of different expression genes and metabolites in KEGG classi cation in POS and NEG mode were shown in Fig. 4C and corresponding quantity of different expression genes and metabolites in each KEGG classi cation were shown in Fig. 4D-E. The results showed that in POS mode the speci c differences in NB3 versus CK3 and Ld3 versus Ld6 are covering more KEGG pathways than other compares, but in NEG mode it have more similarity.
In general, Global and overview maps, Carbohydrate metabolism and Amino acid metabolism present signi cant difference in both of POS and NEG mode, but the amount of difference expression genes and metabolites in Ld3 versus Ld6 were smaller than others, especially the CK3 versus Ld3. Strikingly, the difference of energy metabolism were appeared at POS mode in NB3 versus CK3 and Ld3 versus Ld6, and at NEG mode in NB3 versus Ld3, NB3 versus Ld6 and CK3 versus Ld3.
The interaction and enrichment frequency of covering KEGG pathways were visualized in Fig. 4F to dig out more details. The gure showed that Carbohydrate metabolism, Amino acid metabolism and Lipid metabolism were enriched to most pathways, as Carbohydrate metabolism have 15 pathways, Amino acid metabolism have 14 pathways, and Lipid metabolism have 11 pathways. The KEGG classi cations of Global and overview maps and Nucleotide metabolism were violently uctuation. What's more, the Alanine, aspartate and glutamate metabolism pathway, Lysine degradation pathway and Tyrosine metabolism pathway in Amino acid metabolism, Glycerophospholipid metabolism pathway in Lipid metabolism, ABC transporters pathway in Translation and neuroactive ligand-receptor interaction pathway in Membrane transport were also drastic uctuation. Further analysis showed that the speci c enriched KEGG pathways were arachidonic acid metabolism, propanoate metabolism, mTOR signaling pathway, parkinson disease and alcoholism in NB3 versus CK3; oxidative phosphorylation in NB3 versus Ld3; D-Glutamine and D-glutamate metabolism, nicotinate and nicotinamide metabolism, one carbon pool by folate and FoxO signaling pathway in Ld3 versus Ld6; and glycerolipid metabolismin CK3 versus Ld3.
Except the 54 pathways of each metabolism classi cations, other enrichment only covered 10 pathways, including aminoacyl-tRNA biosynthesis pathway, ABC transporters pathway, FoxO signaling pathway, phosphatidylinositol signaling system, neuroactive ligandreceptor interaction pathway, mTOR signaling pathway, protein digestion and absorption pathway, bile secretion pathway, parkinson disease pathway and alcoholism pathway. Among them, the ABC transporters pathway only keeps invariant in NB3 versus CK3 in NEG mode, which showed severe change. The aminoacyl-tRNA biosynthesis pathway was enriched in NB3 versus Ld3, NB3 versus Ld6 and CK3 versus Ld3 in both POS and NEG mode. The protein digestion and absorption pathway was enriched in NB3 versus Ld3 and NB3 versus Ld6 in both POS and NEG mode. The bile secretion pathway shows slight variation in NB3 versus CK3 and NB3 versus Ld6 in POS mode, and in NB3 versus Ld6 in NEG mode. The mTOR signaling pathway, parkinson disease and alcoholism were only changed in NB3 versus CK3,and the FoxO signaling pathway was only changed in Ld3 versus Ld6. It is remarkable that the changes of neuroactive ligandreceptor interaction pathway and phosphatidylinositol signaling system pathway, which showed incompatibility in POS and NEG mode.
The change of neuroactive ligand-receptor interaction pathway appears at NB3 versus CK3, NB3 versus Ld3, Ld3 versus Ld6 and CK3 versus Ld3 in POS mode, and at NB3 versus Ld6 and CK3 versus Ld3 in NEG mode. And the change of phosphatidylinositol signaling system pathway occurs at CK3 versus Ld3 in POS mode, and at NB3 versus Ld3, Ld3 versus Ld6 and CK3 versus Ld3 in NEG mode. By comparing the different metabolites, it is reveal that in neuroactive ligand-receptor interaction pathway the incompatibility causing by detection of dopamine, acetylcholine, adenosine and trace amine were different in POS and NEG mode, and in phosphatidylinositol signaling system pathway was phosphatidyl-1D-myo-inositol (PI).
Based on the interaction of KEGG pathways and dynamic change between each groups, it is lter out that the change of neurotransmitter, cyclohexanehexol and secondary metabolism may associate with host-parasites interaction. The neurotransmitter like dopamine and acetylcholine, increase in NB3 versus CK3 but remain stable in other compares. The caffeine metabolism pathway was enriched in NB3 versus Ld3 and CK3 versus Ld3. In CK3 versus Ld3, the caffeine was increase, and xanthine dehydrogenase/oxidase (XDH) gene was decrease, which relevant to purine metabolism by oxidizing xanthine. In C.elegans and Drosophila melanogaster cyclohexanehexol was reported to extend lifespan 23,24 . The inositol phosphate metabolism pathway and phosphatidylinositol signaling system pathway were both enriched in CK3 versus Ld3 and Ld3 versus Ld6, which two pathways link by cyclohexanehexol. In CK3 versus Ld3, the myo-inositol (MI) was increase and PI was decrease. The inositol polyphosphate 5-phosphatase (INPP5A) and phosphatidylinositol-4-phosphate 3kinase (PI3K) gene were upregulated and inositol 1,4,5-triphosphate receptor (IP3R) gene was downregulated. And in Ld3 versus Ld6, the PI was increase and MI was decrease. Meanwhile, the phosphatidylinositol-3,4,5-trisphosphate 3-phosphatase (PTEN/MMAC1), multiple inositol-polyphosphate phosphatase/2,3-bisphosphoglycerate 3-phosphatase (MINPP1) and ip3r gene were up-regulated, and the pi3k gene was down-regulated. Other metabolisms like pyruvate, N-acetylglucosamine (GluNAc) and NAD + also relevant to life prolonging. Among them, only GluNAc was increase, which was enriched in Ld3 versus Ld6. In Bombyx mori, after BmNPV infection the trehalose, malate, fumarate, citrate, succinate, α-ketoglutarate, phenylalanine and tyrosine showed accumulate in the resistant strain YeA, but in susceptible strain YeB the gluconate, glutamate and aspartate were upregulation 25 . The results showed the phenylalanine and tyrosine were kept increasing in CK3 versus Ld3 and Ld3 versus Ld6. In Arginine biosynthesis pathway the different metabolites were glutamine, glutamate and citruline, which show changing consistency. In NB3 versus CK3 this three all decreases, with glutamine synthetase (GS) gene, glutamate dehydrogenase (GDH) gene, acylamide amidohydrolase (ASL) gene and amido acid deacylase gene were downregulated, and NADPH-diaphorase gene was upregulated. In CK3 versus Ld3 and Ld3 versus Ld6, the glutamine, glutamate and citruline keep accumulate, but genes expression were altered. In CK3 versus Ld3, gs gene, gdh gene, asl gene and arginosuccinate synthetase (ASS1) gene were upregulated, amido acid deacylase gene was downregulated. In Ld3 versus Ld6, gs gene, and NAD(P) + oxidoreductase gene were upregulated, NADPH-diaphorase gene was downregulated. Moreover, in Ld3 versus Ld6 glutamate as the principal excitatory amino acid neurotransmitters, accumulate indicate targeted utilization by glutamine synthetase (GLUL) gene upregulated, which product was glutamine showed increase. And the glutamate decarboxylase gene was downregulated keeping the content of it product, gamma-Aminobutyric acid (GABA) working as a chief inhibitory neurotransmitter in reducing neuronal excitability, remains stable. The carbohydrate metabolism also revealed noteworthy change after infection. In CK3 versus Ld3, the alpha-D-glucoside glucohydrolase gene and alpha-D-galactoside galactohydrolase gene were upregulated, causing accumulate of galactinol, ra nose, sucrose, stachyose, mannose and MI. Than in Ld3 versus Ld6, the expression level of alpha-D-glucoside glucohydrolase gene and alpha-D-galactoside galactohydrolase gene remained stable, but the galactinol, ra nose, sucrose, mannose and MI present lack.

Discussion
Infection with a wide variety of parasites, including baculovirus, can result in behavior changes in the host. For survival and replication, the parasites have to conquer host defenses barriers of immunity, even manipulate the host's physiology [1][2][3] . Some reports have provided insight into how baculovirus producing in uences on host, depending on some core genes or even just one gene 9, 13-16 . Within the last decade, the primary host behavior manipulation gene of NPV had been brought to light 5, 9, 13-18, 26-29 . The details on host reacting and its mechanism to virus-induced activity rarely had been addressed. In this study, through the integrative analysis of metabolome and transcriptome between the key time points, which were chosen by critical moment of replication characteration and infected host larval behavior changes, and we sum up the alteration of MI, PI and relevant genes involved in larval extended instar duration, and the changes of caffeine, glutamate and GABA involved in LdMNPV inducing behavior change.
Parasites are capable of inducing host behavior alterations in many aspects to increase transmission rates [1][2][3]5 . The NPV infection triggered abnormal physiological status, like increasing food intake, ELA, extended instar duration and climbing to elevated location before death foraging 5 , which leading by alteration of host cellular metabolism and manipulation of signaling transduction. The infected L.
dispar larvae had prolonged their 3th instar duration, which was 4 to 5 days in healthy larvae, to 7 to 8 days 10 . The survival curve in LD90 of LdMNPV-infection (Supplemental Figure S1) was shown that the 3 and 6 dpi were two in ection points. Baculoviuses specially have a biphasic replication cycle, which including two individual virion phenotypes, the ODV and the BV 5, 7 . The ODV and BV were genetically identical but morphologically different, which also were related-time produnction. Combined with the biphasic infections and the early infection mainly happening in midgut which only regional impact, the time point Ld3 and Ld6 were chosen to indicate the proceeding of systemic infection at early and imminent stages of massive replication.
In the four groups (NB3, CK3, Ld3 and Ld6) of L. dispar larvae, twelve transcriptomic samples were sequenced. On the whole, the CK3 and Ld3 were more similary than others groups (Fig. 1A). Based on DEGs annotation, it is showed that the pattern of NB3 versus CK3 and Ld3 versus Ld6 were approximate, and the alteration in cellular mainly concentrates on bene ting virus DNA replication, the energy supply and detoxication relevant metabolism were visibly affected, and immune function relevant pathways, like the JAK/STAT signaling pathway, NF/kB signaling pathway and Toll-like receptors pathway, were not signi cantly affected (Fig. 1B). After being infenction, the DEGs analysis showed regulatory strategy change in different infection stages, involving affect material supply of virion production and antiapoptosis to prolong infected cell survival. Furthermore, by KEGG annotations, the DEGs enrichment showed that between CK3 and Ld3, the violently upregulated pathways were metabolism of xenobiotics by cytochrome P450, drug metabolism-cytochrome P450 and galactose metabolism; and the drastically downregulated pathways were ribosome, carbon metabolism and phagosome, which reveals alteration provide the condition of viral replication, alteration of energy supply and avoid cell damage, toxin and reactive oxygen species (ROS) stress causing by virus repilication. Between Ld3 and Ld6, the violently upregulated pathways were carbon metabolism, phenylalanine metabolism and biosynthesis of amino acids; and the drastically downregulated pathways were lysosome, mTOR signaling pathway and longevity regulating pathway-multiple species. The upregulated of carbon metabolism, phenylalanine metabolism and biosynthesis of amino acids con rmed infection causing globally alteration of material utilization. The downregulated of lysosome, mTOR signaling pathway and longevity regulating pathway were bene t in keeping host cell survival. The genes, which were enrichment in lysosome pathway involving in metabolism disorder, and the upregulated gene could promote glycerolipid collection 30 . The drastic changes of longevity regulating pathway may relate to prolonging instar duration of infected larvae. It is also reported that dietary restriction and genetic down-regulation of nutrient-sensing pathways, like Insulin-like signaling pathway and mTOR signaling pathway can substantially increase life span 31 .
Generally, the GESA result ( Fig. 2A-E) is shown that NB3 versus Ld3 and NB3 versus Ld6 the gene sets which were enriched in have more volume, reveals the larger variations in these two comparisons. From NB3 to CK3 ( Fig. 2A), the enrichment was involved in material metabolism and surroundings signal receptor, which re ecting the larval physiological change of growing process. Between CK3 and Ld3 (Fig. 2E), only the Citrate cycle (TCA cycle) and tyrosine metabolism showed the upregulation, and lots of gene sets downregulated mainly covered maintenance and development signal and genetic information processing, which was accords with the viral manipulation in early replication. The transcriptomic analysis shown that pathways of metabolism, signal transductions and detoxication relevant metabolism were visibly affected in infection L. dispar larvae. Further more, with the infection proceeding, it is revealed that in Ld3 versus Ld6 (Fig. 2D) the alteration mainly consists of material metabolism and protein export, especially in fatty acid metabolism, and the apoptosis was continuous downregulation.
Twenty-four metabolomics samples obtained 1,930 metabolins in POS mode and 1,831 metabolins in NEG mode. In contrast between groups, a large variety of metabolites show signi cant reduction in the virion production stages (Fig. 3A-B). Results of the amount of different metabolites between groups, the variability of different metabolins and the Wayne diagrams (Fig. 3C-F) re ected that the metabolite changes in NB3 versus CK3 and Ld3 versus Ld6 were easement, but in CK3 versus Ld3 was erce. The metabolin enrichment showed the metabolic pathway, ABC-transporters, 2-Oxocarboxylic acid metabolism, biosynthesis of amino acid and protein digestion and absorption merged to produce a violent change in all comparisons. Meanwhile, the alteration of central carbon metabolism in cancer, which was accredited in lots of human source virus infections, only emerged in NB3 versus Ld3, NB3 versus Ld6 and CK3 versus Ld3. In rapidly replication viruses, the hyper-activation of sustained metabolic pathways was su ced to increasing viral titers 32,33 . The results of different metabolites KEGG enrichment revealed that the different metabolites were concentrated in some key pathways. From NB3 to CK3, as normal physiology status, the different metabolites have been enriched evenly distribute in amino acid, lipids, and other small molecular metabolistes. After be infected, from CK3 to Ld3, the enriched pathways covering amino sugar and nucleotide sugar metabolism, biosynthesis of unsaturated fatty acid biotin metabolism, D-glutamine and D-glutamate metabolism, nitrogen metabolism and galactose metabolism were enriched dramatic change. The nutrients uptake and utilization alteration tend to provide enough ATP to meet the energy costs of virus DNA replication. The increasing of lipid material was induced by the needed of additional membrane material for the viral particle envelopment. Though different viruses appear to activate different metabolic altering pattern, the reprogramming of host cell is determinant to viral growth. Since Munger et al reporting the rst eukaryotic virus tampering cellular metabolomics, different studies shown variously manipulate multiple major metabolic pathways to reprogram host cellular nutrients uptake and utilization to meet the increasing material and energy costs to guarantee the optimal microenvironment [34][35][36][37][38] . The research also shown that the viruses inducing aerobic glycolysis as known as the "Warburg effect" 32,33 . The alterations of nutrients uptake and utilization could provide ATP in a rapid fashion to meet the high energy costs and avoid cell ROS stress of virus replication, and the increasing of lipids by altering the fatty acid synthesis, is needed to afford the additional membrane material for the viral particle envelopment.
In comprehensive consideration of existing research and the data of transcription and metabolism, more detail about LdMNPV and host interaction was uncovered. The integrative analysis of metabolome and transcriptome totally enriched 52 correlational pathways in the POS mode and 56 correlational pathways in the NEG mode. The results showed that the speci c differences in NB3 versus CK3 (physiological status) and Ld3 versus Ld6 were covering more KEGG pathways with less enrichment amount of different express genes and metabolites, between which have more similarity than others compares ( Fig. 4A-C). The Global and overview maps, Carbohydrate metabolism and Amino acid metabolism present signi cant difference in both of POS and NEG mode (Fig. 4C), but the amount of difference expression genes and metabolites in Ld3 versus Ld6 were smaller than others, especially the CK3 versus Ld3 (Fig. 4D-E). Comparing with the normal physiological status, the alteration of infected status refers to signal transduction, nucleotide metabolism, fatty acid metabolism, amino acid metabolism and cofactors metabolism.
The gure.4F and table.2 were showed that Carbohydrate metabolism, Amino acid metabolism and Lipid metabolism were enriched to most pathways, Carbohydrate metabolism have 15 pathways, Amino acid metabolism have 14 pathways, and Lipid metabolism have 11 pathways. In metabolism relevant pathway the dynamic change of the Alanine, aspartate and glutamate metabolism pathway, Lysine degradation pathway and Tyrosine metabolism pathway in Amino acid metabolism, Glycerophospholipid metabolism pathway in Lipid metabolism, ABC transporters pathway in Translation and neuroactive ligand-receptor interaction pathway in Membrane transport were striking. It is showed that the speci c enriched KEGG pathways were arachidonic acid metabolism, propanoate metabolism, mTOR signaling pathway, parkinson disease and alcoholism in NB3 versus CK3; oxidative phosphorylation in NB3 versus Ld3; D-Glutamine and D-glutamate metabolism, nicotinate and nicotinamide metabolism, one carbon pool by folate and FoxO signaling pathway in Ld3 versus Ld6; and glycerolipid metabolismin CK3 versus Ld3. For another, the enrichments excepting metabolism relevant only covered 10 pathways, including aminoacyl-tRNA biosynthesis pathway, neuroactive ligand-receptor interaction pathway, ABC transporters pathway, phosphatidylinositol signaling system, FoxO signaling pathway, mTOR signaling pathway, protein digestion and absorption pathway, bile secretion pathway, parkinson disease pathway and alcoholism pathway. Except the aminoacyl-tRNA biosynthesis pathway, neuroactive ligand-receptor interaction pathway, ABC transporters pathway and phosphatidylinositol signaling system, the others just change a few in limited compares. The aminoacyl-tRNA biosynthesis pathway was violent alteration in CK3 versus Ld3, re ecting that the reprogram of host cell transcription preferences which leading by LdMNPV. There have not enrichment any difference in Ld3 versus Ld6, which means alteration of aminoacyl-tRNAs were keeping stable with infection processing. Moreover, it is reported that after infected by BmNPV in Bombyx mori larvae the trehalose, malate, fumarate, citrate, succinate, α-ketoglutarate, phenylalanine and tyrosine were accumulated in the resistant strain YeA, but in susceptible strain YeB the gluconate, glutamate and aspartate showed upregulation 25 . The integrative analysis showed in L. dispar larvae the phenylalanine, tyrosine, gluconate and glutamate were continued accumulation in CK3 versus Ld3 and Ld3 versus Ld6, which revealed different patterns of host-virus interaction. The gluconate and glutamate were gived further attention, in Ld3 versus Ld6, the accumulation glutamate indicate targeted utilization by glul gene upregulated leading glutamine showed increase and the glutamate decarboxylase gene was downregulated causing GABA, a chief inhibitory neurotransmitter, remains stable. In the insect central nervous system (CNS), GABA and glutamate working as antagonists and agonists in neuronal excitability control [39][40][41][42] . It is showed that GABA and glutamate are involved in the mediation of aggressive behavior in invertebrates 43 . In honeybee, the connection between scouting behaviors and catecholamine, glutamate and GABA re ect effects of neurotransmitter to behaviors, and glutamate treatments increased the probability of foragers scouting behaviors 42 . Moreover, in CK3 versus Ld3 the caffeine was increase. It is reasonable to consider that content patterns of caffeine, glutamate and GABA involved in LdMNPV inducing behavior change.
Metabolisms like spermidine, pyruvate, GluNAc and NAD + were relevanted to life prolonging [44][45][46][47] . Among them, only GluNAc showed increasing, which was enriched in Ld3 versus Ld6. In C.elegans and Drosophila melanogaster cyclohexanehexol was reported to extend lifespan 23,24,49 . The inositol phosphate metabolism pathway and phosphatidylinositol signaling system pathway, which linking by cyclohexanehexol, were both enriched after LdMNPV infection. In CK3 versus Ld3, the MI was increase and PI was decrease, and relevent genes, the inn5a and pi3k were upregulated, and pi3r gene was downregulated. And in Ld3 versus Ld6, the PI was increase and MI was decrease, and the pten/mmac1, minpp1 and ip3r gene were upregulated, the pi3k gene was downregulated. As important second messenger, the phosphatidylinositol-3,4,5-triphosphate (PIP3) mediated lots of intracellular signaling molecule which was activating types of growth factors, cellular stimulus or toxic insults and regulates fundamental functions affected transcription, cell cycle progression, proliferation and survival 50,51 , that keeping stable in all groups. The accumulation of PIP3 in cell will trigger stress resistant and apoptosis. Moreover, intracellular phospholipid signal involved in managing endocytic and exocytic processes, regulating iron channaels, pumps and transporters, coordinating vesicular tra cking and modulating lipid distribution and metabolism 52 . Numbers of pathogens, especially virus, has been recognized reprogram cellular processes by hijack host cell phospholipid signal 53 . Refer to alteraction of the inositol phosphate metabolism pathway and phosphatidylinositol signaling system pathway in infected larvae, it is hypothesis that MI, PI and relevant genes involving in larval extended instar duration and apoptosis avoiding.
Considering all the results, LdMNPV biphasic infection mode causing host manipulation changes continuously in different time points of infection. Comparing with the normal physiological status, alterations of nutrients uptake and utilization could provide ATP in a rapid fashion to meet the high energy costs and avoid cell stress resistant and apoptosis of virus replication, and the increasing of lipids by altering the fatty acid synthesis, is needed to afford the additional membrane material for the viral particle envelopment, and immune function relevant pathways were not signi cantly affected. Moreover, maintaining infected cell stability is important to virion production and assembly at terminal of infection and the extended instar duration of host was bene tted to virus production.
In summary, we explored the complicated interactions between LdMNPV and 3th-instar L. dispar larvae at the metabolome and transcriptome levels, conducted the preliminary analysis, uncovered that the content patterns of caffeine, glutamate and GABA involved in LdMNPV inducing behavior change, and the MI, PI and relevant genes involving in larval extended instar duration and apoptosis avoiding. It is anticipated that further research in relevent pathways will reveal in the complex reactions in con ict between L. dispar and LdMNPV.

Materials And Methods
Insect and virus L. dispar eggs were provided of Professor Liang-Jian Qu from the Chinese Academy of Forestry, Beijing, China. Hatched larvae were raised on arti cial diet in a sterile phytotron with conditions as follows: temperature 26 ± 1°C, Relative Humidity 60 ± 5% and the photoperiod of 14: 10 (Light: Dark).
LdMNPV was stored in the laboratory. The ODV s was obtained from inoculated larvae, and prepared by the following steps. Inoculated larvae were washed with 0.25% sodium dodecyl sulfate ltered through lter paper and cheesecloth. Filtrate is washed twice with doubledistilled water and centrifuged at 3000rpm for 30mins. The precipitate was re-suspended with double-distilled water and quanti ed using a bacterial counting chamber, and diluted to the nal concentration of 10 9 occlusion bodies (ODVs)/mL. The virus suspensions were stored at 4°C. Sample preparation Newly molted third-instar larvae were used as the experimental subjects and were randomly divided into two groups: infection group and control group. Larvae of the infection group were fed with 10 6 OBs per larva of LdMNPV. For the control group, the isopyknic doubledistilled water was provided using the same method. After having been fed, larvae were reared in a 24-well insect rearing box in individual wells with ample arti cial diet in each well. Samples for transcriptome and metabolome analysis were obtained at the newly-molted thirdinstar larvae (NB3) stage,at 72 hours after treatment for the control group (CK3), at 72 hours after treatment for the LdMNPV infected group (Ld3), and 144 hours after treatment for the LdMNPV infected group (Ld6). The control group did not have any samples at 144 hours after treatment because mock-treated larvae molted into the next instar before 144 hours. Meanwhile surplus larvae from each group were kept to con rm the treatment was working. Death for infection was committed by liquefaction, and another group was observed till pupation. Sampled larvae were snap frozen in liquid nitrogen more than 24 hours and stored at − 80°C until further steps. Transcriptome sequencing and data analysis Total RNA was extracted from the intact samples by using TRIzol Reagent (Invitrogen) following the protocol of the manufacturer. Quality of RNA-samples was checked on 1% agarose-gel and the concentration were measured using NanoDrop 2000 (Thermo). Sequencing libraries were generated using NEBNext®Ultra™ RNA Library Prep Kit for Illumina® (NEB, USA) and puri ed by AMPure XP system (Beckman Coulter, Beverly, USA). All cDNA libraries sequencing were performed at Biomarker Company (Beijing, China; http://www.biomarker.com.cn).
The raw data were processed through in-house perl scripts and the clean data Q20, Q30, and GC content were calculated. The clean data using for RNA de novo assembly by Trinity 54 which setting like min_kmer_cov as 2 and all others by default. Assembled transcript functions were annotated based on the NCBI non-redundant protein sequences (Nr), Protein family (Pfam), Swiss-Prot, Gene Ontology (GO), Clusters of Orthologous Groups of proteins (KOG/COG/eggNOG) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases using BLASTX to identify the proteins.
The DESeq R package (1.10.1) was used to provide statistical routines for determining differential expression. Genes with adjusting by the Benjamini and Hochberg approach, which P-value < 0.05 were assigned as differentially expressed. The KOBAS software was used to test the statistical enrichment of differential expression genes in KEGG pathways 55 . The Gene Set Enrichment Analysis (GSEA) was using to show the variability of DESeq clursters within KEGG pathway 56 between the combination as NB3 vs CK3, NB3 vs Ld3, NB3 vs Ld6, CK3 vs Ld3, CK3 vs Ld6, Ld3 vs Ld6. The results of GESA were visualized with Cytoscape.

Metabolite extraction and analysis
Each group has six biological replicates, and the samples were prepared with following steps. 400 µL of iced Methanol/Acetonitrile (1: 1, v: v)(rank of LC-MS, CNW Technologies) was added into 100µL sample in 4°C and vortexed for 30sec, at − 20°C incubated for 20min, centrifuged at 14,000 rpm for 15 min at 4°C. Supernatant was transferred to a new tube and dried in a vacuum concentrator than using 100 µL of Acetonitrile/Water (1: 1, v: v) (rank of LC-MS, CNW Technologies) to reconstitute. 60µL aliquot of the supernatant was used to analysis via UHPLC-QTOF-MS.
The LC-MS/MS analyses were performed using 1290 UHPLC system (Agilent Technologies) with ACQUITY UPLC BEH amide column (1.7µm, 2.1 mm × 100 mm, Waters) coupled to a Triple TOF 5600 (Q-TOF, AB Sciex) instrument 57 . Analyst TF 1.7 software (AB Sciex) which was used to control the AB5600 TripleTOF(AB Sciex) continuously evaluates the full scan-surveyed MS data. And the mobile phase consisted of A and B which was passed through the elution gradient as 0 min, 95%B; 0.5 min, 95%B; 7 min, 65%B; 8 min, 40%B; 9 min, 40%B; 9.1 min, 95%B; and 12 min, 95% B, which was delivered at 500µL/min, and the injection volume was 2µL. A was 25mM Ammonium acelate (rank of LC-MS, CNW Technologies) and 25mM Ammonium hydroxide (rank of LC-MS, CNW Technologies) in water (pH = 9.75), and B was Acetonitrile (rank of LC-MS, CNW Technologies). The electrospray ionization (ESI) source conditions were as follows: ion source gas 1 as 60psi; ion source gas 2 as 60; curtain gas as 30; source temperature of 550°C; and ion spray voltage oating (ISVF) at 5000 V or -4000 V in positive or negative mode, respectively. By using the SIMCA14 software package (Umetrics, Umea, Sweden), PCA (principal component analysis) and OPLS-DA (orthogonal projections to latent structures-discriminate analysis) were implemented by with the results including peak numbers, sample names, and normalized peak areas. OPLS-DA was used to identify signi cantly different metabolites between the two comparison groups. The VIP values (variable importance for the projection values) were obtained to clarify the analysis, which exceeding 1.0 was rst selected as different metabolites. Student's T test was used to evaluate the remaining variables. The variables between groups were discarded with P > 0.05. And each metabolite's FC value (fold change value) was calculated by comparing the mean value of the peak area. The KEGG Online databases were used to search relate pathways, and differential metabolites were cross-linked to KEGG pathways.
Integrative analysis of metabolome and transcriptome The PCC (Pearson correlation coe cients) and CCA (canonical correlation analysis) were used to assess the integration of the metabolome and transcriptome. Transcriptome and metabolome datasets were combined to create a general metabolite-transcript network centred on genes according to the KEGG identi cation. The results of KEGG identi cation were visualized with Cytoscape.