RNAseq Reveals Differential Gene Expression Contributing to Phytophthora nicotianae Adaptation to Partial Resistance in Tobacco

: Phytophthora nicotianae is a devastating oomycete plant pathogen with a wide host range. On tobacco, it causes black shank, a disease that can result in severe economic losses. Deployment of host resistance is one of the most effective means of controlling tobacco black shank, but adaptation to complete and partial resistance by P. nicotianae can limit the long-term effectiveness of the resistance. The molecular basis of adaptation to partial resistance is largely unknown. RNAseq was performed on two isolates of P. nicotianae (adapted to either the susceptible tobacco genotype Hicks or the partially resistant genotype K 326 Wz/Wz ) to identify differentially expressed genes (DEGs) during their pathogenic interactions with K 326 Wz/Wz and Hicks. Approximately 69% of the up-regulated DEGs were associated with pathogenicity in the K 326 Wz/Wz -adapted isolate when sampled following infection of its adapted host K 326 Wz/Wz . Thirty-one percent of the up-regulated DEGs were associated with pathogenicity in the Hicks-adapted isolate on K 326 Wz/Wz . A broad spectrum of over-represented gene ontology (GO) terms were assigned to down-regulated genes in the Hicks-adapted isolate. In the host, a series of GO terms involved in nuclear biosynthesis processes were assigned to the down-regulated genes in K 326 Wz/Wz inoculated with K 326 Wz/Wz -adapted isolate. This study enhances our understanding of the molecular mechanisms of P. nicotianae adaptation to partial resistance in tobacco by elucidating how the pathogen recruits pathogenicity-associated genes that impact host biological activities.


Introduction
Plant diseases are estimated to cause crop losses of 13% annually, imposing a major constraint on global crop production [1]. Deployment of complete and partial resistance in host plants is one of the most effective means of managing plant diseases and is an integral part of sustainable disease management that reduces the use of fungicides and other management inputs [2]. However, wide distribution of cultivars with complete resistance places strong selection pressure on pathogen populations to overcome that resistance [3]. Partial resistance selects for isolates that are more aggressive than isolates produced on susceptible cultivars, which can erode the effectiveness of partial resistance over time [4][5][6].
Various mechanisms utilized by plant pathogens to overcome complete resistance have been recognized, including loss of avirulence (Avr) gene products that trigger plant immunity, transposon insertions or mutations to the Avr gene sequence, acquisition of additional epistatic effectors that suppress the plant immune system without disrupting the original Avr gene [7], and endogenous small RNAs silencing Avr genes [8]. Despite our rapid improvement in understanding the molecular basis underlying complete resistance Petri dishes containing oatmeal agar (Difco, Detroit, MI, USA). Petri dishes were incubated in the dark at 28 • C for approximately 2 weeks until dense hyphal mats formed. Hyphal mats were peeled from the oatmeal agar surface and placed into Petri dishes containing 20 mL of sterile 5% sandy-loam soil extract. Five percent soil extract was prepared by mixing 50 g of soil with 1 L of deionized water and letting it sit at room temperature for 48 h. The suspension was filtered through Fisher Brand Qualitative P8 filter paper and Celite 545 (Fisher Scientific, Fair Lawn, NJ, USA), and sterilized by autoclaving for two consecutive days at 121 • C for 60 min. Petri dishes were placed under constant light at room temperature in laboratory for about 5 days at which time numerous sporangia had produced. Sterile 5% soil extract was replaced daily during incubation. Zoospores were released by incubating hyphal mats at 4 • C for 1 h, followed by incubation at 28 • C for 30 min. The concentration of zoospore suspension was determined and adjusted to a concentration of 1 × 10 5 zoospores/mL using a hemocytometer.
Tobacco seeds of K 326 Wz/Wz and Hicks were seeded in potting mix (Fafard 2 Mix; Conrad Fafard, Inc., Agawam, MA, USA) in plastic pots in a greenhouse with a 35 • C/26 • C day/night temperature regime and a 14 h photoperiod supplemented with high intensity lights. After two weeks, germinated seedlings were transplanted to cell trays (cell size 3.8 cm × 3.8 cm × 5.7 cm) containing calcined clay (TURFACE ® All Sport™, PROFILE Products LLC, Buffalo Grove, IL, USA) and grown for about two weeks. Six seedlings of each genotype were removed from calcined clay, washed gently with sterile deionized water, and inoculated by immersing the roots for 3 h in 60 mL of zoospore suspension of one of the two isolates in a Petri dish (25 × 100 mm). After inoculation, seedlings were moved to a new Petri dish (25 × 100 mm) containing 60 mL of 5% soil extract and incubated under constant light at room temperature. Forty-eight hours post inoculation (hpi), roots of individual seedlings were flash frozen in liquid nitrogen in separate 1.7 mL centrifuge tubes and subjected to RNA extraction.

