Comparative transcriptome analysis of pinus densiflora following inoculation with pathogenic (bursaphelenchus xylophilus) or non-pathogenic nematodes (b. Thailandae)

ABSTRACT _Pinus densiflora_ (Korean red pine) is a species of evergreen conifer that is distributed in Korea, Japan, and China, and of economic, scientific, and ecological importance. Korean

red pines suffer from pine wilt disease (PWD) caused by _Bursaphelenchus xylophilus_, the pinewood nematode (PWN). To facilitate diagnosis and prevention of PWD, studies have been conducted

on the PWN and its beetle vectors. However, transcriptional responses of _P. densiflora_ to PWN have received less attention. Here, we inoculated Korean red pines with pathogenic _B.

xylophilus_, or non-pathogenic _B. thailandae_, and collected cambium layers 4 weeks after inoculation for RNA sequencing analysis. We obtained 72,864 unigenes with an average length of 869 

bp (N50 = 1,403) from a Trinity assembly, and identified 991 differentially expressed genes (DEGs). Biological processes related to phenylpropanoid biosynthesis, flavonoid biosynthesis,

oxidation–reduction, and plant-type hypersensitive response were significantly enriched in DEGs found in trees inoculated with _B. xylophilus_. Several transcription factor families were

found to be involved in the response to _B. xylophilus_ inoculation. Our study provides the first evidence of transcriptomic differences in Korean red pines inoculated with _B. xylophilus_

and _B. thailandae_, and might facilitate early diagnosis of PWD and selection of PWD-tolerant Korean red pines. SIMILAR CONTENT BEING VIEWED BY OTHERS TRANSCRIPTOME ANALYSIS REVEALS THE

MOLECULAR MECHANISMS OF RESPONSE TO AN EMERGENT YELLOW-FLOWER DISEASE IN GREEN CHINESE PRICKLY ASH (_ZANTHOXYLUM SCHINIFOLIUM_) Article Open access 23 September 2021 IDENTIFICATION,

CLASSIFICATION, AND CHARACTERIZATION OF AP2/ERF SUPERFAMILY GENES IN MASSON PINE (_PINUS MASSONIANA_ LAMB.) Article Open access 08 March 2021 TRANSCRIPTOMIC PROFILING OF _POA PRATENSIS_ L.

UNDER TREATMENT OF VARIOUS PHYTOHORMONES Article Open access 15 March 2024 INTRODUCTION Pines are conifers in the genus _Pinus_ that are found in the Northern Hemisphere1. They have

economic, ecological, and scientific importance, as they provide timber for construction, furniture, paneling, and flooring2,3, habitats and food for wildlife, and (in their needles) agents

with anticancer, antioxidant, and antimutagenic properties4. However, pine trees (such as _Pinus densiflora_, _P. thunbergii_, and _P. koraiensis_) are threatened by pine-wilt disease (PWD),

a devastating disease that kills trees within a few weeks to a few months from infection, with symptoms that are characterized by wilted and brown-colored needles5,6. PWD is caused by the

pinewood nematode (PWN) _Bursaphelenchus xylophilus_, which is carried by long-horned beetles (_Monochamus_ spp.) and spread when the beetles feed on the trees5. _B. thailandae_ is a

nematode species that was first isolated from pine trees in Thailand, and which has subsequently been detected in Korea. _B. thailandae_ differs morphologically from _B. xylophilus_, and is

not pathogenic, so pine trees inoculated with _B. thailandae_ can survive without showing any visible symptoms of PWD7,8. Therefore, examination of the transcriptome differences between _B.

xylophilus_ and _B. thailandae_ inoculated pines might identify novel resistance mechanisms and candidate genes that play important roles in resistance specifically against _B. xylophilus_.

Furthermore, we also might identify the pathogenesis mechanisms of _B. xylophilus_ to cause PWD. However, comparative transcriptome analysis between trees inoculated with nematodes having

different pathogenicity has not been conducted. Many studies have been carried out to find ways to diagnose and block PWD, and most of them have been focused on _B. xylophilus_ and its

beetle vectors9,10,11. Currently, to diagnose PWD, _B. xylophilus_ or DNA fragments from _B. xylophilus_ must be detected in tree samples12, or PWD symptoms must be observed. However, it is

difficult to detect _B. xylophilus_ and its DNA in pine trees at an early stage of infection, and by the time PWD symptoms can be observed, _B. xylophilus_ has generally already spread

throughout the forest. The physiological symptoms of PWD have been well characterized, but there are only few studies that attempt to understand the comprehensive transcriptome of pine trees

in response to _B. xylophilus_ infection. Previously, we examined the transcriptome differences between trees with and without PWD symptoms in natural forest. However, we didn’t know

whether trees showing PWD symptoms are indeed infected by PWN or not13. In this report, we inoculated trees with PWN and observed the PWD symptoms after inoculation of PWN. Therefore, we can

assure that PWD symptoms are indeed caused by PWN inoculation, and examination of the transcriptome differences among the trees injected by water, non-pathogenic nematode, and pathogenic

pine wood nematode give better understanding of the transcriptional responses against PWD in Korean red pines. In the thale cress _Arabidopsis thaliana_, thousands of genes have been shown

to be differentially expressed upon nematode infection, and important regulators of defense responses against nematode infection have been discovered by transcriptomic analysis14,15. By

extension, it might also be possible to identify genes that are differentially expressed in pine trees upon infection with _B. xylophilus_, and thereby to diagnose PWD infection by

gene-expression analysis. In addition, it might be possible to develop PWD-resistant trees by selection of particular alleles of genes that are differentially expressed after inoculation

with _B. xylophilus_. Despite the importance of developing a comprehensive understanding of the transcriptome of pines in conditions of _B. xylophilus_ infection, there is a lack of

transcriptome analysis and most of them were conducted using saplings under artificial experimental conditions instead of using adult trees in a forest environment16,17. Next-generation

sequencing (NGS) technology has been rapidly developed and widely used for research in plant biology, to enhance our understanding of plant responses under various conditions18. In addition,

