Whole-organism Transcriptomic Insight Into the Effects of Cytco, A Novel Nematotoxic Protein on the Pine Wood Nematode Bursaphelenchus Xylophilus

Background: The pine wood nematode Bursaphelenchus xylophilus is a worldwide destructive pest on Pinus trees and lacks effective control measures. Screening nematotoxic protein toxins has been conducted to develop new strategies for nematode control. Results: The present study provided initial insights into the responses of B. xylophilus exposed to a nematocidal cytolytic-like protein (CytCo) based on the transcriptome proling. A large set of differentially expressed genes (1266 DEGs) were found related to nematode development, reproduction, metabolism, motion, and immune system. In response to the toxic protein, B. xylophilus upregulated DEGs encoding cuticle collagens, transporters, and cytochrome P450. In addition, many DEGs related to cell death, lipid metabolism, major sperm proteins, proteinases/peptidases, phosphatases, kinases, virulence factors, and transthyretin-like proteins were downregulated. And Gene Ontology enrichment analysis showed that CytCo treatment signicantly affecting DEGs functioning in muscle contraction, lipid localization, MAPK cascade. The pathway richness of Kyoto Encyclopedia of Genes and Genomes showed that the DEGs were concentrated in lysosome and fatty acid degradation. The weight co-expression network analysis indicated that the hub genes affected by CytCo were associated with the nematode cuticular collagen. Conclusions: These results showed that the CytCo protein toxin could interference gene expression to produce multiple nematotoxic effects, providing initial insight into its control potential of pine wood nematode.


Background
The pine wood nematode (PWN) Bursaphelenchus xylophilus (Steiner & Buhrer) Nickle (Tylenchida: Aphenlenchoididae), a serious invasive species and the main cause of pine wilt disease (PWD), is listed as a quarantine pest in the legislation of more than 40 countries [1,2]. It has caused severe disaster to coniferous forest ecosystems with timber losses of 10 6 m 3 annually [2,3]. PWN spreads widely by beetle vectors (Monochamus spp.), invades pine trees by secreting various virulence factors such as expansinlike and venom allergen-like proteins, explosively proliferates, and ultimately kills host plants [2,4,5].
Several control measures for PWD are available, including fumigation with methyl bromide for the phytosanitary treatment of log exports, removing and burning dead wood in the infected areas, trunk injection of nematocidal compounds (e.g., emamectin benzoate and abamectin), monitoring and controlling PWN vectors, and breeding of resistant trees [2,6−8].
Pore-forming toxins (PFTs), speci cally the crystal proteins (Cry) and cytolytic (Cyt) δ-endotoxins, have been widely applied in insect pest control [9−11]. Recently, several PFTs, such as Bacillus thuringiensis crystal proteins Cry5B, Cry6A, Cry1E, and Cry 55A, were found to have many nematocidal characteristics in bioassays, holding potential of developing new strategies for nematode control [12−15]. For example, the Cry6Aa2 toxin has been found to suppress the hatching of the root-knot nematode Meloidogyne hapla eggs, with hatch rates being reduced by 29.82-97.37%, and also inhibits motility and penetration into the host plant [12]. However, the large molecular masses and dimensions of Cry proteins decrease their control e cacy on plant-parasitic nematodes. A Cyt-like protein from an entomopathogenic fungus Conidiobolus obscurus (named as CytCo) has been found to hold high nematicide activity against PWN and an inhibitory effect on egg hatching, and also signi cantly suppresses the reproductive and thrashing behaviors of PWN in bioassays [16]. CytCo (less than 22 kDa) characterized with a single domain of the center of β-strands being wrapped with a layer of α-helices can be easily taken up by PWN, presenting a potential of nematode control [16,17].
The modes of action of PFTs are attributed to their cytotoxicity, causing cell lysis by toxin oligomers assembling pores in cell membranes or eliciting lethal reactions in cells through signal transduction, metabolism and immune responses [10,18,19]. Speci cally, in Caenorhabditis elegans, Cry6Aa was found to trigger cell death after binding to the receptor of a CUB-like-domain containing protein RBT-1, and Cry5Ba was found to use an invertebrate-speci c glycolipid as its receptors for triggering cell lysis [20−22]. The different structures of these two Cry toxins contribute to their different modes of action [13].
The structure of a single α/β domain of most Cyt-like proteins is unlike that of the three-domain Cry toxins [17,23]. This unique molecular architecture implies that CytCo utilize a different mode of action on nematode pests. The present study aims to demonstrate PWN responses to CytCo by using transcriptomic pro ling in order to understand its potential mechanism against nematodes.

