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

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.


Results
Histological characterizations of resistant and susceptible clones in masson pine and transcriptome assembly. To understand the molecular basis of nematode response in masson pine, 2 clones of masson pine that displayed distinctive phenotypes (resistant and susceptible) after PWN inoculation were identified. Before inoculating PWN, vertical resin duct characteristics of the shoot were compared between resistant and susceptible phenotypes (Fig. 1A). The result showed that resin duct size, area and number per unit area in cross section were all significantly different (P < 0.05), and the resistant phenotype had a larger resin duct size, area and number (Fig. 1B). We performed next-generation sequencing to generate the transcriptome of a stem of masson pine. After quality assessment and data filtering, we obtained 48,316,242 high quality reads (9.7 gigabase pairs), with 85.25% bases over Q30. Clean reads were assembled into 132,648 transcripts, with a mean length of 904 bp by Trinity (Supplementary Table S1). Finally, 80,340 unigenes were obtained through the clustering of overlapped transcripts with a mean length of 675 bp. Among these unigenes, 14,994 unigenes (18.66%) were larger than 1 kb. The length distributions of transcripts and unigenes (Supplementary Table S1) have a similar tendency, with N50 values of 1,634 bp and 1,245 bp, respectively. The coefficients of correlation between the 36 samples ranged from 0.74 to 1.00 (Supplementary Table S2).
Functional annotation of the unigenes was carried out by aligning the sequence similarity with several public databases, and the annotation results are listed in Supplementary Table S3. A total of 50,227 (62.5%) unigenes were annotated in the public databases (Supplementary Table S4). Based on the annotation, 32,611 unigenes were distributed into one or more Gene Ontology (GO) terms ( Supplementary Fig. S1). In addition, 16,489 sequences could be assigned to the Cluster of Orthologous Groups of proteins (COG) database for functional prediction and classification ( Supplementary Fig. S2), and 10,259 of 80,340 unigenes were annotated into 119 pathways according to the Kyoto Encyclopedia of Genes and Genomes (KEGG). These results suggested that the high-quality transcriptome was appropriate for further analyses.
Transcriptomics profiling in response to nematode inoculation of masson pine. To isolate genes responsive to PWN infection in susceptible and resistant phenotypes of masson pines, 36 sequencing libraries were constructed under two inoculated conditions (nematodes and sterile water) at three time points of 1, 15 and 30 days post-inoculation (dpi), and each sample contained three biological replicates (Fig. 1C). At 30 dpi, a characterization of the resistant trees showed almost no obvious change, but the needles began to turn yellow at 30 dpi, and the trees died at 50 dpi in the susceptible masson pine (Supplementary Fig. S3). In total, we generated 9.00-12.82 million high-quality reads per library. Subsequently, the reads were mapped to the transcriptome database described above. Approximately 80. .87% of the reads matched the sequences in the transcriptome (Supplementary Table S5).
To explore the global difference in the defense responses of susceptible and resistant phenotypes, trees inoculated with sterile water were used as controls in each phenotype, and significant DEGs were compared by applying a cutoff of a false discovery rate (FDR) less than 0.01 and a fold change no less than 2. The obtained DEGs were 260, 371 and 152 in resistant trees and 756, 2179 and 398 in susceptible trees at 1, 15 and 30 dpi, respectively ( Fig. 2). At each time point, the number of DEGs in susceptible trees was consistently higher than in resistant trees, suggesting that more DEGs were induced by PWN in susceptible trees than in resistant trees. We also found that the number of DEGs was highest at 15 dpi and least at 30 dpi in both resistant and susceptible trees, and more DEGs were down-regulated in the susceptible phenotype at the three time points.
Between different phenotypes at the same time point, or different time points for each phenotype, there were small numbers of common DEGs (Fig. 2). For example, only 21 common DEGs were found between resistant and susceptible trees at 1 dpi, and 5 common DEGs were found between 15 and 30 dpi in resistant trees, indicating that the defense responses of resistant and susceptible trees might vary as the phenotype and disease progresses.

