B. xylophilus was isolated from a wilted Masson pine in Guangdong province, China, purified using a single fertilized adult female, and cultured in the fungus Pestalotiopsis sp. grown on potato dextrose agar (PDA) medium at 25°C under 16 h light/8 h dark (16/8 LD) conditions for 6 days.

Approximately 3-year-old P. massoniana seedlings were obtained from Yangshan County (24°17′1.6″N, 112°33′42.4″E), Qingyuan City, Guangdong province, China, and grown at temperatures ranging from 25°C to 30°C. Nicotiana benthamiana LAB were grown in a glasshouse at 25°C with a relative humidity of 60% under 16/8 LD conditions.

B. xylophilus was inoculated into a small, approximately 1 cm, wound in P. massoniana stems. After 12 h, mixed-life-stage nematodes were isolated from the pine stems using a Baermann funnel technique (Hu et al., 2019), and approximately 20,000 nematodes were used for RNA extraction. The same number of B. xylophilus nematodes were also isolated from Pestalotiopsis sp. grown on PDA plates for RNA extraction. Each treatment consisted of two replicates. Total RNA was extracted using the RNAprep Pure Micro Kit (DP420, TianGen Biotech), and first-strand cDNA was synthesized from 1 µg of total RNA using a TransScript One-Step gDNA Removal and the cDNA Synthesis SuperMix kits (AT311, TransGen Biotech, Beijing, China) according to the manufacturer’s instruction.

Sequencing libraries were constructed using the NEBNext Ultra RNA Library Prep Kit (E7530, NEB, USA) according to the manufacturer’s instructions. RNA integrity was determined by the RNA Nano 6000 Assay Kit of the Bioanalyzer 5400 system (Agilent Technologies, CA, USA). Total RNA was used as input material for the RNA library preparations. The library preparations were sequenced on an Illumina Novaseq platform, generating 6.49-6.82 GB of data of 150-base long paired reads per sample (Modi et al., 2021). The raw RNA-Seq data in this study are available through the National Center for Biotechnology Information under the accession number PRJNA1027064.

