Transcriptomics and Co-expression Network Pro ling of Effects of Levamisole Hydrochloride on Bursaphelenchus Xylophilus


 Background: Bursaphelenchus xylophilus is one of the most dangerous forest pathogens in the world, resulting in the devastating death of pine forests and causing great economic losses and forest ecological damage to the affected areas. To explore the effect of levamisole hydrochloride (LH) as an effective control agent on nematode.Results: The results indicated that 2.5 mg/mL and 3.5 mg/mL LH had the toxicological effect on B. xylophilus, and the mortality increased significantly with the increase of concentration (P<0.05). We found and studied Sodium / Chloride and other ion genes belonging nervous system were up-regulated and the ion signal transmitted to the muscle protein and cause disorders, producing body-wall muscle twitchin, paralysis, and ultimately death. In the other hand, 2 Senecionine N-oxygenase and 2 Alcohol dehydrogenase as the central genes in drug metabolism-cytochrome P450 (ko00982) and metabolism of xenobiotics by cytochrome P450 (ko00980), and 8 highly related genes including Cuticle collagen, Cystathionine beta-synthase, Endochitinase, Pyruvate dehydrogenase E1 component subunit beta, Aldehyde dehydrogenase, Lipase, Zinc metalloproteinase. Their expressions up-regulation significantly at low concentrations and were significantly related to the resistance of B. xylophilus to LH. Conclusion: This study has shown that B. xylophilus gene family expansions occurred in xenobiotic detoxification pathways through the degree of genes expression and potential horizontal correlated gene transfer with LH, and shed light on LH lethality and evolutionary mechanisms behind adaptations of B. xylophilus to the environment. These findings contribute in several ways to our understanding of B. xylophilus under LH and provide a basis for control of it.


Background
Pine wilt disease (PWD) caused by Bursaphelenchus xylophilus has caused serious damage to pine trees, resulting in signi cant economic and ecological loss [1,2]. B. xylophilus is native to North America, where it does not cause mortality in endemic hosts, occurring sporadically [3,4]. When it was exported to Japan, South Korea, and China in East Asia with diseased trees, it caused devastating damage to the local pine ecosystem [2,5,6].The problem of B. xylophilus caused widespread concern. Nematocide injection into tree trunks is one of the most widely used comprehensive control measures at present, but the effect is not ideal because of its hidden harm [7][8][9]. So far, there have been many reports about natural nematicidal phytochemicals from various genera. Avermectin injection is often used in South Korea to protect pine [10], but its high prices limit the feasibility of its widespread use. Chemical nematicidal agents, such as fosthiazate carbofuran [11,12], are mostly broad-spectrum and non-speci c. While controlling nematodes, they are also very easy to harm non-harrowing organisms, including natural enemy insects, such as Scleroderma guani Ichneumonidae. It is urgent to nd speci c, e cient, and cheap nematicidal agents.
Levamisole hydrochloride is a white crystalline powder with a broad-spectrum, highly effective and low toxic nematicide. It is effective against gastrointestinal nematodes in many animals [13], but harmless to mammalian hosts, and even has can improve immunity to individuals with low immunity [14]. It was rst proposed that acted on the drug-resistant mutant of acetylcholine receptors in Caenorhabditis elegans [15]. Later, it was con rmed that LH selectively acted on nicotinic acetylcholine receptor (nAChRs) on the somatic muscle and nerve of C. elegans, resulting in contraction and spastic paralysis of nematode [16,17].In addition, nAChRs was found in the non-excited tissues. Among them, the L-AChR is the target of nematocide agents, such as levamisole and pyrantel, which, by acting as potent agonists, produce body-wall muscle hypercontraction, paralysis, and ultimately death [18]. Recent studies have shown that in addition to neuroreceptors, it has been found that intestinal nAChRs, activated by levamisole in Ascaris lumbricoides can promote the uptake and digestion of intestinal contents of nematodes [19]. But the targets of plant-parasitic nematodes have not been studied.
There is an urgent need to address the speci c, e cient, and inexpensive control of these problems caused by B. xylophilus. Little is known about B. xylophilus xenobiotic detoxi cation and immunity, and it is not clear what factors of LH metabolism. To explore the genetic basis of adaptation of B. xylophilus to LH, the transcriptome data were compared with those of B. xylophilus [20]. We studied pathogenesis, xenobiotic detoxi cation pathways, and candidate detoxi cation genes were identi ed. We propose that B. xylophilus response to LH was varied and complex, which facilitated its survival in LH environment. This is the rst time to study the effect of LH to B. xylophilus and to explore the response mechanism of it.