Functional enrichment analysis of DEGs.
To gain a functional perspective of the genes related to resistance to PWN, we performed a GO enrichment analysis of these DEGs in resistant and susceptible phenotypes at three time points. Approximately 41.45% to 52.56% of DEGs were annotated by the GO database (Supplementary  Table S6). The 6 DEG sets underwent GO enrichment analysis to identify over-representative GO terms. We found that in the biological processes category, the fewer enriched GO terms were in the resistant phenotype than in the susceptible phenotype at 1 dpi and 15 dpi, indicating that more biological processes were interfered with in the susceptible phenotype (Table 1). Notably, the enriched GO term "syncytium formation" (a process that was related to the colonization of PWN in plant species) was identified in the resistant phenotype but not in the susceptible phenotype at 1 dpi and 15 dpi, although the unigenes contributing to the GO term were not the same at these two time points (Table 2). However, the most significantly enriched GO terms were "response to stress" for the susceptible phenotype at 1 dpi and 15 dpi. The results indicated that the PWN responses in the two phenotypes were fundamentally different. In biological processes, significantly enriched GO terms involved in the biosynthesis of terpenoid were also observed in the susceptible phenotype at 15 dpi, such as "geranyl diphosphate metabolic process, " "alpha-pinene biosynthetic process" and "monoterpene biosynthetic process, " suggesting that terpenoids might play an important role in the susceptible phenotype. At 30 dpi, the most significant GO terms were "amino acid catabolic process to alcohol via Ehrlich pathway," "ethanol biosynthetic process involved in (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.  glucose fermentation to ethanol" and "NADH oxidation" in the resistant phenotype. However, the only significant GO term was "cellular glucan metabolic process" in the susceptible phenotype at 30 dpi.
Characterization of gene expression related to syncytium formation. To better understand whether the GO process "syncytium formation" contributed to PWN responses, we identified related genes based on the functional annotation. A total of 7 DEGs involved in the GO term "syncytium formation" were identified in the resistant phenotype at 1 dpi and at 15 dpi. At 1 dpi, three up-regulated DEGs encoding expansin (c72778. graph_c0, c81022.graph_c0 and c72881.graph_c0) and one up-regulated DEG encoding pectate lyase precursor (c69842.graph_c0) were found in the resistant phenotype (Table 2). In the susceptible phenotype, these unigenes were not found to be significantly different between inoculation with PWN or water and were down-regulated after inoculation with PWN. Simultaneously, the DEGs encoding two expansins (c72778.graph_c0 and c81022. graph_c0) had higher expression in the resistant phenotype than in the susceptible phenotype after inoculation with PWN at 1 dpi ( Fig. 3). At 15 dpi, two up-regulated DEGs encoding expansins (c68276.graph_c0 and c72881.graph_c0) and one up-regulated DEG encoding RALF-like protein (c49518.graph_c0) were found in the resistant phenotype ( Table 2). In the susceptible phenotype, the expression level of these unigenes was not found to be significantly different between inoculation with PWN or water. Among these unigenes, the unigene encoding expansin (c72881.graph_c0) had a higher expression level in the susceptible phenotype than in the resistant phenotype after inoculation with PWN at 15 dpi (Fig. 3). The results indicated that the expression of expansin was required to reach certain thresholds to ensure syncytium formation and PWN invasion.
At 30 dpi, only delta-selinene synthase (c32505.graph_c0) was found to be down-regulated in the susceptible phenotype (Fig. 4). However, the expression level of the unigene was not significantly different between the resistant phenotype and the susceptible phenotype after PWN inoculation at 30 dpi.
In addition, three unigenes encoding geranylgeranyl pyrophosphate synthase (GGPPS) (c58549.graph_c0, c33360.graph_c0 and c33943.graph_c0) have higher expression in the resistant phenotype at 1, 15 and 30 dpi  Supplementary Fig. S3), although the expression levels of these unigenes did not change significantly in the resistant phenotype or the susceptible phenotype after inoculation with PWN at every time point.
Stress responsive pathways were involved in the resistance to PWN. We found that a total of 38 and 55 DEGs were involved in the significantly enriched GO term "response to stress" in the susceptible phenotype at 1 dpi and 15 dpi, respectively (Table 3). At 1 dpi, 37 down-regulated unigenes and one up-regulated unigene were found. Among the 37 down-regulated unigenes, there were 22 well characterized Hsps (i.e., 9 Hsp70, 4 Hsp 82, 4Hsp ST1, 2 Hsp 17.8, 1 Hsp 80, 1Hsp 83 and 1 Hsp90), 2 heat shock cognate proteins and 6 chaperone proteins. The only up-regulated unigene encoded a heat shock cognate protein (Supplementary Table S7). At 15 dpi, 52 down-regulated DEGs and 3 up-regulated DEGs were found. Among the 52 down-regulated DEGs,   most of the DEGs encoded HSP and chaperone proteins. Simultaneously, there was one down-regulated unigene of ethylene-responsive transcription factor (ERF), which was also involved in the biological process of signal transduction and response to hormone. The up-regulated unigenes included 1 Hsp70 and 2 unknown proteins (Supplementary Table S7). Among the unigenes involved in the GO term "response to stress" in the susceptible phenotype at 1 dpi and 15 dpi, only one unigene encoding Hsp70 (c70394.graph_c0) was down-regulated significantly in the resistant phenotype at 15 dpi. These data highly suggested that substantial changes in the stress defense pathways were activated by PWN in the susceptible phenotype, while the resistant phenotype was more tolerate to PWN inoculation.
The ROS scavenging pathway was required for PWN resistance. To understand the physiological changes involved in PWN infection, we investigated the accumulation of ROS in PWN-inoculated masson pines. We showed that in both phenotypes, the concentration of superoxide radical (O 2 − ) was significantly increased after inoculation of PWN at 1, 15 and 30 dpi (Fig. 6A). Furthermore, the concentration of O 2 − was higher in the susceptible phenotype than in the resistant phenotype when inoculating PWN at 30 dpi. Although the concentration of hydrogen peroxide (H 2 O 2 ) was significantly increased in the resistant phenotype at 1 dpi, it decreased gradually, and no significant difference was found between trees inoculated with PWN at 15 and 30 dpi. However, in the susceptible phenotype, the concentration of H 2 O 2 increased quickly after inoculating with PWN, and it was much higher than in the control inoculated with water and the resistant phenotype inoculated with PWN (Fig. 6B). The results of the physiological experiment showed that more O 2 − and H 2 O 2 accumulated in the susceptible phenotype after inoculation with PWN.
From the transcriptome sequencing, a total of 33 DEGs related to ROS scavenging were identified in two phenotypes at three time points (Fig. 6C). Nine down-regulated DEGs encoding superoxide dismutase, glutathione reductase and L-ascorbate peroxidase were only found in the susceptible phenotype. Fifteen down-regulated DEGs encoding catalase or catalase isozyme were expressed in both phenotypes. The DEGs encoding peroxidase were down-regulated at 1 dpi and up-regulated at 15 and 30 dpi in the resistant phenotype, but no obvious trend was found in the susceptible phenotype. Among these DEGs, the expression levels of catalase (c77348.graph_c0), superoxide dismutase (c68079.graph_c0), glutathione reductase (c60353.graph_c0) and peroxidase (c65380. graph_c0) were involved in the biological process of defense response to nematode according to the annotation by the GO database (Fig. 6D), almost remaining higher in the resistant phenotype than in the susceptible phenotype at the three time points. These results suggested that the resistant phenotype was more competent to scavenge ROS, which was in good agreement with the physiological result.

