Comprehensive Analysis and Functional Verification of the Pinus massoniana NBS-LRR Gene Family Involved in the Resistance to Bursaphelenchus xylophilus

Pinus massoniana Lamb. is a crucial timber and resin conifer in China, but its plantation industry is threatened by outbreaks of pine wilt disease (PWD) caused by Bursaphelenchus xylophilus (pinewood nematode; PWN). However, as of yet, there is no comprehensive analysis of NBS-LRR genes in P. massoniana involved in its defense against PWN. In this study, 507 NBS genes were identified in the transcriptome of resistant and susceptible P. masoniana inoculated with the PWN. The phylogenetic analysis and expression profiles of resistant and susceptible P. massoniana revealed that the up-regulated PmNBS-LRR97 gene was involved in conferring resistance to PWN. The results of real-time quantitative PCR (qRT-PCR) showed that PmNBS-LRR97 was significantly up-regulated after PWN infection, especially in the stems. Subcellular localization indicated that PmNBS-LRR97 located to the cell membrane. PmNBS-LRR97 significantly activated the expression of reactive oxygen species (ROS)-related genes in P. massoniana. In addition, the overexpression of PmNBS-LRR97 was capable of promoting the production of ROS, aiding in plant growth and development. In summary, PmNBS-LRR97 participates in the defense response to PWN and plays an active role in conferring resistance in P. massoniana. This finding provides new insight into the regulatory mechanism of the R gene in P. massoniana.


Introduction
Plants are subjected to external biotic and abiotic stresses throughout their life cycle. Over the course of long-term evolution, a complete immune system has evolved to resist such external threats to plant fitness. Two defense mechanisms, pathogen-associated molecular pattern (PAMP)-triggered Immunity (PTI) and effector-triggered immunity (ETI), constitute the immune system of plants. When plants are attacked by pathogens, pattern recognition receptors (PRRs) on the cell surface identify PAMPs from attacking species and trigger successive defense responses. However, in the process of host-pathogen interaction, the effector proteins secreted by the invading pathogen can inhibit PTI and destroy the first line of defense against entry into plant cells. At this time, the intracellular immune receptor recognizes the effector proteins that enter the cells and triggers a defense response [1][2][3][4]. Pattern recognition receptors in PTI and intracellular immune receptors in ETI are collectively referred to as R genes, which can recognize pathogens, carry out signal transduction, and cause downstream defense responses, collectively rendering plants resistant to their enemies [5]. Up to now, more than 400 disease R genes have been cloned from different plants, and nine molecular mechanisms of R gene-mediated resistance have been identified [6,7].

Identification of NBS Genes in P. massoniana
A total of 507 NBS proteins were obtained from the P. massoniana transcriptome (Tables 1 and S1). Based on the domain predictions, the sequences were divided into eight subclasses: N (NBS), NL (NBS-LRR), TNL (TIR-NBS-LRR), TN (TIR-NBS), CNL (CC-NBS-LRR), CN (CC-NBS), RNL (RPW8-NBS-LRR), and RN (RPW8-NBS). The TNL subclass had the largest number of genes, 120, accounting for 23.67% of the total, whereas the CN subclass had the fewest, only 6, or 1.18% of the total. Four types of genes contained an LRR domain, totaling 292 genes (57.59% of the total), while four types of genes lacked an LRR domain, for a total of 215. The total number and subclass number of NBS-LRR genes varied greatly among different plants. There were 402 NBS genes in P. trichocarpa, 105 fewer than those in P. massoniana, and 209 in A. thaliana, the latter accounting for just 41.22% of the total number present in P. massoniana. Like P. massoniana, TNL also had the most genes in A. thaliana, while in P. trichocarpa, the most abundant subclasses were NL and CNL; likewise, the least abundant subclasses among the three species also differed. They

