Transcriptome Analyses of Pre-parasitic and Parasitic Meloidogyne Chitwoodi Race 1 to Identify Putative Effector Genes

The Meloidogyne chitwoodi isolate race 1 egg inoculum was initially provided by Dr. Charles Brown (USDA-ARS), and M. chitwoodi race 1 was maintained on the susceptible tomato Solanum lycopersicum cv. Rutgers under greenhouse conditions. To obtain nematode eggs from inoculated tomato plants, tomato roots were cut into small pieces and rinsed in a 0.6% sodium hypochlorite solution for 3 min. The solution was then passed through a set of sieves (pore size of 125, 45, 25 μm). Nematode eggs were collected on the 25 μm sieve and then further purified using sucrose floatation (Jacob and van Bezooijen, 1984). The eggs were suspended in 0.1% Plant Preservative Mixture (Plant Cell Technology) and incubated in a modified Baermann beaker (Figure S1) in darkness at room temperature. Hatched infective J2s were collected after 3 days of incubation.

Figure S1:

For inoculation of potato plants with M. chitwoodi, three-week-old tissue culture Russet Burbank plantlets were transferred to 500-mL cone-tainers filled with sand and grown in growth chambers (14-hr light/10-hr darkness, 23^oC) for 14 days. The Russet Burbank plants were then inoculated with 500 freshly hatched M. chitwoodi race 1 J2s. Root samples were harvested at 4-, 8- and 16-day-post-inoculation (dpi) and were stained with acid fuchsin solution (0.35% acid fuchsin, 25% acetic acid) to visualize nematodes within roots (Bybd et al., 1983) on a Nikkon Eclipse Ti inverted microscope at × 40 magnification.

For RNA-seq experiment, newly hatched M. chitwoodi race 1 J2s and galled potato root tissues at 8 dpi were collected, respectively. Three biological samples of M. chitwoodi race 1 J2s were also collected. Total nematode RNA was extracted using PureLink™ RNA Mini Kit (ThermoFisher). Following the inoculation of potato plants with 500 M. chitwoodi J2s (procedure described above), Russet Burbank roots were harvested at 8 dpi, and the infected roots were checked under a stereomicroscope (Zeiss zoom stereomicrope) to collect galled root tissues where nematodes resided. Galled root tissues from two individual potato plants were combined as one biological replicate, and three biological replicates were collected for 8 dpi samples. Total RNA was extracted from root samples using RNeasy Plant Mini Kit (QIAGEN, Germany). High quality total RNA (1-2 µg, 260/280 = 1.8-2.2, with RIN value of RNA of at least 8.5) was packed on dry ice and sent to Novogene Co. Ltd for library preparation and sequencing (30M reads per sample with PE150) using the Illumina platform.

