Five-year-old P. sibirica seedlings used in the experiment were cultivated in the greenhouse of Northeast Forestry University, maintained at 23–28°C and a relative humidity of 65%–75%. PWNs (collected from Fushun City, Liaoning Province, China) were cultured in the dark at 25°C using Botrytis cinerea, which was cultured on potato dextrose agar medium. Nematodes (males: females: juveniles = 1:1:2) were isolated and collected using Baermann funnels to prepare a nematode suspension (3000 nematodes/100 μl). Following the method described in our previous research (Chen et al., 2021), PWNs were inoculated into P. sibirica seedlings (100 μl per plant) as the treatment group and isovolumes of ddH₂O were inoculated into P. sibirica seedlings as the control check (CK) group.

All the treated and CK P. sibirica seedlings were divided into three groups. For the first group, five treated P. sibirica seedlings were randomly selected every day, and nematodes in the plants were collected using a Baermann funnel, then changes in the nematode population were calculated. For the second group, five treated and five CK P. sibirica seedlings were randomly selected every day, and stem segments were cut for plant tissue observation. For the third group, five treated and five CK P. sibirica seedlings were randomly selected every day, and 3 cm long stem segments were extracted from 1 cm below the inoculation site of each seedling, then all the segments were promptly placed in liquid nitrogen and maintained at −80°C for transcriptome sequencing and metabolome analysis.

Total RNA extraction from both inoculated and control tree samples was conducted utilizing the RN38 EASYSpin Plus Plant RNA Kit (Aidlab Biotech, Beijing, China). Quantification and assessment of purity for total RNA were performed using the Bioanalyzer 2100 and RNA 1000 Nano LabChip Kit (Agilent, CA, USA), ensuring a RIN number > 7.0. Poly(A) RNA was isolated from total RNA (5 μg) through two rounds of purification using Poly-T oligomer magnetic beads. The purified mRNA underwent decomposition into smaller fragments under high-temperature conditions in the presence of divalent cations. Then, according to the mRNASeq sample preparation kit (Illumina, San Diego, CA, USA), the cleaved RNA fragments were reverse-transcribed to obtain the final cDNA library. The average insertion length of the paired library was 300 bp ( ± 50 bp). Sequencing was carried out using Illumina Novaseq™ 6000 (LC Sciences, Houston, TX, USA) following the recommended procedures for double c-terminal sequencing.