software developments now enable the _de novo_ assembly of the transcriptome of an organism (such as _P. densiflora_) that does not have a reference transcriptome19. Here, we report the _de

novo_ assembly of the _P. densiflora_ transcriptome and quantification of transcript expression in response to inoculation with either _B. xylophilus_ or _B. thaliandae_ at felling age in a

natural forest environment. Identification of differentially expressed transcripts in pathogenic PWN-inoculated _P. densiflora_ could spur the development of a diagnostic method for PWD

infection and aid in the selective breeding of PWD-resistant _P. densiflora_. MATERIALS AND METHODS PLANT MATERIALS AND INOCULATION WITH PWN Pathogenic PWN (_B. xylophilus_) and

non-pathogenic nematode (_B. thailandae_) were originally isolated from Korean red pines and reared on fungal hyphae of _Botrytis cinerea_ (de Bary) Whetzel grown on potato dextrose agar

medium at 25 °C for 2 weeks. Nematodes were re-isolated from the medium by the Baermann funnel method20. Nine Korean red pines of 11–13 m height and 15–20 cm diameter at breast height in a

forest in Jinju-si, Gyeongsangnam-do province, South Korea, were selected. Water, _B. xylophilus_, and _B. thailandae_ were injected into three independent trees each. The stem of each tree

at breast height was wounded mechanically in three places, and 1 ml sterile water or 1 ml sterile water containing 20,000 nematodes was injected into each wound site (for a total of 60,000

PWNs per tree). Cambium samples of the trees were collected 4 weeks after inoculation and subjected to RNA-Seq analysis. Briefly, hard outer bark was removed and soft cambium layers were

taken from the main stem at breast height by using chisel. RNA EXTRACTION, CDNA LIBRARY PREPARATION, AND SEQUENCING Cambium samples were taken from the main stem at breast height by using

chisel. Total RNA was isolated from cambium samples with an RNA isolation kit (TAESIN Bio Co., Seoul, South Korea). The assessment of RNA integrity (RIN), library construction, and

sequencing were performed as described previously13. Briefly, RNA quality was determined with a 2100 Bioanalyzer (Agilent, Santa Clara, CA, USA), and only samples with an RNA integrity

number >8 were used for library preparation. Preparation of each paired-end non-directional cDNA library (2 × 101 bp) was conducted according to the TruSeq RNA Sample Preparation Guide

(Illumina, San Diego, CA, USA). Sequencing of cDNA libraries was performed on an Illumina HiSeq. 2000 sequencer. _DE NOVO_ TRANSCRIPTOME ASSEMBLY OF NEMATODE INOCULATED _P. DENSIFLORA_

PRINSEQ-lite v0.20.4 was used for read cleaning as described in Lee _et al_. (2018) with minor modification (filtering of sequences <50 bp length and eliminating exact duplicates or

reverse-complement exact duplicates caused by library PCR amplification with –derep 14 option)13,21. Reference transcriptome were generated from all clean reads using Trinity v2.5.1 with

default parameters19 and _de novo_ assembled transcriptome is available as Supplementary Data S1. Candidate coding regions in the all assembled transcripts were identified using TransDecoder

v5.3.0 with default parameters (two steps such as extraction of the long ORFs and prediction of the likely coding regions)22. Clustering of transcripts was performed using CD-HIT-EST v4.6.1

with default parameters23,24, and the longest transcripts in each cluster were used for subsequent analysis. Transcriptome completeness was assessed by using Benchmarking Universal

Single-Copy Orthologs (BUSCO) v3 with the Embryophyta_(odb10) database25. QUANTIFICATION OF THE EXPRESSION OF TRANSCRIPTS AND IDENTIFICATION OF DIFFERENTIALLY EXPRESSED TRANSCRIPTS Clean

paired end reads were mapped to the reference transcriptome using Bowtie software26. Read counts were obtained using RSEM v1.3.027. The raw counts were normalized as trimmed mean of M-values

(TMM)-normalized transcripts per kilobase million (TPM) values for each transcript and differentially expressed genes (DEGs) showing more than 2-fold expression change with a

false-discovery rate (FDR)-adjusted _P_-value ≤ 0.05 among all pairwise sample comparisons were obtained using EdgeR v3.16.5, a component of Trinity19,28. GENE ANNOTATION, GENE ONTOLOGY

ENRICHMENT, AND MAPMAN ANALYSIS OF DEGS Translated protein sequences corresponding to assembled unigenes were compared with those of _A. thaliana_, with an E-value ≤ 1E-7 using the Basic

Local Alignment Search Tool for proteins (BLASTX) with default parameters except for max_target_seqs 129 and annotation result was presented as Supplementary Data S2. Putative _P.

densiflora_ transcription factors (TFs) that aligned and annotated with _A. thaliana_ TFs were classified into TF families based on the Plant Transcription Factor Database v4.0