Results
Effects of LH on the activity of B. xylophilus By toxicity test, the higher the concentration of LH, the higher the mortality of B. xylophilus with signi cant difference (P < 0.05) (Fig. 1). The highest mortality rate of LH (3.5 mg/mL) was 55% (Additional le 1).
RNA sequencing and quality evaluation.
After removing rRNA and low-quality readings, 515,955,240 high-quality clean readings were obtained from nine samples (Table 1). In all data, the mass value ≥ 30 base ratios were more than 94.13%, and the proportion of total readings located to the reference genome was between 68.15% and 75.28%. The conversion of the Counts to FPKM. Dendrogram clustering of the samples for each RNA sequencing showed a conservative comparison between samples. At the same time, the Pearson correlation value of repeated samples was signi cantly higher than that of inter-processing correlation (Fig. 2).  Predictive functional analysis of the highly expressed genes for different types of pro le 487 DGEs were obtained from nine groups (control, 2.5 mg/mL and 3.5 mg/mL). According to the different trends of genes expression under different conditions, a set of unique model expression trends (similarity > 0.7) were identi ed. The trend of FPKM of different concentrations of DGEs was analyzed by STEM software, and eight differential mRNA expression trends were calculated by using the signi cant expression of gene array, of which three pro les were drawn in the color part (pro le 1,6 and 7 DEGs p < 0.01) (Fig. 4A). The signi cant trend of Pro le1 (including 172 DEGs) showed that the expression level in DEPC-H2O (CK) was signi cantly higher than that in 2.5 mg/mL LH (YD). However, there was no signi cant change from YD group to 3.5 mg / ml LH (YG) group (Fig. 4B). The signi cant trend of pro le 6 and pro le 7 containing 161and 58 DEGs, from pro le 6 and 7 showed a sharp up-regulation from CK group to YD group, and there was no signi cant change from YD group to YG group in pro le 6 ( Fig. 4C), but the DEGs expression level increased signi cantly in pro le 7 (Fig. 4D).

Construction and analysis of the co-expression network
According to the sequencing data, DEGs were selected as the centre genes in the xenobiotic detoxi cation pathways, including drug metabolism-cytochromeP450 (ko00982) and metabolism of xenobiotics by cytochrome P450 (ko00980  . 8) (Additional le 8). And the expression trends of these genes were always up-regulated, consistent with RNA sequencing analysis (Fig. 9).