Sequence, Phylogenetic, and Conserved Motif Analyses of NBS-LRR Resistance Genes in P. massoniana
The basic physicochemical properties, subcellular localization, presence or absence of signal peptides, and transmembrane domains of the protein sequences encoded by 292 NBS-LRR genes of NL, TNL, CNL, and RNL were analyzed (Table S2).
The NBS domain region sequences were extracted from 292 NBS-LRR genes, and MEGA7.0 was used for their phylogenetic analysis. The resulting phylogenetic tree shows that the 292 NBS-LRR disease resistance proteins could be clearly divided into two large branches ( Figure 1): TNL and non-TNL, with the non-TNL branch further divided into five small branches (non-TNL-A, non-TNL-B, non-TNL-C, non-TNL-D, and non-TNL-E). TNL subclass genes are basically homogenous in the large TNL branch. The NL subtype genes were distributed in the two large branches, with more genes in the non-TNL large branch than in the TNL large branch, indicating that the origin of these subtypes' genes was more diverse. The genes of each subclass in the non-TNL large branch were not clustered on one branch as clearly as the TNL subclass genes were. The NL subclass genes were mainly clustered in the non-TNL-C and non-TNL-E small branches. The genes of RNL subtypes were mainly clustered in the non-TNL-B and non-TNL-D branches, and the CNL subtype genes were mainly clustered in the non-TNL-A branch. Apart from most genes of each subtype being clustered in one branch, there were a few genes scattered in other small branches.
The conserved motifs of 292 NBS-LRR genes were analyzed using the MEME online software. This revealed that the composition and arrangement of the conserved domains and motifs of genes in the same clade were basically the same, whereas the differences between different clades were pronounced (Table S3 and Figure S1). There was little difference between TNL subclass genes, with 4, 10, and 7 motifs found in the TIR, NBS, and LRR domains, respectively. The seven motifs of the LRR domain were named L1-7, respectively. Compared with the TIR and NBS domains, the structure of the LRR domain was relatively less conserved. Due to the large sequence differences among NL subtype genes, they had multiple origins; hence, the composition and arrangement of motifs differed markedly between different small branches. A special motif (motif 24), named RPW8-1, was found in the RPW8 domain of RNL subclass genes. No motifs were found in the CC domain of CNL subclass genes, probably because the N-terminal sequence was not highly conserved. There was a motif (motif 19), named NBS-4, in the NBS domain of CNL subclass genes, which was not present in TNL subclass genes. Compared with the CC domain, the NBS domain was evidently highly conserved. The conserved motifs of 292 NBS-LRR genes were analyzed using the MEME online software. This revealed that the composition and arrangement of the conserved domains and motifs of genes in the same clade were basically the same, whereas the differences between different clades were pronounced (Table S3 and Figure S1). There was little difference between TNL subclass genes, with 4, 10, and 7 motifs found in the TIR, NBS, and LRR domains, respectively. The seven motifs of the LRR domain were named L1-7, respectively. Compared with the TIR and NBS domains, the structure of the LRR domain was relatively less conserved. Due to the large sequence differences among NL subtype genes, they had multiple origins; hence, the composition and arrangement of motifs differed markedly between different small branches. A special motif (motif 24), named RPW8-1, was found in the RPW8 domain of RNL subclass genes. No motifs were found in the CC domain of CNL subclass genes, probably because the N-terminal sequence was not highly conserved. There was a motif (motif 19), named NBS-4, in the NBS domain of CNL subclass genes, which was not present in TNL subclass genes. Compared with the CC domain, the NBS domain was evidently highly conserved.

