Transcriptome Profiling of Potato (Solanum tuberosum L.) Responses to Root-Knot Nematode (Meloidogyne javanica) Infestation during A Compatible Interaction

Root-knot nematode (RKN) Meloidogyne javanica presents a great challenge to Solanaceae crops, including potato. In this study, we investigated transcriptional responses of potato roots during a compatible interaction with M. javanica. In this respect, differential gene expression of Solanum tuberosum cultivar (cv.) Mondial challenged with M. javanica at 0, 3 and 7 days post-inoculation (dpi) was profiled. In total, 4948 and 4484 genes were detected, respectively, as differentially expressed genes (DEGs) at 3 and 7 dpi. Functional annotation revealed that genes associated with metabolic processes were enriched, suggesting they might have an important role in M. javanica disease development. MapMan analysis revealed down-regulation of genes associated with pathogen perception and signaling suggesting interference with plant immunity system. Notably, delayed activation of pathogenesis-related genes, down-regulation of disease resistance genes, and activation of host antioxidant system contributed to a susceptible response. Nematode infestation suppressed ethylene (ET) and jasmonic acid (JA) signaling pathway hindering JA/ET responsive genes associated with defense. Genes related to cell wall modification were differentially regulated while transport-related genes were up-regulated, facilitating the formation of nematode feeding sites (NFSs). Several families of transcription factors (TFs) were differentially regulated by M. javanica infestation. Suggesting that TFs play an indispensable role in physiological adaptation for successful M. javanica disease development. This genome-wide analysis reveals the molecular regulatory networks in potato roots which are potentially manipulated by M. javanica. Being the first study analyzing transcriptome profiling of M. javanica-diseased potato, it provides unparalleled insight into the mechanism underlying disease development.


Introduction
Potato, Solanum tuberosum (L), belongs to the Solanaceae family, which comprises several economically important crops such as tomato, pepper, aubergine, and tobacco. Plant parasitic nematodes, particularly root-knot nematodes (RKNs) are among the most destructive and economically important pests of potatoes worldwide [1,2]. In this context, Meloidogyne spp. are obligate and highly polyphagous pests that form an intricate relationship with their host, causing drastic morphological and physiological changes in plant cell architecture [3].
Taylor et al. [23] ranking scale was used to determine susceptibility while Sasser et al. [24] reproduction factor (RF) formula was used to assess the host status of potato cultivars. For RNA experiments, whole root tissues of a compatible potato cultivar were collected at 0, 3, and 7 days post-inoculation (dpi) with two biological replicates per time point. Samples were washed and immediately frozen in liquid nitrogen to prevent RNA degradation and later stored at −80 • C until RNA extraction.

RNA Extraction, Library Preparation, and Sequencing
RNA extraction, library preparation, and sequencing were carried out at Novogene (HK) Company Limited. Total RNA for individual time course and replicates was extracted using TiaGen extraction kit (Biotech Beijing Co., Ltd., Beijing, China) and treated with sigma DNase1 (D5025). RNA degradation and contamination was measured on 1% agarose gel while RNA purity was assessed using the NanoPhotometer ® spectrophotometer (IMPLEN, Westlake Village, CA, USA). RNA concentration and integrity were assessed using Qubit ® RNA Assay kit in Qubit ® 2.0 Fluorometer (Life Technologies, Carlsbad, CA, USA) and RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA), respectively. Three micrograms of RNA samples were used as input for library construction. Libraries were constructed using NEBNext ® Ultra™ RNA Library Prep Kit for Illumina ® (NEB, Ipswich, MA, USA) according to the manufacturer's instructions and index codes were added to attribute sequences to each sample. Finally, PCR products were purified using AMPure XP system and quality of the library assessed using the Agilent Bioanalyzer 2100 system. A cBot Cluster Generation System using HiSeq PE Cluster Kit cBot-HS (Illumina) was used to cluster the index-coded samples. After cluster generation, the library preparations were sequenced on an Illumina Hiseq platform 2500 generating 150-bp paired-end reads.

