Identifying Virulence-Associated Genes Using Transcriptomic and Proteomic Association Analyses of the Plant Parasitic Nematode Bursaphelenchus mucronatus

Bursaphelenchus mucronatus (B. mucronatus) isolates that originate from different regions may vary in their virulence, but their virulence-associated genes and proteins are poorly understood. Thus, we conducted an integrated study coupling RNA-Seq and isobaric tags for relative and absolute quantitation (iTRAQ) to analyse transcriptomic and proteomic data of highly and weakly virulent B. mucronatus isolates during the pathogenic processes. Approximately 40,000 annotated unigenes and 5000 proteins were gained from the isolates. When we matched all of the proteins with their detected transcripts, a low correlation coefficient of r = 0.138 was found, indicating probable post-transcriptional gene regulation involved in the pathogenic processes. A functional analysis showed that five differentially expressed proteins which were all highly expressed in the highly virulent isolate were involved in the pathogenic processes of nematodes. Peroxiredoxin, fatty acid- and retinol-binding protein, and glutathione peroxidase relate to resistance against plant defence responses, while β-1,4-endoglucanase and expansin are associated with the breakdown of plant cell walls. Thus, the pathogenesis of B. mucronatus depends on its successful survival in host plants. Our work adds to the understanding of B. mucronatus’ pathogenesis, and will aid in controlling B. mucronatus and other pinewood nematode species complexes in the future.


Introduction
The genus Bursaphelenchus is currently comprised of~100 recognized species, most of which are mycophagous or plant parasitic [1,2]. A number of Bursaphelenchus species, including Bursaphelenchus xylophilus (B. xylophilus) and Bursaphelenchus mucronatus (B. mucronatus), are found in damaged, dying or dead coniferous trees. These species, along with a few others, which share several morphological features and biological characteristics, are known as pinewood nematode species complexes (PWNSCs) [3]. PWNSCs that originate from different regions of the world vary in their pathogenicity, but little is known about their virulence-associated gene components or their regulation in pathogenic processes [4][5][6]. However, most studies have focused on the identification and virulence testing of PWNSCs because of the limitations of traditional biotechniques. Although a number of Int. J. Mol. Sci. 2016, 17,1492 2 of 17 studies have examined the factors contributing to virulence at the molecular level [7][8][9][10][11][12][13], the possible molecular mechanisms of PWNSC-associated virulence are still not very clear.
B. mucronatus is a PWNSC species that is a sister species to B. xylophilus, and it is widely distributed in natural pine forests throughout the northern hemisphere [14][15][16]. Both B. mucronatus and B. xylophilus parasitic on the living pine trees. When the hosts die, they begin to feed on fungi [16,17]. The nematodes could not leave their hosts independently and, thus, they are transferred to new hosts by insect vectors in natural conditions [18][19][20]. B. xylophilus complete one generation much faster than B. mucronatus, theirs life cycles having completed in 4-5 and 6-7 days, respectively [21]. Therefore, the virulence of B. mucronatus has attracted increasing attention. B. mucronatus was considered as nonpathogenic, or as having a very low virulence level only under stress conditions, although it had killed some pine seedlings in the past [16,[22][23][24]. In recent years, research showed that the virulence of B. mucronatus varies to a certain degree, such as in some isolates from Eurasia that were highly virulent to pine trees under greenhouse or outdoor conditions [25][26][27][28][29]. To our knowledge, no studies have investigated the molecular events occurring during the pathogenic processes of B. mucronatus, and thus, determining the genes and proteins involved in the pathogenic processes of B. mucronatus will help us better understand the molecular events during the pathogenic processes of B. mucronatus and other PWNSCs.
To obtain a better understanding of the genes and proteins involved in the pathogenic processes of PWNSCs, various multilevel hierarchical approaches can be applied. Current advances in molecular biological technology, especially novel high-throughput sequencing technologies, such as RNA-Seq and isobaric tags for relative and absolute quantitation (iTRAQ), may provide valid methods for detecting the genes and proteins associated with PWNSCs' virulence [30,31]. In iTRAQ, a multiplexed set of isobaric reagents are utilized to label the amine of peptides for the qualitative and quantitative analysis of different samples synchronously [32]. Furthermore, the isobaric reagent labelling increases tandem mass spectrometric (MS/MS) fragmentation, resulting in more reliable results than former techniques [32]. Additionally, RNA-Seq transcriptomic analyses were restricted to species having available reference genomes until the debut of de novo RNA-Seq technology [33,34]. The novel technology has been an extremely effective in rapidly and inexpensively obtaining abundant gene resources for many organisms, however, few researches on PWNSCs have been reported to this day [35]. Transcriptome analyses have been used to identify the differentially expressed genes (DEGs) involved in many physiologic processes [36][37][38]. Nevertheless, the analyses were unable to reveal the main proteins involved in the pathogenic process. Hence, in addition to a transcriptome analysis, we also surveyed proteins involved in the pathogenic processes of B. mucronatus by iTRAQ, which has been employed for proteomics analysis of many species [39][40][41].
The aim of the present study was to investigate the differences in both the transcriptome and proteome that occur during the pathogenic processes of two B. mucronatus (a low-virulence nematode species compared with B. xylophilus) isolates. These isolates originated from different regions and vary in their pathogenicity (Bm5, high-virulence isolate; Bm7, low-virulence compared with Bm5) [29]. The comparisons at two "omic" levels allow for a more comprehensive analysis of processes that are occurring during the pathogenesis of B. mucronatus.

Identification of the Re-Isolated Nematodes
All of the re-isolated nematodes had distinct mucro of~3-5 µm (Supplementary Figure S1) that were morphologically consistent with of the terminal mucro of B. mucronatus and previous research results [16]. Furthermore, a dendrogram was constructed based on internal transcribed spacer (ITS) sequences indicated that the re-isolated nematodes clustered with B. mucronatus (Supplementary Figure S2). Thus, the morphological characteristics combined with the phylogenetic analysis confirmed that the re-isolated nematodes in this study were B. mucronatus.

Quantitative Transcript and Protein Profiling
Bm5 (high virulence) and Bm7 (low virulence) isolates of B. mucronatus were used for RNA-Seq and iTRAQ analysis . The transcriptome sequencing resulted in 114,154,604 raw reads (55,168,960 and  58,985,644 reads for Bm5 and Bm7, respectively), and 107,097,176 clean reads were generated after filtering. The raw reads were deposited into the NCBI Sequence Read Archive (Study accession number SRP075786). After clustering the clean reads, we finally obtained 40,355 unigenes, including 17,851 clusters and 22,504 singletons. The unigenes had an average length of 1220 bp and N50 value of 2007 bp ( Table 1). All of the 40,355 unigenes were annotated by BLASTX and BLASTN algorithms to search protein databases and nucleotide databases (E-value < 0.00001). In Nr, Swiss-Prot, KEGG, COG, GO, and Nt databases, 26,623,22,507,19,776,12,180,17,041, and 9615 unigenes were aligned, respectively. For gene expression determination, we used FPKM (fragments per kilobase of transcript per million mapped reads) values to calculate the gene expression levels in the different B. mucronatus isolates, and DEGs were defined as having a false discovery rate < 0.001 and log 2 (ratio) > 1. We determined that 16,527 genes were differentially expressed between the two isolates (Supplementary Table S1). Additionally, compared with Bm7, 4919 and 11,608 genes were up-and down-regulated in Bm5 (Supplementary Table S1). Based on the iTRAQ proteomic analysis of the two isolates, 5092 proteins were identified from the 329,868 spectra, which included 16,205 peptides, of which 12,452 were unique. Of the 5092 proteins, 1463 were statistically significantly differentially expressed proteins (DEPs) with fold-changes > 1.5 and p-values < 0.05 (Supplementary Tables S2 and S3). When Bm7 served as the control, 835 of the DEPs were up-regulated and the surplus 628 were down-regulated in Bm5 (Supplementary Tables S2 and S3).

Validation of Gene Expression Data Using Real-Time Quantitative PCR (qRT-PCR)
To validate the RNA-Seq data, 12 genes were selected from the differentially abundant sequences in the transcriptome analyses and analysed by qRT-PCR. The results were consistent (r = 0.96) with the expression profiles estimated from RNA-Seq data (Figure 1), indicating the dependability of the RNA-Seq data.

Gene Ontology (GO) Analysis of Differentially Expressed Transcripts and Proteins
Gene ontology (GO) was employed to categorize the functions of the DEGs and DEPs in the two B. mucronatus isolates. All of the DEGs and DEPs were categorized into three groups: biological process, cellular component, and molecular function. In the transcriptomic analysis, 17,515 GO terms were assigned to 7454 transcripts, and a single transcript might be annotated to several GO terms. Among the 17,515 GO terms, 69% fell into "biological process" group, 11% into "cellular component" group, and 20% into "molecular function" group. "Single-organism" was the largest category in the "biological processes" group, occupying 11%. "Cell" and "cell parts" occupied nearly half the "cellular components" group. In the group "molecular function", 41% were involved in "binding" and 35% in "catalytic activity". The percentages of the GO categories in each of the three groups are shown for up-regulated transcripts ( Figure 2) and down-regulated transcripts ( Figure 3). In the proteomic analysis, 586 up-regulated proteins were annotated to 3719 GO terms while 444 down-regulated proteins were annotated to 2831 GO terms. As with transcripts, a single protein might be annotated to several GO terms. The "biological process" group of up-regulated proteins contained 19 categories; "single-organism process" was the largest, accounting for 10%, followed by 9% each for "cellular process", "developmental process", "metabolic process", and "multicellular organismal process". In the group "cellular component", "cell", "cell part" and "organelle" accounted for 23%, 23%, and 17%, respectively. Moreover, "molecular function" had two major categories, "catalytic activity" and "binding", which accounted for 45% and 43%, respectively ( Figure 2). Meanwhile, the component and percentages of the categories in the three groups for down-regulated proteins were like that of the up-regulated proteins ( Figure 3). while 444 down-regulated proteins were annotated to 2831 GO terms. As with transcripts, a single protein might be annotated to several GO terms. The "biological process" group of up-regulated proteins contained 19 categories; "single-organism process" was the largest, accounting for 10%, followed by 9% each for "cellular process", "developmental process", "metabolic process", and "multicellular organismal process". In the group "cellular component", "cell", "cell part" and "organelle" accounted for 23%, 23%, and 17%, respectively. Moreover, "molecular function" had two major categories, "catalytic activity" and "binding", which accounted for 45% and 43%, respectively ( Figure 2). Meanwhile, the component and percentages of the categories in the three groups for down-regulated proteins were like that of the up-regulated proteins ( Figure 3).    while 444 down-regulated proteins were annotated to 2831 GO terms. As with transcripts, a single protein might be annotated to several GO terms. The "biological process" group of up-regulated proteins contained 19 categories; "single-organism process" was the largest, accounting for 10%, followed by 9% each for "cellular process", "developmental process", "metabolic process", and "multicellular organismal process". In the group "cellular component", "cell", "cell part" and "organelle" accounted for 23%, 23%, and 17%, respectively. Moreover, "molecular function" had two major categories, "catalytic activity" and "binding", which accounted for 45% and 43%, respectively ( Figure 2). Meanwhile, the component and percentages of the categories in the three groups for down-regulated proteins were like that of the up-regulated proteins ( Figure 3).

KEGG Pathway Enrichment Analysis of Differentially Expressed Genes and Proteins
To investigate the activated biological pathways between the two B. mucronatus isolates, DEGs and DEPs were assigned to the reference pathways in KEGG. As a result, all of the 16,527 DEGs were mapped to 254 KEGG pathways. Sixty-three of these were significantly enriched (p-value < 0.05), including "calcium signalling pathway", "salivary secretion" and "spliceosome" (Table 2). In the proteome pathway analysis, 1463 DEPs were assigned to 217 pathways. Among them, 38 pathways were significantly enriched (p-value < 0.05), such as "fatty acid metabolism", "drug metabolism-cytochrome P450", and "α-linolenic acid metabolism" ( Table 3). The KEGG pathway enrichment analysis also found some enriched peptides-such as cytochrome P450s, glucuronosyltransferases, glutathione S-transferases, glucosyltransferases, and uridine 5′-diphosph-glucuronosyltransferases (UDP-glucuronosyltransferases, UGT)-that were involved in xenobiotic metabolism and helped nematodes adapt to life on pine hosts.

KEGG Pathway Enrichment Analysis of Differentially Expressed Genes and Proteins
To investigate the activated biological pathways between the two B. mucronatus isolates, DEGs and DEPs were assigned to the reference pathways in KEGG. As a result, all of the 16,527 DEGs were mapped to 254 KEGG pathways. Sixty-three of these were significantly enriched (p-value < 0.05), including "calcium signalling pathway", "salivary secretion" and "spliceosome" (Table 2). In the proteome pathway analysis, 1463 DEPs were assigned to 217 pathways. Among them, 38 pathways were significantly enriched (p-value < 0.05), such as "fatty acid metabolism", "drug metabolism-cytochrome P450", and "α-linolenic acid metabolism" ( Table 3). The KEGG pathway enrichment analysis also found some enriched peptides-such as cytochrome P450s, glucuronosyltransferases, glutathione S-transferases, glucosyltransferases, and uridine 5 -diphosph-glucuronosyltransferases (UDP-glucuronosyltransferases, UGT)-that were involved in xenobiotic metabolism and helped nematodes adapt to life on pine hosts.

Correlation Analysis of Protein and RNA Expression
To examine the correlations between protein and mRNA expression levels, the proteins obtained from the B. mucronatus iTRAQ analysis were compared with the corresponding transcripts from the RNA-Seq analysis. When we matched all of the proteins with their detected transcripts, a low correlation coefficient of r = 0.138 was found, indicating that the transcriptomic and proteomic data shared a low identity (Supplementary Figure S3). In the proteomic analysis, 1463 DEPs were identified, 795 of which (54.3%) had corresponding genes in the transcriptomic analysis. Then, we detected whether the changes between proteins and transcripts were in the same direction, and found that only 135 DEPs (17%) were accordant (Supplementary Table S4). The present results indicate a low correlation between the transcript and protein abundances, probably because post-transcriptional regulation occurred in the pathological process.

Analysis of Pathogenesis-Related Genes
To date, many genes involved in the pathogenic processes of plant parasitic nematodes have been identified. Little is known as to whether these virulence-associated genes are related to the pathogenic processes of PWNSCs, such as B. mucronatus. We, therefore, used the BLAST algorithm to query the DEPs against the NCBI protein database and Uniprot, and five proteins related to pathogenic processes in other plant parasitic nematode species were obtained. Three of them, peroxiredoxin (TPX), fatty acid-and retinol-binding protein (FAR) and glutathione peroxidase (GPX), are related to resistance against plant defence responses, while the other two, β-1,4-endoglucanase (ENG) and expansin (EXP), are associated with the nematode breakdown of plant cell walls. In addition, all of them were detected to have corresponding genes in iTRAQ analysis based on their amino acid sequences and regulatory changes (Table 4). However, most important, all five virulence-associated proteins have higher expression levels in the high-virulence isolate Bm5 than in the low-virulence isolate Bm7.

Characterization and Analysis of RNAi in B. mucronatu
To confirm their involvement in B. mucronatus' virulence, FAR (CL4567.Contig1) and ENG (CL5080.Contig2) were selected for RNAi construction because they were highly expressed at both mRNA and protein levels in the highly virulent B. mucronatus isolate (Bm5). Once the RNAi were constructed and introduced into the nematode, qPCR was used to assess the FAR and ENG mRNA expression levels in different treatments. By comparing the expression levels of the FAR gene in the nematodes exposed to the target dsRNA-soaking solution and controls, it was shown that the RNAi procedure caused a significant reduction in the target mRNA's expression level (Figure 4a). Taking the mRNA expression level in M9 medium (non-dsRNA treated) as 100%, the average expression levels of samples with the green fluorescent protein (GFP) dsRNA treatment was 99.8%, while with the target dsRNA treatment it was only 26.9%. The expression levels of the ENG gene produced similar results after soaking in an ENG dsRNA solution (Figure 4b). The number of dsRNA-treated nematode offspring was used to assess the influence of RNAi on the development and propagation of B. mucronatus. The number of offspring from the dsRNA-treated nematodes was less than that of non-dsRNA-treated nematodes (Figure 4c). In the non-dsRNA treatment, an average of 1239 individuals in the F1 generation were produced by 20 pairs of females and males, but an average of only 945 and 985 offspring were produced in the dsRNA treatments of FAR and ENG, respectively. The difference between the target dsRNA treatment and M9 was highly significant (p-value < 0.001). Treatment with GFP dsRNA soaking resulted in an average of 1133 individuals in the F1 generation, which is a little less than that of non-dsRNA treatment, but not highly significant. However, it is highly significantly different (p-value < 0.001) from the target dsRNA treatment. Thus, the knockdowns of the two target genes had marked effects on the development and propagation of B. mucronatus. Additionally, a possible dsRNA toxicity to the nematodes was not obvious at a concentration of 800 ng/µL. To assess the influence of RNAi on the virulence of B. mucronatus, the dsRNA-treated nematodes were inoculated on two-year-old Pinus thunbergii (P. thunbergii) seedlings under greenhouse conditions. The first wilted seedling infected with ENG dsRNA-, GFP dsRNA-, and M9-treated nematodes were observed 5 weeks after inoculation, and the first wilted FAR dsRNA-treated seedling was observed 6 weeks after inoculation. All of the treatments of B. mucronatus caused some seedlings to wilt (ranging from 4 to 10 seedlings); while the FAR dsRNA-treated nematodes did not cause as many wilted seedlings as the others, and the ENG dsRNA-treated nematodes caused a slower wilting of seedlings than GFP dsRNA-and M9-treated nematodes (Figure 4d). Not a single seedling was wilted in the uninoculated control group. In addition, nematodes were re-isolated from each wilted seedling.
dsRNA-treated nematodes were inoculated on two-year-old Pinus thunbergii (P. thunbergii) seedlings under greenhouse conditions. The first wilted seedling infected with ENG dsRNA-, GFP dsRNA-, and M9-treated nematodes were observed 5 weeks after inoculation, and the first wilted FAR dsRNA-treated seedling was observed 6 weeks after inoculation. All of the treatments of B. mucronatus caused some seedlings to wilt (ranging from 4 to 10 seedlings); while the FAR dsRNA-treated nematodes did not cause as many wilted seedlings as the others, and the ENG dsRNA-treated nematodes caused a slower wilting of seedlings than GFP dsRNA-and M9-treated nematodes (Figure 4d). Not a single seedling was wilted in the uninoculated control group. In addition, nematodes were re-isolated from each wilted seedling.

Discussion
PWNSC species are found in wilted coniferous trees, including B. xylophilus, B. mucronatus and a few other Bursaphelenchus species, causing enormous ecological and economic losses worldwide [3,6,27,42]. However, the molecular mechanisms behind the pathogenic processes of these PWNSC species are still poorly understood. B. mucronatus isolates are known to be vary in their level of pathogenicity [25][26][27][28][29]. With the recent advancements in RNA-Seq approach, we report the first study to investigate transcriptomic data from B. mucronatus isolates gained from wilted trees and we present proteomic data obtained by the iTRAQ technology. Furthermore, when the transcriptome and proteome data were coupled, we were capable of making, for the first time,

Discussion
PWNSC species are found in wilted coniferous trees, including B. xylophilus, B. mucronatus and a few other Bursaphelenchus species, causing enormous ecological and economic losses worldwide [3,6,27,42]. However, the molecular mechanisms behind the pathogenic processes of these PWNSC species are still poorly understood. B. mucronatus isolates are known to be vary in their level of pathogenicity [25][26][27][28][29]. With the recent advancements in RNA-Seq approach, we report the first study to investigate transcriptomic data from B. mucronatus isolates gained from wilted trees and we present proteomic data obtained by the iTRAQ technology. Furthermore, when the transcriptome and proteome data were coupled, we were capable of making, for the first time, integrated and precision measurements of gene and protein expression patterns in B. mucronatus isolates from wilted trees in the field.
Constructing an RNA-Seq transcriptome is a novel approach we have used to identify virulence-associated genes in B. mucronatus isolates having different virulence levels during the pathogenic processes. RNA-Seq approach has been successfully applied to detect common and different molecular idiosyncrasies between B. xylophilus and B. mucronatus, when they were cultured on fungal colonies [35]. Similar to that report, RNA-Seq data revealed a remarkable change in gene expression levels between the highly and lowly virulent isolates, indicating that the RNA-Seq is a high-efficiency technique for detecting different molecular idiosyncrasies among nematode isolates. Many proteins referred to xenobiotic biodegradation pathways were enriched in highly virulent B. mucronatus isolates' transcriptomes based on the KEGG annotation (Supplementary Table S5), which was similar to the results of a previous study [35]. Recently, transcriptome and genome sequencing have been conducted on strongly virulent and weakly virulent B. xylophilus isolates to investigate the molecular differences between them. More than 3 Mb of genetic variations were investigated as single-nucleotide polymorphisms (SNPs) or small indels between different Japanese B. xylophilus isolates, and the genetic diversity was involved with their virulence and ecological characteristic differences [12]. Meanwhile, a number of transcripts and exons revealed remarkable differences between the strongly and weakly virulent Chinese B. xylophilus isolates, and functional study of those differences suggested that various virulence B. xylophilus isolates showed differences in their development, propagation, and oxidoreductase activities [43]. Similarly, there were considerable differences in strongly and weakly virulent isolates of both nematode species. It is known that B. mucronatus is a low-virulence nematode species compared with B. xylophilus, but the two Bursaphelenchus species have evolved similar parasitic mechanisms for living on host trees [35], which could perhaps be explained by the fact that the developmental and movement speed of B. xylophilus is much faster than B. mucronatus [21,44,45].
Meanwhile, in the proteomic profiling, we also found some enriched proteins that are involved in xenobiotic biodegradation, such as cytochrome P450s, glucuronosyltransferases, glutathione S-transferases and UDP-glucuronosyltransferases, based on KEGG pathway analysis (Supplementary Table S6). A global correlation analysis between iTRAQ and RNA-Seq data revealed a low correlation in B. mucronatus during the pathogenic processes (Supplementary Figure S3), indicating that post-transcriptional regulation may be involved in the pathogenic processes. Interestingly, in the vast majority of organisms that have been examined in recent years, transcript abundance only partially predicted protein abundance, suggesting that after experimental errors have been removed, other modes of regulation must be employed to determine the protein abundances within cells [46]. First, mRNAs are less stable, having an average half-life of less than 7 h, while that of proteins is approximate 46 h [47]. Second, experiments in living cells found that variations in protein abundance levels are primarily determined by the regulation of translation and protein degradation [47,48]. Therefore, these post-transcriptional regulations contribute to the variation in protein abundance, leading to a low correlation between the iTRAQ and RNA-Seq data.
Proteins direct the work of living cells, therefore, we used the BLAST algorithm to query the DEPs against the NCBI protein and Uniprot databases. Five DEPs were homologous to proteins with previously identified roles in general parasitic nematode pathogenic processes (Table 4). More importantly, all five of these proteins have higher expression levels in the highly virulent isolate than in the weakly virulent isolate. In addition, all five were detected to have corresponding transcripts in the RNA-Seq analysis, based on their amino acid sequences and read direction (Table 4). They were divided into two categories based on their functions: proteins related to resistance against plant defence responses and those that were associated with breaking down plant cell walls. To overcome host defence responses and survive in different tissues of the host, parasitic nematodes have evolved an ingenious strategy to manipulate the host's metabolism to their own advantage. Three DEPs related to resistance against host defence responses were found to be more highly expressed in the highly virulent isolates than in the weakly virulent isolate of B. mucronatus. First, the FAR protein facilitates parasitic nematode infections by removing and transporting fatty acids, which are necessary for the development and reproduction of the nematodes [49]. Furthermore, the FAR protein also takes part in disturbing intra-and inter-cellular lipid signals which are related to host defence responses [50]. Second, GPX functions in parasitic nematodes by either removing hydrogen peroxide or transforming peroxidised fatty acids, which related to inhibition of lipid peroxidation [51]. It may take part in active oxygen from the in vivo metabolism, whereas the protein protects nematodes from host defence responses [52]. Finally, TPX particularly metabolises hydrogen peroxide and may also take part in protecting nematodes from host defences [53]. For a plant parasitic nematode that completes most of its life cycle within the plant, the plant cell wall is supposed to be its main obstacle, and the ability to get through the wall is, thus, quite important. ENG and EXP proteins, which can modify plant cell walls, were also more highly expressed in the highly virulent isolates of B. mucronatus. ENG, which degrades polysaccharides possessing β-1,4-glucan backbones (for instance, cellulose and xyloglucan) presumably facilitates nematode migration through plant tissue, and has been found among various plant parasitic nematodes [8,54,55]. EXP rapidly induces extensions of plant cell walls by melting the network of wall polysaccharides to allow turgor-driven cell enlargement [56]. However, plant parasitic nematodes can also produce functional expansin that is used to loosen cell walls during host plant invasion [57].
The ability to inactivate a target gene transiently by RNAi was first found in the nematode Caenorhabditis elegans [58]. Since then, RNAi has been a fairly effective strategy to provoke sequence-specific silencing in most eukaryotic cells [59]. In the present study, the RNAi strategy was employed to assess the influence of RNAi knockdowns of FAR and ENG genes in B. mucronatus. The inhibition of the seedling death rate may be explained by the ENG knockdown resulting in a reduced ability to locate and invade plant cells, and a decreased number of established nematodes [60,61]. The results of ENG gene RNAi in B. mucronatus are in general agreement with the previous report of its knockdown in B. xylophilus [59]. FAR is a member of the nematode-specific FAR family of proteins, and it exists in animal-and plant-parasitic nematodes [49,50]. After B. mucronatus was soaked in the FAR dsRNA solution, the cumulative mortality rate of its infected seedlings was sharply reduced compared with the other treatments (Figure 4d). Similar results were also found for the tomato root-knot nematode Meloidogyne javanica [50]. FAR facilitates invasion by manipulating host lipid-based defences, and this is a key component for the nematodes' successful parasitism [62].

Nematode Sample Preparation
Bm5 and Bm7 isolates used in this study were originally isolated from wilted Pinus massoniana trees in different towns of China by the Baermann funnel technique [63] and have been maintained for over 6 years in the Institute of Forest Protection Lab in Nanjing Forestry University (Nanjing, China), have significantly different in virulence levels according to our previous studies [29]. Both isolates are from inbred lines that were established by sister-brother mating (full-sib mating repeated 10 times) and were used to generate materials for subsequent experiments. Nematodes isolates were cultured on the Botrytis cinerea colonies of potato-dextrose agar plates at 25 • C for 7 days and then separated for inoculation. Then, each 12-year-old P. thunbergii tree was inoculated with~1 mL of nematode suspension (~10,000 of mixed-stage nematodes). The treated trees were grown under natural conditions, observed once a week, and cut when wilted. The mixed-stage nematodes were extracted from the small pieces of the wilted trees using Baermann funnel technique (place at 25 • C for 3 h) and cleaned by sucrose flotation, then rinsed three times in 1× poly(butylene succinate-co-butylene terephthalate) (PBST) [8,64], and promptly immersed in liquid nitrogen and deposited at −80 • C. In addition, these re-isolated nematodes were identified based on morphological and molecular analyses [65,66].

RNA Preparation and Transcriptome Sequencing
The RNA-Seq experiments were conducted with the assistance of the Beijing Genomics Institute (BGI, Shenzhen, China). For RNA-seq, total RNAs were extracted from the frozen samples using the RNAprep kit (Tiangen, Beijing, China) and then purified further with the RNA clean kit (Tiangen, Beijing, China) according to the manufacturer's instructions. RNA was quantified with an Agilent 2100 Bioanalyzer RNA Nanochip (Agilent Technologies, Inc., Waldbronn, Germany). Finally, the SolexaHiSeq™ 2000 platform was employed for sequencing according to the manufacturer's instructions (Illumina, San Diego, CA, USA).

Redundant Data Filtering and de Novo Assembly
Raw reads obtained from the RNA sequencing machines contained redundant reads that comprised of adapters, and unknown or inferior quality bases. These data have an adverse impact on the bioinformatics analysis and, thus, they were dumped using the following steps: (I) removed reads with adaptors; (II) removed reads with unknown nucleotides larger than 5%; (III) removed low-quality reads; and (IV) isolated the clean reads. Unigene de novo assembly was conducted with the short clean reads-assembling program Trinity [67].

Unigene Annotation and Classification
The unigenes annotation and classification methods and software were conducted as reported in the previous research [68]. For each isolate, the homologous gene expression ratios were counted based on the RPKM (reads per kilobase of transcript per million mapped read). The unigenes with two-fold changes between Bm5 and Bm7 were identified as DEGs.

qRT-PCR
To ensure the reliability of the RNA-Seq data, we selected 12 genes having different abundances in the transcriptome analysis, and validated the data by performing qRT-PCR. The β-tubulin gene (β-tub, AB500151) was used as the reference gene, and the Primer Premier 5 software (Premier Biosoft International, Palo Alto, CA, USA) was employed to design gene-specific primers for the target and reference genes (Supplementary Table S7). The total RNAs were extracted as described above and subjected to reverse transcription using the TransScript ® One-Step gDNA Removal and cDNA Synthesis SuperMix Kit (TransGen Biotech, Beijing, China). qRT-PCR was performed with the SYBR Green Master Mix (Vazyme, Nanjing, China) using the StepOne Plus Real-time PCR System (Applied Biosystems, Foster City, CA, USA). All of the qPCR experiments were performed with three replicates.

Protein Preparation and iTRAQ Labelling
The iTRAQ experiments were conducted with the help of the Beijing Genomics Institute (BGI, Shenzhen, China). The frozen nematodes were suspended in the Lysis buffer and sonicated in the ice. The proteins were reduced with 10 mM DL-dithiothreitol solution (DTT) at 56 • C for 60 min and then alkylated by 55 mM iodoacetamide under dark conditions for 60 min. After that, the reduced and alkylated protein mixtures were precipitated by adding 4× volume of chilled acetone at −20 • C overnight. After centrifugation at 30,000× g at 4 • C, the pellet was dissolved in 0.5 M tetraethylammonium bromide (Applied Biosystems, Milan, Italy) and sonicated in ice. After centrifuging at 30,000× g at 4 • C, an aliquot of the supernatant was taken to measure the protein concentration by Thermo NanoDrop 2000C (Thermo Fisher Scientific Inc., Wilmington, DE, USA). Total protein (100 µg) samples were removed, and then the proteins were digested with Trypsin Gold (Promega, Madison, WI, USA) at the protein: trypsin ratio of 30:1 at 37 • C for 16 h. After the digestion, peptides were dried and reconstituted in 0.5 M tetraethylammonium bromide. Finally, the reconstituted samples were labelled by the 8-plex iTRAQ reagent (Applied Biosystems, Foster City, CA, USA) according to the manufacturer's instruction.

Strong Cation Exchange (SCX) Fractionation and LC-ESI-MS/MS Analysis
Strong Cation Exchange (SCX) fractionation and LC-ESI-MS/MS analysis were performed as reported in the previous articles [69,70].

Proteomic Data Analysis
The iTRAQ data analyses were conducted as reported in the previous research [71]. Proteins with a 1.5-fold change between Bm5 and Bm7 samples and p-value <0.05 were identified as DEPs.

Construction and Analysis of RNAi
In the present study, the RNAi strategy was used to evaluate the influence of RNAi knockdowns of two genes, FAR (CL4567.Contig1) and ENG (CL5080.Contig2), on B. mucronatus. GFP was used as non-endogenous control gene, and gene-specific primers for the target and control genes were designed using Primer Premier 5 software (Supplementary Table S8). In addition, M9 (42.3 mM Na 2 HPO 4 , 22 mM KH 2 PO 4 , 85.6 mM NaCl, and 1 mM MgSO 4 , pH 7.0) was used as a non-dsRNA-treated control. The dsRNAs for RNAi were synthesized using the MEGscript ® RNAi Kit (Ambion Inc., Austin, TX, USA) according to the manufacturer's instructions (FAR dsRNA is 588 bp, ENG dsRNA is 441 bp, and GFP dsRNA is 670 bp). The RNAi-soaking approach utilized in the present study was conducted according to previous reports [10,59]. The RNAi efficacy was detected by qRT-PCR, and gene function was assessed by comparative fecundity and virulence assessments between the wild-type and the RNAi-treated nematodes. All of the experiments were conducted with three replicates.

Conclusions
In conclusion, we conducted transcriptomic and proteomic analyses as a complementary approach to investigate B. mucronatus isolates during the pathogenic processes. Lots of DEGs and DEPs were identified which are potentially involved in the pathogenic processes of B. mucronatus. Interestingly, the highly virulent isolate had higher expression levels of TPX, FAR, GPX, ENG, and EXP compared with the weakly virulent isolate. These genes are critical factors in B. mucronatus that are associated with breaking down plant cell walls and resisting plant defence responses in pine trees, therefore, we hypothesize that the pathogenesis of B. mucronatus is a result of their survival in the host tree. It is in accordance with the hypothesis which has been proposed in B. xylophilus researches [8,9,13]. This study increases our understanding of B. mucronatus pathogenesis, and may be useful in the future for controlling B. mucronatus and other PWNSC species.