Paired-end reads were aligned to the reference genome (https://parasite.wormbase.org/Bursaphelenchus_xylophilus_prjea64437/Info/Index) using Hisat2 v2.0.5 (Kim et al., 2019). The reads mapped to each gene were measured using FeatureCounts v1.5.0-p3 (Liao et al., 2014). The FPKM values of each gene were calculated based on the length of the gene and the reads count mapped to this gene (Trapnell et al., 2010).

Differential expression analysis was performed using the DESeq2 R package (1.20.0). The resulting P-values were adjusted using the Benjamini and Hochberg’s approach for controlling the false discovery rate (Benjamini and Hochberg, 1995). Genes identified by DESeq with an adjusted P-value<0.05 and |Log₂(fold change)|>1 were considered differentially expressed genes.

Real-time quantitative PCR (RT-qPCR) was performed using Green qPCR SuperMix (AQ101, TransGen Biotech, Beijing, China) on a Thermal Cycler Dice Real Time System (Takara, Beijing, China). The RT-qPCR reaction mixture consisted of 1.0 μL cDNA, 0.4 μL of each of the 10 pM primers, 10 μL 2 × Green qPCR SuperMix, and 8.2 μL RNase-free ddH₂O. The ubiquitin-conjugating enzyme gene BxUBI2 (GenBank accession number CAD5208953.1) and NbEF1α (Heese et al., 2007) were used as endogenous controls for gene expression normalization in the different species. All the gene-specific primers used in this study are listed in Supplementary Table 1 . The RT-qPCR cycling conditions were as follows: 94°C for 30 sec, followed by 40 cycles of 94°C for 10 s, 55°C annealing for 15 s, 72°C elongation for 10 s. The relative changes in gene expression were determined using the 2^−ΔΔCT method (Livak and Schmittgen, 2001). Three independent experiments were conducted.

Putative secreted proteins encoded by differentially expressed genes (DEGs) were predicted using both SignalP 5.0 (https://services.healthtech.dtu.dk/service.-php?SignalP-5.0) and TMHMM-2.0 (https://services.healthtech.dtu.dk/service.php?TMHMM-2.0) for the presence of an N-terminal signal peptide (SP) and for the absence of a transmembrane domain, respectively. These steps were implemented as a ‘secreted protein prediction’ workflow (Petitot et al., 2016). Protein functional domains were predicted by EGGNOG-MAPPER (http://eggnog-mapper.embl.de/). The phylogenetic tree was constructed using the NJ method and visualized using MEGA-X (Kumar et al., 2018).

The full-length cDNA sequences of 10 candidate genes were amplified by using the corresponding gene-specific primers. The PCR reaction was performed in a total volume of 50 μL containing 1.0 U KOD-FX (KFX-101, Toyobo, Osaka, Japan), 25 μL 2 × PCR Buffer, 10 μL 2 mM dNTPs, 1.5 μL of each of the 10 pM primers, 200 ng cDNA, and sterilized distilled water up to 50 μL. The PCR cycling conditions were as follows: 95°C for 5 min, followed by 30 cycles of 98°C for 10 s, annealing at 55°C for 30 s, elongation at 68°C for 30 s, and final elongation at 68°C for 5 min.

The CaMV35S driven β-glucuronidase (GUS) gene of the pCAMBIA1305.1 vector was removed and replaced with the eGFP gene to generate the binary vector pCambia1305.1^ΔGUS and pCambia1305.1:eGFP. To construct the overexpression vector, the 10 candidate genes mentioned above, with or without the native SP sequence, were amplified and ligated to the linearized pCambia1305.1:eGFP digested by NcoI and PmlI (FD0574 and FD0364, Thermo Fisher Scientific, MA, USA) to generate pCambia1305.1:gene:eGFP, following the instructions of the ClonExpress II One Step Cloning Kit (C112-01, Vazyme, Nanjing, China). BxICD1:Flag was amplified and ligated to the linearized pCambia1305.1^ΔGUS to generate pCambia1305.1:BxICD1:Flag. For pSUC2, vectors were digested by EcoRI and XhoI (FD0274 and FD0694, Thermo Fisher Scientific, MA, USA) in the appropriate conditions. The SP sequences of BxICD1 and Avr1b were linked to linearized pSUC2 to generate pSUC2:SP^BxICD1 and pSUC2:SP^Avr1b following the instructions of the ClonExpress II One Step Cloning Kit.

Approximately 10,000 mixed-life-stage nematodes were collected from pine stems. The primer pairs ISH-BxICD1-F/ISH-BxICD1-R were employed to synthesize digoxygenin (DIG)-labeled sense and antisense DNA probes based on the BxICD1 fragments of 96–281 bp using a PCR DIG Probe Synthesis Kit (1636090910, Roche Applied Science, Rotkreuz, Switzerland). The primer pair ISH-BxICD2-F/ISH-BxICD2-R was employed to synthesize sense and antisense DNA probes based on the BxICD2 fragments of 103-311 bp in size. The in situ hybridization and staining of the nematodes was performed as described previously (De Boer et al., 1998). Nematodes were then examined using a Nikon ECLIPSE Ni microscope (Nikon, Tokyo, Japan).

The cell death assay was performed as described previously (Chen et al., 2017). Briefly, the pCambia1305.1:gene:eGFP constructs mentioned above were individually introduced into Agrobacterium tumefaciens GV3101. Then the transformed bacteria were suspended in a buffer containing 10 mM 2-(N-morpholino) ethanesulfonic acid (MES) (pH 5.5) and 200 μM acetosyringone until a 600 nm absorbance (OD₆₀₀) of 0.6 was reached. Then, the A. tumefaciens cell suspensions were infiltrated into 4-week-old N. benthamiana leaves. After 3 days, the cell-death phenotype was observed and photographed. eGFP and the P. infestans elicitor gene INF1, cloned individually into the pCambia1305.1 vector, were used as the negative and positive control, respectively. The infiltration assay was performed three times, and five different plants with two inoculated leaves were used in each assay.

The secretion function of the BxICD1 SP was verified using the yeast signal sequence trap system described previously (Yin et al., 2018). Briefly, the SP sequence of BxICD1 was cloned into the pSUC2 vector containing the invertase gene but lacking Methionine (Met) and SP sequence to generate pSUC2:SP^BxICD1-Invertase. Additionally, the SP of Phytophthora sojae RXLR effector Avr1b has been demonstrated to have secretion function in yeast (Dou et al., 2008) and was also cloned into pSUC2 to generate pSUC2:SP^Avr1b-Invertase as a positive control. pSUC2:SP^BxICD1-Invertase and pSUC2:SP^Avr1b-Invertase were transformed into the yeast strain YTK12 by the lithium acetate method (Gietz and Schiestl, 2007; Yin et al., 2018). The yeast strains YTK12 and YTK12 carrying the pSUC2 empty vector were used as negative controls. YTK12 and YTK12 containing pSUC2-derived plasmids were grown on a CMD-W medium (0.67% yeast N base without amino acids, 0.075% W dropout supplement, 2% sucrose, 0.1% glucose, and 2% agar) to detect the expression of pSUC2-derived plasmids. The enzymatic activity of invertase was confirmed by the reduction of the dye 2,3,5-triphenyltetrazolium chloride (TTC) to the insoluble red-colored 1,3,5-triphenylformazan (TPF).

The subcellular localization assays were performed as described previously (Chen et al., 2017). Briefly, BxICD1 sequences with or without the native SP were cloned into the pCambia1305.1 vector to generate BxICD1:eGFP and BxICD1^Δsp:eGFP, respectively. eGFP alone was used as the control. The constructs were used for the transformation of 4-week-old N. benthamiana leaves by agroinfiltration. Then the N. benthamiana plants were cultured for 3 days at 25°C under 16/8 LD. N. benthamiana leaves expressing the constructs were infiltrated with 30% glycerol for 20 min to induce plasmolysis for the observation of extracellular localization (Chen et al., 2021). The fluorescence was observed using an SP5 Leica confocal microscope (Nikon, Tokyo, Japan).

Western blotting was performed to verify the production of intact BxICD1:eGFP and BxICD1^Δsp:eGFP fusion proteins in N. benthamiana, as described previously (Chen et al., 2021). Briefly, the total proteins from N. benthamiana leaves were extracted using RIPA lysis buffer. The proteins were denatured, separated on an SDS-PAGE gel, and transferred to a nitrocellulose membrane (PALL, Washington, NY, United States). After blocking with 5% (w/v) non-fat milk for 2 h at room temperature, the membranes were incubated with a primary mouse anti-GFP antibody (1:5000 dilution) (HT801, TransGen Biotech, Beijing, China) in a blocking solution for 2 h. Then membranes were incubated with an anti-mouse horseradish peroxidase-conjugated secondary antibody at a 1:5000 dilution (HS201, TransGen Biotech, Beijing, China). Proteins were visualized using the Immobilon Western Chemiluminescent system with Pierce ECL Western Blotting Substrate (Thermo Fisher Scientific, MA, USA). Ponceau S staining was used to assess equal loading (Goldman et al., 2016).

Electrolyte leakage assays were performed in N. benthamiana, as described previously (Yu et al., 2012). The A. tumefaciens harboring BxICD1:eGFP, BxICD1^Δsp:eGFP, eGFP and INF1 were infiltrated into N. benthamiana leaves as described above. 36 h after agroinfiltration, five leaf discs (9 mm diameter) were floated on 5 mL distilled water for each sample and shaken at 25°C for 3 h. Then the conductivity of the solution was measured with a conductivity meter (DDS-307A, Shanghai INESA Scientific Instrument) to obtain the value A. In addition, the conductivity of the solution containing the leaf discs was measured after boiling for 30 min to obtain the value B. Ion leakage was calculated according to the formula (value A/value B) × 100%. The experiment was performed three times. Western blotting was performed to detect protein expression as described above.

Virus-induced gene silencing targeting NbBAK1 and NbSOBIR1 in N. benthamiana was performed as described previously (Lin et al., 2016). Briefly, the vectors pTRV1, pTRV2:NbBAK1, pTRV2:NbSOBIR1, and pTRV2:eGFP were individually transformed into A. tumefaciens GV3101. A. tumefaciens GV3101 harboring pTRV1 was mixed in a 1:1 ratio with those harboring pTRV2:NbBAK1, pTRV2:NbSOBIR1 or pTRV2:eGFP, respectively, and infiltrated into N. benthamiana leaves. The silencing efficiency of NbBAK1 and NbSOBIR1 was validated by RT-qPCR. RT-qPCR was performed as described above. The infiltration experiment was performed thrice and three different plants, with three inoculated leaves in each, were used in each assay.

The BxICD1 gene was silenced using in vitro RNAi, as described previously (Liu et al., 2022b). Briefly, the T7-promoter sequence was introduced in the sense and antisense direction of a 186-bp BxICD1 fragment using the primer pairs RNAi-T7BxICD1-F/RNAi-BxICD1-R and RNAi-BxICD1-F/RNAi-T7BxICD1-R. The PCR was performed as described above. The obtained PCR products were used to synthesize double-stranded RNA (dsRNA) using the Thermo T7 Transcription Kit (TSK-101, Toyobo, Shanghai, China) according to the manufacturer’s instructions. The same approach was used to synthesize dsRNA for eGFP using the primer pairs RNAi-T7eGFP-F/RNAi-eGFP-R and RNAi-eGFP-F/RNAi-T7eGFP-R. Nematodes were then inoculated in a dsRNA and a non-dsRNA solution at 25°C for 24 h. The treated nematodes were washed with ddH₂O three times to remove external dsRNA. Subsequently, approximately 5,000 nematodes were collected for RNA extraction. The remaining nematodes were used for nematode virulence, migration, and reproduction analysis. The extent of BxICD1 silencing was assessed using RT-qPCR. RNA extraction and RT-qPCR were performed as described above, except for the different primer pairs that were used in the RT-qPCR.

Each 3-year-old P. massoniana seedling was inoculated with 1,000 B. xylophilus treated with a BxICD1 dsRNA, eGFP dsRNA, and non-dsRNA solution. Nine to ten pine trees were inoculated with each treatment. The experiment was repeated 3 times. Following the classification by Hu et al. (2022a), the morbidity degree of the P. massoniana seedlings was categorized into five different grades: 0, all needles were green; I, a few needles have turned yellow; II, approximately half of the needles have turned yellow or brown; III, most of the needles turned brown; and IV, the entire seedling was withered. The following formulas were used to calculate morbidity and the disease severity index.

The migration ability of B. xylophilus in P. massoniana was determined as described previously (Kusumoto et al., 2014). Briefly, approximately 800 nematodes treated with BxICD1 dsRNA, eGFP dsRNA, and non-dsRNA solutions were used to inoculate 3-year-old P. massoniana seedlings. 12 h after inoculation, the pine stem segments 1-2, 2-3, 3-4 and 4-5 cm below the inoculation point were used for nematode extraction using the Baermann funnel method. Nematodes isolated from these stem segments were counted. The experiment was repeated twice with four biological replicates each time. Statistically significant differences between treatments and controls were determined by Duncan’s multiple-range test.

The reproduction and feeding rate of B. xylophilus in fungi was assessed as described previously (Hu et al., 2022a). Approximately 100 nematodes treated with BxICD1 dsRNA, eGFP dsRNA, and non-dsRNA solutions were inoculated to PDA plates covered with Pestalotiopsis sp. After culturing at 25°C for 8 d, the nematodes were collected from each PDA plate using the Baermann funnel method and counted to calculate their reproduction. The amount of Pestalotiopsis sp. remaining on PDA plates was measured to determine the feeding rate of B. xylophilus. Six PDA plates were used for each treatment, and the experiment was repeated three times. Statistically significant differences between treatments and controls were determined by Duncan’s multiple-range test.

RNA was extracted from B. xylophilus during the mycetophagous and phytophagous stages to construct respective RNA-Seq libraries. In total, 181,144,358 raw reads were obtained from all samples. After filtering out low-quality sequences, 177,125,722 clean reads were retained. The GC content of individual libraries was between 48%-49%. On average, 92% of the clean reads mapped to the B. xylophilus reference genome ( Supplementary Table 2 ). The high genome coverage of our RNA-Seq data indicated that the transcriptome data were reliable for the further bioinformatics analyses.

Differential expression analysis showed that 1,233 genes were differentially expressed between the mycetophagous and phytophagous stages of B. xylophilus. Compared with the mycetophagous stage, 765 genes were significantly up-regulated, and 468 genes were significantly down-regulated in the phytophagous stage ( Figure 1A ). To validate the identified DEGs based on the RNA-seq data, seven genes were selected and validated by RT-qPCR. Similar expression levels were observed for all seven genes in both RNA-seq and RT-qPCR datasets, further confirming the reliability of the RNA-seq transcriptome data ( Figure 1B ).

Among the 765 up-regulated genes, 394 encoded proteins contained an SP and lacked transmembrane domains, as predicted by SignalP 5.0 and TMHMM-2.0. Of these, 210 genes had no functional annotations and were not identified in the transcriptome analysis of B. xylophilus infecting Pinus thunbergii. 31 of these 210 genes were further considered to be candidate genes encoding for effectors ( Supplementary Table 3 ) based on the following criteria: Log₂(fold change) > 3 and FPKM > 10 with regards to their differential expression, and a gene open reading frame (ORF) is smaller than 600 bp. From this set, 10 genes were randomly selected and were transiently expressed in N. benthamiana. The results showed that the genes BXY_0304900 and BXY_1398900 can induce cell death (ICD) in the presence of the native SP. They were named BxICD1 and BxICD2, respectively ( Figure 2 ).

The nematode tissue localization of BxICD1 and BxICD2 in the nematodes was further determined by in situ hybridization. Hybridization signals were detected in the esophageal gland cells with the DIG-labeled antisense cDNA probe of BxICD1 and in the intestinal terminus with the antisense cDNA probe of BxICD2. No hybridization signals were detected in B. xylophilus tissues with the sense cDNA probe of BxICD1 and BxICD2 ( Figure 3A , Supplementary Figure 1 ). The results indicated that BxICD1 should be a nematode effector because it was expressed in the typical secretory organs of nematodes. Therefore, we selected the effector BxICD1 for further analysis.

The BxICD1 gene has a 396-bp ORF, which encodes a 131-amino acid protein. The protein possesses a 19-amino acid SP, as predicted by the SignalP program, and has no putative transmembrane domain according to TMHMM. It also shows similarities to the proteins of B. xylophilus by a protein BLAST search, including CAD5219553 (62.3% identity), CAD5219557 (44.6% identity), CAD5219551 (40.2% identity) etc., while no sequences were matched to other organisms. A neighbor-joining (NJ) tree was constructed to examine the phylogenetic relationships among these proteins, with BxICD1 and CAD5219553 being located in the same branch ( Supplementary Figure 2B ).

To determine whether the SP of BxICD1 is functional, the predicted SP was cloned into the yeast vector pSUC2. The growth of yeast YTK12 yeast carrying pSUC2, pSUC2:SP^BxICD1 and pSUC2:SP^Avr1b constructs on CMD-W demonstrated the successful transformation with the plasmids. The yeast strain YTK12 carrying pSUC2:SP^BxICD1-Invertase or the positive control pSUC2:SP^Avr1b-Invertase could reduce TTC to the red-colored TPF. By contrast, no color change was observed in the YTK12 strain used as a negative control, and YTK12 carried the empty pSUC2 vector ( Figure 3B ). These results suggested that BxICD1 carries a functional secretory SP.

Based on the above results, BxICD1 could induce PCD in the presence of the native SP, therefore the subcellular localization of BxICD1 was determined. Three constructs, BxICD1:eGFP, BxICD1^Δsp:eGFP and eGFP were generated and transiently expressed in N. benthamiana leaves. Based on the results, the BxICD1:eGFP fusion protein was accumulated in the edges of the tobacco cell, while BxICD1^Δsp:eGFP and eGFP were localized in the cytoplasm ( Figure 4A ). Plasmolysis was further performed to distinguish the plasma membrane from the apoplast. After the plasmolysis treatment, BxICD1:eGFP fusion proteins were observed in the extracellular space, while BxICD1^Δsp:eGFP and eGFP were still localized in the intracellular space ( Figure 4B ). Western blot analysis using a GFP antibody revealed bands at approximately 42, 40 and 27 kDa, corresponding to BxICD1:eGFP, BxICD1^Δsp:eGFP, and eGFP, respectively ( Figure 4C ), indicating that BxICD1:eGFP and BxICD1^Δsp:eGFP fusion proteins were expressed intact. At the same time, BxICD1:eGFP could induce tobacco cell death, while BxICD1^Δsp:eGFP and eGFP could not ( Figure 4D ). In addition, the electrolyte leakage induced by INF1 or BxICD1:eGFP was significantly greater than that induced by BxICD1^Δsp:eGFP or eGFP ( Figure 4E ). These results suggested that the presence of BxICD1 in N. benthamiana apoplasts triggers cell death.

Many studies have demonstrated that BAK1 and SOBIR1 are indispensable for cell death triggered by many apoplastic effectors (Liebrand et al., 2014). To determine whether NbBAK1 and NbSOBIR1 are involved in BxICD1-triggered cell death, they were knocked down using virus-induced gene silencing (VIGS) in N. benthamiana. Two weeks after viral inoculation, N. benthamiana plants were agroinfiltrated with the pCAMBIA1305.1:BxICD1:Flag, pCAMBIA1305.1:INF1, and pCAMBIA1305.1:eGFP : Flag-expressing constructs, respectively. Based on the results, INF1 failed to induce cell death in NbBAK1- or NbSOBIR1-silenced plants. Notably, although BxICD1 was unable to induce cell death in NbBAK1-silenced plants, it retained the capacity to induce cell death in NbSOBIR1-silenced plants ( Figure 5A ). Neither eGFP nor the EV, used as negative controls, triggered cell death in pTRV2:NbBAK1-, pTRV2:NbSOBIR1- and pTRV2:eGFP- infiltrated plants ( Figure 5A ). Western blot analysis was used to confirm BxICD1 protein expression in NbBAK1-, NbSOBIR1- and eGFP-silenced plants ( Figure 5B ), and RT-qPCR confirmed the successful silencing of NbBAK1 and NbSOBIR1 in plants ( Figures 5C, D ). Taken together, these results indicated that BxICD1-triggered cell death is dependent on NbBAK1.

To assess the role of BxICD1 in nematode parasitism, an RNAi assay was performed by soaking B. xylophilus in BxICD1 dsRNA. The transcript levels of BxICD1 were significantly decreased after incubation with the BxICD1 dsRNA, compared to the treatments with eGFP dsRNA and non-dsRNA solution ( Figure 6A ), demonstrating that RNAi was successful.

Subsequently, inoculation experiments were performed using B. xylophilus that underwent different pre-treatment. Compared with the eGFP dsRNA-treated and H₂O-treated nematodes at 12 dpi, BxICD1 dsRNA-treated nematodes caused to a significantly lower morbidity and disease severity index (DSI) in P. massoniana. The morbidity was calculated at 36.7%, 34.4%, and 13.7% after inoculation with eGFP dsRNA-, H₂O- and BxICD1 dsRNA-treated nematodes, respectively ( Figure 6B ), while corresponding DSI was 11.7, 9.5, and 3.4 ( Figure 6C ). At 20 dpi, the DSI of P. massoniana inoculated with BxICD1 dsRNA-treated nematodes was still significantly lower compared with P. massoniana inoculated with eGFP dsRNA and non-dsRNA solution-treated nematodes. However, no significant differences in morbidity were observed ( Figures 6C, D ). Moreover, at 12 h post-inoculation, B. xylophilus was isolated from pine stem segments at different distances from the inoculation site. No significant differences were observed in the number of BxICD1 dsRNA, eGFP dsRNA, and the non-dsRNA solution-treated nematodes at pine stem segments 1-2 and 2-3 cm below the inoculation site. However, the number of BxICD1 dsRNA-treated nematodes was significantly lower compared with eGFP dsRNA and non-dsRNA solution-treated nematodes isolated from pine stem segments 3-4 cm below the inoculation site. No nematodes were isolated from pine stem segments 4-5 cm below the inoculation point in all treatments ( Figures 6E, F ).

In addition, BxICD1 dsRNA, eGFP dsRNA, and the non-dsRNA solution-treated nematodes were cultured on a PDA medium covered with Pestalotiopsis sp. After eight days, no significant differences were observed in the amount of nematodes and fungi mycelium in different treatments ( Supplementary Figure 3 ). The above results showed that BxICD1 affected B. xylophilus parasitism and migration in P. massoniana, but did not affect the feeding rate and reproduction of B. xylophilus in fungi.