Validation of RNA-seq expression data by qRT-PCR.
To validate the reliability of the RNA-seq results, 43 unigenes from the DEGs described above were selected for qRT-PCR analysis (Supplementary Table S8). The expression of these unigenes was significantly different between the inoculated PWN and water trees for each phenotype or between the resistant and susceptible phenotype after PWN infection at 1, 15 or 30 dpi (Fig. 7A). The unigenes selected for qRT-PCR analysis were those involved in syncytium formation, terpenoid biosynthesis, pathogenesis-related genes, cell wall-related genes, heat shock protein and ROS responsive genes. The expression pattern of almost all unigenes indicated by qRT-PCR agreed well with RNA-seq except for one unigene (c68240. graph_c0). Simultaneously, a solid correlation was observed between the data for these unigenes from qRT-PCR and RNA-seq, with a correlation coefficient of 0.88 (Fig. 7B).

Discussion
More DEGs were obtained between the resistant and susceptible masson pines using next-generation sequencing in our study than in Japanese black pines using suppression subtractive hybridization 22 . To identify high confidence DEGs in response to PWN, it is important to set appropriate controls for comparison in large-scale gene expression analysis. In this study, the control (inoculation with water) was set at each sampled time point in both resistant and susceptible masson pines, which enhanced the reliability of the obtained DEGs. In the resistant phenotype, the number of DEGs was always less than for the susceptible phenotype at each time point, which implied that the resistant masson pine is not easily interfered with by PWN compared with the susceptible masson pine. It has been shown that the resistance and susceptibility of plant species depend upon qualitative and quantitative differences in the activated genes in response to insects and pathogens 27 . According to the GO classification of DEGs, the difference in the significant GO terms in the resistant phenotype and susceptible phenotype indicated a qualitative difference in activated defense genes between the two phenotypes to PWN infection.
We observed that the most significant GO terms were "syncytium formation" for the DEG between inoculation with PWN or water in the resistant phenotype at 1 dpi and 15 dpi (Table 1). Syncytium formation is considered to be evoked by the nematode releasing secretions into initial syncytial cells; then, the partial cell wall between the initial syncytial cells and neighboring cells is dissolved, and the protoplasts are fused 28 . In this study, the DEGs involved in the biological process of syncytium formation included putative expansin, pectate lyase precursor and RALF-like protein. Expansins, which are cell wall-loosening proteins, were involved in the growth regulation of various cell types in response to different stimuli at different plant developmental stages 29 . In tomato and soybean, expansin was identified after inoculation with cyst nematodes 30,31 , indicating the importance of cell wall modification during the plant defense response. Hirao et al. 22 also reported that cell wall-related genes played a role in the effective defense response against PWN infection in P. thunbergii. RALF-like proteins, forming a large family, are small polypeptides 32 . Pearce et al. 33 found that RALF was involved in the cell-cell signaling biological process and can activate intracellular MAPKs. Although previous studies have thus far revealed that expansins participated in cell wall disassembly and cell growth during syncytium formation when plants were infected with root-knot nematodes in roots 34 , syncytium was not found in the stem of masson pine after inoculation with PWN. However, it is well known that new traumatic resin ducts are induced in conifers under biotic or abiotic stress 35 , which are differentiated by clusters of swollen and irregularly shaped cells in the cambial zone. Nagy et al. 36 also indicated that resistant trees induced traumatic resin ducts earlier than susceptible trees in response to fungal infection. During the development of traumatic resin ducts, the cell wall was changed greatly. Therefore, the function of expansin might be associated with the formation of traumatic resin ducts. In our study, we found that the expression of expansin was higher in resistant trees than susceptible trees at 1 dpi, but expansin was induced much more in susceptible trees at 15 dpi.
Oleoresin, a mixture of monoterpenes, sesquiterpenes and diterpene resin acids, plays an important role in defense against pathogens and insects in conifers 37,38 . The biosynthesis of the terpenoid backbone is usually completed through two separate pathways of methyl-erythritol 4-phosphate (MEP) and mevalonate (MVA) in conifers. Monoterpenes and diterpenes are biosynthesized by the MEP pathway 39  To date, only two P450 genes, CYP720B1 from loblolly pine (P. taeda) and CYP720B4 from white spruce (Picea sitchensis), which take part in the formation of the diterpene resin acids of oleoresin, have been functionally characterized 40,41 .
After the trunk of a conifer suffers an insect attack, pathogen invasion or mechanical wounding, inducible oleoresin can be synthesized in addition to constitutive oleoresin, and the volume and composition of oleoresin varies depending by herbivores and pathogens 36,42,43 . The turpentine (monoterpenes and sesquiterpenes) can directly affect herbivores through emitting toxic volatiles, and diterpene resin acids present physical barriers at the wound 38,44,45 . Among the components of oleoresin, limonene is considered the most abundant toxic substance to beetles 46 . Zhao et al. 47 reported that the concentration of longifolene was changed after inoculation with PWN in masson pine. As the products of (+)-and (−)-alpha-pinene synthase, (+)-alpha-pinene and (−)-alpha-pinene present simultaneously in masson pine. Alpha-pinene is a mixture of the isomers (+)-alpha-pinene and (−)-alpha-pinene. The content of (−)-alpha-pinene is approximately 10-to 20-fold higher than (+)-alpha-pinene in masson pines 41 . The mono-TPS and sesqui-TPS might play important roles in controlling PWN infestation. The candidate genes of (−)-beta-pinene synthase, (+)-alpha-pinene synthase, (−)-alpha/beta-pinene synthase, and (−)-limonene synthase might play positive regulatory roles, and the candidate genes of longifolene synthase and delta-selinene synthase might play negative regulatory roles. In this study, we found that the expression levels of many genes involved in terpenoid biosynthesis were changed after PWN inoculation. However, among these unigenes, the expressed difference of DXS, HDS, MECPS, HMGR, and (−)-alpha-pinene synthase was not significant after PWN inoculation between the resistant phenotype and the susceptible phenotype. The arrangement of the resin ducts might affect PWN activities. Kuroda 48 reported that intricate and winding structure of resin ducts at the joints between the branch and the main stem might be effective in retarding PWN. In P. massoniana, we found that the resistant phenotype had a larger resin duct size, area and number than the susceptible phenotype when the trees were in the original state (Fig. 1A). The area of the axial resin canals was positively related to the oleoresin yield in P. pinaster 49 . The trees having a large amount of oleoresin had stronger resistance in loblolly pine 50 . In this study, although the relationship between GGPPS and the defense response to PWN infestation is not clear, it is interesting that the expression of GGPPS has a higher level in resistant than in susceptible trees. Our previous study demonstrated that GGPPS was expressed at higher levels in the high oleoresin-yielding phenotype than in the low oleoresin-yielding phenotype of masson pine 51 . Therefore, resistant masson pines might have a higher oleoresin yield than susceptible ones. However, except for the genes encoding CYP450, no genes involved in terpenoid biosynthesis had differential expression between resistant and susceptible trees in P. thunbergii 22 , which suggested that the defense mechanism against PWN may vary by tree species.
In addition to oleoresin, releasing ROS as signal molecules or toxic molecules plays an important role in pathogen defense 52 . Hirao et al. 22 found that the genes involved in the ROS scavenging were rapidly induced in P. thunbergii infected with PWN. The major ROS-scavenging enzymes of plants include superoxide dismutase, ascorbate peroxidase, catalase, monodehydroascorbate, dehydroascorbate reductase, glutathione reductase, and glutathione peroxidase 53 . As the H 2 O 2 scavenging enzyme, ascorbate peroxidase might be responsible for the fine modulation of ROS for signaling purposes, but catalase might be responsible for removing excess ROS during stress. In this study, we found that the content of ROS (O 2 − and H 2 O 2 ) was increased after inoculating PWN in both the resistant and susceptible phenotypes at 1 dpi, and the content of the ROS in the resistant phenotype then decreased gradually at 15 and 30 dpi (Fig. 6A,B). Especially for H 2 O 2 , the content was almost maintained at the normal level in the resistant phenotype at 30 dpi. From the results of the transcriptome, the expression levels of ROS scavengers, such as catalase, superoxide dismutase, peroxidase and glutathione reductase, were higher in the resistant phenotype than in the susceptible phenotype at almost every time point, which indicated that the resistant phenotype had a stronger ROS-scavenging capability. The enhanced H 2 O 2 can induce the expression of PRs and phytoalexin in plants 54 . PR-1, PR-2 and PR-5 are salicylic acid-responsive genes, and PR-6 is a jasmonic acid and ethylene responsive gene 55 . In our study, the expression of PR-1, PR-3, PR-5, PR-9 and PR-10 was changed significantly in the resistant or susceptible phenotype after PWN inoculation. However, no obvious trend was observed for PR-1, PR-3, PR-5 and PR-10. Therefore, PRs (PR-1, PR-3, PR-5 and PR-10) might not be effective in controlling PWN infection, which is consistent with the result for P. thunbergii 22 .