Discussion
B. xylophilus is highly destructive to pines. Until now, there is a lack of effective prevention and control of it. So, genes in the xenobiotic metabolism were very important to speculate the target for the control of B. xylophilus [21]. In eukaryotes, genes in the xenobiotic metabolism can remove drugs from cells, resulting in a decrease in intracellular drug concentration, contributing for drug resistance [22]. Several studies have already shown the multilayer detoxi cation of B. xylophilus against plant defences, which involves the detoxi cation and xenobiotic metabolic pathways of the nematode. Even there is a study on the effect of emamectin benzoate in B. xylophilus with the same approach. So, it is necessary to study the function of xenobiotic metabolism in order to provide theoretical basis for the development of LH control methods and targets. LH is a broad-spectrum, high-e ciency, and low-toxic nematicide, which has a good killing effect on chicken, sheep, and human gastrointestinal nematodes [13]. The effect of LH, xenobiotic detoxi cation and pathogenic mechanism of B. xylophilus has important value and provides a theoretical basis for prevention B. xylophilus.
Nine RNA libraries directly re ected the effect of LH on the genetic activity of B. xylophilus. Among the 32 common Gene Ontology (Fig. 3A), antioxidant activity, catalytic activity, detoxi cation, immune system process, molecular transducer activity, synapse, synapse part, transporter activity, signal transducer activity showed that LH treatment had a signi cant effect on the neural transmission and immune system of B. xylophilus [23]. The single-organism process was only in the group of YDvsYG, includes a 3'5'-cyclic nucleotide phosphodiesterase and 3 Cuticle collagen, which were signi cantly up-regulated from low concentration to high concentration. The Cuticle collagen, an extracellular collagen matrix that surrounds C. elegans, is the rst barrier against environmental damage, protecting nematodes from pathogens, dryness, and other stress [24][25][26]. The up-regulation of the high concentration of Cuticle collagen indicated that the integument of B. xylophilus thickened, the resistance of nematodes was enhanced. Xenobiotic detoxi cation pathways with drug metabolism -cytochrome P450 (ko00982) and metabolism of xenobiotics by cytochrome P450 (ko00980) were found, including most of the upregulation FMO, UDT, GST, ADH and a Juvenile hormone epoxide hydrolase down-regulations (Fig. 3B) (Additional le 4). Compared with C. elegans, B. xylophilus has speci c gene up-regulation in FMOs branch (FMO-5, FMO-10), which has a function similar to cytochrome P450 and adds molecular oxygen to insoluble xenobiotic organisms. Many insects and microorganisms can metabolize chemicals into less toxic or more soluble compounds. They use cytochrome P450 or FMO to oxidize it to alcohols, then ADHs to oxidize alcohols to aldehydes/ketones [22], and nally ALDHs to oxidize these products to acids for subsequent metabolism or transport. ADH is one of the key genes regulating the dynamic behavior of B. xylophilus [27]. P450, UGT, glutathione -S transferase, Redox, dehydrogenase, and reductase of C. elegans under ve different Benzimidazole drugs increased. UGT are enzymes that are glucuronidated the substrates include drugs, environmental toxins, and endogenous compounds. And the products are watersoluble and readily excreted[28]. UGT inhibitors resulted in a decrease of metabolite production in nematode and the decrease of drug resistance [29]. UGT inhibitor increased the sensitivity of wild type H. contortus to naphtalophos (anthelmintic) [30]. GST can catalyze the binding of glutathione with electrophilic groups of chemical substances, and eventually form mercapturic acid to be excreted in vitro [31]. And 3 GST and P450 were upregulated contributed to mitigating tetrabromobisphenol A-induced toxicity on the C. elegans [32]. Thus, based on the results and previous studies, UDT, GST, FMO, ADH, Juvenile hormone epoxide hydrolase were considered as one of the main factors of B. xylophilus detoxi cation ability.
Sets of pro les were prede ned which were expected in the process of genetic changes. These sets of pro les cover all of the possible gene expression changes, and each represents a single drug concentration expression pattern. Pro le signi cance levels were determined by gene enrichment, and signi cant pro les indicate the most common functions of co-expressed genes, with the main biological character [33][34][35]. Xenobiotic detoxi cation pathways, such as Glycolysis/ Gluconeogenesis, Glutathione metabolism, Metabolism of xenobiotics by cytochrome P450, Drug metabolism-cytochrome P450 were signi cant enrichment. Also, Lysosome, Phagosome, Peroxisome, and other immune-related pathways were signi cantly enriched (Fig. 5B Fig. 6B), indicating that the innate immune defense of B. xylophilus increased, such as Lysosome, Cysteine proteinase, and transthyretin-like proteins [36]. On the other hand, it has been reported that LH has the function of enhancing immune response by enhancing phagocytosis [37][38][39]. LH is a typical deworming drug that activates nematode nAChR, which acts as a powerful agonist in muscle contraction and movement and is a key target [36,40]. NAChR transfers the ion signal to disorganized muscle protein through Sodium/Chloride and causes gene up-regulation, produce body-wall muscle twitchin, paralysis, and ultimately death in B. xylophilus (Additional le 9). This is one of the direct causes of the death of B. xylophilus, but how LH acts on the nervous system signal transduction process of B. xylophilus awaits further exploration (Fig. 10) (Additional le 10 A).
Genes in the same biological pathway are more likely to be co-expressed to synchronize an array of biochemical reactions [41,42]. We studied the connectivity patterns of drug metabolism-cytochromeP450 (ko00982) and metabolism of xenobiotics by cytochrome P450 (ko00980) pathways in pro le1and pro le 6 gene expression data. In Pro les 6, we observed that 54 genes were associated with the correction of Pearson, that is, the corresponding P-value was < 0.05. FMO (BUX. Gene. S00600.44, BUX. Gene. S00713.800) and ADH (BUX. gene. S01281.195, BUX. gene. S01147.198) were signi cant correlation in the center gene (Fig. 7), which were highly consistent with the previous conjecture that its play important role in xenobiotics metabolism. Among them, highly correlated Cuticle provides cysteine, which is important for the formation of the stratum corneum in C. elegans. As well as, it produces neuromodulator and smooth muscle relaxant hydrogen sul de in muscle and subcutaneous tissue [45][46][47].
In conclusion, from the transcriptomes of B. xylophilus treated with 2.5 mg /mL,3.5 mg/mL LH, and DEPC-H 2 O, with STEM and Spearman correlation, pathogenesis, xenobiotic detoxi cation pathways, and candidate detoxi cation genes were identi ed. We propose that B. xylophilus response to LH was varied and complex which facilitated its survival in LH environment (Fig. 10) (Additional le 10 B). However, this is the rst time to study the effect of LH on B. xylophilus and to explore the resistance mechanism of B. xylophilus and these ndings contribute in several ways to our understanding of B. xylophilus under LH. Lays the groundwork for future research into LH could also be conducted to determine the effectiveness of control of B. xylophilus. xylophilus was added to the nal concentrations YD group (2.5 mg/mL), the YG group (3.5 mg/mL), and the CK group (DEPC-H 2 O) respectively, and mixed well. After 12 hours, the number of nematodes was counted by the dilution concentration method on the 24-well plate, under the anatomical microscope. The phenotype of the nematodes was observed, and the mortality rate was calculated (the body was stiff and the nematodes that did not move when touching the needle were dead). The biologicalexperiment was repeated 6 times. Nematodes were washed with double distilled water and DEPC-H 2 O, removing by 4000 rpm in 5 min. Then nematodes were frozen in liquid nitrogen and stored at -80℃ for RNA extraction. The One-way ANOVA (Duncan's) of the data was performed using IBM SPSS Statistics 19.

Materials And
RNA extraction, library construction, RNA sequencing and assembly From 9 samples, TRIzol reagent (Invitrogen, USA) was used to extract total RNA. RNA degradation and contamination was monitored on 1% agarose gels. RNA purity was checked using the NanoPhotometer spectrophotometer (IMPLEN, CA, USA). RNA concentration was measured using Qubit RNA Assay Kit in Qubit 2.0 Flurometer (Life Technologies, CA, USA). RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA).
A total amount of 1 µg RNA per sample was used as input material for the RNA sample preparations. Sequencing libraries were generated using NEBNext UltraTM RNA Library Prep Kit for Illumina (NEB, USA) following manufacturer's recommendations and index codes were added to attribute sequences to each sample. Brie y, mRNA was puri ed from total RNA using poly-T oligo-attached magnetic beads. Fragmentation was carried out using divalent cations under elevated temperature in NEBNext First Strand Synthesis Reaction Buffer(5X). First strand cDNA was synthesized using random hexamer primer and M-MuLV Reverse Transcriptase. Second strand cDNA synthesis was subsequently performed using DNA Polymerase I and RNase H. Remaining overhangs were converted into blunt ends via exonuclease/polymerase activities. After adenylation of 3' ends of DNA fragments, NEBNext Adaptor with hairpin loop structure were ligated to prepare for hybridization. In order to select cDNA fragments of preferentially 200-250 bp in length, the library fragments were puri ed with AMPure XP system (Beckman Coulter, Beverly, USA). Then 3 µl USER Enzyme (NEB, USA) was used with size-selected, adaptor-ligated cDNA at 37 °C for 15 min followed by 5 min at 95 °C before PCR. Then PCR was performed with Phusion High-Fidelity DNA polymerase, Universal PCR primers and Index (X) Primer. At last, PCR products were puri ed (AMPure XP system) and library quality was assessed on the Agilent Bioanalyzer 2100 system The clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v4-cBot-HS (Illumia) according to the manufacturer's instructions. After cluster generation, the library preparations were sequenced on an Illumina Hiseq 2500 platform and paired-end reads were generated by BIOMARKER Biotechnology Co., Ltd. (Beijing, China). Raw data (raw reads) of fastq format were rstly processed through in-house perl scripts. In this step, clean data(clean reads) were obtained by removing reads containing adapter, reads containing ploy-N and low quality reads from raw data. At the same time, Q20, Q30, GC-content and sequence duplication level of the clean data were calculated. All the downstream analyses were based on clean data with high quality.
The adaptor sequences and low-quality sequence reads were removed from the data sets. Raw sequences were transformed into clean reads after data processing. These clean reads were then mapped to the reference genome sequence. Only reads with a perfect match or one mismatch were further analyzed and annotated based on the reference genome. Then use TopHat2 (version 2.0. 3.12)[48] to map the clean readings of each sample to the reference genome [20]. Based on the Cu inktware, the FPKM method quanti ed the gene expression level by comparing the location information on the genome [49].