Differential Gene Expression Pattern of NBS-LRR and Cloning of PmNBS-LRR97 in P. massoniana
A heat map was drawn using the expression data of differentially expressed NBS-LRR genes in the transcriptome of resistant versus susceptible P. massoniana at different days since the inoculation with the PWN (Figure 2). These results showed that the  PmNBS-LRR184, and PmNBS-LRR203 were down-regulated genes whose expression levels in susceptible P. massoniana surpassed those in resistant P. massoniana. Further, the expression levels of PmNBS-LRR25 and PmNBS-LRR62 were decreased in both resistant and susceptible P. massoniana after the inoculation with the PWN, while those of PmNBS-LRR97, PmNBS-LRR16, and PmNBS-LRR64 increased by 2.17-fold, 1.66-fold, and 1.21-fold in resistant P. massoniana. Accordingly, PmNBS-LRR97 was regarded as the dominant gene involved in PWN infection. Sequence analysis showed that the PmNBS-LRR97 gene sequences in resistant and susceptible P. massoniana were entirely identical. The full-length cDNA of the sequence was 2373 bp, encoding a total of 790 aa ( Figure 3A), with a molecular weight of 90.213 kDa and an isoelectric point of 6.18. It is an unstable protein given its instability index of 98.04; it is a thermophilic protein whose grand average for hydrophilicity is less than 0, meaning that it is also a hydrophilic protein; it lacks both a signal peptide and a transmembrane domain; and it contains 61 phosphorylation sites.
The PmNBS-LRR97 sequence contained TIR, NBS, and LRR domains and was evidently a TNL subclass gene. An analysis of its secondary structure showed that the PmNBS-LRR97-encoded protein had four structural elements: an alpha helix, an extended strand, a beta turn, and a random coil ( Figure 3B). Among them, the alpha helix is composed of 448 aa, accounting for the most amino acids (56.71%), while the beta turn consisted of 28 aa, the fewest (3.54%). Sequence analysis showed that the PmNBS-LRR97 gene sequences in resistant and susceptible P. massoniana were entirely identical. The full-length cDNA of the sequence was 2373 bp, encoding a total of 790 aa ( Figure 3A), with a molecular weight of 90.213 kDa and an isoelectric point of 6.18. It is an unstable protein given its instability index of 98.04; it is a thermophilic protein whose grand average for hydrophilicity is less than 0, meaning that it is also a hydrophilic protein; it lacks both a signal peptide and a transmembrane domain; and it contains 61 phosphorylation sites.

Expression Pattern of PmNBS-LRR97 in Resistant and Susceptible P. massoniana
Taking the susceptible P. massoniana inoculated with water at 0 days post-inoculation (dpi) as the control, the expression pattern of PmNBS-LRR97 in the stems of resistant and The PmNBS-LRR97 sequence contained TIR, NBS, and LRR domains and was evidently a TNL subclass gene. An analysis of its secondary structure showed that the PmNBS-LRR97encoded protein had four structural elements: an alpha helix, an extended strand, a beta turn, and a random coil ( Figure 3B). Among them, the alpha helix is composed of 448 aa, accounting for the most amino acids (56.71%), while the beta turn consisted of 28 aa, the fewest (3.54%).

Expression Pattern of PmNBS-LRR97 in Resistant and Susceptible P. massoniana
Taking the susceptible P. massoniana inoculated with water at 0 days post-inoculation (dpi) as the control, the expression pattern of PmNBS-LRR97 in the stems of resistant and susceptible P. massoniana inoculated with the PWN was analyzed by qRT-PCR at 0, 1, 15, and 30 dpi ( Figure 4A). At each time point, after inoculating the resistant and susceptible P. massoniana with water, the gene's expression level basically did not exceed that at 0 dpi; in stark contrast, after resistant and susceptible P. massoniana were inoculated with PWN, the gene featured an 'up-down-down' expression trend. Compared with the respective expression levels at 0 dpi, the expression level of resistant P. massoniana after inoculation was six times higher than that of susceptible P. massoniana. Next, the expression pattern of PmNBS-LRR97 in different tissues was explored by comparing the roots of resistant and susceptible P. massoniana exposed to PWN at different times post-inoculation ( Figure 4B,C). The gene's expression level in stems was higher than that in roots or needles across all time points, regardless of whether P. massoniana was resistant or susceptible.   The subcellular localization results showed that PmNBS-LRR97 located to the cell membrane ( Figure 4D), which is consistent with the function of the R gene of recognizing effectors.