To ensure the acquisition of high-quality reads, the in-house Cutadapt (Kechin et al., 2017) and perl (https://registry.hub.docker.com/_/perl/) scripts were employed to eliminate adapter contamination, low-quality base, and reads with undetermined base. The quality of sequences was subsequently assessed using FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/), encompassing Q20, Q30, and GC content of the clean data. All downstream analyses were exclusively conducted based on this high-quality clean data. De novo transcriptome assembly was carried out using Trinity 2.4.0 (Grabherr et al., 2011). In Trinity, transcripts were grouped into clusters based on shared sequence content, often loosely referred to as “genes”. The longest transcript within each cluster was designated as the representative “gene” sequence, also referred to as a single gene.

All assembled unigenes were associated with the non-redundant (Nr) protein database (http://www.ncbi.nlm.nih.gov/), the Gene Ontology (GO, http://www.geneontology.org), the SwissProt (http://www.expasy.ch/sprot/), the Kyoto Encyclopedia of Genes and Genomes (KEGG, http://www.genome.jp/kegg/) and the eggnog (http://eggnogdb.embl.de/) databases through DIAMOND (Buchfink et al., 2015), employing a threshold of E value < 0.00001.

Expression levels for unigenes were quantified using Salmon (Patro et al., 2017) by calculating Transcripts Per Kilobase of exon model per Million mapped reads (TPM) (Mortazavi et al., 2008). The differentially expressed Unigenes were selected with log₂ (fold change) > 1 or log₂ (fold change) < -1 and with statistical significance (P value < 0.05) by R package edgeR (Robinson et al., 2010).

Advanced Mfuzz analysis was performed using the OmicStudio tools at https://www.omicstudio.cn/tool.

After thawing the collected samples on ice, metabolites were isolated using an 80% methanol buffer. The Liquid Chromatograph Mass Spectrometer (LC-MS) system acquired all samples, following machine orders. Initially, all chromatographic separations were executed utilizing the UltiMate 3000 UPLC System (Thermo Scientific, Bremen, Germany). Reversed-phase separation occurred on an ACQUITY UPLC T3 column (100 mm × 2.1 mm, 1.8 μm, Waters, Milford, USA). The detection of metabolites eluted from the column was accomplished using the high-resolution tandem mass spectrometer Q-Exactive (Thermo Scientific, Bremen, Germany).

The MS data underwent various pretreatments, encompassing peak picking, peak grouping, retention time correction, second peak grouping, and annotation of isotopes and adducts utilizing XCMS (http://metlin.scripps.edu/download/) (Smith et al., 2006). LC-MS raw data files were transformed into mzXML format and processed through the XCMS, CAMERA (Kuhl et al., 2012), and metaX (http://metax.genomics.cn/) (Wen et al., 2017) toolbox integrated with the R software. The online KEGG and Human Metabolome Database (HMBD, https://hmdb.ca) were utilized for metabolite annotation by aligning the exact molecular mass data (m/z) of samples with those presented in the database. Additionally, an in-house fragment spectrum library of metabolites was employed to validate metabolite identification.

Student t-tests were executed to identify variations in metabolite concentrations between 2 phenotypes. The P value was adjusted utilizing a false discovery rate (FDR) for multiple tests. Supervised partial least squares-discriminant analysis (PLS-DA) was executed utilizing metaX to discern distinct variables between groups. The variable important for the projection (VIP) value was computed. A VIP cut-off value of 1.0 was utilized for the selection of the significant features. Using the criteria of a fold change ≥ 1.5 or ≤ 1/1.5, P value ≤ 0.05, and VIP ≥ 1 to identify differentially abundant metabolites.

RT-qPCR was performed with the GoTaq 2-Step RT-qPCR System Kit (Promega, Madison, WI, USA, catalogue number: A6010) and the Stratagene Mx3000P qPCR system (Agilent Technologies, Santa Clara, CA, USA) to validate the transcript levels of the genes (Chen et al., 2021). Actin was used as the internal control. All primers used in this study are listed in Supplementary Table S1 . The normalization of the data was performed according to the instructions for the GoTaq 2-Step RT-qPCR System Kit and by the 2^−ΔΔC _T method (Livak and Schmittgen, 2001). The experiments were repeated three times. Significance was determined by Student’s t-test.

Symptoms in P. sibirica seedlings inoculated with PWNs became apparent at 21 days post inoculation (dpi), with approximately 25% of the needles displaying chlorosis. All inoculated P. sibirica seedlings showed symptoms by 27 dpi, with three-quarters of the needles exhibiting substantial yellowing. At 79 dpi, the inoculated P. sibirica seedlings had withered and died. The duration from the onset of symptoms to complete wilting and death spanned 58 days. Conversely, the CK of P. sibirica seedlings displayed no notable changes throughout the entire observation period ( Figure 1A ).

Observations from plant tissue sections of P. sibirica seedlings revealed that stem segments near the point of nematode inoculation initially became desiccated and shriveled. Symptoms of wilting appeared 1 cm below the inoculation site at 9 dpi ( Figure 1B ). At this stage, there were no apparent external changes in the P. sibirica seedlings ( Supplementary Figure S1 ). Thus, internal symptoms in P. sibirica infected by PWNs appeared more quickly than external symptoms. Examination of the nematode population within the inoculated P. sibirica seedlings showed a decrease in nematode numbers at 5 dpi, possibly due to the mass mortality of adult nematodes following completion of their life cycle. Subsequently, nematode numbers gradually increased at 7 dpi, indicating the hatching of new larvae and signifying that PWNs could grow and reproduce normally within P. sibirica seedlings ( Figure 1C ).

Based on the onset of internal stem segment symptoms occurring 9 days after PWN inoculation in P. sibirica seedlings, an analysis was conducted to examine changes in gene expression from nematode inoculation to the onset of symptoms. Samples were selected corresponding to time points when the change in nematode population trend shifted (3 dpi, 5 dpi, 7 dpi, 9 dpi), along with control samples for each time point, leading to a total of eight samples. For each selected time point, five randomly collected pine seedlings were mixed in equal amounts and then subjected to transcriptome sequencing. An average of 5.8425 GB of valid data was obtained from each sample ( Supplementary Table S2 ). Upon assembling all samples, 136,125 transcripts, and 71,938 genes were acquired. These sequences had average lengths of 558 bp and 516 bp, N50 values of 1,489 bp and 1,526 bp, and GC ratios of 43.97% and 44.77%, respectively. All genes were matched against six databases, including KEGG, eggNOG, Swiss-Prot, Pfam, Nr, and GO, utilizing BLASTX to obtain comprehensive genetic information ( Supplementary Table S3 ). Functional annotations were acquired for 41,287 genes with respect to GO (57.39%), 27,097 genes for KEGG (37.67%), 37,611 genes for Pfam (52.28%), 34,353 genes for Swiss-Prot (47.75%), 45,181 genes for eggNOG (62.81%), and 41,914 genes for NR protein database (58.26%).

In order to explore the dynamic changes of genes in P. sibirica seedlings after PWN infection, the expression trends of all genes were analyzed and all genes were assigned to 12 clusters ( Figure 2A ; Supplementary Table S4 ). Among them, cluster 10 had the largest number of genes (8552 genes), and cluster 8 had the smallest number of genes (4611 genes). There was no obvious rule in the changes of gene expression patterns in different clusters. KEGG enrichment analysis was then performed on the genes in each cluster to further explore the differences in gene functions and genes significantly enriched in plant defense-related pathways were screened in each cluster ( Figure 2B ; Supplementary Tables S5-16 ). The results showed that all clusters had genes significantly enriched in plant defense-related pathways with high rich factor, indicating that genes related to plant defense-related pathways in P. sibirica seedlings might continue to respond after PWN infection.

Gene expression trend analysis showed that genes related to plant defense-related pathways in P. sibirica seedlings exhibited different changing trends after PWN infection. To determine whether the differential expression changes of these genes were significant and to screen for genes associated with the PWN inoculation in P. sibirica seedlings, genes differentially expressed after inoculation were analyzed. The results showed that there were 4679 genes exhibiting differential expression at 3 dpi, with 2531 genes up-regulated and 2148 genes down-regulated. At 5 dpi, 6962 genes showed differential expression, with 3280 genes up-regulated and 3682 genes down-regulated. At 7 dpi, 6098 genes exhibited differential expression, including 3068 genes up-regulated and 3030 genes down-regulated. At 9 dpi, 8100 genes showed differential expression, with 3981 genes up-regulated and 4119 genes down-regulated ( Figure 3A ; Supplementary Table S17 ).

KEGG enrichment analysis was then performed on differentially expressed genes at different time points. Each time point had differentially expressed genes enriched in multiple plant defense-related pathways ( Supplementary Tables S18-25 ). However, it cannot be ignored that each time point had a large number of up-regulated genes enriched in flavonoid biosynthesis (ko00941) with high rich factor ( Figure 3B ), indicating that flavonoid biosynthesis continued to respond in P. sibirica seedlings after PWN inoculation. And each time point had a large number of down-regulated genes enriched in photosynthesis (ko00195), photosynthesis-antenna proteins (ko00196) and carotenoid biosynthesis (ko00906) with high rich factors ( Figure 3C ), indicating that photosynthesis continued to be inhibited in P. sibirica seedlings after PWN inoculation.

The KEGG enrichment results of differentially expressed genes indicated that the flavonoid biosynthesis pathway might respond quickly and sustainably in P. sibirica seedlings after PWN inoculation. Therefore, the expression patterns of up-regulated genes that enriched in the flavonoid biosynthesis pathway at each time point were further analyzed to explore whether these genes were associated with the inoculation time. After removing duplicate genes, a total of 63 genes enriched in the flavonoid biosynthesis pathway were found to be up-regulated at any time point ( Supplementary Table S26 ). Expression pattern analysis showed that most of these genes were up-regulated at multiple time points, and only a few were down-regulated at some time points ( Figure 3D ). Among them, 16 genes were significantly up-regulated at all time points ( Figures 3E, F ; Supplementary Table S27 ). Ten genes were randomly selected from these 16 genes for RT-qPCR verification, including three chalcone-related genes (No. 11, 19 and 20), two caffeoyl-CoA O-methyltransferase genes (No. 6 and 10), a pinosylvin synthase gene (No. 7), an anthocyanidin reductase gene (No. 13), a flavonol synthase gene (No. 16), a shikimate O-hydroxycinnamoyltransferase gene (No. 4), and a bifunctional dihydroflavonol 4-reductase/flavanone 4-reductase gene (No. 21). The trends of RT-qPCR results and RNA-seq results of the 10 genes were consistent ( Supplementary Figure S2 ), indicating the reliability of the transcriptome sequencing and the importance of these genes in P. sibirica seedlings after PWN inoculation.

In order to further explore the role of flavonoid biosynthesis in P. sibirica seedlings after PWN inoculation, non-targeted metabolomic analysis was executed to identify the differentially abundant metabolites at 3, 5, 7, and 9 dpi. A total of 1937 metabolites were identified and classified into 15 categories, with lipids and lipid molecules being the predominant metabolites ( Figure 4A ; Supplementary Table S28 ). For a comprehensive understanding of the metabolic changes during the infection by PWNs at the four different time points, a quantitative analysis was performed on the metabolome data. Principal component analysis showed that, except for at 9 dpi, the other samples could not be effectively distinguished from the control sample ( Figure 4B ). At 9 dpi, the samples intersected with those at 7 dpi and 5 dpi while remaining distinct from 3 dpi ( Figure 4B ). This indicated that only the metabolome data at 9 dpi had a high degree of differentiation from the control. Furthermore, the metabolome data at 9 dpi exhibited a higher degree of differentiation from the metabolome data at 3 dpi. These results indicated that the metabolic impact on P. sibirica seedlings was relatively minor at 3 dpi. In contrast, the impact of PWNs on the metabolism of P. sibirica seedlings was notably significant during the early stages of disease at 9 dpi.

Subsequently, using the criteria of the methods, differentially abundant metabolites were identified when comparing 3, 5, 7, and 9 dpi with the control. A comprehensive analysis of these differentially abundant metabolites was conducted across the different comparisons ( Figure 4C ; Supplementary Table S29 ). The statistical examination revealed 347 significantly differentially abundant metabolites between 3 dpi and the control, encompassing 138 upregulated and 209 downregulated metabolites. When comparing metabolites from 5 dpi sample with the control, 210 significantly differentially abundant metabolites were observed, including 116 increased and 94 decreased metabolites. Approximately 237 significantly differentially abundant metabolites were observed between metabolites from 7 dpi sample and the control, with 142 increased and 95 decreased metabolites. Finally, the comparison of metabolites from 9 dpi sample with the control unveiled 266 differentially abundant metabolites, comprising 168 increased and 98 decreased metabolites.

Following this, KEGG enrichment analysis was executed on the differentially abundant metabolites within these comparative groups. Subsequently, all identified differentially abundant metabolites were assigned to distinct metabolic pathways ( Supplementary Figure S3 ). At 3 dpi, the pathways implicated included glycerolipid metabolism (ko00561), pantothenate and CoA biosynthesis (ko00770), glycosylphosphatidylinositol (GPI)-anchor biosynthesis (ko00563), plant hormone signal transduction (ko04075), and flavone and flavonol biosynthesis (ko00944) ( Supplementary Figure S3 , Supplementary Table S30 ). At 5 dpi, the enriched pathways for differentially abundant metabolites included glycosylphosphatidylinositol (GPI)-anchor biosynthesis (ko00563), flavonoid biosynthesis, biosynthesis of phenylpropanoids (ko01061), alpha-linolenic acid metabolism (ko00592), and galactose metabolism (ko00052) ( Supplementary Figure S3 , Supplementary Table S31 ). At 7 dpi, the observed differentially abundant metabolites were allocated to ABC transporters (ko02010), phenylpropanoid biosynthesis (ko00940), alpha-linolenic acid metabolism (ko00592), arginine biosynthesis (ko00220), and biosynthesis of phenylpropanoids ( Supplementary Figure S3 , Supplementary Table S32 ). Finally, at 9 dpi, the differentially abundant metabolites were linked to flavonoid biosynthesis, ABC transporters, arginine biosynthesis, biosynthesis of phenylpropanoids, and alanine, aspartate and glutamate metabolism (ko00250) ( Supplementary Figure S3 , Supplementary Table S33 ).

The transcriptome analysis results showed that differentially expressed genes significantly enriched in flavonoid biosynthesis, flavone and flavonol biosynthesis, and phenylpropanoid biosynthesis pathways, suggesting that flavonoid-related pathways might play an important role in P. sibirica seedlings after PWN inoculation. Similarly, as differentially abundant metabolites were significantly enriched in flavonoid biosynthesis at 5 dpi and 9 dpi samples, in flavone and flavonol biosynthesis at 3 dpi, and in phenylpropanoid biosynthesis at 7 dpi, differentially abundant metabolites enriched in flavonoid-related pathways were selected at various time points for further research ( Figure 4D ). Upon screening, it became evident that all differentially abundant flavonoid-related metabolites were increased. Specifically, the increased flavonoid metabolites reached the highest count of 9 at 9 dpi, followed by 5 dpi and 3 dpi with counts of 5, and the lowest count of 3 at 7 dpi. Notably, pinocembrin was increased at 3, 5, 7, and 9 dpi, while naringenin was increased at 5, 7, and 9 dpi, and apigenin was increased at 7 and 9 dpi. In addition, the levels of these three metabolites steadily increased each day after inoculation. It was suggested that the sustained significant increases of these metabolites in the late stages of PWN inoculation should have contributed to the resistance of P. sibirica seedlings to PWN infection.

By analyzing the synthesis pathways of these three metabolites, it became apparent that naringenin and pinocembrin shared consistent prerequisites and enzymes required for their synthesis, suggesting a competitive relationship between them. In contrast, apigenin was identified as a downstream metabolite of naringenin ( Figure 5A ). Through RT-qPCR and metabolome analysis, the three genes [two chalcone synthase (chs) genes (No. 19 and 20) and a chalcone isomerase (chi) gene (No.11)] and three metabolites (naringenin, pinocembrin and apigenin) screened out in the pathway map were analyzed. Analysis of expression levels showed that the expression of chs1, chs2 and chi continued to increase in P. sibirica seedlings after PWN inoculation, and with the up-regulation of the three genes, the contents of naringenin, pinocembrin and apigenin also increased significantly ( Figure 5B ). Next, we conducted a correlation analysis between expression levels of chs1, chs2 and chi, and abundance of naringenin, pinocembrin and apigenin, with the inoculation time and the number of nematodes in the P. sibirica seedlings after PWN inoculation. The results showed that the three genes and the three metabolites were positively correlated with the inoculation time and the number of nematodes in the seedlings ( Figure 5C ).