Results
General features of B. xylophilus transcriptome after treatment with CytCo Following quality assessment and data ltering, the resultant transcriptome contained 358,738,780 clean reads (Table S1). In total, 18,034 genes were predicted by mapping to the PWN reference genome, and 1266 DEGs (fold change ≥2, 379 upregulated and 886 downregulated) were ltered out between the RNAseq libraries from the PBS and CytCo treatments. After consolidating similar sequences, 726 (57.3%) DEGs were annotated in Wormbase and 541 (42.7%) DEGs were annotated in Swiss-Prot. One hundred and ninety-six (15.5%) DEGs were subjected to GO classi cation and 163 (12.9%) DEGs were annotated to 91 pathways according to KEGG.
The 196 GO-annotated DEGs were divided into 38 classes (level 2 subcategories) in the three ontologies of molecular function (18 classes), cellular component (7 classes), and biological process (13 classes).
As shown in Fig. 1, the largest class of DEGs was single-organism process (50 DEGs upregulated and 63 downregulated), followed by the classes of developmental process (36 upregulated and 52 downregulated) and cellular process (24 upregulated and 57 downregulated) in the biological process ontology. Most DEGs in the classes of transporter activity, membrane, and membrane part were upregulated, and those in the classes of response to stimulus, multi-organism process, reproduction, metabolic process, and biological regulation were downregulated. Functional enrichment analysis of GO terms (Fig. S1) further showed that CytCo treatment signi cantly affected PWN genes related to lipid localization (GO: 0010876), smooth muscle contraction (GO: 0006939), and MAPK (mitogen-activated protein kinase) cascade (GO: 0000165). KEGG richness showed that the DEGs were concentrated in lysosome, fatty acid degradation, transporters, and drug metabolism by cytochrome P450 (Fig. 2 and Table S2). Ten DEGs were found that are putatively involved in 14 signaling pathways (Table S3).
Differentially expressed genes re ect nematode responses and CytCo toxicity Many of the upregulated genes demonstrate potential self-protection of PWNs in response to the nematode toxicity of the membrane pore-forming toxin CytCo (Fig. 3), including 19 DEGs related to nematode cuticular collagen and epidermal growth factor (Table S4), 23 DEGs related to transporter (Table S5), and six DEGs encoding cytochrome P450 (Table S6). In addition, 10 out of 13 DEGs related to programmed cell death were downregulated (Table S7).

Co-expression network analysis
To further shed light on the key genes of the PWN responses to CytCo toxin, a weight gene co-expression network analysis (WGCNA) was built based on the pairwise correlation of genes across all of the samples. Highly interconnected genes were grouped into the same module, and 11 modules were obtained ( Fig. 4A). Among them, the MEturquoise module contained the largest number of genes ( Fig. 4B). As shown in Fig. 4A, MEyellow had the strongest correlation with the CytCo treatments, respectively. Most DEGs in the MEyelllow module were grouped into GO terms related to cuticle development ( Fig. 4C), implying the nematode response to the damage caused by CytCo.

Validation of RNA-seq expression data by RT-qPCR
The relative expression levels of the selected DEGs in PWNs treated with CytCo protein were evaluated by RT-qPCR (Fig. 5). The expression levels of all selected DEGs were consistent with those shown in the corresponding transcriptome data, speci cally that the genes encoding cuticular collagen, proteinase inhibitor, and transporter were upregulated and the genes related to proteinase, kinase, major sperm protein, and lipid metabolism were downregulated.