RNA Isolation and Transcriptome Sequencing
A total of 24 infected root samples (2 isolates × 2 tobacco genotypes × 6 biological replicates) were subjected to RNA isolation. Total RNA of infected roots of each seedling was extracted using Qiagen Plant RNeasy Kits (Qiagen, Hilden, Germany) following the manufacturer's instructions. Genomic DNA was removed by on-column digestion with DNase I (Zymo Research Corporation, Irvine, CA, USA). The concentration and quality of total RNA were determined by BioAnalyzer RNA analysis. Three of the six biological replicates of each treatment with the highest RNA quality were selected for RNA sequencing.
RNA library preparation and sequencing were conducted at the Genomic Science Library at North Carolina State University (Raleigh, NC, USA). Briefly, the RNAseq library was constructed using a NEBNext ® Ultra™ Directional RNA library Prep Kit for Illumina (New England BioLabs, Ipswish, MA, USA) using 1 µg of each of the RNA samples followed by a 350-500 bp final library size selection. Libraries of the three biological replicates for each of the four treatments (2 isolates × 2 tobacco genotypes) were multiplexed and sequenced in a single Illumina NextSeq 500 lane, generating 75 bp paired-end reads.

Analyses of RNAseq Data from P. nicotianae 2.2.1. Detection of Differentially Expressed Genes (DEGs) in P. nicotianae
Sequence quality was assessed using FastQC v0.11.8 [18]. No trimming was performed since the Phred quality score of each sequenced base was above 30 for all samples. Reads of all samples were aligned to the P. nicotianae genome (phytophthora_parasitica_inra_ 310. 3.scaffolds.fasta) using Hisat2 v2.1.0 [19] with default parameters and maximum intron length of 5000 bp.
Reads mapped to coding sequences (CDS) of annotated genes were counted using featureCounts [20] with default settings. DEGs in a given isolate were identified by comparing infected K 326 Wz/Wz samples to infected Hicks samples using edgeR [21] with Agronomy 2021, 11, 656 4 of 28 TMM normalization, a generalized linear model, and false discovery rate (FDR) calculations based on the Benjamini-Hochberg method. Genes with a false FDR < 0.05 were considered to be DEGs. DEGs were divided into up-and down-regulated groups for further analyses.