Conclusion
The gene expression profiling between resistant and susceptible masson pines after inoculation with PWN is studied for the first time in this work. We systematically identified genes that were potentially related to PWN resistance using a combination of DEG analyses. We showed that most of the DEGs were down-regulated, possibly leading to the physiologic disorders in the susceptible phenotype. We revealed that the resistant and susceptible phenotypes had a different defense mechanism in response to PWN. The GO processes of "syncytium formation, " enriched in the resistant phenotype, and "response to stress, " enriched in the susceptible phenotype, highlighted the different roles of the underlying pathways that were related to PWN tolerance. Moreover, detailed gene expression analysis suggested that terpenoids were prominent defense compounds against PWN. The mono-TPS of (−)-beta-pinene synthase, (+)-alpha-pinene synthase, (−)-alpha/beta-pinene synthase, and (−)-limonene synthase might play positive regulatory roles, and the sesqui-TPS of longifolene synthase and delta-selinene synthase might play negative regulatory roles in resistance against PWN infection. In addition, the higher activity of ROS-scavenging enzymes was effective to inhibit the death of resistant masson pines after PWN inoculation. This study provides a starting point for understanding a successful defense mechanism against PWN in masson pine. In future studies, more detailed functional analysis of thekey genes obtained in this study should be carried out to unravel the complex defense mechanisms.