The raw sequencing data were trimmed using Trimmomatic-v0.38 to remove adapters and low-quality sequences (Bolger et al., 2014). The quality of trimmed sequences was checked using FastQC-v0.11.7 (http://www.bioinformatics.babraham.ac.uk/projects/fastqc). The trimmed sequences were evaluated according to the summary output of FastQC, and threshold of per base sequence quality was set as 20 by default. The trimmed data was analyzed using Tophat/2.1.1 to map to the M. chitwoodi race 1 draft genome (Bali et al., 2021). Mapped read counts were calculated using Cuffdiff of the Cufflinksv2.2.1 package, and the gene read counts were converted to tab-delimited text file (Trapnell et al., 2012). The differentially expressed genes (DEGs, adjusted p-value < 0.05) between M. chitwoodi J2s and at 8 dpi (n = 3 biological samples for each time point) were analyzed using the Cuffdiff.

We predicted secreted proteins from the M. chitwoodi race 1 genes up-regulated at 8 dpi based on the criteria: presences of N-terminus signal peptide and absence of transmembrane domains. First, the coding DNA sequences (CDS) of M. chitwoodi race 1 genes were predicted using the TransDecoder v5.1.0 with minimal length of deduced proteins set as 30 amino acids (Haas et al., 2013). The deduced protein sequences of up-regulated genes were used to predict presence of N-terminus signal peptide for secretion using the software SignalP version 5.0 (Almagro Armenteros et al., 2019; Nielsen, 2017). Proteins predicted to have Signal Peptide (Sec/SPI) likelihood above 0.6 were retained for further analyses. These protein sequences were then subjected to TMHMM v. 2.0 analysis, which predicts membrane proteins (Krogh et al., 2001). Proteins were classified as secreted proteins if they either had no predicted transmembrane domains or only had predicted transmembrane domains in the secretion signal sequence.

For gene expression analysis, potato galled root tissues inoculated with 500 M. chitwoodi race 1 J2s were collected at 4-, 8- and 16-day-post-inoculation. Simultaneously, freshly hatched M. chitwoodi race 1 pre-parasitic J2s were also collected. Total RNA was extracted from galled root tissues and J2s separately and was treated with DNA-free™ DNA Removal Kit (ThermoFisher) to remove residual DNA. cDNA was synthesized using ProtoScript II First Strand cDNA Synthesis Kit (New England Biolabs, USA) with oligo-dT primer. qRT-PCR was performed using SsoAdvanced™ Universal SYBR^® Green Supermix on a CFX96 Real-Time PCR Detection System (Bio-Rad, USA). The qRT-PCR conditions were: 95^oC for 3 min, 40 cycles of 95^oC for 15 sec, 53^oC for 15 sec and 72^oC for 20 sec, and followed by a melting curve analysis from 65^oC to 95^oC with 0.5^oC increased incrementally at 5 sec increments. Expression levels of genes of interest in nematode were normalized to the expression of M. chitwoodi housekeeping gene Internal transcribed spacer 2 (ITS2) rRNA. The relative expression levels of each target gene at 4-, 8- or 16-dpi were calculated by comparing with those in M. chitwoodi race 1 J2s using the 2–ΔΔCt method (Livak and Schmittgen, 2001). The primer pairs used for RT-qPCR are summarized in Table 1. Each qRT-PCR experiment consisted of two technical replicates per time-course of each gene, and the average of three biological replicates. Statistical tests were performed using GraphPad Prism version 8.00 for Windows (GraphPad Software, La Jolla California USA, www.graphpad.com).

Root-knot nematode effector genes critical for parasitism are typically up-regulated during early infection stages, when the nematodes enter the root, evade plant defenses, and form their giant cells. To identify M. chitwoodi genes expressed during early infection stages, we first established infection timepoints in our experimental system. We monitored nematode infections at 4-, 8- and 16-days post-inoculation (dpi) by staining the roots with acid fuchsin, allowing us to observe nematodes within the roots. At 4 dpi, J2 could be seen within the roots, indicating that root penetration had occurred. At 8 dpi, most nematodes within the roots were swollen, but still in the J2 stage. The swollen bodies indicated that the nematodes had begun to feed but had not yet molted to the next life stage. At 16 dpi, galls on the roots were clearly visible and the nematodes were in the globose J4 or young female stage (Fig. 1). Based on our observations and experimental setup, at 4 dpi, the juveniles were still migrating through the roots, and by 8 dpi, the nematodes had settled and started to feed. Because we wanted to capture the J2 stage when initial feeding and giant cell formation occurs, we chose 8 dpi as the timepoint for the RNA-seq experiment.

Figure 1:

Our next goal was to compare gene expression profiles between pre-parasitic and parasitic J2 to determine which nematode genes are differentially regulated during parasitism. We collected freshly hatched M. chitwoodi Race 1 J2 as the pre-parasitic J2 samples. To obtain the parasitic nematodes, Russet Burbank potato plants were inoculated with M. chitwoodi Race 1 J2 and galled root tissues were collected at 8 dpi. We collected whole galls to avoid isolating nematodes from the roots, which meant we had a mixture of plant and nematode tissue in the sample. Six cDNA libraries were constructed – three libraries of pre-parasitic J2 and three libraries of galls isolated at 8 dpi. The libraries were sequenced using an Illumina platform for paired end 150-bp reads, generating between 43 and 55 million raw read pairs for each library. The RNA-seq raw reads were pre-processed to remove low-quality bases and adapter sequences, resulting in between 40 and 52 million clean read pairs (Table 2). Interestingly, the GC content of individual libraries were different among samples: an average of 34% GC content for the M. chitwoodi J2 samples and 41% of GC content for the 8 dpi samples. The low GC content of the M. chitwoodi J2 transcripts corresponds to the low GC content of the M. chitwoodi and other root-knot nematode genomes (Koutsovoulos et al., 2020; Paganini et al., 2012). The higher GC content of the 8 dpi samples may be due to the mixture of plant and nematode reads at 8 dpi.

The differentially expressed genes (DEGs) of M. chitwoodi between J2 and 8 dpi were analyzed using the TopHat and Cufflinks/Cuffdiff pipeline (Trapnell et al., 2012). First, the RNA-seq read pairs were mapped to the M. chitwoodi race 1 genome using TopHat2. Mapping the reads from the J2 libraries generated 21 to 25 million mapped read pairs. Meanwhile, mapping reads from the 8dpi gall library resulted in 1.2 to 1.5 million mapped read pairs. By mapping to the M. chitwoodi race 1 genome, we removed the plant sequences from our dual RNA-seq experiment (Table 2). The reduction in the read pairs suggest that the majority of the reads from the dual sequencing were from potato. Next, we used Cufflinks to merge the information from all six libraries to assemble a M. chitwoodi transcriptome. DEGs between pre-parasitic and parasitic (8 dpi) J2 were determined by Cuffdiff. In the parasitic nematodes, we found 693 M. chitwoodi genes up-regulated (fold-change ≥ 4 and adjusted P < 0.05) and 1727 M. chitwoodi genes down-regulated (at least a 4-fold decrease and adjusted P < 0.05) (Supplementary Table 1, https://figshare.com/articles/dataset/Zhang_and_Gleason/16746286). To limit the number of genes in our subsequent analyses, we focused on genes encoding proteins of 30 or more amino acids. This left us with 679 genes up-regulated and 1570 genes down-regulated. A BLASTx-fast search of the NCBI non-redundant (nr) database found that of the 679 upregulated genes, 341 (50%) did not have any Blast hits (E value ≤ 1.0E-3) (Supplementary Table 2, https://figshare.com/articles/dataset/Zhang_and_Gleason/16746286). Second level GO term analysis of the up-regulated genes indicate that sequences related to cellular, metabolic, and developmental processes were enriched during nematode parasitism (Fig. 2A). Of the 1570 down regulated genes, 408 (26%) did not have any Blast hits (E value ≤ 1.0E-3) (Supplementary Table 3, https://figshare.com/articles/dataset/Zhang_and_Gleason/16746286). Second level GO term analysis of the down-regulated genes indicate that sequences related to cellular, signaling, and response to stimulus were enriched in pre-parasitic nematodes (Fig. 2B). When comparing genes expressed in the pre-parasitic J2 versus genes expressed in the parasitic nematode, we comprised a list of the top 20 most differentially expressed genes and their annotations from the BLASTx homology search (Table 3).

Figure 2:

We hypothesized that nematode genes induced at early infection stages may be important for successful infection of host plants, therefore, we focused on the 679 genes that were significantly up-regulated at 8 dpi (Fig. 3, Supplementary Table 2, https://figshare.com/articles/dataset/Zhang_and_Gleason/16746286). We performed a bioinformatic analysis on the 679 M. chitwoodi genes up-regulated at 8 dpi to predict if the genes encode secreted proteins. The prediction of secretion was based on two criteria: the presence of a signal peptide at the N-terminus of the peptide sequence and the lack of any transmembrane domains. The SignalP 5 program predicted that 133 of the proteins contained a signal peptide at the N-terminus. Among these 133 sequences, 6 were excluded from further analysis because the TMHMM2 predictions indicated they contained one or more transmembrane domains within the protein sequence, minus the signal peptide sequence. Therefore, in silico analyses identified 127 M. chitwoodi genes (18.7% of the 679 up-regulated genes) encoding predicted secreted proteins (Supplemental Table 4, https://figshare.com/articles/dataset/Zhang_and_Gleason/16746286).

Figure 3:

Using the sequence information for the 127 identified putative secreted proteins, we performed a Blastx search of the NCBI nr database. We found that 62.9% (80/127) did not have any BLAST hits (E value ≤ 1.0E-3) (Supplemental Table 4, https://figshare.com/articles/dataset/Zhang_and_Gleason/16746286). Twelve of the genes (12/127) were annotated as “putative esophageal gland cell secretory protein” (Huang et al., 2003). Three known effectors were found in the annotations of the 127 M. chitwoodi genes, including VAP1, which was originally discovered in cyst nematodes as an effector that suppresses basal plant defense responses (Lozano-Torres et al., 2014). A gene encoding a chorismate mutase domain-containing protein was also identified in our in-silico M. chitwoodi secretome analysis. Root-knot nematode chorismate mutases may be altering salicylic acid levels within the host cells to suppress plant immune responses (Huang et al., 2005; Lambert et al., 1999; Wang et al., 2018). The third known effector found in our list of M. chitwoodi genes is homologous to the candidate secreted effector Minc11888, which was previously identified in M. incognita and is highly expressed in the parasitic juvenile stages (Nguyen et al., 2018).

Similarly, to assess protein homology, a BLASTp analysis of the 127 proteins using the nr database was performed. The analysis revealed that 66% (84/127) had significant hits (60% query coverage and E value ≤ 1.0E-3), to plant parasitic nematodes (Supplemental Table 5, https://figshare.com/articles/dataset/Zhang_and_Gleason/16746286). Interestingly, 53.5% (45/84) of those genes were hypothetical or unnamed proteins from the rice root-knot nematode (M. graminicola), which is a close relative to M. chitwoodi in the Meloidogynidae (Bert et al., 2008; Phan et al., 2020). We also found that 35/84 genes in our BLASTp analysis had top hits to genes encoding unnamed proteins in M. enterolobii, perhaps owing to the fact that M. enterolobii has one of the most complete genomes available for root-knot nematodes at this time (Koutsovoulos et al., 2020).

To validate our list of up-regulated genes, we selected 11 genes from the list of 127 that were highly up-regulated at 8 dpi (FC > 10 and FPKM ≥ 2, Table 4). We measured the expression of these 11 genes by quantitative reverse transcription polymerase chain reaction (RT-qPCR) in different nematode life stages: pre-parasitic J2 and parasitic nematodes in potato roots at 4, 8, and 16 dpi. Expression of 6 genes (Mc03439, Mc03992, Mc04677, Mc04921, Mc08895 and Mc10400) peaked at the early infection stages (4 or 8 dpi) and decreased at a later life stage (16 dpi) (Fig. 4). For the other genes, the expression levels were either similar across all infection stages (Mc01565, Mc06869, Mc08299) or highly expressed at the later infection stage, 16 dpi (Mc07118 and Mc10102) (Fig. 4). All 11 genes were confirmed to be expressed at higher levels at 8 dpi by RT-qPCR compared to the J2 stage, confirming our transcriptome data (Fig. 4). These results highlight the fidelity of dual RNA-seq analyses to identify differentially expressing genes.

Figure 4: