Transcriptome Analysis of Bursaphelenchus xylophilus Uncovers the Impact of Stenotrophomonas maltophilia on Nematode and Pine Wilt Disease

Stenotrophomonas maltophilia influences the reproduction, pathogenicity, and gene expression of aseptic Bursaphelenchus xylophilus after inoculation of aseptic Pinus massoniana. Pine wilt disease is a destructive pine forest disease caused by B. xylophilus, and its pathogenesis is unclear. The role of bacteria associated with B. xylophilus in pine wilt disease has attracted widespread attention. S. maltophilia is one of the most dominant bacteria in B. xylophilus, and its effect is ambiguous. This study aims to explore the role of S. maltophilia in pine wilt disease. The reproduction and virulence of aseptic B. xylophilus and B. xylophilus containing S. maltophilia were examined by inoculating aseptic P. massoniana seedlings. The gene expressions of two nematode treatments were identified by transcriptome sequencing. The reproduction and virulence of B. xylophilus containing S. maltophilia were stronger than that of aseptic nematodes. There were 4240 differentially expressed genes between aseptic B. xylophilus and B. xylophilus containing S. maltophilia after inoculation of aseptic P. massoniana, including 1147 upregulated genes and 2763 downregulated genes. These differentially expressed genes were significantly enriched in some immune-related gene ontology (GO) categories, such as membrane, transporter activity, metabolic processes, and many immune-related pathways, such as the wnt, rap1, PI3K-Akt, cAMP, cGMP-PKG, MAPK, ECM-receptor interaction, and calcium signaling pathways. The polyubiquitin-rich gene, leucine-rich repeat serine/threonine-protein kinase gene, glyceraldehyde-3-phosphate dehydrogenase (GAPDH) gene, acetyl-CoA carboxylase gene, and heat shock protein genes were the key genes associated with immune resistance. Moreover, there were four cell wall hydrolase genes, thirty-six detoxification- and pathogenesis-related protein genes, one effector gene and ten cathepsin L-like cysteine proteinase genes that were differentially expressed. After inoculation of the host pine, S. maltophilia could affect the virulence and reproduction of B. xylophilus by regulating the expression of parasitic, immune, and pathogenicity genes of B. xylophilus.


Introduction
Pine wilt disease (PWD) is a devastating disease in pines caused by the pine wood nematode (PWN), Bursaphelenchus xylophilus. It is native to North America [1], and epidemic in East Asian (China, Japan, and Korea) [2][3][4], causing serious ecological damage and economic loss. Since PWD was reported in 1982 in Nanjing, China, it has spread to 18 provinces with an area of more than 1.1 million hm 2 , and more 1 = 0-25% of needles turned yellow; 2 = 25-50% of needles turned yellow; 3 = 50-75% of needles turned yellow; and 4 = 75-100% of needles turned yellow. The calculation formula was as follows: Infection rate (%) = Number of infected plants with symptoms Total number of plants

×100
Disease severity index (DSI) = Number of diseased plants × symptom stage Total number of plants × highest symptom stage ×100 Seven days post inoculation (dpi), the PWNs of each treatment were separated under aseptic conditions by the Berman funnel method, and aseptic PWNs (PmBx_a) and PWNs with S. maltophilia (PmBx_b) were obtained and their reproduction were counted. The ratio of the final population to the initial reproduction of nematodes were presented as means ± standard deviation (Mean ± S.D.). All parameters were calculated using Microsoft Excel and GraphPad Prism 5 (GraphPad Software, San Diego, CA, USA). The statistical significance was determined using SPSS Statistics 24.0 software (IBM China Company Ltd., Beijing, China). A Student's t-test was used to compare two samples. The level of significance was p < 0.05. A small amount of nematode solution from each treatment was taken for the sterility test and morphological observation. The rest of the nematodes were rinsed with sterile water 3-5 times and centrifuged, and the supernatant was removed, frozen in liquid nitrogen, and stored in a refrigerator at −80 • C for subsequent analyses.