Materials and Methods
Plant materials. Plant materials and treatments. A resistant variety of masson pine, 'Xiuning 5,' which has the higher resistance to PWN (mortality rates of its clones were all 0% after inoculating with PWN at three sites in the previous test) and is selected from 324 resistant families, was planted in the germplasm nursery of Anhui Academy of Forestry, Anhui Province, China. A susceptible variety 'Huangshan 1' (mortality rates of its clones were all 100% after inoculation with PWN at three sites in the previous test), selected as the control, was also planted in the germplasm nursery of the Anhui Academy of Forestry. Both clones were obtained from the original trees using the pith-cambium pairing grafting method onto 2-year-old seedlings at the Anhui Academy of Forestry in 2010. The PWN used in our study was the highly virulent Guangzhu-3B isolate.
The experiment was conducted using eighteen four-year-old ramets for each clone on July 20, 2014. The average heights of the resistant and susceptible clones were 2.43 m and 2.51 m, respectively. The tree branch approximately 2 cm below the node of the annual shoot was cut at a slant with a knife so that it was approximately 3 cm long and 1 mm deep into the xylem. Then, the section was scraped using a saw to imitate bites from cerambycid beetles. Subsequently, 10,000 nematodes suspended in 20 µL sterile water were injected into the longitudinal wound of nine ramets of each clone. For the other nine ramets of each clone, only sterile water was injected into the wound as a control. At 1, 15 and 30 dpi, the stem tissue of inoculated PWN trees and control trees was collected 5 cm above the inoculation point. Three ramets for each clone were selected as three biological replicates. Then, a 2 cm-long segment of stem was cut off and put into liquid nitrogen immediately in the field. Finally, these samples were stored at −80 °C for RNA extraction. Simultaneously, the needles per sample above the 1 cm inoculation point were also taken and stored at −20 °C for ROS measurement. The last sampled time was according to previous results that the needles of susceptible trees turned yellow after 30 dpi, which was observed in this experiment.
Total RNA isolation and cDNA synthesis. Total RNA from each sample of inoculated PWN trees and control trees at 1, 15 and 30 dpi was separately extracted using the Plant RNA kit RN38 EASYSpin plus (Aidlab Biotech, Beijing, China), containing DNA digestive enzyme. The concentration of the total RNA was detected using an Ultraspec TM 2100 Pro UV/visible spectrophotometer, and the integrity of the total RNA was detected using 1% agarose gel electrophoresis. Then, the integrity of the RNA was also tested using an Agilent 2100 Bioanalyzer before proceeding. High-quality RNA was defined as having a concentration of more than 400 ng/μL, OD260/280 ranging from 1.8 to 2.2, OD260/230 over 2.0, 28S/18S no less than 1.5 and RIN over 8.0. A total of thirty-six RNA samples were obtained, including three biological replicates for inoculated PWN trees and control trees at three time points.
Transcriptome sequencing and assembly of masson pine. To obtain a comprehensive list of transcripts, equal amounts of high-quality RNA from thirty-six RNA samples were pooled together. Then, the total RNA was delivered to the Biomarker Technology Company (Beijing, China) to construct a cDNA library and perform sequencing. The cDNA library was paired-end sequenced on the Illumina HiSeqTM 2000 sequencing platform, and 2× 101 bp reads were produced from either end of a cDNA fragment.
After filtering the raw reads, we assembled the high-quality reads into contigs using the Trinity method 56 . Then, the transcripts were assembled with the paired-end information of the contigs. Finally, the longest transcripts from the potential alternative splicing transcripts were selected as the unigenes.
Differentially expressed genes and GO enrichment. All high-quality reads of 36 samples were aligned to the assembled transcripts above using Bowtie, allowing no more than one nucleotide mismatch. The data has been submitted to NCBI, and the accession number is SRP103562. To compare the expression abundance among 36 samples, RPKM was used to represent the expression abundance of the unigenes. The differential expression of genes was analyzed by the edgeR package (Version 3.2.4) in BioConductor (release 2.12, R version 2.15.0) 57 . The right-sided hypergeometric enrichment test was performed at a medium network specificity selection, and the p-value correction was performed using the Benjamini-Hochberg method. A GO enrichment analysis of the DEGs was performed using the GOseq R package, and GO terms with a corrected p-value below 0.05 were considered to have significant enrichment. Heatmaps and hierarchical clustering analysis were performed using the HemI software (Heatmap Illustrator, version 1.0.3.3) 58 , with the data normalized as in Equation (1): where y ij represents the value of the jth replicate of the ith phenotype used in the heatmap and hierarchical clustering (i = 1 to 2 and j = 1 to 3), x ij represents the value of the jth replicate of the ith phenotype obtained by RNAseq, and u is the overall mean.
Quantitative RT-PCR analysis. The RNA samples used in the qRT-PCR and transcriptome sequencing were identical. The primer pairs (Supplementary Table S9) for the representative candidate genes were designed using the Primer 3.0 software with the optimal Tm at 58-62 °C, primer length of 19-21 nucleotides, PCR product size of 120-200 bp and GC content of 45-55%. Quantitative RT-PCR was run in a 7300 Real Time PCR System (Applied Biosystems, CA, USA) with the SYBR Green detection method to verify the transcriptome sequencing results. The cDNA was amplified in a total volume of 20 μL, including 0.4 μL of ROX reference dye, 2 μL of cDNA, 0.4 μL of primers and 10 μL of SYBR Premix ExTaqTM. The PCR program was 95 °C for 10 s and 40 cycles of 95 °C for 5 s, 58 °C for 30 s and 72 °C for 30 s. Three technical replicates for each of the three biological replicates were performed. The transcript profiles were normalized using the reference gene of elongation factor 1-alpha, and the relative expression levels of candidate genes were calculated with the 2 −ΔΔCt method 59 .