Overexpression of PmNBS-LRR97 Enhances Plant Growth and Increases the Content of Defense Signaling Molecules
PmNBS-LRR97 was transformed in N. benthamiana to verify its function, because the genetic transformation system of P. massoniana has not been successfully established. The qRT-PCR results demonstrated that, compared with WT, the PmNBS-LRR97 expression levels of all transgenic tobacco lines were higher, with those of PmNBS-LRR97-OE5, PmNBS-LRR97-OE4, and PmNBS-LRR97-OE3 being the greatest; hence, these three lines were used for further verification ( Figure 5A).
To explore whether PmNBS-LRR97 can affect plant growth, the phenotypes of WT and transgenic tobacco were observed ( Figure 5B-D). The transgenic tobacco plants were taller than WT counterparts, irrespective of their age (i.e., 30, 60, and 90 d) (p < 0.01). It was also found that the transgenic tobacco bloomed earlier than WT and had more flowers and greater branching at 90 days, indicating that the overexpression of PmNBS-LRR97 can effectively promote the growth and development of plants.
Plant disease resistance genes can recognize effectors secreted by nematodes, activate signaling pathways, and cause a series of defense responses. The burst of reactive oxygen species (ROS) is the initial defense response. The production and scavenging of ROS in plants are normally maintained in a dynamic equilibrium [34], but once stimulated by foreign organisms, that dynamic balance will falter and cause damage to plants. The expression patterns of PmSOD, PmPOD, and PmCAT were analyzed in the stems of resistant and susceptible P. massoniana ( Figure 6A-C) to explore the relationship between R genes and ROS-related genes. These results uncovered an 'up-down-down' expression trend for PmSOD, PmPOD, and PmCAT after inoculating the resistant and susceptible P. massoniana with the PWN. Moreover, the up-regulated fold-changes of PmSOD, PmPOD, and PmCAT were far greater in resistant than susceptible P. massoniana. This result was similar to the expression pattern of PmNBS-LRR97. Next, the expression of NbSOD, NbPOD, and NbCAT was analyzed in WT and transgenic tobacco ( Figure 6D-F), and the contents of known plant defense signaling molecules (H 2 O 2 , NO) and enzymes related to the ROS scavenging system (SOD, POD, CAT) were detected in tobacco ( Figure 6G-K). The expression levels of NbSOD, NbPOD, and NbCAT were all markedly higher in transgenic tobacco than WT, and the contents of H 2 O 2 , NO, SOD, POD, and CAT were also higher. This indicated that the overexpression of PmNBS-LRR97 can increase the content of defense signal molecules and related enzymes, which is vital for the resistance of P. massoniana to the PWN. Plant disease resistance genes can recognize effectors secreted by nematodes, activate signaling pathways, and cause a series of defense responses. The burst of reactive oxygen species (ROS) is the initial defense response. The production and scavenging of ROS in plants are normally maintained in a dynamic equilibrium [34], but once stimulated by foreign organisms, that dynamic balance will falter and cause damage to plants. The expression patterns of PmSOD, PmPOD, and PmCAT were analyzed in the stems of resistant and susceptible P. massoniana ( Figure 6A-C) to explore the relationship between R genes and ROS-related genes. These results uncovered an 'up-down-down' expression trend for PmSOD, PmPOD, and PmCAT after inoculating the resistant and susceptible P. massoniana with the PWN. Moreover, the up-regulated fold-changes of PmSOD, PmPOD, and PmCAT were far greater in resistant than susceptible P. massoniana. This result was

Discussion
Plant R genes can specifically recognize effector proteins secreted by pathogens and initiate a cascade of downstream defense mechanisms, which is indispensable for plant resistance to pathogens [35][36][37]. The NBS-LRR gene family, the largest category of R genes, has been identified and systematically analyzed for its involvement in the defense processes of many plants. P. massoniana, as a pioneer conifer species, is vulnerable to the PWN. However, the defense strategy of R genes in P. massoniana facing an attack from the PWN remains unknown. Therefore, the identification, characterization, and functional analysis of R genes are crucial steps for understanding how P. massoniana is able to achieve resistance against PWN.
The members of the NBS-LRR gene family were identified in P. massoniana to gain insight into how they conferred resistance to the PWN. Given the differing genome sizes of plant species, their numbers of NBS-LRR gene families and subfamilies are also different. The results showed that 507 NBS genes were identified in the P. massoniana transcriptome by bioinformatic methods, more than those in A. thaliana [11] and P. trichocarpa [12] but fewer than those in E. grandis [13] or O. sativa [14]. The NBS-LRR gene family of P. massoniana can be separated into eight subclasses, of which the number of TNL subclass genes was 3.87 times that of CNL subclass genes; the number of TNL subclass genes in A. thaliana was 1.8 times that of CNL; and the number of CNL subclass genes in P. trichocarpa was 1.5 times that of TNL (Table 1). Different subtype genes occupy the main position in different plants, which may be caused by the long-term immune evolution of plants against pathogens. In short, different plants face different pathogens and thus likely rely on different resistance genes.
One phylogenetic tree of 292 NBS-LRR genes was built in P. massoniana (Figure 1). These 292 NBS-LRR genes divided neatly into two prominent branches (TNL and non-TNL), not unlike the classification reported for other plants [11,12,[38][39][40][41]. It can be deduced from the phylogenetic tree results that the TNL-type genes evolved earlier than the non-TNLtype genes, which further supports the inference that TIR-type genes originated earlier than non-TIR-type genes [38,42]. Furthermore, the types of non-TNL genes are more diverse: CNL and RNL subclass genes clustered on the non-TNL branch, and NL subclass genes are distributed on each branch, indicating that NL-type genes have more diversified origins. The results of the gene conserved domain and conserved motif analyses (Table S3 and Figure S1) also demonstrated that genes within the same branch had similar conserved domains and conserved motifs, which strongly suggests that those genes have similar structures and functions. The NBS domain is the most conserved, while the N-terminal TIR/CC/RPW8 and C-terminal LRR domains are less conserved. Besides the conserved motifs of TIR-1-4, P-loop, Kinase2, RNBS-B, RNBS-C, GLPL, RNBS-D, and MHDV in the TNL subclass, the NBS domain lacks the RNBS-A motif identified in previous studies [43] but harbors three new conserved motifs, which may have been lost during the evolution of P. massoniana [44]. Compared with TNL subclass genes, the NL, RNL, and CNL subclass genes have a higher diversity in the arrangement and composition of domains and conserved motifs. Motif 19 and motif 24 are unique motifs in CNL and RNL subclass genes, respectively. This finding could be explained by the differential resistance genetic variation in P. massoniana in the face of different biotic or abiotic stresses. In summary, phylogenetic tree and conserved motif analyses provide compelling evidence for the complex evolutionary relationships of NBS-LRR gene family members and further explain the diversity of resistance in P. massoniana.
The expression profile of PmNBS-LRR97 in the stem was analyzed to further clarify this gene's involvement in conferring PWN resistance ( Figure 4A). When a plant is not infected by one or more pathogens, its expression of the R gene is at a low level. Once infected by pathogens, the R gene will respond quickly, and its level of expression will increase rapidly [45][46][47]. Compared with the control group (inoculated with water), the expression of PmNBS-LRR97 in resistant and susceptible P. massoniana was always greater at different times after the inoculation with the PWN, indicating that the inoculated wounds did not affect the expression of the PmNBS-LRR97. The expression of PmNBS-LRR97 in resistant and susceptible P. massoniana increased after their inoculation with the PWN, but the fold increase was significantly higher in the former than in the latter, indicating that PmNBS-LRR97 was linked to the resistance in P. massoniana to the PWN. At 15 dpi, this gene's expression level in resistant P. massoniana basically reverted to its level before the inoculation, whereas in susceptible P. massoniana, it remained significantly higher. At 30 dpi, the expression of PmNBS-LRR97 in resistant and susceptible P. massoniana was on par with that before their inoculation with the PWN.
We speculate that the PmNBS-LRR97 of resistant P. massoniana can respond quickly to PWN invasion and significantly increase its expression level to generate enough defense substances to resist PWN-induced damage at 15 dpi. Although the susceptible P. massoniana can also respond quickly to invasion, the augmented expression of its PmNBS-LRR97 is not enough to eliminate the damage incurred from the PWN, which can continue to impact P. massoniana, leaving this gene still expressed at a higher level than before the inoculation. At 30 dpi, there were almost no PWNs in the resistant P. massoniana due to its early resistance response to the PWN; in contrast, because the susceptible P. massoniana did not eliminate the PWN, resulting in wilting or even tissue death at the inoculation site, its expression level was lower than that before the inoculation. Zhang et al. [48] also found similar results in the response of grape (Vitis quinquangularis Rehd.) to Plasmoparaviticola infection. VqCN, a CC-NBS-LRR-type R gene, was significantly induced, and its expression level increased rapidly after the resistant variety was infected by P. viticola. The expression profiles of R genes in Pennisetum glaucum (L.) R.Br. in response to Sclerosporagraminicola [47], Solanum aculeatissimum in response to Meloidogyne incognita [49], and P. monticola in response to Cronartiumribicola [50] were similar to those found here for P. massoniana. Further, for both resistant and susceptible P. massoniana, PmNBS-LRR97 was expressed more in their stems than in their needles and roots at any post-inoculation time ( Figure 4B,C). This may be due to the fact that the Monochamusalternatus vector carrying the PWN mainly ate the stems of P. massoniana, enabling PWN's entry into the host via this tissue.
Knowing the subcellular localization of genes can further help in understanding their biological functions in plants [51]. Some reported R genes locate to the nucleus [52][53][54], and others locate to the cell wall, cytoplasm [55], or cell membrane [48,56]. The PmNBS-LRR97 located exclusively to the cell membrane ( Figure 4D). It could have recognized the effector proteins secreted by the PWN on the cell membrane, thus causing a defense response in the attacked P. massoniana.
The generation of ROS is an early event in the immune response caused by the recognition of effector proteins by the R gene, which is pivotal for inducing plant resistance responses, among which H 2 O 2 is quite important [57][58][59][60]. NO is also an important signaling molecule that regulates plant resistance, and it is known that an increased NO concentration can enhance resistance to pathogens [61][62][63]. The contents of H 2 O 2 and NO will both increase rapidly after a plant is attacked by pathogens. Appropriate amounts of H 2 O 2 and NO can participate in the transmission of disease resistance signals and thereby enhance resistance. However, too much H 2 O 2 and NO will upset the redox balance, impairing it. At this time, active oxygen-scavenging enzymes, such as SOD, POD, and CAT, are activated to remove the excess H 2 O 2 and NO to protect plant cells from damage [64,65]. For PmSOD, PmPOD, and PmCAT, their respective expression patterns in the resistant and susceptible P. massoniana were characterized by an 'up-down-down' trend ( Figure 6A-C). After the inoculation with the PWN, their expression levels were significantly increased at 1 dpi, but the up-regulated fold-change of resistant P. massoniana significantly surpassed that of susceptible P. massoniana. Crucially, the expression pattern of these genes was basically similar to that of PmNBS-LRR97, which indicates that PmNBS-LRR97 may cause a significant increase in ROS-related genes. Next, the expression of NbSOD, NbPOD, and NbCAT was analyzed in the WT and transgenic tobacco, and the contents of defense signal molecules (H 2 O 2 , NO) and ROS scavenging system-related enzymes (SOD, POD, CAT) were also detected ( Figure 6D-K). Evidently, the expression levels of these genes and these indicators' contents were higher in transgenic tobacco than in WT, suggesting that the overexpression of PmNBS-LRR97 can cause ROS production, activate enzymes that scavenge ROS, and render plants resistant to the PWN. Interestingly, its overexpression also caused tobacco's early flowering and promoted its growth and development ( Figure 5B-D). The findings can further explain the resistance conferred by PmNBS-LRR97 to P. massoniana against the PWN.

Plant Materials and Treatment
The test materials were procured from a 5-year-old plantation of P. massoniana clones at the Linhai Nursery (latitude: 120 • 12 55.836 E, 30 • 15 11.088 N), located in Linhai City, Zhejiang Province, China. The current-year shoots of three resistant P. massoniana clones (all Huang 13-1) and three susceptible clones (all Guang 27-2) were selected for inoculation with the PWN in August [66]. The resistant clones were clone individuals whose growth was unaffected at 30 days post-inoculation (dpi), whereas the susceptible clones had withered and died at about 30 dpi. Each shoot was inoculated with 10,000 heads/200 µL of PWN; the counterpart control groups consisted of three selected strains of resistant and susceptible clones which received the same amount of sterile water. At 0, 1, 15, and 30 dpi with PWN, the roots, stems, and needles of the resistant and susceptible clones of P. massoniana in the infected and control groups were collected as three biological replicates and stored in a refrigerator at −80 • C for further study.

Phylogenetic and Conserved Motif Analysis of NBS-LRR Genes
The protein sequences containing both NBS and LRR domains were subjected to multiple sequence alignment using ClustaW [75]. Then, a phylogenetic tree was constructed using the maximum likelihood method in MEGA 7.0 [76] with these param-eters: ML tree method, JTT+G+I+F model, partial deletion, and 1000 bootstrap replicates. The evolutionary tree was visualized and polished using EvolView v3 software (http://www.evolgenius.info/evolview (accessed on 13 August 2022)) [77].
Both the motif search and analysis of NBS-LRR were performed by MEME (https: //meme-suite.org/meme/tools/meme (accessed on 13 August 2022)) [78] under these parameter settings: site distribution; any number of repetitions; number of motifs, 25. Finally, the protein domains and conserved motifs were visualized using TBtools v1.098765 [79].

Heatmap of Differentially Expressed NBS-LRR Genes
The NBS-LRR genes from resistant and susceptible P. massoniana responsive to PWN inoculation were obtained from the transcriptome data. An FDR ≤0.05 and a fold-change ≥2 were used for the identification of the differentially expressed genes (DEGs). A hierarchical Cluster heatmap based on NBS-LRR DEGs analysis was conducted with TBtools.

Gene Cloning and Sequence Analysis
The total RNA was extracted from the resistant and susceptible P. massoniana clones using an EASY38 Spin Plus plant RNA Kit (AidLab, Beijing, China), and then each RNA sample was reverse-transcribed into cDNA with the PC18-TRUEscript 1st Strand cDNA Synthesis Kit (OneStep gDNA Removal) (AidLab, Beijing, China). According to the reference sequence of PmNBS-LRR97 in the transcriptome, Primer3 Input (https: //bioinfo.ut.ee/primer3-0.4.0/ (accessed on 25 December 2022)) [80,81] was used to design the cds-specific primers (Table S1), whose PCR amplification was carried out with Phanta Max Super-Fidelity DNA Polymerase (Vazyme, Nanjing, China). The amplified sequence was ligated to the pMD20-T vector (TaKaRa, Beijing, China), transformed into DH5α Escherichia coli, and sequenced externally by Zhejiang Shangya Biotechnology Co., Ltd. after undergoing colony PCR.

Real-Time Quantitative PCR (qRT-PCR) Analysis
The total RNA, which was extracted from the roots, stems, and needles of P. massoniana clones inoculated with the PWN and sterile water at different time points (0, 1, 15, and 30 dpi), was reverse-transcribed using PrimeScript TM RT Master Mix (TaKaRa, Beijing, China), with Primer3 Input used to design the fluorescent quantitative primers for PmNBS-LRR97, PmSOD, PmPOD, and PmCAT (Table S4) in the non-conserved domain region of the above RNA. The EF2 of P. massoniana served as an internal reference, and the ABI PRISM 7300 Real-Time PCR System (Foster City, CA, USA) was operated according to the instructions of the TB Green ® Premix Ex Taq TM II Kit (TaKaRa, Beijing, China). Three technical replicates were performed, and the relative expression levels were calculated by applying the 2 − CT method [86] to the data. The total RNA in the leaves of wild-type (WT) and transgenic tobacco was extracted and reverse-transcribed. Using NbACTIN as the internal reference gene, Primer3 Input was used to design fluorescent quantitative primers for NbSOD, NbPOD, and NbCAT (Table S1), whose relative expression levels were calculated (as described above).

Subcellular Localization of PmNBS-LRR97
To clarify the subcellular localization of the PmNBS-LRR97, the full-length sequence of the gene was ligated to the pCambia1300 expression vector with GFP at its C-terminus, and this then transformed into Agrobacterium GV3101. Using the empty vector as a control, the Agrobacterium was suspended to OD 600 = 1.2, incubated for 2 h, and then injected into the lower epidermis of tobacco leaves. Two days later, the transient expression of the GFP fusion protein was observed by the LSM900 confocal microscope imaging system (Zeiss, Oberkochen, Germany). Nicotiana benthamiana Domin was cultured in a light incubator at 25 • C, under a 16 h light/8 h dark cycle.

Generation of PmNBS-LRR97 Transgenic Lines in Tobacco
The pCambia1300-GFP expression vector containing the target gene was transformed into Agrobacterium GV3101, and the leaf disc method was used to transform N. benthamiana. In this way, the medium containing hygromycin B and PCR amplifications were performed to select the positive transgenic lines. The WT and transgenic tobacco plants' RNA was extracted and reverse-transcribed, and the expression level of each transgenic line was analyzed by qRT-PCR, using NbACTIN as the internal reference gene.

Determination of H 2 O 2 , NO, SOD, POD, and CAT in Tobacco
The contents of H 2 O 2 , NO, superoxide dismutase (SOD), peroxidase (POD), and catalase (CAT) in 1-month-old WT and transgenic tobacco leaves were respectively determined using a hydrogen peroxide content kit (enzyme-labeled method), nitric oxide content determination kit (enzyme-labeled method), superoxide dismutase kit (WST-8 method) (enzymelabeled method), peroxidase determination kit instruction (enzyme-labeled method), and CAT kit instruction (enzyme-labeled method), according to Suzhou Mengxi Biomedical Technology Co., Ltd.'s instructions, with three biological replicates used.

Conclusions
In this study, 507 NBS genes were identified and comprehensively analyzed in the P. massoniana transcriptome. The 292 NBS-LRR genes were divided into two major branches: TNL and non-TNL. PmNBS-LRR97, one of the important differentially expressed genes, is rapidly activated after the infection with the PWN, and its expression levels are higher in the stems than in the roots and needles of P. massoniana. According to its subcellular localization, PmNBS-LRR97 may recognize effector proteins on the cell membrane to elicit defense responses. PmNBS-LRR97 may cause the significant up-regulation of ROS-related genes in P. massoniana. The overexpression of PmNBS-LRR97 increases the contents of active oxygen and active oxygen scavenging enzymes in N. benthamiana, which may foster resistance in P. massoniana. Overall, these findings provide a fundamental insight into the regulation of NBS-LRR resistance genes in the defense response of P. massoniana to PWN.