Discussion
The present study showed that the nematocidal protein toxin, CytCo, displayed diverse effects on PWN development, reproduction, infectivity, motility, and immune defenses based on whole-organism transcriptome pro ling, in consistent with the reported bioassay results [16].
There is a set of DEGs putatively involved in the interaction of CytCo with PWNs, affecting the immune system, metabolism, and proteostasis. It should be noted that the functions of genes discussed here have been studied in C. elegans and there have been limited studies in PWN. Dafachronic acids play multiple functions in suppressing dauer formation, lipid metabolism, and inducing reproductive growth in host trees, by binding directly to the nuclear hormone receptor daf-12 [24,25]. The DEG encoding DHS16 (3 βhydroxysteroid dehydrogenase), which is involve in the production of Δ 7 -dafachronic acid from cholesterol, was found to be downregulated, in agreement with the lowering lipid droplet formation in CytCo-treated PWN [16]. Nipi-3, a tribbles-like kinase is required in the response of C. elegans to fungal infection [26]. Here, a potential pathogen-associated DEG encoding TRIB3 (Tribbles homolog 3) was upregulated and might be involved in the PWN response to the fungus-originating CytCo toxicity.
Tyrosinases that may be evolutionarily related to the prophenoloxidases are responsible for melanin biosynthesis in response to wounding [27]. Here, a DEG encoding TYRO1 (Tyrosinase-like protein 1) was found upregulated and may be involved in responses to the cell lysis caused by PFT. Transthyretin-like proteins are involved in nematode innate immune defenses during the interaction of C. elegans with B. thuringiensis [28]; here, the downregulation of most of these proteins suggests a reduced immune defense during CytCo attack.
Many DEGs encoding serine carboxypeptidases, aspartyl proteases, cysteine proteases, and zinc metallopeptidases were downregulated in this study, which may be linked to impaired digestive capabilities [28]. Two upregulated DEGs encoding proteinase inhibitors (BXY_0886000 and BXY_0357700) may be related to the downregulation of three proteasome-related DEGs (BXY_1164700, BXY_1737700, and BXY_0260400), which could be linked to a dysfunction in protein homeostasis [29]. Many DEGs annotated as "kinase and phosphatase" were downregulated (Tables S9 and S10), which may cause multiple effects on development, reproduction, metabolism, and immune responses. For example, serine/threonine-protein phosphatase 1 (PP1) is essential for sperm meiosis and motility in C. elegans [30], in agreement with the previous report that the reduced fecundity of the PWN under CytCo treatment [16]. The abovementioned genes may favor developing new molecular targets to control PWN.
In response to protein toxin stress, the rst molecular PFT defense pathways identi ed in nematodes were two mitogen-activated protein kinase (MAPK) pathways (p38 and c-Jun N-terminal kinase (JNK)like) in C. elegans in response to Cry5B toxin [31]. The MAPK cascades are central signaling pathways that regulate a wide variety of stimulated cellular processes, including proliferation, differentiation, apoptosis and stress response [32]. Here, a DEG (BXY_0768000) encoding CRE-HSP-70 was found upregulated, which may inhibit apoptosis through a JNK-like MAPK pathway and involve defense against CytCo in PWNs (Fig. S2). In this study, many DEGs related to programmed cell death were downregulated (Table S7), which may be attributed to this activated pathway.
Activation of the necrosis signaling pathway by Cry6Aa has been shown to play an important role in cell death in C. elegans [13]. Necrosis is characterized by the loss of plasma membrane integrity and the ensuing cell death can contribute to in ammation [33]. Two necrosis-related DEGs encoding TFIP8 (tumor necrosis factor α-induced protein 8-like protein) and LITAF (lipopolysaccharide-induced tumor necrosis factor-α factor-like protein) were downregulated in this study and may act as negative mediators of cell death. This implies that there are differences in modes of action between Cry and Cyt toxins.
Additionally, the upregulated DEG encoding CREB binding protein isoform X1 (a cyclic AMP-responsive element binding protein) may in uence the homeostasis of lipids and proteins in PWNs via the Jak-STAT signaling pathway, which is also involved in the immune system [34]. Downregulated GSK3 (a serine/threonine-protein kinase) may cause an increase in β-catenin and activate Wnt signaling, which links to metabolism and stem cell self-renewal [35,36]. Downregulated ADT2 (an adenine nucleotide translocator) may in uence PWN development and body size by modulating TGF-β signaling activity and organizing cuticle collagen brils in C. elegans [37].
In our study, the upregulated DEGs related to sodium/sulphate symporter and potassium channel proteins (Table S5) might also cause ionic imbalance besides that produce consequences similar to the pore-forming effects of CytCo on cell membrane. Moreover, the effects of CytCo in the bioassays were analogous to the adverse effects of the chemical nematicide, emamectin benzoate (EB), including lowered fecundity, hatching rate, and thrashing frequency [8,16]. However, some of the observed DEGs were unique, and even the shared DEGs had different expression patterns (Fig. S3). For example, many cuticular collagen-related DEGs were upregulated and programmed cell death-related DEGs were downregulated with CytCo treatment, but the opposite was reported with EB treatment [8]. Thus, the compound of protein toxin and chemical agents could be developed for the nematode control.