Gene Ontology Analysis, KEGG Pathway Enrichment Analysis and PHIB-Blast
Gene Ontology (GO) ID and protein sequences were linked to individual DEGs using the UniProt website (https://www.uniprot.org/. Accessed on 21 May 2019) [22]. GO term enrichment analysis was performed using the BiNGO plugin [23] in Cytoscape v3.7.1 [24]. Over-represented GO terms were evaluated against the P. nicotianae genome in the categories "biological process", "molecular function", and "cellular component". The DEGs were subjected to the Kyoto Encyclopedia of Genes and Genomes (KEGG) [25] pathway enrichment analysis to understand their roles in biological pathways using KOBAS [26] with background species set to Phytophthora infestans, statistical method set to Hypergeometric test/Fisher's exact test, and FDR correction method set to Benjamini and Hochberg. Protein sequences of the DEGs were subjected to a blast search in PHIB-Pathogen Host Interactions base [27] to identify DEGs associated with pathogenicity. A DEG was considered pathogenicity-associated if it or its ortholog was verified in association with "loss of pathogenicity", "reduced virulence", "lethal", or "effector" in pathogens with an Evalue cutoff of 1.0 × 10 −5 .

Detection of Differential Transcript Usage
To detect differential transcript usage, transcripts in each sample were assembled and quantified using StringTie [28] by comparing the BAM file (aligned using Hisat2) to the annotated reference genome (phytophthora_parasitica_inra_310.3.genes.gtf). Analysis of alternative splicing and isoform switches was conducted using IsoformSwitchAnalyzeR [29] package in R. Isoforms in each isolate of P. nicotianae found by comparing infected samples of K 326 Wz/Wz to infected samples of Hicks were considered differentially switched if difference in isoform fraction (dIF) > 0.1 and FDR corrected q-value < 0.05. Genes with differential transcript usage were subject to GO, KEGG, and PHIB blast analyses as well.

Identification of Single-Nucleotide Polymorphisms (SNPs)
To identify SNPs in the two isolates of P. nicotianae, sorted bam files of individual samples were subjected to variant calling using samtools mpileup and filtered using bcftools [30]. SNPs identified between the two isolates were located in genes. Gene sequences and corresponding protein sequences were blasted in NCBI to further identify their potential roles in aggressiveness in each of the two isolates of P. nicotianae.  [19]. Reads mapped to coding sequences (CDS) of annotated genes were counted using featureCounts [20] with default settings. DEGs in K 326 Wz/Wz inoculated with a given isolate were identified by comparing it to Hicks inoculated with the same isolate by edgeR [21] using TMM normalization, a generalized linear model, and false discovery rate (FDR) calculations based on the Benjamini-Hochberg method. Genes with a false FDR < 0.05 were considered DEGs. DEGs were divided into up-and down-regulated datasets for further analysis.

Gene Ontology and KEGG Enrichment Analyses
DEGs were subjected to GO enrichment analysis using AgriGO v2 [32] against Nitab4.5 ID (solgenomics) as background with default settings. Corresponding protein sequences of DEGs were extracted and subjected to KEGG enrichment analysis using KOBAS [26] with background species set to Nicotiana tabacum, statistical method set to Hypergeometric test/Fisher's exact test, and FDR correction method set to Benjamini and Hochberg.

Quantitative Real-Time PCR (qRT-PCR) Validations
Two up-regulated and two down-regulated DEGs detected in both of the two P. nicotianae isolates with largest fold changes (Table 1) were chosen for qRT-PCR quantification to validate the DEGs called in RNAseq analysis. RNA samples extracted from two biological samples of each treatment were used as templates in qRT-PCR validations. First strand cDNA synthesis was initiated using ProtoScript®II reverse transcriptase (New England Biolabs, Beverly, MA, USA) following first strand cDNA synthesis standard protocol NEB#M0277. The ubiquitin-conjugating enzyme (Ubc) and the 40S ribosomal protein S3A (WS21), constitutively expressed throughout P. nicotianae development stages, were used as internal control genes [33] (Table 1). A total of 20 µL of reaction solution, including 1 µL of cDNA, 10 µL of iTaq Universal SYBR Green SuperMix (BioRad, Hercules, CA, USA), 0.6 µL of forward and 0.6 µL of reverse primers (10 µM), and 7.8 µL of molecular grade water was used for qRT-PCR. qRT-PCR was performed on 96-well plates using the Applied Biosystems QuantStudioTM 6 Flex Real-Time PCR system with the following settings: one cycle of 95 • C for 20 s (hold stage), followed by 40 cycles of 95 • C for 15 s, 60 • C for 30 s (PCR stage), with a final melt curve stage: 95 • C for 15 s, 60 • C for 1 min and 95 • C for 15 s. Three technical replicates were performed for each sample and primer set combination.

Infected Tobacco Root Samples for RNAseq
Root tissue colonized by P. nicotianae was obtained for RNAseq by inoculating and harvesting the roots of seedlings 48 h post inoculation (hpi). At 48 hpi, slight browning of the roots was present and abundant sporangia were present around roots.

RNA Sequencing and Sequence Mapping to the Reference Genome
Approximately 33 million reads were obtained per sample. An average of 27% of the reads were mapped to the P. nicotianae genome and an average of 6% of the reads were mapped to the N. tabacum genome ( Table 2).

DEGs Identified in P. nicotianae
For each of the two isolates, DEGs were identified by comparing infected K 326 Wz/Wz samples to infected Hicks samples. The DEGs identified in the two isolates were compared to view the dynamics in gene expression in the two isolates after infecting their adapted and their non-adapted tobacco host genotypes.
Forty-six genes in Wz-Wz and 50 genes in Wz-H isolates were differentially expressed. Specifically, 16 up-regulated and 30 down-regulated genes were identified in the Wz-Wz isolate (Table 3), and 29 up-regulated and 21 down-regulated genes were detected in the Wz-H isolate (Table 4). qRT-PCR of the four selected DEGs indicated expression pattern consistent with those captured in RNAseq where PPTG_08585 and PPTG_20266 had a lower Agronomy 2021, 11, 656 6 of 28 expression, and PPTG_10666 and PPTG_19949 had a higher expression in the two isolates when infecting K 326 Wz/Wz comparing to when they were infecting Hicks ( Figure 1). Both isolates up-regulated genes PPTG_19949, PPTG_06767, PPTG_10666, PPTG_05470 that encode uncharacterized proteins in P. nicotianae, and down-regulated 8 genes including genes that encode 60S ribosomal protein L38, phosphoadenosine phosphosulfate reductase, and NAD(P)H:quinone oxidoreductase.
Thirty-four DEGs were detected exclusively in the Wz-Wz isolate. Most of these genes encoded for uncharacterized proteins in P. nicotianae. Up-regulated genes with known function included PPTG_02121, PPTG_08145, and PPTG_12158 that encode Hsp70like protein, 4-aminobutyrate transaminase, and ULK/ULK protein kinase. Similarly, the majority of the 38 DEGs identified only in the Wz-H isolate encoded for uncharacterized proteins. Genes with known function included an up-regulated gene, PPTG_17442, that encodes protein-S-isoprenylcysteine O-methyltransferase and down-regulated genes, PPTG_17561, PPTG_00501, PPTG_21942, and PPTG_15084 predicted to encode for a glycine cleavage system H protein, homoserine O-acetyltransferase, phosphate acetyltransferase, and TKL/DRK protein kinase.   qRT-PCR validation of the DEGs identified in P. nicotianae by RNAseq analysis. Difference in relative gene expression between groups was determined using t-test in Prism 8. Significance is indicated using * for p value < 0.05, **** for p value < 0.0001. Genes PPTG_08585 and PPTG_20266 showed a lower relative expression in P. nicotianae when infecting K 326 Wz/Wz compared to infecting Hicks. Genes PPTG_10666 and PPTG_19949 showed a higher relative expression in P. nicotianae when infecting K 326 Wz/Wz compared to infecting Hicks. The expression of the genes was normalized to Ubc and WS21 as two internal controls.

Over-Represented Gene Ontology Analysis
To obtain insight into the types of differentially expressed genes in the two isolates, GO enrichment analysis was performed using BiNGO to test over-representation for the DEGs against the annotated genes in P. nicotianae. Down-regulated genes in the Wz-Wz isolate were significantly enriched into three GO terms for the functional classes "sulfate reduction" (GO: 0019419), "sulfate assimilation, phosphoadenylyl sulfate reduction by phosphoadenylyl-sulfate reductase (thioredoxin)" (GO: 0019379) and "sulfate assimilation" (GO: 0000103). Up-regulated genes in the Wz-Wz isolate were enriched in "monooxygenase activity" (GO: 0004497) ( Figure 2). In the Wz-H isolate, down-regulated genes were enriched into a broader spectrum of GO terms. In addition to the ones found in the Wz-Wz isolate, down-regulated genes in the Wz-H isolate were also enriched in "glycine cleavage complex" (GO: 0019419), "NAD(P)H dehydrogenase (quinone) activity" (GO: 0004604), "FMN binding" (GO: 0050662), and "oxidoreductase activity" (GO: 0016651, GO: 0000103, GO: 0009071; GO: 0016671). No GO terms were particularly enriched for up-regulated genes in the Wz-H isolate (Figure 3).

KEGG Analysis of DEGs
A KEGG pathway enrichment analysis was conducted on DEGs identified in the two isolates to help understand the interaction between isolates and their adapted and non-adapted tobacco host genotypes.
For up-regulated DEGs, "spliceosome" was enriched in the Wz-Wz isolate. No enriched pathway was identified for up-regulated DEGs in the Wz-H isolate. For downregulated DEGs, "ubiquinone and other terpenoid-quinone biosynthesis" and "biosynthesis of secondary metabolites" were enriched in the Wz-Wz isolate. In addition to the two pathways enriched, another eight enriched pathways were identified in the Wz-H isolate including "biosynthesis of antibiotics", "carbon metabolism", and "sulfur metabolism" (Table 5).

PHIB-Database Blast
The pathogen-host interactions database (PHI-base) stores curated molecular and biological information on genes experimentally proven to alter the outcome of pathogen-host interactions. Given that a majority of the DEGs were identified to encode uncharacterized proteins in the annotated Phytophthora nicotianae genome, the protein sequences of the identified DEGs were subjected to a blast search in the PHI-base to determine their role or the role of their orthologs to have a better knowledge of how they could potentially contribute to aggressiveness in individual isolates.
Among the up-regulated genes in the Wz-Wz isolate, 68.75% (11/16) were associated with pathogenicity. In contrast, 31.03% (9/29) of the up-regulated genes in the Wz-H isolate were found to have a role in pathogenicity. The two isolates showed a similar percentage of down-regulated genes involved in pathogenicity, with 46.67% (14/30) in the Wz-Wz isolate and 47.62% (10/21) in the Wz-H isolate (Table S1).
Specifically, the Wz-Wz isolate up regulated pathogenicity-associated genes including, but not limited to, PPTG_10595 encoding a protein belonging to ABC transporter superfamily, PPTG_12158 encoding ULK/ULK protein kinase, PPTG_02121 encoding a heat shock protein 70, and PPTG_06886 encoding a protein within solute carrier family. Down-regulated pathogenicity-associated genes included PPTG_19261 and PPTG_00236 encoding WRKY transcription factor, PPTG_02595 encoding aldehyde dehydrogenase, and PPTG_02595 encoding a sugar transport protein (Table S1).
The Wz-H isolate up regulated pathogenicity-associated genes including PPTG_22560 which matched to a gene encoding an effector protein in P. infestans and PPTG_15982 matched to a gene encoding a glycoside hydrolase in P. palmivora. The down regulated genes in association with pathogenicity in the Wz-H isolate included PPTG_15084 encoding a TKL/DRK protein kinase and PPTG_13205 matched to a gene encoding an elicitin-like protein in P. infestans (Table S1).

Genes with Differential Transcript Usage in P. nicotianae
Differential transcript usage allows a single gene to produce multiple transcript isoforms. To explore possible molecular mechanisms other than selectively expressing genes in a specific pathogen isolate-host genotype interaction, genes with differential transcript usage were identified and analyzed in P. nicotianae. Sixty-six and 128 annotated genes in the Wz-Wz and Wz-H isolates, respectively, were differentially transcribed when the given isolate was infecting K 326 Wz/Wz compared to its infecting Hicks. Intronretention (IR), alternative transcription start sites (ATSS), and alternative transcription termination sites (ATTS) were the three most common alternative splicing types in the two isolates.
Twenty-seven and 60 differentially transcribed genes were predicted to have functional consequences in the Wz-Wz and Wz-H isolates, respectively (Table S2). These genes were subjected to GO and PHIB blast analyses. No GO terms were enriched for either of the two sets of genes identified in the two isolates. The PHIB blast results showed that 9 out of 27 (33.3%) genes with differential transcript usage in the Wz-Wz isolate were involved in pathogenicity, while 13 (21.7%) pathogenicity-associated genes were identified from the 60 genes in the Wz-H isolate (Table S3). For example, the Wz-Wz isolate alternatively spliced genes PPTG_10075 encoding a serine/threonine protein kinase and PPTG_00215 encoding eukaryotic translation initiation factor. The Wz-H isolate differentially transcribed genes such as PPTG_06129 encoding pre-mRNA 3 end processing protein and PPTG_03522 encoding ankyrin repeat protein.

SNPs Identified in the Wz-Wz and Wz-H Isolates of P. nicotianae
A total of 8 SNPs were identified between the two isolates of P. nicotianae (Table 6)   DEGs were identified in K 326 Wz/Wz by comparing it to Hicks inoculated with a given P. nicotianae isolate. When inoculated with the Wz-Wz isolate, K 326 Wz/Wz had 305 up-regulated and 303 down-regulated genes compared to 174 up-regulated and 393 down-regulated genes when inoculated with the Wz-H isolate. The DEGs identified in K 326 Wz/Wz inoculated with the two isolates were further analyzed for commonalities and differences. There were 94 up-regulated and 163 down-regulated genes in common from samples of K 326 Wz/Wz inoculated with the two isolates (Table 7; Table S4). Among the commonly up-regulated genes were Nitab4.5_0007488g0040.1 and Nitab4.5_ 0001477g0080.1 that encode pathogenesis-related (PR) protein 1a. In addition, five genes, Nitab4.5_0003154g0030.1, Nitab4.5_0000754g0140.1, Nitab4.5_0003324g0100.1, Nitab4.5_ 0014015g0010.1, Nitab4.5_0013087g0020.1, were predicted to encode proteinase inhibitors. Of the 163 commonly down-regulated genes, 13 genes were found to encode Glutathione S-transferase or Glutathione S-transferase-like protein.
No significant GO terms were enriched for exclusively up-regulated genes in K 326 Wz/Wz inoculated with the Wz-H isolate, and 8 GO terms were enriched for downregulated genes after inoculation with the Wz-H isolate (Table 9).  To investigate the potential resistance mechanism in K 326 Wz/Wz, commonly upregulated and down-regulated DEGs found across samples of K 326 Wz/Wz inoculated with each of the two isolates were analyzed for enriched GO terms. In particular, GO terms in the functional categories "serine-type endopeptidase inhibitor activity" (GO: 0004867), "peptidase inhibitor activity" (GO: 0030414), "peptidase regulator activity" (GO: 0061134), "endopeptidase inhibitor activity" (GO: 0004866), "endopeptidase regulator activity" (GO: 0061135) were enriched for the up-regulated genes (Table 10). GO terms enriched for the down-regulated genes in the functional categories included "transferase activity, transferring hexosyl groups" (GO: 0016758), "transferase activity, transferring glycosyl groups" (GO: 0016757), "ATPase activity, coupled to transmembrane movement of substances" (GO: 0042626), "ATPase activity, coupled to movement of substances" (GO: 0043492), "ATPase activity" (GO: 0016887), "ATPase activity, coupled" (GO: 0042623) ( Table 10).

KEGG Pathway Enrichment Analysis of DEGs
DEGs were subjected to KEGG pathway enrichment analysis to identify their biological roles in K 326 Wz/Wz in response to the two P. nicotianae isolates. Up-regulated DEGs in the Wz-Wzand Wz-H-inoculated samples were enriched in 35 and 18 KEGG pathways, respectively (Figures 4 and 5). Sixteen KEGG pathways were commonly enriched for up-regulated DEGs identified in K 326 Wz/Wz inoculated with each of the two isolates, with the most significantly enriched pathway being "valine, leucine and isoleucine degradation" in both Wz-Wzand Wz-H-inoculated samples. Down-regulated DEGs in the Wz-Wzand Wz-H-inoculated samples were enriched in 9 and 12 KEGG pathways, respectively (Figures 6 and 7). KEGG pathways that included "sulfur metabolism", "metabolic pathways", "glutathione metabolism", "ascorbate and aldarate metabolism", and "ABC transporters" were commonly enriched for down-regulated DEGs in K 326 Wz/Wz inoculated with each of the two isolates.

Discussion
Black shank is one of the most devastating tobacco diseases globally. The use of host resistance is the most important and effective way to manage the disease worldwide, but with reduced value of complete resistance, an integrated approach using partial resistance is needed to effectively manage the disease [10]. Partial resistance is vital to the management of black shank, like many other root diseases, due to the absence of complete resistance genes or the loss of complete resistance due to pathogen new race development. Pathogen adaptation to partial resistance has been observed in various pathosystems, including P. nicotianae and tobacco, resulting in a pathogen population that is more aggressive than wild type populations [12][13][14]. It is urgent to better understand how pathogens overcome partial resistance in host plants so that resistance deployment strategies can be optimized to preserve the durability of the resistance.
The major goal of this study was to explore the molecular mechanisms underlying adaptation by P. nicotianae to partial resistance in tobacco [14,15]. Adaptation to partial resistance involves many genes and is generally considered to be more complex than overcoming complete resistance, which can result from a single nucleotide mutation in an Avr gene [34]. To obtain a holistic view of genetic differences that occur during adaptation to partial resistance and between pathogen isolates with distinctly different aggressiveness levels on a single source of partial resistance, we kept DEGs with a false FDR < 0.05 without a specified fold change of gene expression. Aggressiveness and partial resistance are two quantitative traits involving various biological activities supported by a broad spectrum of genes in the pathogen and host plant. The cumulative effect of slight changes in multiple genes could potentially influence the outcome of the interaction between a specific pathogen isolate and host genotype.
The DEGs detected involve a broad spectrum of biological activities related to pathogenicity factors in P. nicotianae. The PHIB blast showed a much higher percentage (68.75%) of up-regulated DEGs involved in pathogenicity of the Wz-Wz isolate compared to the Wz-H isolate (31.03%), indicating that the Wz-Wz isolate more efficiently recruited pathogenicity-associated genes when infecting partially resistant K 326 Wz/Wz. Particularly, some genes with essential roles in pathogenicity were only found up-regulated in the Wz-Wz isolate. These genes included PPTG_10595 that encodes a protein belonging to the ABC transporter superfamily, and PPTG_12158 that encodes ULK/ULK, a protein kinase, which has an important role in autophagy.
ABC transporters, also known as ATP Binding Cassette transporters, are of significant importance in regulating ion transport, chromosome condensation and DNA repair, mRNA processing in eukaryotes [35]. Studies have demonstrated ABC transporters are also involved in virulence [36] and toxicant efflux [37,38] in plant pathogenic fungi. It was speculated that ABC transporters export toxic phytoalexins in pathogens, therefore, contributing to pathogenicity. Comparing to P. infestans, ABC transporter gene family in P. nicotianae was significantly expanded, which suggested their crucial roles in evolutionary host adaptation [39]. The specific function of PPTG_10595 remains unclear, and given its versatility in biological processes, it is likely a higher expression of this gene could contribute to a higher aggressiveness in P. nicotianae.
ULK/ULK protein kinase (autophagy related protein 1, Atg1) is localized at the autophagy initiation site and initiates autophagy, which is critical in cell differentiation, secondary metabolism, and programmed cell death in eukaryotes [40]. In plant pathogens, autophagy has a vital role in pathogenicity. Silencing of Atg1 highly reduced conidiation and led to a reduction or loss of pathogenicity in Magnaporthe oryzae [41], Botrytis cinerea [42], Fusarium graminearum [43]. In P. sojae, expression of multiple autophagy related protein genes was increased during infection, and autophagy was highly induced during sporangium formation and cyst germination. Silencing autophagy related genes in P. sojae significantly reduced sporulation and pathogenicity, and in some cases led to defective haustorium formation, suggesting a central role of autophagy in both the development and pathogenicity in P. sojae [44]. Little is known of Atg1 in P. nicotianae, but it is possible that the up-regulated Atg1 expression in the Wz-Wz isolate could be a major contributor to aggressiveness on the resistant host K 326 Wz/Wz.
Three significantly enriched GO terms were assigned to down-regulated genes and one to up-regulated genes in the Wz-Wz isolate. Sixteen enriched GO terms were attributed to down-regulated genes in Wz-H isolate, however, no enriched GO terms were linked to up-regulated genes. These observations indicate that a broader spectrum of biological functions were affected because of the down-regulation of the genes in the Wz-H isolate.
The three enriched GO terms, "sulfate reduction" (GO: 0019419), "sulfate assimilation, phosphoadenylyl sulfate reduction by phosphoadenylyl-sulfate reductase (thioredoxin)" (GO: 0019379), and "sulfate assimilation" (GO: 0000103) that were assigned to downregulated genes in the Wz-Wz isolate, were also enriched for the down-regulated genes in Wz-H isolate. Sulfate reduction and sulfate assimilation are two important biological processes changing sulfate into sulfide, which is then used for the synthesis of methionine, cysteine, and other metabolites [45]. The biological significance of methionine is predominantly because the methionine codon AUG is the most common start codon that initiates protein synthesis [46]. Cysteine is a strong antioxidant with the potential to trap reactive oxygen species (ROS), stabilizes the high-order structures of proteins, and serves as an active center for the bioactivity of proteins [47]. When genes involved in sulfate assimilation and sulfate reduction processes are down-regulated, their effects on the synthesis and bioactivity of proteins can be profound. The fact that these 3 GO terms were enriched for the down-regulated genes in both Wz-Wz and Wz-H isolates seems to suggest that Wz resistance in tobacco was disrupting the protein synthesis in P. nicotianae as a defense mechanism.
In addition to the GO terms mentioned above, another 13 enriched GO terms were assigned to the down-regulated genes in Wz-H isolate, including "oxidoreductase activity" (GO: 0016651, GO: 0000103, GO: 0009071; GO: 0016671). The "oxidation reduction process" catalyzed by oxidoreductases was found to be one of the two enriched GO terms distinguished the fungal pathogen Colletotrichum kahawae on coffee compared to its nonpathogenic sibling species [48], suggesting a substantial contribution of oxidoreductases in general pathogenicity in plant pathogens. Prospectively, a vital role of oxidoreductases in P. nicotianae aggressiveness is speculated. Down-regulation of genes annotated as oxidoreductases in the Wz-H isolate could potentially hamper its performance on K 326 Wz/Wz compared to Wz-Wz isolate.
Genes with differential transcript usage (DTU) also was investigated in the two isolates of P. nicotianae. Differential transcript usage is primarily the result of alternative splicing events that regulate translational processes. This allows for the formation of protein variants (isoforms) with different cellular functions or properties originating from a single gene, tremendously diversifying proteins encoded by genomes. Detection of genes with DTU has been widely used in RNAseq data analysis in human research [49]. The importance of DTU in genes in plant pathogens is largely unexplored, but a study with Pseudoperonospora cubensis showed that a functional effector protein PscRXLR1 was generated from the alternative splicing of the Psc_781.4 gene that encodes a putative multi-drug transporter [50]. This finding prompted our interest in genes with DTU in P. nicotianae. A total of 27 genes in Wz-Wz and 60 genes in Wz-H with DTU were identified with predicted functional consequences. Among those genes, 33.3% in Wz-Wz and 23.3% in Wz-H were pathogenicity-associated genes found in PHIB. How DTU in the genes could alter the aggressiveness in P. nicotianae needs to be further researched, but DTU in genes such as PPTG_00215 encoding eukaryotic translation initiation factor 1A, PPTG_06129 encoding pre-mRNA 3 end processing protein WDR33, and PPTG_17135 encoding CCR4-NOT transcription complex subunit 1 signifies its role in translational biological processes in response to Wz resistance. Evidence of the importance of transcription factors in pathogen aggressiveness was confirmed in the genome-wide-association study (GWAS) in C. kahawae [48]. A slight change in pathogen genes or transcription factors associated in gene regulatory networks has been recognized to have a profound evolutionary impact [51].
to partial resistance in tobacco were able to interrupt biological processes in host nuclei to facilitate infection.

Conclusions
Few molecular studies have been completed to identify factors that might determine the overall aggressiveness of plant pathogens. Quantitative traits such as aggressiveness are challenging to dissect. A GWAS study on aggressiveness in C. kahawae strongly suggested that aggressiveness is associated with some small effect SNPs and is not regulated by causal mutations, which would indicate that aggressiveness might be a variable and very complex trait regulated by differential gene expression and corresponding regulatory mechanisms [48]. Our study was designed to enhance our understanding of the genetic mechanisms underlying P. nicotianae adaptation to partial resistance in tobacco using dual RNAseq. Overall, results from this study suggest that isolates of P. nicotianae adapted to partial resistance are able to recruit a high percentage of pathogenicity-associated genes when infecting a partially resistant genotype, and are more tolerant to the defenses expressed by Wz resistance. Finally, isolates adapted to the source of partial resistance used were able to severely hinder nuclear synthesis processes in K 326 Wz/Wz.
Wz resistance in K 326 Wz/Wz potentially disrupts protein synthesis in P. nicotianae as a defense mechanism. A broad spectrum and a high level of expression of defense related genes were recruited by K 326 Wz/Wz to inhibit non-adapted isolates of P. nicotianae compared to the adapted isolate. This was confirmed by the observation of a wide range of biological activities were affected by the down-regulated DEGs in the non-adapted isolate on K 326 Wz/Wz.
It would be beneficial if additional studies involving more isolates with distinct aggressiveness levels could be used to confirm our findings. Notwithstanding, sets of differentially expressed genes and genes with differential transcript usage were generated and can be researched via additional functional analyses to substantiate their roles in P. nicotianae aggressiveness. These findings provide a foundation for further investigation of the molecular mechanisms underlying pathogen adaptation to partial resistance in host plants.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/agronomy11040656/s1, Table S1: PHIB blast results for DEGs identified in the Wz-Wz and Wz-H isolates of Phytophthora nicotianae, Table S2: DEGs identified in the Wz-Wz and Wz-H isolates of Phytophthora nicotianae with differential transcript usage (DTU), Table S3: PHIB blast results for DEGs identified in Wz-Wz and Wz-H isolates of Phytophthora nicotianae with differential transcript usage (DTU), Table S4: DEGs identified in K326 Wz/Wz inoculated with either Wz-Wz or Wz-H isolate of Phytophthora nicotianae.  Data Availability Statement: The datasets generated during the current study are available in the GEO repository with GEO accession GSE168516.