Transcriptomic Profiling Reveals Differentially Expressed Genes Associated with Pine Wood Nematode Resistance in Masson Pine (Pinus massoniana Lamb.)

Abstract

Pine wilt disease caused by pine wood nematode (Bursaphelenchus xylophilus, PWN) is a severe forest disease of the genus Pinus. Masson pine as an important timber and oleoresin resource in South China, is the major species infected by pine wilt disease. However, the underlying mechanism of pine resistance is still unclear. Here, we performed a transcriptomics analysis to identify differentially expressed genes associated with resistance to PWN infection. By comparing the expression profiles of resistant and susceptible trees inoculated with PWN at 1, 15, or 30 days post-inoculation (dpi), 260, 371 and 152 differentially expressed genes (DEGs) in resistant trees and 756, 2179 and 398 DEGs in susceptible trees were obtained. Gene Ontology enrichment analysis of DEGs revealed that the most significant biological processes were "syncytium formation" in the resistant phenotype and "response to stress" and "terpenoid biosynthesis" in the susceptible phenotype at 1 and 15 dpi, respectively. Furthermore, some key DEGs with potential regulatory roles to PWN infection, including expansins, pinene synthases and reactive oxidation species (ROS)-related genes were evaluated in detail. Finally, we propose that the biosynthesis of oleoresin and capability of ROS scavenging are pivotal to the high resistance of PWN.

Figure 1

Resistant and susceptible phenotypes used in this study and their symptoms after PWN inoculation. (A) Vertical resin ducts in the cross section of the stem for resistant and susceptible phenotypes; the bar is 100 μm. (B) Difference in the parameters of the resin ducts between resistant and susceptible phenotypes in which resistant phenotypes had a larger AC size (resin duct size), AC area (resin duct area) and AC freq (resin duct frequency) than the susceptible phenotype (AC size comparison: n = 30, p = 0.001; AC area comparison: n = 30, p = 0.003; AC freq comparison: n = 30, p = 0.009). The mean + SE is represented by each bar. Different letters “a” and “b” on top of bars indicate significant differences (t-test, α = 0.05). (C) The materials used for transcriptome sequencing and each treatment contained three biological replicates. R: resin canal; X: xylem; P: pith; nX: new xylem.

Figure 2

Comparisons of DEGs between two phenotypes at three points. Number of DEGs obtained in resistant and susceptible phenotypes at 1, 15 and 30 dpi and Venn diagram depicting the number and overlapping relationships of DEGs between different phenotypes at the same time point or same phenotype at different time points. R: resistant phenotype; S: susceptible phenotype; I: inoculating PWN; W: inoculating water.

Figure 3

Difference in expression of six DEGs between resistant and susceptible phenotypes after PWN inoculation. These DEGs are involved in the biological process of syncytium formation in the resistant phenotype. *0.01 < P < 0.05 and **P < 0.01. Error bars represent the SE.

Figure 4

DEGs involved in terpenoid biosynthesis in stem of masson pine. I represents PWN inoculation, and W represents water control. The numbers following I and W represent the biological replicate. GT represents phenotype and time point. Enzymes involved in each step are shown in blue. DXS: 1-deoxy-D-xylulose-5-phosphate synthase; DXR: 1-deoxy-D-xylulose-5-phosphate reductoisomerase; MCT: 4-diphosphocytidyl-2C-methyl-D-erythritol synthase; CMK: 4-diphosphocytidyl-2C–methyl-D-erythritol kinase; MECPS: 2C-methyl-D-erythritol 4-phosphate cytidylyltransferase; HDS: 1-hydroxy-2-methyl-2-(E)-butenyl-4-diphosphate synthase; HDR: 1-hydroxy-2-methyl-2-(E)-butenyl-4-diphosphate reductase; AACT: acetoacetyl-CoA thiolase; HMGS: 3-hydroxy-3-methylglutaryl-CoA synthase; HMGR: 3-hydroxy-3-methylglutaryl-CoA reductase; MVK: mevalonate kinase; PMVK: phosphomevalonate kinase; MVD: mevalonate diphosphate decarboxylase; IPPI: isopentenyl-diphosphate isomerase; FPPS: farnesyl diphosphate synthase; GPPS: geranyl diphosphate synthase; SesquiTPS: sesquiterpene synthase; MonoTPS: monoterpene synthase; DiTPS: diterpene synthase; CYP720B: abietadienol/abietadienal oxidase PtAO.

Figure 5

Difference in expression of terpene synthases and GGPS involved in terpenoid biosynthesis between resistant and susceptible phenotypes after PWN inoculation. Expression of these terpene synthases was significantly changed in resistant or susceptible phenotypes after PWN inoculation. Error bars represent SE. (−)βpin syn: (−)-β-pinene synthase; (−)Lim syn: (−)-Limonene synthase; (+)αpin syn: (+)-α-pinene synthase; Lon syn: longifolene synthase; Selsyn: delta-selinene synthase; GPPS: geranyl diphosphate synthase. Error bars represent the SE.

Figure 6

The content of ROS and the expression of genes encoding ROS-scavenging enzymes in resistant and susceptible phenotypes at three points. (A) Difference in O₂ ⁻ content among resistant and susceptible phenotypes inoculating PWN and water. The results of multiple comparisons among resistant and susceptible phenotypes inoculated with PWN or water at the same time point are shown with lowercase letters. (B) Difference in H₂O₂ content among resistant and susceptible phenotypes inoculated with PWN and water. (C) Cluster analysis of DEGs encoding ROS-scavenging enzymes. The number following I and W represents the biological replicate. (D) Difference in expression of unigenes encoding catalase, superoxide dismutase, glutathione reductase and peroxidase between resistant and susceptible phenotypes after PWN inoculation. *0.01 < P < 0.05 and **P < 0.01. Error bars represent the SE.

Figure 7

qRT-PCR validation of differentially expressed genes. (A) qRT-PCR validation of unigenes associated with resistance to PWN. (B) Correlation of 43 unigene expression results obtained from qPCR and RNA-seq. Relative expression levels of qRT-PCR calculated using Elongation factor 1-alpha as the internal control. The data are expressed as the mean (±SE) of three replicates. The expression data are presented as expression values of genes in the stigma sample relative to their expression. Error bars represent the SE.