Conclusion
The present study provided initial insights into multiple effects of the nematocidal CytCo on PWN on the transcriptomic level. We hypothesized two modes of action of CytCo nematotoxicity (Fig. 6). One mode based on the mainstream recognition of pore-formation protein toxins is that the CytCo targets phospholipid to cause apoptosis and necrosis of intestinal cells, and then digestive enzymes and ions enter the coelom to kill the nematode. It associates with the nematode responses to CytCo via upregulating the genes encoding cuticle collagens and transporters and downregulating the genes related to programed cell death and hydrolases. Compared to the mode for the acute toxicity, we also hypothesized the other lethal mechanism via energy metabolism. It associated with lipid metabolismrelated genes affected by CytCo, warranting further studies.

Preparation Of Cytco
According to the methods described by Zhou et al. [16], CytCo protein was expressed and puri ed. Brie y, Escherichia coli Arctic-Express™ cells (Agilent Technologies, USA) with the recombinant plasmid (pCzn1-CytCo) was inoculated for heterologous expression. The CytCo-expressing cells were harvested by centrifugation and lysed by sonication in an ice-water bath. The CytCo was eluted from the a nity chromatography by loading the cleared bacterial lysate onto a 1-mL Ni-IDA-Sepharose Cl-6B a nity column (Novagen). The protein was extensively dialyzed overnight with PBS (pH 7.4), the nal protein concentration was assessed using Bradford Protein Assay Kit (Takara) and bovine serum albumin as a standard.

Preparation of nematodes
PWNs (isolate NB-6) were collected from PWD-outbreak forests in Ningbo city, Zhejiang, China, and fed on 7-d-cultivated Botrytis cinerea Pers. using potato dextrose agar (PDA) plates at 25°C. Newly emerged stage larvae (L2) were collected and inoculated on B. cinerea plates in batches. In 3 days, the larvae developed into adults. The Baermann funnel method was used to separate the nematodes from each PDA plate, and the nematode samples (10,000 nematodes/ml) were collected by centrifugation (4000 g) for 4 min [8]. PWN adults (2000 nematodes/sample) were collected after being treated with 20 µg/ml puri ed CytCo or phosphate buffered saline (PBS, pH 7.4) in the dark for 24 h at 25°C.

RNA sampling
There were three biological replicates for total RNA extraction per treatment. The TRIzol Max Bacterial RNA Isolation Kit (Thermo Fisher Scienti c, New York, USA) was used following the manufacturer's protocol. The RNA concentration and purity were measured on a NanoDrop2000 (Thermo Fisher Scienti c, USA) and the integrity was veri ed by 1% agarose gel electrophoresis and on the Agilent 2100 BioAnalyzer (Agilent, USA). The extracted RNA samples (~ 3 µg of total RNA per sample) were stored at − 80°C and then sent to Woosen Co. (Hangzhou, China) for sequencing.

Transcriptome analysis
Nematode mRNA was enriched from each total RNA sample using Oligo(dT) magnetic beads. Paired-end RNA-seq libraries of different treatments were prepared following Illumina's library construction protocol, and the libraries were then sequenced on the Illumina HiSeq 2500 platform (Illumina, USA). FASTQ les were produced and sorted. To obtain high-quality clean reads, raw reads were removed if they contained adapters, had ≥10% unknown nucleotides, or possessed other low-quality indicators. The leftover clean reads were mapped to a reference genome (PRJEA64437, www.wormbase.org) using TopHat2 (version 2.0.3.12). The raw sequencing data were deposited in the China National GeneBank DataBase (CNGBdb, https://db.cngb.org/) with accession number CNP0001233.
The R package RSEM was used to calculate the fragments per kilobase of exon per million fragments mapped (FPKM) value [38]. Differentially expressed genes (DEGs) between libraries were ltered out using the R package DEGsEq. DEGs were identi ed with a fold change ≥2 and a false discovery rate <0.001 [39]. In order to provide further insight into the DEGs involved in the modes of CytCo working on PWNs, the functions of DEGs were predicted by annotation using several databases. Speci cally, the genes were annotated by BLASTx search (E-value < 10 −5 ) against Wormbase, Swiss-Prot (www.uniprot.org), Gene Ontology (GO, www.geneontology.org), and Kyoto Encyclopedia of Genes and Genomes (KEGG; www.genome.jp/kegg/kegg2.html) databases. Functional enrichment analysis for GO terms and KEGG pathways was carried out using the R package clusterPro ler [40].