Transcriptomic Data Analysis
Quality analysis of sequenced reads were initially analyzed using FASTQC package (https://www. bioinformatics.babraham.ac.uk/projects/fastqc). Clean reads were obtained by removing reads containing adapter reads with poly-N and low-quality reads from raw data. Trimming of low-quality regions was performed using Trimmomatic v 0.36 [25]. All the subsequent downstream analyses were based on high-quality data. Solanum tuberosum genome v4.03 [26] was used for reference-guided mapping of RNA-seq reads. Paired-end clean reads were aligned to the potato genome using hisat2 v 2.1.0 software [27]. Unmapped reads were progressively trimmed at the 3 end and re-mapped to the genome. Next, featureCounts package [28] was used to perform raw-reads counts in R environment (https://www.r-project.org/). The read counts were then used for differential expression analysis using edgeR package [29]. Further, to investigate the responses at different time points (3 dpi and 7 dpi), the expression profiles were compared to mock-inoculated (0 dpi) data sets. The transcripts were then classified as differentially expressed genes (DEGs) based on both (a) false discovery rate (FDR) [30] cut-off of 0.05 and (b) log 2 fold change ≥ 1 or ≤−1 for induced and repressed genes, respectively.

Gene Ontology (GO) and Enrichment Analysis
The GO and enrichment analysis were performed using agriGO v.2.0 [31] and categorized by WEGO v 2.0 tool [32]. Parametric gene set enrichment analysis based on differential expression levels (log 2 fold change) was performed and FDR correction was performed using the default parameters to adjust the p-value. Functional annotations and pathway analyses were obtained through sequence search performed on eggNOG database utilizing eggmapper [33]. Annotations from eggNOG were then integrated with Kyoto Encyclopedia of Genes and Genomes (KEGG) database to reach pathway annotation level. In addition, we used MapMan software to visualize the biotic stress pathway to identify the genes that are known to be part of the cascade of defense signals in response to pathogen or pest invasion [34].

Validation for DEGs by qRT-PCR
For qRT-PCR, first-strand cDNA was done from total RNA using Superscript IV First-Strand cDNA Synthesis SuperMix Kit (Invitrogen, Carlsbad, CA, USA) following manufacturer's protocol. Quantitative real-time PCR was performed using SYBR Green Master Mix in the QuantStudio 12k Flex Real-Time PCR system (Life Technologies, Carlsbad, CA, USA) to validate DEGs. Two micrograms of the sample (concentration) were added to 5 µL of Applied Biosystems SYBR Green Master Mix and primers at a concentration of 0.4 µM. The amplification cycle consisted of the following: initial denaturation at 50 • C for 5 min and 95 • C for 2 min followed by 45 cycles of 95 • C for 15 s and 60 • C for 60 s. Each sample was run in triplicates. Specific qRT-PCR primers for six target genes were designed using an online tool Prime-Blast (http://www.ncbi.nlm.nih.gov/tools/primer-blast) (Table S1). Each sample was run in triplicates. The 18S rRNA and elongation factor 1-α (PGSC0003DMG400020772,ef1α), [35] were used as the reference genes for normalization, and the mock-treated samples used as calibrators. The comparative 2 −∆∆Ct method was used to determine the relative fold change according to Schmittgen and Livak [36]. Despite the two techniques (RNA-seq and qRT-PCR) being different, the expression patterns of selected genes upon nematode infestation was consistent between the two procedures ( Figure S1).

Data Access
Both raw and processed sequencing data have been deposited to the Gene Expression Omnibus (GEO) repository at the National Center for Biotechnology Information (NCBI) with accession no. GSE134790.

Statistical Analyses
To identify differentially expressed genes with significant differences between the control (0 dpi) and 3 and 7 dpi infestation stages, statistical analyses were performed in R environment using edgeR package. The p-values were corrected using the Benjamini and Yekutieli [30] approach.

Analysis of Meloidogyne javanica Infestation and Functional Annotation of Differentially Expressed Genes
To establish whether Solanum tuberosum cv. Mondial is susceptible or tolerant to M. javanica, five-week-old potato plants were inoculated with 1000 J2s eggs and the control mock-inoculated with sterile water. At eight weeks post-inoculation, the results showed that S. tuberosum cv. Mondial plants infested with M. javanica have a reproduction factor (Rf > 1) of 3.11 and galling index of 5 ( Figure 1A).
Based on the GI (>100 per root system) S. tuberosum cv. Mondial potato cultivar was classified as highly susceptible. Figure 1B,C illustrate mature galls at eight weeks after inoculation exhibited either as single or egg masses, confirming M. javanica infectivity and ability to reproduce and complete their life cycle on this potato cultivar.
To understand the molecular basis of this compatible interaction between M. javanica and S. tuberosum cv. Mondial, RNA sequencing was performed at two infestation stages: 3 and 7 dpi. These time points correspond to nematode stages of induction of feeding sites at 3 dpi and nutrient acquisition stage that starts from 7 dpi to 8 weeks after infestation [5]. Approximately 1.3 billion paired-end reads were generated yielding an average of 23 million high-quality reads for individual samples. Successfully mapped reads onto the S. tuberosum reference genome (v4.03) [26] accounted for 78-86% of the total number of reads generated per sample (Table S2). Log 2 fold change ≥ ±1 and adjusted p-value (FDR) < 0.05 were used as cut off values to obtain DEGs through pairwise comparison between the mock-inoculated and infested samples at 3 and 7 dpi. Overall, 4948 genes were differentially expressed at 3 dpi. Of these, 2867 were down-regulated and 2081 up-regulated. At 7 dpi, fewer genes (4484) were differentially expressed compared to 3 dpi. The number of down-regulated genes at 7 dpi (2871) was similar to that observed at 3 dpi (2867); however, there were 22% fewer genes up-regulated at 7 dpi compared to those at 3 dpi ( Figure 2A and Table S3). Based on the GI (>100 per root system) S. tuberosum cv. Mondial potato cultivar was classified as highly susceptible. Figure 1 B,C illustrate mature galls at eight weeks after inoculation exhibited either as single or egg masses, confirming M. javanica infectivity and ability to reproduce and complete their life cycle on this potato cultivar.
To understand the molecular basis of this compatible interaction between M. javanica and S. tuberosum cv. Mondial, RNA sequencing was performed at two infestation stages: 3 and 7 dpi. These time points correspond to nematode stages of induction of feeding sites at 3 dpi and nutrient acquisition stage that starts from 7 dpi to 8 weeks after infestation [5]. Approximately 1.3 billion paired-end reads were generated yielding an average of 23 million high-quality reads for individual samples. Successfully mapped reads onto the S. tuberosum reference genome (v4.03) [26] accounted for 78-86% of the total number of reads generated per sample (Table S2). Log2 fold change ≥ ± 1 and adjusted p-value (FDR) < 0.05 were used as cut off values to obtain DEGs through pairwise comparison between the mock-inoculated and infested samples at 3 and 7 dpi. Overall, 4948 genes were differentially expressed at 3 dpi. Of these, 2867 were down-regulated and 2081 up-regulated. At 7 dpi, fewer genes (4484) were differentially expressed compared to 3 dpi. The number of downregulated genes at 7 dpi (2871) was similar to that observed at 3 dpi (2867); however, there were 22% fewer genes up-regulated at 7 dpi compared to those at 3 dpi ( Figure 2A and Table S3). To understand the molecular basis of this compatible interaction between M. javanica and S. tuberosum cv. Mondial, RNA sequencing was performed at two infestation stages: 3 and 7 dpi. These time points correspond to nematode stages of induction of feeding sites at 3 dpi and nutrient acquisition stage that starts from 7 dpi to 8 weeks after infestation [5]. Approximately 1.3 billion paired-end reads were generated yielding an average of 23 million high-quality reads for individual samples. Successfully mapped reads onto the S. tuberosum reference genome (v4.03) [26] accounted for 78-86% of the total number of reads generated per sample (Table S2). Log2 fold change ≥ ± 1 and adjusted p-value (FDR) < 0.05 were used as cut off values to obtain DEGs through pairwise comparison between the mock-inoculated and infested samples at 3 and 7 dpi. Overall, 4948 genes were differentially expressed at 3 dpi. Of these, 2867 were down-regulated and 2081 up-regulated. At 7 dpi, fewer genes (4484) were differentially expressed compared to 3 dpi. The number of downregulated genes at 7 dpi (2871) was similar to that observed at 3 dpi (2867); however, there were 22% fewer genes up-regulated at 7 dpi compared to those at 3 dpi ( Figure 2A and Table S3).  Table S3). One possible explanation of this scenario is that down-regulation of genes might be essential for proper formation of galls induced by RKN as reported previously by Jammes et al. [37]. Secondly, as an obligate biotroph, M. javanica establishes an intricate relationship with its host; therefore, a larger suppression of defense-related genes is also plausible.  Table S3).
One possible explanation of this scenario is that down-regulation of genes might be essential for proper formation of galls induced by RKN as reported previously by Jammes et al. [37]. Secondly, as an obligate biotroph, M. javanica establishes an intricate relationship with its host; therefore, a larger suppression of defense-related genes is also plausible.
GO enrichment analyses revealed the main regulatory trends in root tissues in response to M. javanica infestation. The GO terms were grouped into three main functional categories at adjusted p-value < 0.05 and categorized using WEGO software [32]. Within the biological process class, the highest percentage of the DEGs was down-regulated and fell under metabolic process category. Within this category, we found the following sub-categories: Primary and cellular metabolic process, biosynthetic and oxidation-reduction process as well as regulation of metabolic processes. Other significant GO terms in this category include response to stimulus, cellular process, localization and signalling processes, and regulation of biological processes ( Figure 3 and Table S4).  Further, MapMan analysis using the biotic stress overview classified potato genes putatively involved in mediating responses to nematode attack. A total of 1471 and 1408 genes at 3 and 7 dpi, respectively, belonging to hormone signaling, cell wall organization, proteolysis, redox homeostasis, signaling, transcriptional regulation as well as defense-related genes (secondary metabolites, PR proteins, heat shock proteins, and proteinase inhibitors) were found to be under differential regulation following M. javanica infestation ( Figure 4A, B). Further, MapMan analysis using the biotic stress overview classified potato genes putatively involved in mediating responses to nematode attack. A total of 1471 and 1408 genes at 3 and 7 dpi, respectively, belonging to hormone signaling, cell wall organization, proteolysis, redox homeostasis, signaling, transcriptional regulation as well as defense-related genes (secondary metabolites, PR proteins, heat shock proteins, and proteinase inhibitors) were found to be under differential regulation following M. javanica infestation ( Figure 4A,B).  Genes related to redox state, signaling, beta-glucanase and cell wall displayed similar expression patterns at both stages 3 and 7 dpi. Genes encoding for peroxidases, ERF, and WRKY were downregulated mostly at 7 dpi compared to the 3-dpi infestation stage. Genes related to secondary metabolism and proteolysis exhibited a general up-regulation at 3 dpi in comparison to 7dpi ( Figure  4A, B). Genes related to redox state, signaling, beta-glucanase and cell wall displayed similar expression patterns at both stages 3 and 7 dpi. Genes encoding for peroxidases, ERF, and WRKY were down-regulated mostly at 7 dpi compared to the 3-dpi infestation stage. Genes related to secondary metabolism and proteolysis exhibited a general up-regulation at 3 dpi in comparison to 7dpi ( Figure 4A,B).

Pathogen Perception and Regulation of Defense Response Genes by Meloidogyne javanica Infestation
Plant recognition of M. javanica presence, migration, and feeding is the initial step in triggering activation of downstream signaling cascades to induce defense responses. Large arsenals of plant receptors are responsible for recognition of phytonematodes and subsequent induction of basal plant defenses [38,39]. Eight genes encoding for membrane-localized proteins, mildew resistance locus o (MLO) family protein were down-regulated at 3 and/or 7 dpi ( Figure 5, Table S5). Plant recognition of M. javanica presence, migration, and feeding is the initial step in triggering activation of downstream signaling cascades to induce defense responses. Large arsenals of plant receptors are responsible for recognition of phytonematodes and subsequent induction of basal plant defenses [38,39]. Eight genes encoding for membrane-localized proteins, mildew resistance locus o (MLO) family protein were down-regulated at 3 and/or 7 dpi ( Figure 5, Table S5). Additionally, an important group of pattern recognition receptors (PRRs) was differentially regulated with the highest proportion down-regulated (62.86% of 175 genes). These include leucinerich repeat receptor-like protein kinase (LRR-RLK), receptor-like protein kinases (RLP), and wallassociated receptor kinases (WAKs) (Table S5) indicating the importance of PRRs in mediating plant immunity to M. javanica. Among the down-regulated LRR-RLK genes were two genes encoding for BAK1-interacting receptor-like kinase 1. Further, we found an LRR-RLK encoding gene, Additionally, an important group of pattern recognition receptors (PRRs) was differentially regulated with the highest proportion down-regulated (62.86% of 175 genes). These include leucine-rich repeat receptor-like protein kinase (LRR-RLK), receptor-like protein kinases (RLP), and wall-associated receptor kinases (WAKs) ( M. javanica. Among the down-regulated LRR-RLK genes were two genes encoding for BAK1-interacting receptor-like kinase 1. Further, we found an LRR-RLK encoding gene, polygalacturonase inhibiting protein 2 (PGIP 2) down-regulated 7-fold at 3 dpi and 4-fold at 7 dpi ( Figure 5 and Table S5). In addition to suppression of a cell wall-associated kinase, plant elicitor peptides (PEP1 receptor, three genes) encoding genes that perceive damage-associated molecular proteins (DAMPs) were down-regulated at both 3 and 7 dpi ( Figure 5 and Table S5). These results suggest that M. javanica likely interferes with DAMP-mediated immunity subduing activation defense response.
Transmission of perceived signals from the PRRs is mediated through the MAPK cascade and calcium (Ca 2+ ) signaling pathway which transfers downstream components of plant immunity. In this regard, we found the expression of MAPKs genes was largely repressed (9 out of 11 genes) by nematode infestation ( Figure 5 and Table S5). Furthermore, M. javanica repressed 88% of the genes involved in Ca 2+ signaling pathway (66 out of 75 genes) ( Table S5).
The intracellular nucleotide-binding domain leucine-rich repeat (NBS-LRR) proteins are important in recognizing PPNs effectors leading to effector-triggered immunity. The majority of these NBS-LRR proteins are a major class of disease resistance (R) proteins encoding resistance R-genes [39]. In the current study, we identified 32 genes encoding for either toll-interleukin 1 receptor (TIR) domain or coiled-coil (CC), which are the two subgroups of NBS-LRR proteins. Of these, 20 genes were down-regulated at 3 and/or 7 dpi following nematode infestation interfering with the host immune response ( Figure 5 and Table S5).
PTI activation is associated with expression of pathogenesis-related (PR) proteins which defense-related genes expressed upon infection by diverse pathogens including nematode pests. These PR proteins contribute to systemic acquired resistance (SAR) and are induced through the action of signaling compounds such as SA, JA, and ET [40]. Generally, we detected a delayed activation of the PR genes as the majority of the PR encoding genes were up-regulated at 7 dpi ( Figure 5 and Table S5). This could reflect a strategy adopted by the M. javanica to suppress defense response associated with PR proteins at early stages of invasion to ensure successful nematode infection. However, the induction of these PR proteins at 7dpi could have other roles in mediating a successful M. javanica infestation. For instance, we identified nine genes coding osmotin34 (PR-5), where 7 genes were up-regulated at 7 dpi ( Figure 5 and Table S5). Additionally, 11 genes encoding basic chitinase and chitinase-like family protein (PR-3 and PR-4) were detected with nine genes up-regulated at 7dpi ( Figure 5 and Table S5).
Rapid generation of reactive oxygen species (ROS) is one of the early PTI cellular events that trigger several defense responses such as activation of several defense genes and cell wall reinforcement [39]. In this study, respiratory burst homolog (RBOHs) and Riboflavin synthase-like superfamily which are important players in production of ROS in plants, were down-regulated both at 3 and/or 7 dpi ( Figure 5 and Table S5). In a resistant tomato genotype, ROS was strongly induced in response to RKN infestation mediating the R-gene resistant response [41]. Moreover, genes coding for peroxidases (50 genes) were detected mostly down-regulated (70%) during disease development. Emerging evidence shows that RKN can utilize the host ROS scavenging system to reduce the damaging effects of oxygen species [42,43]. Here, we detected one gene encoding for peroxiredoxin (PGSC0003DMG401002721), the main detoxifying antioxidant enzyme in the plant-nematode interface [44], being up-regulated at both time points. Glutathione S transferase (GSTs) encoding genes were also identified in this study. The majority of these were up-regulated at both time points (of these, 19 genes out of 25 were up-regulated ( Figure 5 and Table S5).
The production of anti-nematode compounds in the form of secondary metabolites play a critical role in induced plant immune responses against PPNs [39]. Analysis of genes related to secondary metabolism alterations was revealed in phenlyprononaids, lignin, flavonoids, isoprenoid, and phenols biosynthetic pathways ( Figure S2). The three genes encoding for key enzymes (phenylalanine ammonia-lyase (PAL), cinnamate 4-hydroxylase (CH4), and 4-coumarate CoA-ligase 2 (4CL) essential for the synthesis of phenylpropanoids and downstream synthesis of other metabolites such as flavonoids, lignin, and SA [45] were down-regulated in this study (Table S5). Other downstream genes detected to be differentially regulated include chalcone synthase, chalcone-flavone, and cinnamyl alcohols dehydrogenase involved in flavonoids and lignin pathway (Table S5). This can imply that nematode interacts with phenylpropanoid metabolism to interfere with defense compounds that originate from this biosynthetic pathway.
Further, genes encoding for enzymes involved in isoprenoid biosynthetic pathway such as hydroxymethylglutaryl-CoA synthase and hydroxy methylglutaryl CoA reductase 1 were identified with more genes detected at 7 dpi than 3 dpi. A higher proportion of these genes were largely up-regulated, 34 out of 48 genes (Table S5). Studies show that isoprenoids have various key roles in plant physiology including synthesis of key phytohormones (cytokinins, abscisic acid, gibberellins, and brassinosteroids) as well as plant defense compounds [46]. Taken together, we can hypothesize that M. javanica alters isoprenoid biosynthetic pathway to modulate hormone signaling which determines the outcome of plant-nematode interaction as reviewed by Gheysen and Mitchum [47]. Furthermore, up-regulation of genes related to secondary metabolism suggests mounting of a defense response as a result of damage caused by M. javanica feeding.

Nematode Responsive Phytohormones and Transcription Factors
Plant hormones are key players in mediating plant defense following nematode attack. Similarly, in this study M. javanica infestation influenced the expression of several genes associated with the synthesis and signaling of JA, SA, ET, auxin, abscisic acid (ABA), and brassinosteroids (BR) (Figures 4 and 6). This can imply that nematode interacts with phenylpropanoid metabolism to interfere with defense compounds that originate from this biosynthetic pathway. Further, genes encoding for enzymes involved in isoprenoid biosynthetic pathway such as hydroxymethylglutaryl-CoA synthase and hydroxy methylglutaryl CoA reductase 1 were identified with more genes detected at 7 dpi than 3 dpi. A higher proportion of these genes were largely upregulated, 34 out of 48 genes (Table S5). Studies show that isoprenoids have various key roles in plant physiology including synthesis of key phytohormones (cytokinins, abscisic acid, gibberellins, and brassinosteroids) as well as plant defense compounds [46]. Taken together, we can hypothesize that M. javanica alters isoprenoid biosynthetic pathway to modulate hormone signaling which determines the outcome of plant-nematode interaction as reviewed by Gheysen and Mitchum [47]. Furthermore, up-regulation of genes related to secondary metabolism suggests mounting of a defense response as a result of damage caused by M. javanica feeding.

Nematode Responsive Phytohormones and Transcription Factors
Plant hormones are key players in mediating plant defense following nematode attack. Similarly, in this study M. javanica infestation influenced the expression of several genes associated with the synthesis and signaling of JA, SA, ET, auxin, abscisic acid (ABA), and brassinosteroids (BR) (Figures 4, 6). The salicylic acid pathway is effective against biotrophic parasites while JA and ET pathway synergy is associated with enhanced resistance against necrotrophic microbes and herbivorous insects [48]. In this study, 9 out of 12 genes associated with SA biosynthesis were up-regulated ( Figure  6A). Moreover, genes encoding for PR-1 (five genes) a robust marker of SA responsive genes expression were up-regulated at 3 and/or 7 dpi ( Figure 5D and Table S5), indicating activation of systemic induced resistance. Genes encoding for enzymes involved in the JA synthesis allene oxide synthase (AOS 2 genes), allene oxide cyclase (AOC, one gene), lipoxygenase (LOX, four genes), and The salicylic acid pathway is effective against biotrophic parasites while JA and ET pathway synergy is associated with enhanced resistance against necrotrophic microbes and herbivorous insects [48]. In this study, 9 out of 12 genes associated with SA biosynthesis were up-regulated ( Figure 6A). Moreover, genes encoding for PR-1 (five genes) a robust marker of SA responsive genes expression were up-regulated at 3 and/or 7 dpi ( Figure 5D and Table S5), indicating activation of systemic induced resistance. Genes encoding for enzymes involved in the JA synthesis allene oxide synthase (AOS 2 genes), allene oxide cyclase (AOC, one gene), lipoxygenase (LOX, four genes), and 12-oxophytodienoate (12-OPR, three genes) were down-regulated at both time points ( Figure 6B and Table S5). Consequently, down-regulation of LOX and 12-OPR enzymes might have played a role in initiating a susceptible interaction through interfering with JA-mediated defense pathway. The LOX pathway mediates resistance against phytopathogens including nematodes [49]. Plants incapable of producing JA or 12-oxo-phytodieonoic acid (OPDA) are more susceptible to phytonematodes [50]. Moreover, four protease inhibitors coding genes were down-regulated in this study at 3 and 7 dpi ( Figure 5D and Table S5). Jasmonic acid enhances expression of protease inhibitors which prevent proteolytic activity of the insect's enzymes to debilitate their growth and reproduction. Nematodes essentially, depend on proteases to acquire nutrients as their source of food [47]. Hence, it is likely that JA plays a role in mediating defense against M. javanica. The antagonistic effect of SA on the JA pathway is often expressed at the gene transcription level, where SA targets the JA signaling. For instance, WRKY50 plays a role in SA mediated suppression [48] was up-regulated in the current study ( Figure 7C and Table S5). Furthermore, TGA, sub-class of basic leucine zipper (bZIP) family of TFs are important molecular players in mediate JA suppression through SA. In this study, two genes coding for TGA4 were up-regulated (Table S5) further signifying suppression of JA pathway through induction of SA signaling pathway. Taken together, we can deduce that M. javanica adopts the strategy of targeting the SA pathway which concurrently suppresses JA-dependent defense responses resulting in potato susceptibility.
In addition to suppression of JA defense pathway, genes involved in ET synthesis, 1-aminocyclopropane-1-carboxylate synthase, and 1-aminocyclopropane-1-carboxylate were differentially regulated in this study with 4 and 5 genes down-regulated, respectively ( Figure 7C and Table S5). We also detected a gene encoding for ETHYLENE RESPONSE SENSOR 1 up-regulated ( Figure 7C and Table S5) which acts as negative regulators of ET responses suggesting ET suppression. Previous research demonstrated that ET inhibits RKN infestation, probably through a decrease in nematode attraction to the roots [47]. In the same way, rice resistant plants show more up-regulation of ET biosynthesis and response genes than susceptible plants infested with RKN [14,51].
Transcriptional reprogramming is a key hallmark of plant innate immunity. Members of several families of TFs (WRKY, bZIP, bHLH, NAC, MYB, and AP2/ERF) are well-known in regulating plant immune response [52]. In this study, we detected differential regulation of genes encoding for ERF, WRKY, MYB, bZIP, and DOF family of TFs among others ( Figure 7A). Classification and identification of the differentially expressed TFs was attained from the Plant Transcription Factor Database (http://planttfdb.cbi.pku.edu.cn/ v.4.0) [53]. Most of the differentially expressed TFs in the current data set were down-regulated (298/532) after nematode infestation (Table S6). The importance of TFs in mediating potato disease resistance to Pectobacterium brasiliense 1692 and Phytophthora infestans infestation has been reported previously [17,18]. This is an indication that TFs are key regulators of potato immune response to various biotic stresses.
In this study, most of the genes (70.37%) encoding for AP2/ERF TF family were down-regulated at both 3 and/or 7 dpi. In addition, 15 genes were suppressed to a greater extent at 7 than 3 dpi ( Figure 7B and Table S6).
This could be ascribed to the secretion of nematode effectors and subsequent suppression of defense response associated with the activation of AP2/ERF TFs. The ERF subfamily of AP2/ERF TFs are important transcriptional regulators of genes responsive to biotic stress, particularly genes related to the JA and ET hormone signaling pathway [52]. Interestingly, we found three genes encoding for ERF6 up-regulated following M. javanica infestation ( Figure 7B and Table S6). This could imply that ERF6 TFs have a negative role in mediating potato susceptibility to M. javanica. Furthermore, seven genes coding for DREB were down-regulated in this study ( Figure 7B and Table S6).
It is generally accepted that pathogen-directed modulation of WRKY genes in plants is an important aspect that enhances success rates of pathogen infestation. Cyst nematode's successful infestation process in A. thaliana roots was attributed to the nematode's control over the expression of WRKY genes [54]. In agreement with this notion, we found 23 genes down-regulated WRKY-encoding genes, including WRKY40, WRKY23, and WRKY29 both infection stages ( Figure 7C and Table S6). The MYB TFs have been characterized as an important regulator of both biotic and abiotic stresses [55]. Among the 34 down-regulated MYB TFs in this study, we found genes encoding for MYB108 (3 genes), MYB105, and MYB 14 at 3 and 7 dpi ( Figure 7D and Table S6).
We also identified genes encoding for DNA-binding with one finger (DOF) TFs (14 out 17 genes up-regulated ( Figure 7E and Table S6), transcriptional regulators acting on the formation and function of vascular tissues [56]. This further highlights the importance of vascular-related genes in establishment and formation of NFSs, therefore successful M. javanica colonization. Overall, the large number of differentially expressed TFs in this study reflects the complexity associated with plant defense regulation which orchestrates plant immunity leading to physiological alteration of the host to favor M. javanica infestation. This could be ascribed to the secretion of nematode effectors and subsequent suppression of defense response associated with the activation of AP2/ERF TFs. The ERF subfamily of AP2/ERF TFs are important transcriptional regulators of genes responsive to biotic stress, particularly genes related to the JA and ET hormone signaling pathway [52]. Interestingly, we found three genes encoding for ERF6 up-regulated following M. javanica infestation ( Figure 7B and Table S6). This could imply that ERF6 TFs have a negative role in mediating potato susceptibility to M. javanica. Furthermore, seven genes coding for DREB were down-regulated in this study ( Figure 7B and Table S6).
It is generally accepted that pathogen-directed modulation of WRKY genes in plants is an important aspect that enhances success rates of pathogen infestation. Cyst nematode's successful

Cell Wall Organisation and Transport Processes Regulation by Meloidogyne javanica
This study showed that genes related to cell wall synthesis, degradation, modification, and cell wall proteins were differentially expressed following nematode attack ( Figure 8A and Table S6).
to favor M. javanica infestation.

Cell Wall Organisation and Transport Processes Regulation by Meloidogyne javanica
This study showed that genes related to cell wall synthesis, degradation, modification, and cell wall proteins were differentially expressed following nematode attack ( Figure 8A and Table S6). In addition, genes encoding for hydrolytic enzymes involved in pectin degradation such as polygalacturonase (PG), pectate lyases (PL), and pectin esterase (PE) were differentially expressed at 3 and 7 dpi by nematode infestation (Table S5). Further, genes encoding for pectin methylesterase inhibitor (PMEI) were largely down-regulated (10 out of 12 genes) ( Figure 8A and Table S5). More PMEI genes were suppressed at 7 dpi compared to 3 dpi. Repression of PMEI by nematode attack shows that the activity of PME was activated leading to the breakdown of pectin bonds, increasing the vulnerability of the cell wall to microbial pectic enzymes, and other degrading enzymes, culminating to a susceptible response. We also detected xyloglucan endotransglucosylase/hydrolase and expansin encoding genes, cell wall modifying enzymes under differential regulation following nematode attack. Of these, 19 out of 21 genes encoding for xyloglucan endotransglucosylase/hydrolase were suppressed. Two genes (PGSC0003DMG400004670 and PGSC0003DMG400021877) were down-regulated at 3 dpi were up-regulated at 7dpi. Moreover, 4 genes encoding for expansin and expansin-like were precisely up-regulated at 7 dpi ( Figure 8A and Table S5). Xyloglucan endotransglucosylase /hydrolase and expansin enzymes have cell wall loosening properties required for plant cell wall expansion during plant development. Therefore, we can hypothesis that M. Javanica induces cell wall modifying enzymes at 7 dpi for expansion and In addition, genes encoding for hydrolytic enzymes involved in pectin degradation such as polygalacturonase (PG), pectate lyases (PL), and pectin esterase (PE) were differentially expressed at 3 and 7 dpi by nematode infestation (Table S5). Further, genes encoding for pectin methylesterase inhibitor (PMEI) were largely down-regulated (10 out of 12 genes) ( Figure 8A and Table S5). More PMEI genes were suppressed at 7 dpi compared to 3 dpi. Repression of PMEI by nematode attack shows that the activity of PME was activated leading to the breakdown of pectin bonds, increasing the vulnerability of the cell wall to microbial pectic enzymes, and other degrading enzymes, culminating to a susceptible response. We also detected xyloglucan endotransglucosylase/hydrolase and expansin encoding genes, cell wall modifying enzymes under differential regulation following nematode attack. Of these, 19 out of 21 genes encoding for xyloglucan endotransglucosylase/hydrolase were suppressed. Two genes (PGSC0003DMG400004670 and PGSC0003DMG400021877) were down-regulated at 3 dpi were up-regulated at 7dpi. Moreover, 4 genes encoding for expansin and expansin-like were precisely up-regulated at 7 dpi ( Figure 8A and Table S5). Xyloglucan endotransglucosylase/hydrolase and expansin enzymes have cell wall loosening properties required for plant cell wall expansion during plant development. Therefore, we can hypothesis that M. Javanica induces cell wall modifying enzymes at 7 dpi for expansion and enlargement of GCS a crucial event for nematode development. Further, we identified the alternation of genes related to beta-glucanase, 18 out of 29 genes down-regulated during disease progression ( Figure 8A and Table S5). Plant cell wall modification, the deposition of callose a 1,3-β glucan cell wall polymer, is associated with cell wall thickening at the site of pathogen attack acting as a physical barrier derailing further pathogen invasion [57]. Therefore, suppression of beta-glucanase genes is probable to interfere with plant cell wall reinforcement to ward off further nematode invasion, hence successful M. javanica colonization.
With the increased demand for nutrients in nematode feeding cells, nematodes deploy specialized membrane transporters to control the flow of nutrients in and out of the NFS [58]. Several genes encoding for peptide transporters (31 genes), major intrinsic proteins (22 genes), amino acid transporters (32 genes), metabolite transporters (27 genes), and sugar transporters (16 genes) were differentially expressed. Overall, we found that 55.29% of transporter encoding genes in the DEGs were up-regulated following nematode infestation ( Figure 8B and Table S5). The activation of genes encoding amino acid transporters (17 genes) and sugar transporters (7 genes) indicates regulation of amino acid and carbohydrate metabolism, respectively. Similarly, the induction of sugar transporters increases soluble sugar content in RKN infested tomato plants, which is crucial for nematode development [59]. Furthermore, multidrug transporter-encoding genes were differentially expressed in our samples which encompasses ATP-binding cassette (ABC, 42 genes) and multidrug and toxin extrusion proteins (MATE, 17 genes) ( Figure 8B and Table S5). These results may suggest that M. javanica recruits some of these transporters to flush out toxic secondary metabolites or to disperse nematode effectors produced during the infestation process.

Regulation of Proteolysis Processes by Meloidogyne javanica
Gene Ontology enrichment analyses showed that genes involved in primary metabolism were overrepresented among the down-regulated genes (Table S5). Protein degradation was among the significantly altered processes involved in metabolism. In this study, genes coding for various classes of proteolytic enzymes were identified as differentially expressed include cysteine proteases, aspartate proteases, AAA-type, metacaspase, metalloprotease, subtilases, and serine proteases. The immense turnover of cellular proteins during nematode feeding can enhance a compatible nematode interaction due to regulation of plant immunity associated with proteolysis process. Further, we found several genes involved in ubiquitin-proteasome system (UPS), the main protein degradation pathway in a cell. Various classes of E3 ubiquitin ligases enzymes that play a central role during protein ubiquitination process were differentially expressed. This includes SKP1-CUL1-F-box-protein (SCF) and the U-box domain RING finger protein-encoding genes which were co-regulated, 102 down-and 105 genes up-regulated (Table S5). The fact that the UPS regulates degradation of many proteins in the cell affecting various cellular processes such as signal transduction, hormone signaling, and immune responses [60], makes it an attractive target for pathogen virulence factors including nematode effectors [61,62].
Collectively, this study uncovers the molecular networks regulated during compatible interaction between potato and M. javanica. Our results lay a foundation for functional studies in the future for genes highlighted to have a role in mediating susceptibility to M. javanica infestation to reveal their precise role in this interaction. In addition to providing further insights on plant-nematode interactions, further studies in the area development of target-specific control strategies against Meloidogyne species will originate. Supplementary Materials: The following are available online at http://www.mdpi.com/2076-2607/8/9/1443/s1, Table S1: List of primer used for RT-qPCR for validation of RNA-seq results; Table S2: Mapping summary of RNA-Seq reads from Solanum tuberosum samples using the reference genome (v4.03); Table S3: Differentially expressed genes induced by Meloidogyne javanica at 3 and 7 dpi; Table S4: Gene ontology enrichment analysis; Table S5: Differentially expressed biotic stress-related genes induced by Meloidogyne javanica at 3 and 7 dpi as defined by Mapman analysis; Table S6: Classification of differentially expressed transcription factors; Figure S1: The expression profile of genes experimentally validated using real-time qPCR in response to Meloidogyne javanica infestation; Figure S2

Conflicts of Interest:
The authors declare no conflict of interest.