(http://planttfdb.cbi.pku.edu.cn/)30. Gene ontology (GO) enrichment analysis was performed with PANTHER GO classification system (http://www.geneontology.org)31. Gene ontology information

files from The Arabidopsis Information Resource (TAIR) were used as a reference32. Biological process gene ontology (GOBP) with a Fisher’s exact test with FDR corrected _P_-value va0.05 were

considered to be significantly enriched. REVIGO (http://revigo.irb.hr) were used for reducing and visualizing the genes ontology33. For MapMan analysis, homologous _A. thaliana_ IDs and

log2 fold-change values of TMM normalized TPM (_B. xylophilus_ versus _B. thailandae_) were mapped to biotic-stress pathways34. Pictorial representation for the biotic-stress pathways was

downloaded from the MapMan website (https://mapman.gabipd.org/)34. CORRELATION ANALYSIS BETWEEN QUANTITATIVE PCR AND NGS DATA Preparation of total RNA from Korean red pines was performed

using RNA isolation kit (TAESIN), and quantitative reverse-transcription PCR (RT-qPCR) was conducted using the TOPreal One-step RT qPCR Kit (Enzynomics, Daejeon, South Korea). 12 DEGs in

trees inoculated with _B. xylophilus_ relative to those inoculated with _B. thailandae_ were used for RT-qPCR with a CFX96 Touch Real-Time PCR Detection System (Bio-Rad, Hercules, CA, USA).

Expression levels of the genes were calculated by the comparative threshold method, with _EIF4A-2_ as the internal control. The Pearson correlation coefficient between RT-qPCR and RNA-Seq

was analyzed with the R statistical software35. Primer sequences are listed in Supplementary Table S1 and a related script for correlation analysis in R were presented in Supplementary Data 

S3. DATA DEPOSITION All the raw read sequences were deposited in the NCBI sequence read archive under the accession number SRP165817. RESULTS VISIBLE PWD SYMPTOMS IN NEMATODE-INOCULATED

KOREAN RED PINES To examine the relative pathogenicity of two species of the genus _Bursaphelenchus_, we monitored the visible phenotype of selected _P. densiflora_ trees injected with

water, _B. xylophilus_, or _B. thailandae_. None of these trees showed any visible PWD symptoms 4 weeks after treatment. However, 3 months after treatment, some of the green needles turned

to brown and drooped in trees inoculated with _B. xylophilus_, whereas trees treated with water or _B. thailandae_ still displayed no visible PWD symptoms (Fig. 1). This result is consistent

with the previous observation that _B. xylophilus_ is pathogenic, whereas _B. thailandae_ is non-pathogenic. _DE NOVO_ TRANSCRIPTOME ASSEMBLY OF NEMATODE-INOCULATED KOREAN RED PINES To

examine the transcriptional programs in response to inoculation with B. xylophilus or B. thailandae, we sampled cambium layers of trees 4 weeks after treatment and conducted RNA-Seq,

generating 271,466,695 raw paired reads (54,836,272,390 bp). Finally, 213,695,875 cleaned paired reads were used for _de novo_ assembly of _P. densiflora_ transcriptome (Supplementary Table 

S2). In total, 72,864 unigenes (161,111,300 nucleotides) were generated by using the Trinity assembler19 (Table 1). BUSCO analysis showed that 87.4% (1,166 single copy genes and 36

duplicated genes) complete BUSCO genes were presented in the transcriptome. In addition, 6.3% (86 genes) of all BUSCO genes were presented as fragmented form and 6.3% (87 genes) were missing

in the _de novo_ assembled transcriptome, respectively (Supplementary Fig S1). Candidate coding regions were predicted within the transcripts and the number of each type of transcripts was

listed in Supplementary Table S3 and sequencing depth for both transcripts and unigenes were also presented in Supplementary Table S4. IDENTIFICATION OF DEGS IN RESPONSE TO NEMATODE

INOCULATION IN KOREAN RED PINES At first, we examined the correlation of the biological samples to investigate relationships among them. PCA and correlation heatmap showed that water

injected tree - 3 (water - 3), _B. thailandae_ inoculated tree - 1 (_B. thailandae_ - 1) and _B. xylophilus_ inoculated tree - 1 (_B. xylophilus_ - 1) were not closely clustered together

with their biological replicates (Supplementary Fig. S2). In this experiment, we inoculated and sampled the trees in forest. Therefore, experimental condition (light, temperature, watering,

etc.) is not under control. In addition, the trees are not genotypically identical (not clones). These might cause different responsiveness to inoculation of PWN and subsequent transcriptome

differences within the biological replicates. In nature, genetic and phenotypic heterogeneity is commonly observed thus we used all the biological samples to reflect the real ecological

system for further transcriptome analysis. To identify the genes that are related to response against nematode inoculation in Korean red pines, we examined differentially expressed genes

among the comparisons using the three biological replicates each (Supplementary Fig. S3). In total, 991 genes were differentially expressed (with a cutoff of greater than 2-fold change with

a _P_-value for FDR <0.05) in all the pairwise comparisons (Fig. 2A). Among the DEGs, 102 were significantly up-regulated and 126 were down-regulated in trees inoculated with _B.

thailandae_, compared with injection of water only. In trees inoculated with _B. xylophilus_, 373 DEGs were up-regulated and 172 were down-regulated, compared with water-only controls.

Comparison of the transcriptomes of trees inoculated with _B. xylophilus_ or _B. thailandae_ identified 595 transcripts as reliable DEGs. Among them, 422 were up-regulated and 173 were

down-regulated in trees inoculated with _B. xylophilus_ compared with trees inoculated with _B. thailandae_ (Fig. 2B and Supplementary Fig. S4). Intriguingly, no unigenes were identified as

reliable DEGs that commonly observed in all pairwise comparisons. In addition, the largest number of unigenes (270) was commonly identified as DEGs in the comparison between water and _B.

xylophilus_, and _B. thaliandae_ and _B. xylophilus_ injected trees (Fig. 2C). The Venn diagrams for up- and down-regulated genes are separately presented in Supplementary Fig. S5. GO

ENRICHMENT ANALYSIS FOR DEGS To obtain comprehensive functional features associated with transcriptional programs in response to _B. xylophilus_ inoculation, we categorized the DEGs on the

basis of Gene Ontology Biological Processes (GOBPs) (Fig. 3 and Supplementary Fig. S6). GOBP analysis was conducted for the best hits from the _Arabidopsis_ genome. GOBPs related to defense

response and response to stress were enriched in the DEGs identified in the comparison between trees inoculated with _B. thailandae_ or treated with water. In the comparison between trees

inoculated with _B. xylophilus_ and those treated with water, in addition to terms shared with other comparisons (such as the defense response and response to stress, and phenylpropanoid and

flavonoid biosynthetic processes), the GOBP terms catabolic process, response to chemical, polysaccharide catabolic process, and cellular catabolic process were specifically enriched.

Notably, the largest number of GOBPs was enriched in the DEGs identified in the comparison between trees inoculated with _B. xylophilus_ and _B. thailandae_. Among these terms, response to

bacterium, cell communication, oxidation–reduction process, transmembrane receptor protein tyrosine kinase signaling pathway, lignin biosynthetic process, plant-type hypersensitive response,

defense response to nematode, and innate immune response were specific to this comparison. Heatmap of GOBPs enriched in the DEGs in each category of Venn diagram in Fig. 2C was shown in

Supplementary Fig. S7. IDENTIFICATION OF TFS INVOLVED IN THE RESPONSE TO INOCULATION WITH _B. XYLOPHILUS_ TFs govern transcriptional programs through regulation of expression of target

genes. To investigate the regulation of transcriptional programs in Korean red pines in response to _B. xylophilus_ infection, we identified differentially expressed TFs (DETFs) in the

comparison between trees inoculated with _B. xylophilus_ or with _B. thailandae_. We identified 34 TFs as DETFs and the most highly represented DETF families were the WRKY (7 members)

followed by LBD (6 members), bHLH and MYB families (5 members each) (Fig. 4A and Supplementary Table S5), and the expression patterns of these TFs in the different cambium samples are shown

in Fig. 4B. Many of the individual DETFs were up-regulated in trees inoculated with _B. xylophilus_, compared with their expression in trees inoculated with _B. thailandae_. This result

suggested that these DETFs might have important roles in controlling transcriptional programs in response to inoculation with _B. xylophilus_ by activating or repressing expression of target

genes through binding to cis-acting elements. KEY ELEMENTS OF BIOTIC-STRESS PATHWAYS ARE INVOLVED IN THE RESPONSE TO _B. XYLOPHILUS_ To investigate key signaling elements in biotic-stress

pathways, we identified DEGs in the comparison between trees inoculated with _B. xylophilus_ or with _B. thailandae_, and located each gene (and its log2 fold-change value) in biotic-stress

pathway with the MapMan visualization software (Fig. 5). Genes involved in signaling of phytohormones (such as auxin, abscisic acid, ethylene, salicylic acid, and jasmonic acid) were more

highly expressed in response to _B. xylophilus_ than to _B. thailandae_. Genes related to cell-wall modification and proteolysis were differentially expressed, and many of them showed

elevation of expression in response to _B. xylophilus_. In addition, several genes for pathogen recognition, signaling, and defense response (PR proteins) were more highly expressed in

response to _B. xylophilus_ than to _B. thailandae_. Oxidation-related processes make up one of the most important pathways in the control of defense responses in plants36, and genes related

to oxidation–reduction processes were also differentially expressed, as were a number of TFs, including several members of both the WRKY and MYB families. These types of TFs are known to

have important roles in regulation of the defense response against pathogen infection37,38,39,40. The MapMan metabolism and regulation overviews are also shown as Supplementary Fig. S8 and

S9, respectively. These results suggest that the biotic-stress pathways and related components identified here are involved in control of the defense response to _B. xylophilus_ infection.

QRT-PCR CONFORMATION OF EXPRESSION LEVELS OF DEGS To validate transcriptome results, we conducted qRT-PCR analysis with 12 DEGs in trees inoculated with _B. xylophilus_ relative to those

inoculated with _B. thailandae_. The expression changes of the DEGs from qRT-PCR and RNA-seq analyses were highly correlated (Supplementary Fig. S10), and it indicated the reliability of the

RNA-seq results. DISCUSSION In this study, we inoculated mature Korean red pines with _B. xylophilus_ or _B. thailandae_ and performed transcriptome analysis to comprehensively understand

the responses of trees against pathogenic nematodes. We identified novel resistance mechanisms and candidate genes that play important roles in resistance specifically against _B.

xylophilus_. Furthermore, we also tried to identify the pathogenesis mechanisms of _B. xylophilus_ to cause PWD. Overall, 991 DEGs were identified, and _B. xylophilus_ infection resulted in

595 DEGs compared with _B. thailandae_ infection, and 545 DEGs compared with water injection. Notably, _B. thailandae_ inoculation only resulted in 228 DEGs compared with water injection,

which suggests a limited transcriptional response of the trees to the low pathogenicity of _B. thailandae_. The physiological nature of the responses to nematode infection was examined by

GOBP analysis. Notably, the defense response GOBP was enriched in the DEGs observed in all pairwise comparisons of treatments, suggesting that both _B. xylophilus_ and _B. thailandae_

commonly affected expression of the genes involved in the defense response, regardless of their pathogenicity in the inoculated trees. To identify GOBPs related to different extents of PWD

in trees infected with different species of nematodes, we examined the genes that were differentially expressed between trees inoculated with _B. xylophilus_ and _B. thailandae_.

Phenylpropanoid and flavonoid biosynthetic processes were enriched in the DEGs in this comparison. The involvement of phenylpropanoids in the defense response is a well-established

phenomenon41. Flavonoids, isoflavonoids, hydroxycinnamic acids, monolignols, and stilbenes are types of phenylpropanoids that function as defensive molecules, acting as physical barriers and

signaling molecules to induce the defense response against pathogen invasion42. These results suggested that expression of genes related to phenylpropanoid biosynthesis is regulated for

defense against _B. xylophilus_ infection. We found that the lignin biosynthetic process GOBP was enriched in the DEGs in the comparison between trees inoculated with _B. xylophilus_ and _B.

thailandae_, and it is well known that lignin is rapidly deposited after nematode invasion, and serves as a mechanical barrier43. Correlation between the increase of lignin content and

resistance to nematodes, and an influence of lignin composition on nematode resistance, has been observed in several plant species44. Previous report showed that lignin concentration had

significantly increased after infestation of PWN at an early stage of the infestation in _P. abies_ and _C. lusitanica_45. In addition, Ishida _et al_. also showed that inoculation of _B.

xylophilus_ to Japanese black pine caused accumulation of lignin around the resin canals in the cortex46. These results could support the validity of our analysis. We also found that

expression of genes related to the hypersensitive response and oxidation–reduction process was affected by _B. xylophilus_ inoculation. The hypersensitive response is a type of cell death

that is associated with plant resistance to pathogen infection47. This localized cell death blocks the migration of the pathogen to adjacent cells, and hypersensitive response-associated

resistance is also observed in nematode–plant interactions48. Plants produce reactive oxygen species (ROS) upon nematode infection, thereby activating defense responses, and ROS are closely

related to the hypersensitive response36,49. However, ROS also function as pathogenicity factors to facilitate nematode infection in _A. thaliana_50,51. Therefore, it is also possible that

_B. xylophilus_ modulates expression of host oxidation–reduction-related genes to suppress plant defense responses. In our comparison between trees inoculated with _B. xylophilus_ and _B.

thailandae_, transmembrane receptor protein tyrosine kinase signaling pathway GOBP were enriched in the DEGs. An effective plant defense against pathogens is achieved through recognition of

pathogen-associated molecular patterns by surface-localized receptor kinases, and by the consequent downstream signaling cascade52,53. In this analysis, homologues of _Arabidopsis FLS2_ and

_NILR1_ were identified as DEGs in this biological process. Both of them were well known as important regulators in PAMP triggered immunity. In _A. thaliana_, leucine-rich repeat

receptor-like kinase NILR1 is required for innate immunity against parasitic nematodes54. Therefore, interaction between nematode-associated molecular patterns and transmembrane receptor

protein tyrosine kinases might be important for successful resistance to _B. xylophilus_ in Korean red pines. We explored the expression patterns of TFs that mainly govern transcriptional

programs responding to _B. xylophilus_ infection by comparing the expression of TFs between trees inoculated with _B. xylophilus_ and with _B. thailandae_. The family of DETFs with the most

members represented by DEGs was WRKY, followed by theLBD, bHLH, MYB, ERF, and MIKC_MADS families. WRKY TFs are involved in many developmental processes, especially in defense and

senescence37. _AtWRKY6_ homologue was one of the DETF genes in the present study, and _AtWRKY6_ is a positive regulator in the defense response against the beet-cyst nematode _Heterodera

schachtii_, and induces expression of salicylic acid-dependent defense-response genes in _A. thaliana_38,55. In addition, _AtWRKY51_ homologue was also induced after _B. xylophilus_

inoculation, and it is reported that _WRKY51_ mediates the defense response against _Pseudomonas syringae_56. Ethylene response factor (ERF) TFs are regulators of pathogenesis-related genes,

as well as ethylene-, salicylic acid-, and jasmonic acid-inducible genes57. We identified an _ERF9_ homologue as a DETF, with induction by _B. xylophilus_ inoculation. Results from a

previous study indicate that ERF9 is a negative regulator of resistance against necrotrophic fungi and acts as a molecular brake for sustained activation of a defense response58. However, it

is also plausible that _B. xylophilus_ activates expression of _ERF9_ homologue to repress the host defense system. Therefore, the DETF families identified here might play roles in

fine-tuning of the defense response through regulation of expression of downstream signaling components. In accordance with DETF and GOBP analysis, genes in the _WRKY_ and _MYB_ TF families,

cell-wall biosynthesis, proteolysis, hormone signaling, and defense showed altered expression in Korean red pines inoculated with _B. xylophilus_ in a MapMan analysis of the biotic-stress

response. In conclusion, this is the first report that describes the differences in transcriptomes between Korean red pines of felling age inoculated with _B. xylophilus_ or _B. thailandae_.

Our analysis defined biological processes that are involved in regulation of the defense response against inoculation of _P. densiflora_ with _B. xylophilus_. These results might enable the

discovery of the genes for the fast identification of _B. xylophilus_ infected trees, as well as breeding programs to produce PWD-resistant Korean red pines. DATA AVAILABILITY The data

generated and analyzed during this study are available from the corresponding author on request. REFERENCES * Hao, Z. Z., Liu, Y. Y., Nazaire, M., Wei, X. X. & Wang, X. Q. Molecular

phylogenetics and evolutionary history of sect. Quinquefoliae (Pinus): implications for Northern Hemisphere biogeography. _Mol Phylogenet Evol_ 87, 65–79,

https://doi.org/10.1016/j.ympev.2015.03.013 (2015). Article  PubMed  Google Scholar  * Canas, R. A. _et al_. Understanding developmental and adaptive cues in pine through metabolite

profiling and co-expression network analysis. _J Exp Bot_ 66, 3113–3127, https://doi.org/10.1093/jxb/erv118 (2015). Article  CAS  PubMed  PubMed Central  Google Scholar  * Kim, Y. B. _et

al_. Regulation of resin acid synthesis in Pinus densiflora by differential transcription of genes encoding multiple 1-deoxy-D-xylulose 5-phosphate synthase and

1-hydroxy-2-methyl-2-(E)-butenyl 4-diphosphate reductase genes. _Tree Physiol_ 29, 737–749, https://doi.org/10.1093/treephys/tpp002 (2009). Article  CAS  PubMed  Google Scholar  * Kwak, C.

S., Moon, S. C. & Lee, M. S. Antioxidant, antimutagenic, and antitumor effects of pine needles (Pinus densiflora). _Nutr Cancer_ 56, 162–171, https://doi.org/10.1207/s15327914nc5602_7

(2006). Article  CAS  PubMed  Google Scholar  * Futai, K. Pine wood nematode, Bursaphelenchus xylophilus. _Annu Rev Phytopathol_ 51, 61–83,

https://doi.org/10.1146/annurev-phyto-081211-172910 (2013). Article  CAS  PubMed  Google Scholar  * Mamiya, Y. History of pine wilt disease in Japan. _J Nematol_ 20, 219–226 (1988). CAS 

PubMed  PubMed Central  Google Scholar  * Braasch, H. & Braasch-Bidasak, R. _First record of the genus Bursaphelenchus Fuchs, 1937 in Thailand and description of B. thailandae sp. n.

(Nematoda: Parasitaphelenchidae)_. Vol. 4 (2002). * Han, H., Chung, Y. J. & Shin, S. C. First Report of Bursaphelenchus thailandae on Pinus densiflora in Korea. _Plant Disease_ 94,

922–922, https://doi.org/10.1094/PDIS-94-7-0922A (2010). Article  CAS  PubMed  Google Scholar  * Zhou, L. _et al_. Identifying Virulence-Associated Genes Using Transcriptomic and Proteomic

Association Analyses of the Plant Parasitic Nematode Bursaphelenchus mucronatus. _Int J Mol Sci_ 17, https://doi.org/10.3390/ijms17091492 (2016). Article  PubMed Central  Google Scholar  *

Kikuchi, T. _et al_. Genomic insights into the origin of parasitism in the emerging plant pathogen Bursaphelenchus xylophilus. _PLoS Pathog_ 7, e1002219,

https://doi.org/10.1371/journal.ppat.1002219 (2011). Article  CAS  PubMed  PubMed Central  Google Scholar  * Alves, M. _et al_. Bacterial community associated to the pine wilt disease insect

vectors Monochamus galloprovincialis and Monochamus alternatus. _Sci Rep_ 6, 23908, https://doi.org/10.1038/srep23908 (2016). Article  ADS  CAS  PubMed  PubMed Central  Google Scholar  *

Kikuchi, T., Aikawa, T., Oeda, Y., Karim, N. & Kanzaki, N. A rapid and precise diagnostic method for detecting the Pinewood nematode Bursaphelenchus xylophilus by loop-mediated

isothermal amplification. _Phytopathology_ 99, 1365–1369, https://doi.org/10.1094/PHYTO-99-12-1365 (2009). Article  CAS  PubMed  Google Scholar  * Lee, I. H. _et al_. _De novo_ assembly and

transcriptome analysis of the Pinus densiflora response to pine wilt disease in nature. _Plant Biotechnology Reports_ 12, 229–236, https://doi.org/10.1007/s11816-018-0488-5 (2018). Article 

Google Scholar  * Hewezi, T. _et al_. Cyst Nematode Parasitism Induces Dynamic Changes in the Root Epigenome. _Plant Physiol_ 174, 405–420, https://doi.org/10.1104/pp.16.01948 (2017).

Article  CAS  PubMed  PubMed Central  Google Scholar  * Jammes, F. _et al_. Genome-wide expression profiling of the host response to root-knot nematode infection in Arabidopsis. _Plant J_

44, 447–458, https://doi.org/10.1111/j.1365-313X.2005.02532.x (2005). Article  CAS  PubMed  Google Scholar  * Shin, H. _et al_. Identification of genes upregulated by pinewood nematode

inoculation in Japanese red pine. _Tree Physiol_ 29, 411–421, https://doi.org/10.1093/treephys/tpn034 (2009). Article  CAS  PubMed  Google Scholar  * Liu, Q. _et al_. Transcriptomic

Profiling Reveals Differentially Expressed Genes Associated with Pine Wood Nematode Resistance in Masson Pine (Pinus massoniana Lamb.). _Sci Rep_ 7, 4693,

https://doi.org/10.1038/s41598-017-04944-7 (2017). Article  ADS  CAS  PubMed  PubMed Central  Google Scholar  * van Dijk, E. L., Auger, H., Jaszczyszyn, Y. & Thermes, C. Ten years of

next-generation sequencing technology. _Trends Genet_ 30, 418–426, https://doi.org/10.1016/j.tig.2014.07.001 (2014). Article  CAS  PubMed  Google Scholar  * Haas, B. J. _et al_. _De novo_

transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis. _Nat Protoc_ 8, 1494–1512, https://doi.org/10.1038/nprot.2013.084 (2013).

Article  CAS  PubMed  Google Scholar  * Ayoub, S. M. _Plant nematology: an agricultural training aid_. (Nema Aid Publications, 1980). * Schmieder, R. & Edwards, R. Quality control and

preprocessing of metagenomic datasets. _Bioinformatics_ 27, 863–864, https://doi.org/10.1093/bioinformatics/btr026 (2011). Article  CAS  PubMed  PubMed Central  Google Scholar  * Tang, S.,

Lomsadze, A. & Borodovsky, M. Identification of protein coding regions in RNA transcripts. _Nucleic Acids Research_ 43, e78–e78, https://doi.org/10.1093/nar/gkv227 (2015). Article  CAS 

PubMed  PubMed Central  Google Scholar  * Fu, L., Niu, B., Zhu, Z., Wu, S. & Li, W. CD-HIT: accelerated for clustering the next-generation sequencing data. _Bioinformatics_ 28,

3150–3152, https://doi.org/10.1093/bioinformatics/bts565 (2012). Article  CAS  PubMed  PubMed Central  Google Scholar  * Li, W. & Godzik, A. Cd-hit: a fast program for clustering and

comparing large sets of protein or nucleotide sequences. _Bioinformatics_ 22, 1658–1659, https://doi.org/10.1093/bioinformatics/btl158 (2006). Article  CAS  PubMed  Google Scholar  * Simao,

F. A., Waterhouse, R. M., Ioannidis, P., Kriventseva, E. V. & Zdobnov, E. M. BUSCO: assessing genome assembly and annotation completeness with single-copy orthologs. _Bioinformatics_ 31,

3210–3212, https://doi.org/10.1093/bioinformatics/btv351 (2015). Article  CAS  PubMed  Google Scholar  * Langmead, B. Aligning short sequencing reads with Bowtie. _Curr Protoc

Bioinformatics_ CHAPTER 11, Unit 11 17, https://doi.org/10.1002/0471250953.bi1107s32 (2010). Article  Google Scholar  * Li, B. & Dewey, C. N. RSEM: accurate transcript quantification

from RNA-Seq data with or without a reference genome. _BMC Bioinformatics_ 12, 323, https://doi.org/10.1186/1471-2105-12-323 (2011). Article  CAS  PubMed  PubMed Central  Google Scholar  *

Robinson, M. D. & Smyth, G. K. Small-sample estimation of negative binomial dispersion, with applications to SAGE data. _Biostatistics_ 9, 321–332,

https://doi.org/10.1093/biostatistics/kxm030 (2008). Article  PubMed  MATH  Google Scholar  * Altschul, S. F. _et al_. Gapped BLAST and PSI-BLAST: a new generation of protein database search

programs. _Nucleic Acids Res_ 25, 3389–3402 (1997). Article  CAS  PubMed  PubMed Central  Google Scholar  * Jin, J. _et al_. PlantTFDB 4.0: toward a central hub for transcription factors

and regulatory interactions in plants. _Nucleic Acids Res_ 45, D1040–D1045, https://doi.org/10.1093/nar/gkw982 (2017). Article  CAS  PubMed  Google Scholar  * Mi, H. _et al_. PANTHER version

11: expanded annotation data from Gene Ontology and Reactome pathways, and data analysis tool enhancements. _Nucleic Acids Res_ 45, D183–D189, https://doi.org/10.1093/nar/gkw1138 (2017).

Article  CAS  PubMed  Google Scholar  * Vetushko, A. _et al_. Sustainable funding for biocuration: The Arabidopsis Information Resource (TAIR) as a case study of a subscription-based funding

model. _Database_ 2016, https://doi.org/10.1093/database/baw018 (2016). Article  PubMed  PubMed Central  Google Scholar  * Supek, F., Bosnjak, M., Skunca, N. & Smuc, T. REVIGO

summarizes and visualizes long lists of gene ontology terms. _PLoS One_ 6, e21800, https://doi.org/10.1371/journal.pone.0021800 (2011). Article  ADS  CAS  PubMed  PubMed Central  Google

Scholar  * Thimm, O. _et al_. MAPMAN: a user-driven tool to display genomics data sets onto diagrams of metabolic pathways and other biological processes. _Plant J_ 37, 914–939 (2004).

Article  CAS  PubMed  Google Scholar  * Development Core Team, R. _R: A Language and Environment for Statistical Computing_. Vol. 1 (2011). * Shah, S. J. _et al_. Damage-associated responses

of the host contribute to defence against cyst nematodes but not root-knot nematodes. _J Exp Bot_ 68, 5949–5960, https://doi.org/10.1093/jxb/erx374 (2017). Article  CAS  PubMed  PubMed

Central  Google Scholar  * Eulgem, T. & Somssich, I. E. Networks of WRKY transcription factors in defense signaling. _Curr Opin Plant Biol_ 10, 366–371,

https://doi.org/10.1016/j.pbi.2007.04.020 (2007). Article  CAS  PubMed  Google Scholar  * Pandey, S. P. & Somssich, I. E. The role of WRKY transcription factors in plant immunity. _Plant

Physiol_ 150, 1648–1655, https://doi.org/10.1104/pp.109.138990 (2009). Article  CAS  PubMed  PubMed Central  Google Scholar  * Seo, P. J. & Park, C. M. MYB96-mediated abscisic acid

signals induce pathogen resistance response by promoting salicylic acid biosynthesis in Arabidopsis. _New Phytol_ 186, 471–483, https://doi.org/10.1111/j.1469-8137.2010.03183.x (2010).

Article  CAS  PubMed  Google Scholar  * De Vos, M. _et al_. The Arabidopsis thaliana Transcription Factor AtMYB102 Functions in Defense Against the Insect Herbivore Pieris rapae. _Plant

Signal Behav_ 1, 305–311 (2006). Article  PubMed  PubMed Central  Google Scholar  * Vaganan, M. M. _et al_. Phenylpropanoid enzymes, phenolic polymers and metabolites as chemical defenses to

infection of Pratylenchus coffeae in roots of resistant and susceptible bananas (Musa spp.). _Indian J Exp Biol_ 52, 252–260 (2014). CAS  PubMed  Google Scholar  * Dixon, R. A. _et al_. The

phenylpropanoid pathway and plant defence-a genomics perspective. _Mol Plant Pathol_ 3, 371–390, https://doi.org/10.1046/j.1364-3703.2002.00131.x (2002). Article  CAS  PubMed  Google

Scholar  * Cano-Delgado, A., Penfield, S., Smith, C., Catley, M. & Bevan, M. Reduced cellulose synthesis invokes lignification and defense responses in Arabidopsis thaliana. _Plant J_

34, 351–362 (2003). Article  CAS  PubMed  Google Scholar  * Holbein, J., Grundler, F. M. & Siddique, S. Plant basal resistance to nematodes: an update. _J Exp Bot_ 67, 2049–2061,

https://doi.org/10.1093/jxb/erw005 (2016). Article  CAS  PubMed  Google Scholar  * Nunes da Silva, M., Lima, M. R. M. & Vasconcelos, M. W. Susceptibility evaluation of Picea abies and

Cupressus lusitanica to the pine wood nematode (Bursaphelenchus xylophilus). _Plant Pathology_ 62, 1398–1406, https://doi.org/10.1111/ppa.12037 (2013). Article  CAS  Google Scholar  *

Ishida, K., Hogetsu, T., Fukuda, K. & Suzuki, K. _Cortical responses in Japanese black pine attack by the pine wood nematode_. Vol. 71 (2011). * Morel, J. B. & Dangl, J. L. The

hypersensitive response and the induction of cell death in plants. _Cell Death Differ_ 4, 671–683, https://doi.org/10.1038/sj.cdd.4400309 (1997). Article  CAS  PubMed  Google Scholar  *

Branch, C., Hwang, C. F., Navarre, D. A. & Williamson, V. M. Salicylic acid is part of the Mi-1-mediated defense response to root-knot nematode in tomato. _Mol Plant Microbe Interact_

17, 351–356, https://doi.org/10.1094/MPMI.2004.17.4.351 (2004). Article  CAS  PubMed  Google Scholar  * Lozano-Torres, J. L. _et al_. Apoplastic venom allergen-like proteins of cyst

nematodes modulate the activation of basal plant innate immunity by cell surface receptors. _PLoS Pathog_ 10, e1004569, https://doi.org/10.1371/journal.ppat.1004569 (2014). Article  CAS 

PubMed  PubMed Central  Google Scholar  * Siddique, S. _et al_. Parasitic worms stimulate host NADPH oxidases to produce reactive oxygen species that limit plant cell death and promote

infection. _Sci Signal_ 7, ra33, https://doi.org/10.1126/scisignal.2004777 (2014). Article  CAS  PubMed  Google Scholar  * Feng, B. & Shan, L. ROS open roads to roundworm infection. _Sci

Signal_ 7, pe10, https://doi.org/10.1126/scisignal.2005273 (2014). Article  CAS  PubMed  PubMed Central  Google Scholar  * Nurnberger, T., Brunner, F., Kemmerling, B. & Piater, L.

Innate immunity in plants and animals: striking similarities and obvious differences. _Immunol Rev_ 198, 249–266 (2004). Article  PubMed  Google Scholar  * Thomma, B. P., Nurnberger, T.

& Joosten, M. H. Of PAMPs and effectors: the blurred PTI-ETI dichotomy. _Plant Cell_ 23, 4–15, https://doi.org/10.1105/tpc.110.082602 (2011). Article  CAS  PubMed  PubMed Central  Google

Scholar  * Mendy, B. _et al_. Arabidopsis leucine-rich repeat receptor-like kinase NILR1 is required for induction of innate immunity to parasitic nematodes. _PLoS Pathog_ 13, e1006284,

https://doi.org/10.1371/journal.ppat.1006284 (2017). Article  CAS  PubMed  PubMed Central  Google Scholar  * Ali, M. A., Wieczorek, K., Kreil, D. P. & Bohlmann, H. The beet cyst nematode

Heterodera schachtii modulates the expression of WRKY transcription factors in syncytia to favour its development in Arabidopsis roots. _PLoS One_ 9, e102360,

https://doi.org/10.1371/journal.pone.0102360 (2014). Article  ADS  PubMed  PubMed Central  Google Scholar  * Gao, Q. M., Venugopal, S., Navarre, D. & Kachroo, A. Low oleic acid-derived

repression of jasmonic acid-inducible defense responses requires the WRKY50 and WRKY51 proteins. _Plant Physiol_ 155, 464–476, https://doi.org/10.1104/pp.110.166876 (2011). Article  CAS 

PubMed  Google Scholar  * Gutterson, N. & Reuber, T. L. Regulation of disease resistance pathways by AP2/ERF transcription factors. _Curr Opin Plant Biol_ 7, 465–471,

https://doi.org/10.1016/j.pbi.2004.04.007 (2004). Article  CAS  PubMed  Google Scholar  * Maruyama, Y. _et al_. The Arabidopsis transcriptional repressor ERF9 participates in resistance

against necrotrophic fungi. _Plant Sci_ 213, 79–87, https://doi.org/10.1016/j.plantsci.2013.08.008 (2013). Article  CAS  PubMed  Google Scholar  Download references ACKNOWLEDGEMENTS This

work was supported by the National Institute of Forest Science, Republic of Korea (FE0702-2016-01 and FE0702-2017-03). AUTHOR INFORMATION AUTHORS AND AFFILIATIONS * Department of Forest

Bio-Resources, National Institute of Forest Science, Suwon, 16631, Republic of Korea Il Hwan Lee, In Sik Kim, Seok-Woo Lee & Donghwan Shim * Division of Forest Insect Pests and Diseases,

National Institute of Forest Science, Seoul, 02455, Republic of Korea Hyerim Han * Ilsong Institute of Life Science, Hallym University, Anyang, Republic of Korea Young Ho Koh Authors * Il

Hwan Lee View author publications You can also search for this author inPubMed Google Scholar * Hyerim Han View author publications You can also search for this author inPubMed Google

Scholar * Young Ho Koh View author publications You can also search for this author inPubMed Google Scholar * In Sik Kim View author publications You can also search for this author inPubMed

 Google Scholar * Seok-Woo Lee View author publications You can also search for this author inPubMed Google Scholar * Donghwan Shim View author publications You can also search for this

author inPubMed Google Scholar CONTRIBUTIONS I.H.L., D.S. and H.H. conceived and designed the experiments. I.H.L., D.S., H.H. and Y.H.K. performed the experiments. I.H.L. and D.S. analyzed

the data and wrote the manuscript. I.H.L., D.S., I.S.K. and S.W.L. carefully checked and revised the manuscript. CORRESPONDING AUTHOR Correspondence to Donghwan Shim. ETHICS DECLARATIONS

COMPETING INTERESTS The authors declare no competing interests. ADDITIONAL INFORMATION PUBLISHER’S NOTE: Springer Nature remains neutral with regard to jurisdictional claims in published

maps and institutional affiliations. SUPPLEMENTARY INFORMATION SUPPLEMENTARY MATERIALS DATASET 1 DATASET 2 DATASET 3 DATASET 4 RIGHTS AND PERMISSIONS OPEN ACCESS This article is licensed

under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate

credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article

are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and

your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this

license, visit http://creativecommons.org/licenses/by/4.0/. Reprints and permissions ABOUT THIS ARTICLE CITE THIS ARTICLE Lee, I.H., Han, H., Koh, Y.H. _et al._ Comparative Transcriptome

Analysis of _Pinus densiflora_ Following Inoculation with Pathogenic (_Bursaphelenchus xylophilus_) or Non-pathogenic Nematodes (_B. thailandae_). _Sci Rep_ 9, 12180 (2019).

https://doi.org/10.1038/s41598-019-48660-w Download citation * Received: 14 November 2018 * Accepted: 06 August 2019 * Published: 21 August 2019 * DOI:

https://doi.org/10.1038/s41598-019-48660-w SHARE THIS ARTICLE Anyone you share the following link with will be able to read this content: Get shareable link Sorry, a shareable link is not

currently available for this article. Copy to clipboard Provided by the Springer Nature SharedIt content-sharing initiative