Co-expression network construction
To screen the distinct hub genes in uenced by CytCo toxin, gene co-expression networks were established through the WGCNA (v1.69) package of R [41]. Gene expression values were imported into WGCNA to construct co-expression modules using the automatic network construction function block-wise modules with default settings. The expression level of the DEGs was log transformed using log2 (FPKM + 1). Pearson's correlations coe cient was used to measure the co-expression relationship between each pair of genes. The WGCNA network was constructed with a soft thresholding power of β = 17, a minimum module size of 30 genes, and the TOM-Type was unsigned, the merge Cut Height was 0.25. The moduletrait relationship was used to differentiate the hub genes between CytCo and emamectin benzoate (EB) treatments. Original data of the transcriptome of EB treatments are downloaded from the Gene Expression Omnibus (www.ncbi.nlm.nih.gov/geo) under accession number GSE135014.

RT-qPCR analysis of selected DEGs
The RNA samples were extracted from PWNs exposed to 20 µg/ml puri ed CytCo or PBS for 24 h. To validate the CytCo-affected DEGs, RT-qPCR was performed to assess the transcript levels of selected genes, namely those encoding proteinase inhibitor (BXY_0886000), serine carboxypeptidase (BXY_0963400), cathepsin (BXY_0408100), epidermal growth factor (BXY_1430300), major sperm protein (BXY_0820100), arginine kinase (BXY_1237900), ATP-binding cassette transporter (BXY_0203900), elongation of very long chain fatty acids protein (BXY_1705500), serine/threonine-protein kinase (BXY_0372900) and ATP-binding protein permease (BXY_1678600). For RT-qPCR analysis, total RNA (1 µg) was rst reverse-transcribed into cDNA with a PrimeScript™ RT Reagent Kit with gDNA Eraser (TaKaRa, Tokyo, Japan). The qPCR of the cDNA samples was then performed using a SYBR Green PCR kit (SYBR Premix Ex Taq™ II, TaKaRa). Paired primers were designed and are listed in Table S1. PCRs were performed on a real-time PCR thermal cycler (qTOWER 2.2, Analytikjena, Germany) under the following conditions: 95°C for 30 s, followed by 40 cycles at 95°C for 5 s and 60°C for 30 s. The data were analyzed using qPCRsoft 1.1 (Analytikjena). Transcript quanti cation for each gene was performed using at least three independent replicates. The relative fold-change in gene expression was normalized to that of elongation factor 1-alpha (ef-1α) (BXY_0569100).

Declarations Ethics approval
This article does not contain any studies with human participants or animals performed by any of the authors.  Relative expression levels of the selected differentially expressed genes (DEGs) in pine wood nematode (PWN). cDNA samples were derived from PWN treated with phosphate buffered saline (PBS, control), and 20 µg/ml CytCo protein, and fold changes (FC) of relative expression levels of DEGs were calculated based on the analysis of real-time quantitative PCR. The selected DEGs contain the genes encoding proteinase inhibitor (BXY_0886000), serine carboxypeptidase (BXY_0963400), cathepsin (BXY_0408100), epidermal growth factor (BXY_1430300), major sperm protein (BXY_0820100), arginine kinase (BXY_1237900), ATP-binding cassette transporter (BXY_0203900), elongation of very long chain fatty acids protein (BXY_1705500), serine/threonine-protein kinase (BXY_0372900) and ATP-binding protein permease (BXY_1678600). Error bars: SEM from three biological replicates. Primers listed in Table S14.