RNA Extraction and Sequencing
Six samples (PmBx_a_1, PmBx_a_2, PmBx_a_3, PmBx_b_1, PmBx_b_2, and PmBx_b_3, each treatment had three replicates) were used for RNA sequencing. Total RNA was extracted from the nematodes with TRIzol reagent (Invitrogen, Waltham, MA, USA) according to the manufacturer's protocol. The RNA quality was assessed using an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Total RNA with a concentration higher than 20.0 ng/µL, a total mass higher than 1.0 µg, an RNA integrity number (RIN) over 7.0, and 28S/18S no less than 1.0 was used for library preparation. RNA-Seq library construction and paired-end sequencing on an Illumina HiSeq X-Ten was completed at Beijing Genomics Institute (Shenzhen, China). The RNA-Seq experiments have three technical replicates. The data have been submitted to NCBI, and the accession number is SRR11474049-SRR11474054.

RNA-Seq Analysis
The RNA-Seq data were cleaned using SOAPnuke (https://github.com/BGI-flexlab/SOAPnuke) to remove adaptors, reads with unknown bases (N) > 5%, and low-quality sequences (quality score (Q20) < 20). After filtering, high-quality reads were then mapped to the B. xylophilus genome (BioProject PRJEA64437) [27] using Hierarchical Indexing for Spliced Alignment of Transcripts (HISAT) [28]. StringTie was used to do transcript reconstruction of each sample [29]. Cuffmerge was used to integrate the reconstruction information of all samples, and then Cuffcompare was used to compare the integrated transcript with the reference annotation information [30]. Transcripts with class code types of "u", "i", "o", and "j" were selected and defined as the new transcripts. The protein coding potential of the new transcripts was predicted by CPC [31], and the new transcripts with predicted protein coding potential were added to the reference gene sequences to obtain a complete reference sequence information, which would be analyzed in the future. Clean reads were compared to the reference sequence to calculate the gene alignment rate by Bowtie 2 [32]. The expression levels of genes and transcripts were calculated using RSEM (RNA-Seq by Expectation Maximization) [33]. Genes with an absolute value of log 2 fold change ≥ 1 and an adjusted p-value (padj) ≤ 0.05 were considered to be differentially expressed genes (DEGs) [34]. The reference genome of B. xylophilus was used as the background to determine significantly enriched gene ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways enriched within DEGs. A false discovery rate (FDR) ≤ 0.01 was used to indicate statistically significant enrichment. DIAMOND [35] software was used to compare the DEGs to the STRING [36] database, resulting in the identification of DEG-encoded protein interactions based on their homology with known proteins. The whole interaction relationship was visualized by Cytoscape software [37].

qPCR Validation
The real-time quantitative polymerase chain reaction (qPCR) was carried out to confirm the expression of 9 randomly selected DEGs (significantly different expression) using the same RNA samples that were used in sequencing. First-Strand cDNA was synthesized using HiScript ® II Q RT SuperMix for qPCR (+gDNA wiper) (Vazyme, Nanjing, China) following the product manual. The gene-specific primers for q-PCR were designed by Primer Premier 5 (Table S1). The housekeeping actin gene of B. xylophilus was utilized as a reference gene. The cDNA products were amplified by SYBR Green Master Mix (Vazyme, Nanjing, China) on the ABI PRISM 7500 Real-Time PCR System (Applied Biosystems, Foster City, CA, USA). Each q-PCR was conducted in triplicate, and then relative gene expression levels were analyzed by the 2 −∆∆Ct method.

Effect of S. maltophilia on the Virulence of B. xylophilus
At 3 dpi, thirteen P. massoniana seedlings inoculated with aseptic PWNs (Bx_a) began to lose water and turn brown; the infection rate was 54.17%, and the disease severity index (DSI) was 13.54. Twenty-one pine seedlings inoculated with the nematode with S. maltophilia (Bx_b) showed infection; the infection rate was 87.50%, and the DSI was 21.88. At 5 dpi, 23 pine seedlings inoculated with Bx_a showed symptoms of PWD; the infection rate was 95.83%, and the DSI was 23.96. All of the pines inoculated with Bx_b had symptoms of PWD; the infection rate was 100%, and the DSI was 25. At 7 dpi, all the seedlings of P. massoniana inoculated with Bx_a were infected, the infection rate was 100%, the DSI was 37.50; the infection rate of P. massoniana inoculated with Bx_b was 100%, and the DSI was 38.54. The P. massoniana inoculated with S. maltophilia alone did not show any symptoms ( Figure 1, Table 1). The seedlings inoculated with Bx_b showed earlier and more visible symptoms than those inoculated with Bx_a. At this time, the number of PmBx_a was 6.74 times of the initial inoculation (Bx_a), and the number of PmBx_b was 7.91 times of the initial inoculation (Bx_b), but the difference between the two treatments was not significant ( Figure 2). It is suggested that S. maltophilia could enhance the virulence of PWNs.

Transcriptome Sequencing and Differentially Expressed Genes
RNA sequencing of six independent B. xylophilus samples (three PmBx_a and three PmBx_b) was performed on an Illumina HiSeq platform to characterize the transcriptome and measure the differential expression of genes, generating 54.03 to 58.30 Mb of raw reads per sample. After quality

Transcriptome Sequencing and Differentially Expressed Genes
RNA sequencing of six independent B. xylophilus samples (three PmBx_a and three PmBx_b) was performed on an Illumina HiSeq platform to characterize the transcriptome and measure the differential expression of genes, generating 54.03 to 58.30 Mb of raw reads per sample. After quality filtering, 44.55 (76.43%) to 45.31 (83.78%) Mb of clean reads were obtained, with an average of 6.75 Gb bases. Subsequently, the clean reads were mapped to the B. xylophilus reference genome. Approximately 77.12-79.77% of reads matched the sequences in the genome, and 62.13-66.25% of the reads were uniquely mapped. When the clean reads mapped to reference genes, the total mapping ratio was 67.39-69.71%, and the unique ratio was 58.62-60.26% (Table 2). A total of 10,187 new transcripts were detected in the transcripts of the two treatments, of which 8359 were predicted to have coding potential, 7901 were predicted to be new transcripts of known genes with protein coding potential, 458 were new transcripts of new genes with protein coding potential, and 1828 were non-coding transcripts (Table 3). Expression levels of genes and transcripts were calculated using RSEM (RNA-Seq by Expectation Maximization). There were 4240 DEGs (absolute value of log 2 fold change ≥ 1 and adjusted p-value ≤ 0.05) between the bacteria-treated and bacteria-free B. xylophilus, including 2763 downregulated and 1477 upregulated DEGs ( Figure 3). The number of no-DEGs (absolute value of log 2 fold change < 1 or adjusted p-value > 0.05) was 11,942.

qPCR Verification of Differentially Expressed Genes
Nine genes (7 upregulated genes and 2 downregulated genes) with high expression differences were randomly selected from PmBx_a and PmBx_b transcripts for qPCR verification. The expression trends of gene transcripts were consistent with the results of qPCR verification, indicating that the sequencing data were reliable ( Figure 4).

qPCR Verification of Differentially Expressed Genes
Nine genes (7 upregulated genes and 2 downregulated genes) with high expression differences were randomly selected from PmBx_a and PmBx_b transcripts for qPCR verification. The expression trends of gene transcripts were consistent with the results of qPCR verification, indicating that the sequencing data were reliable ( Figure 4).

Characterization of Differentially Expressed Transcripts
All the transcripts were further functionally characterized into three GO categories: biological process, cellular component, and molecular function. Among the set of 4240 DEGs, 356 genes were associated with at least one GO term, and a total of 50 different GO terms were found. Twenty-four GO terms were involved in biological processes, 16 GO terms were associated with cellular components, and 10 GO terms had various molecular functions ( Figure 5). Among the biological process terms, the cellular process, single-organism process, and metabolic process terms were the

Characterization of Differentially Expressed Transcripts
All the transcripts were further functionally characterized into three GO categories: biological process, cellular component, and molecular function. Among the set of 4240 DEGs, 356 genes were associated with at least one GO term, and a total of 50 different GO terms were found. Twenty-four GO terms were involved in biological processes, 16 GO terms were associated with cellular components, and 10 GO terms had various molecular functions ( Figure 5). Among the biological process terms, the cellular process, single-organism process, and metabolic process terms were the three terms containing the greatest number of DEGs. In the cellular component category, the top three GO terms were cell, cell part, and membrane, and for the molecular function category, catalytic activity, binding, and transporter activity were the top three terms ( Figure 5).
Forests 2020, 11, x FOR PEER REVIEW 9 of 20 GO terms were cell, cell part, and membrane, and for the molecular function category, catalytic activity, binding, and transporter activity were the top three terms ( Figure 5). The putative functions of DEGs were also analyzed with KEGG classification. The 871 DEGs could be mapped to 301 different pathways that might be affected by S. maltophilia. KEGG Orthology (KO) analysis showed that most DEGs were assigned to human diseases, metabolism, and organismal systems. In the human diseases class, "cancers: overview", "endocrine and metabolic diseases", and "infectious diseases: bacterial" had many DEGs. "Global and overview maps", "lipid metabolism", "carbohydrate metabolism", and "amino acid metabolism" were the dominant terms in the metabolism class. Few transcripts were assigned to "xenobiotic biodegradation and metabolism" (23 transcripts) and "metabolism of terpenoids and polyketides" (6 transcripts) ( Figure 6). This result showed that S. maltophilia mainly affected the growth and development of nematodes. In the organismal systems category, "immune system", "endocrine system", and "digestive system" had more transcripts. The top five enriched groups among the KEGG categories included signal transduction, cancers: overview, global and overview maps, cellular community-eukaryotes, and immune system. The putative functions of DEGs were also analyzed with KEGG classification. The 871 DEGs could be mapped to 301 different pathways that might be affected by S. maltophilia. KEGG Orthology (KO) analysis showed that most DEGs were assigned to human diseases, metabolism, and organismal systems. In the human diseases class, "cancers: overview", "endocrine and metabolic diseases", and "infectious diseases: bacterial" had many DEGs. "Global and overview maps", "lipid metabolism", "carbohydrate metabolism", and "amino acid metabolism" were the dominant terms in the metabolism class. Few transcripts were assigned to "xenobiotic biodegradation and metabolism" (23 transcripts) and "metabolism of terpenoids and polyketides" (6 transcripts) ( Figure 6). This result showed that S. maltophilia mainly affected the growth and development of nematodes. In the organismal systems category, "immune system", "endocrine system", and "digestive system" had more transcripts. The top five enriched groups among the KEGG categories included signal transduction, cancers: overview, global and overview maps, cellular community-eukaryotes, and immune system. A list of the top 30 pathways and number of transcripts mapped to those pathways are shown in Figure 7. The most highly active pathway was cancer (ko05200), with 72 active transcripts, followed by protein digestion and absorption (ko04974, 59 transcripts). A number of signaling pathways that are known to be involved in development and diverse functions were active in PmBx_a, such as the wnt signaling pathway (ko04310), rap1 signaling pathway (ko04015), PI3K-Akt signaling pathway (ko04151), cAMP signaling pathway (ko04024), cGMP-PKG signaling pathway (ko04022), ECMreceptor interaction, MAPK signaling pathway (ko04010), and calcium signaling pathway (ko04020). Thus, it can be concluded that S. maltophilia affects the immune system of B. xylophilus. A list of the top 30 pathways and number of transcripts mapped to those pathways are shown in Figure 7. The most highly active pathway was cancer (ko05200), with 72 active transcripts, followed by protein digestion and absorption (ko04974, 59 transcripts). A number of signaling pathways that are known to be involved in development and diverse functions were active in PmBx_a, such as the wnt signaling pathway (ko04310), rap1 signaling pathway (ko04015), PI3K-Akt signaling pathway (ko04151), cAMP signaling pathway (ko04024), cGMP-PKG signaling pathway (ko04022), ECM-receptor interaction, MAPK signaling pathway (ko04010), and calcium signaling pathway (ko04020). Thus, it can be concluded that S. maltophilia affects the immune system of B. xylophilus.

Nematode Response Gene Network
The 4240 DEGs were compared to the STRING database, 2196 of which could be mapped to 2007 protein clusters (Table S1). Cytoscape software was used to analyze and visualize these DEGs, resulting in many interaction networks (Figure 8). The genes central to a highly connected network have been shown to be more likely to have essential functions. Thus, we hypothesized that genes central to the identified network would be fundamental to the B. xylophilus response to the host pine seedlings.

Differential Expression of Virulence Genes in B. xylophilus
According to the annotation information, a series of virulence genes were found to be differentially expressed in PmBx_a and PmBx_b (Table 5). Four genes were related to cell wall hydrolysis (including the pectinase gene and cellulase gene), one of which was upregulated, and the other three were downregulated. Thirty-Six detoxification-related genes (including the cytochrome P450 gene, toxic anaphylactic protein gene, glutathione S-transferase gene, ABC transporter gene, retinol dehydrogenase gene, carboxylesterase gene, etc.) were differentially expressed, including ten upregulated and twenty-six downregulated genes. Ten cathepsin L-cysteine protease genes were differentially expressed, including five upregulated and five downregulated genes. In addition, an effector gene (BxSapB2) of PWN was upregulated.

Discussion
S. maltophilia is a kind of Gram-negative bacteria that is widely distributed in nature [38]. This species is also widely found in different PWN strains, and its relative abundance is high [15,17,19].
He et al. found that wild-type three-year-old P. massoniana inoculated with nematodes containing S. maltophilia wilted faster than those inoculated with aseptic pine wood nematodes [20]. Since wild pine carries a variety of microorganisms, four-month-old aseptic P. massoniana seedlings were used for the inoculation experiment in our study. The results provided more direct evidence that S. maltophilia can enhance the reproduction and virulence of B. xylophilus.
To explore the mechanism of the effect of S. maltophilia on PWNs in PWD, transcriptome sequencing of aseptic PWNs and those carrying S. maltophilia after inoculation was carried out. The results showed that 4240 genes were differentially expressed between the two treatments, including 1147 upregulated genes and 2763 downregulated genes. Compared with the transcriptomes of aseptic PWNs and the PWNs carrying S. maltophilia cultured on B. cinerea, which have 891 significantly DEGs, including 61 upregulated genes and 830 downregulated genes [21], the number of DEGs in the comparative transcriptome of aseptic PWNs and PWNs with S. maltophilia in pines was much higher. This suggests that different habitats could affect the gene expression of PWN, which is also reflected in the studies of Qiu et al. [39] and Tsai et al. [40].
Furthermore, the potential functions of DEGs in the comparative transcriptome of PWNs after inoculation were analyzed with gene ontology (GO) enrichment. The GO functional analysis showed that the DEGs were significantly enriched in the cell membrane, transport activity, and metabolic process. These GO classes are believed to be related to the digestion and innate immunity of nematodes [41,42]. White et al. studied the interaction between Caenorhabditis elegans and S. maltophilia and found that the DEGs were significantly enriched in these GO categories, such as components of membrane, transport, oxidation, and reduction and metabolic process [43]. It was suggested that several DEGs in the comparative transcriptome of PWNs after inoculation might be related to the immunity of PWN.
KEGG classification was also performed to estimate the functions of the DEGs. We found that the DEGs between two PWN treatments were enriched in the Wnt, Rap1, PI3K Akt, cAMP, cGMP PKG, ECM-receptor, MAPK, and calcium signaling pathways. It has been reported that ECM-receptor interaction, Wnt, PI3K Akt, and MAPK signaling pathways play immunomodulatory roles when organisms interact with bacteria or fight against other stresses [44][45][46][47]. In the interaction between Meloidogyne incognita and bacteria, the expression of many genes involved in Wnt, PI3K Akt, MAPK, and other immune-related pathways, as well as cAMP, cGMP PKG, calcium, and other pathways changed significantly [48,49]. Compared with aseptic PWNs, it could be considered that the PWNs carrying bacteria changed significantly in various immune regulatory pathways after being inoculated pines. This also meant that S. maltophilia might affect the immune process of B. xylophilus in infected P. massoniana.
According to the analysis of the protein interactions of the DEGs, a number of key genes emerged from our study, including polyubiquitin gene, leucine rich repeated serine/threonine protein kinase gene, glyceraldehyde-3-phosphate dehydrogenase gene (GAPDH), acetyl CoA carboxylase gene, and heat shock protein gene. One of these genes, Acetyl CoA carboxylase gene of B. xylophilus is related to the propagation, migration, and movement behavior of B. xylophilus [50,51]. For the HSP genes, the proteins they encode are conserved and related to biotic and abiotic stresses [52,53]. Following interference with the HSP70 gene, the reproductive ability and amount of oviposition of B. xylophilus decrease [54]. The Hsp90 gene of B. xylophilus is necessary for maintaining body temperature and ensuring reproduction in nematodes [55]. However, little is known about other HSP genes of B. xylophilus, such as the Hsp 75 and Hsp 110 genes. Similarly, the functions of the ubiquitin genes, leucine-rich repetitive serine/threonine protein kinase genes, and GAPDH genes of B. xylophilus still need to be studied and verified. The ubiquitination pathway of nematodes could help them resist the infection of pathogenic factors [56][57][58][59]. Two ubiquitin genes of B. xylophilus were also cloned, but their functions remain to be further studied [60,61]. Leucine-rich repetitive serine/threonine protein kinase genes are often associated with plant resistance [62][63][64]. The glyceraldehyde-3-phosphate dehydrogenase gene (GAPDH) is a gene with a wide range of functions, including gene expression control, membrane transport, cell signal transmission, interaction with RNA and other proteins, and avoidance of host immunity [65]. It was speculated that these key genes might affect the survival, reproduction, immunity, and many other processes of B. xylophilus, although their characterization remains to be explored.
In addition, many cell wall degrading enzyme genes, detoxification-related genes, and virulencerelated genes, for example, pectate lyase genes, cellulase genes, cytochrome p450 genes, and venom allergen-like protein genes, were differentially expressed under the two nematode treatments in our study. Cellulose, hemicellulose, and pectin are important components of the plant cell wall. B. xylophilus can secrete pectinase, and cellulase allows the colonization and migration of PWNs in pine [66][67][68][69]. Cytochrome P450 genes were associated with the vitality, reproduction, virulence, and detoxification of B. xylophilus [70,71]. The venom allergen-like protein gene was related to the migration, parasitism, and virulence of B. xylophilus [72]. The silencing of cathepsin L-cysteine protease genes in B. xylophilus reduces the feeding, reproduction, and virulence of B. xylophilus [73]. The glutathione S-transferase gene, ABC transporter gene, uridine diphosphate glucose transferase related gene, retinol dehydrogenase gene, carboxylesterase gene, etc., are also related to the reproduction and virulence of nematodes [39]. The number of these DEGs in nematodes after inoculation of pine was much greater than that in nematodes cultured on B. cinerea [21]. Therefore, it can be concluded that S. maltophilia may affect the invasion, parasitism, and virulence of B. xylophilus, whether on B. cinerea or in trees.

Conclusions
One of the most dominant bacteria in the B. xylophilus, S. maltophilia, can promote the propagation of nematodes and cause pine wilt disease faster. At the level of transcription, we found that 4240 genes were differentially expressed between aseptic B. xylophilus and B. xylophilus containing S. maltophilia after inoculation of aseptic P. massoniana. These differentially expressed genes were significantly enriched in some immune-related GO categories and KEGG pathways. The polyubiquitin-rich gene, leucine-rich repeat serine/threonine-protein kinase gene, and heat shock protein genes, etc., may play the key roles on immune resistance. In addition, S. maltophilia also make four cell wall hydrolase genes, and several detoxification-and pathogenesis-related protein genes differentially expressed. These results indicate that S. maltophilia can affect the expression of some pathogenesis-and immune-related genes, thus enhancing the virulence of B. xylophilus and aggravating the pine wilt disease.