Identi cation of differentially expressed genes (DEGs), Sequencing Data Analysis
The DESeq2 package (https://bioconductor.org/) [23] was used to identify DEGs across groups with fold changes ≥ 2and a false discovery rate (FDR) < 0.05. Basic annotation of genes includes NCBI non-redundant protein (Nr) database, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway annotation, Clusters of orthologous groups of proteins (COG) functional annotation, Swiss-Prot protein database, and Gene Ontology (GO) annotation. DEGs were had to enrichment analysis by using the ClusterPro ler package (https://bioconductor.org/).

Trend analysis of DEGs
The differential gene expression was was calculated and normalized to FPKM. To evaluate the gene expression pattern of nematode under three different concentrations, DEGs were assigned to 8 different types of pro les by Short Time-series Expression Miner software (STEM) [50]. The maximum Unit Change in model pro les between points was 1 and the maximum output pro le number was 20. The correlation coe cient between the gene expression pattern and the most similar trend > 0.7 was classi ed into the same model. Signi cant pro les were determined and the threshold of signi cance was considered as p < 0.05. Then the signi cant pro les were subjected to GO term enrichment and KEGG enrichment to research functions.
Construction of the co-expression network of the signi cant pro les-Spearman's rank correlation In pro les, DEGs expression data set matrix calculated the Spearman correlation coe cient (rs) between the central gene and the other gene, and corrected (p) by cor.test package (R code). Removing the in uence of gene itself and gene quantity [41], the genes with signi cant connection (p < 0.05) were selected to construct the expression-related network. The networks were visualized using Cytoscape_3.7.1. [51]. d i = the difference between the ranks of the i th observations of the two variables. n = the number of pairs of values [41].

Validation of gene expression by qRT-PCR
Twelve DEGs were selected for qRT-PCR analysis. Primer 6.0 software was used to design the primers and completed by Sangon Biotech (Shanghai, China) (Additional le 9). Bx28S and Actin gene were selected as the housekeeping gene [52]. Agilent Technologies Mx3000P was used for real-time PCR detection, and qPCR was performed following the protocol provided by the TaKaRa One-Step RT-PCR Kit (Baori medical technology (Beijing) co., LTD, Beijing, China) and compared to using the threshold cycle (-ΔΔCt) method of gene expression.

Declarations
Ethics approval Figure 1 The mortality rate of B. xylophilus with different concentrations of LH. (*) It is a signi cant difference between the two groups (p<0.05).     Interaction of gene co-expression patterns in pro les 6, those were used for the construction of the network in Cytoscape 3.7.1. Note: Different colors represent different kinds of genes Interaction of gene co-expression patterns in pro les 7, those were used for the construction of the network in Cytoscape 3.7.1. Note: Different colors represent different kinds of genes.

Figure 10
The hypothesis effects of LH on B. xylophilus.