Meloidogyne incognita nematodes were propagated on 3-week-old pepper (Capsicum annuum, Qiemen) plants in an artificial environment from a single female egg mass after isolated from the Sijiqing farm (Beijing, China). Three lines (Avir-1, Avir-2, and Avir-3) of the nematodes that each from a female egg mass were continually reproduced for 30 generations to form stable populations and then used for experiments. Infective, second stage juveniles were collected by hatching the handpicked egg masses at 26°C.

Nicotiana benthamiana was used for plant hypersensitive response (HR) assay and subcellular localization observation. Arabidopsis thaliana (Col-0 ecotype) was used for transgene, plant-nematode interaction and parasitism analyses. Both plants were cultivated in a growth chamber at 23°C.

Nematode genomes and transcripts used to predict effector proteins containing nuclear localization signals (NLSs) and for sequence homologue analysis in were downloaded from the National Center for Biotechnology Information (NCBI) (Supplementary Table S1). Signal peptide predictions were performed using SignalP 4.1¹, predisi², and phobius³. Subcellular location of proteins was predicted by TargetP⁴, and the proteins located in mitochondria were removed. The prediction of transmembrane domain and NLS were performed using TMHMM⁵ and NucPred (Brameier et al., 2007) based on a sequence score ≥ 0.5. Domains were predicted in local PuTTY based on HMMER 3.1b1 (Wong et al., 2015) and the active conserved domains were blasted to the Pfam database (Finn et al., 2014). A phylogenetic analysis of homologues was carried out using the maximum likelihood method with MEGA 5 (Tamura et al., 2011).

The suppression of cell death elicited by B. glumae (from a gift of Dr. Sun, Chinese Agricultural University) on N. benthamiana leaves was performed as described in previous study (Sharma et al., 2013). Candidate effectors and GFP coding sequences were cloned into the pEDV vector to generate pEDV::effector and pEDV::GFP constructs, respectively. All constructs were introduced into B. glumae by electroporation, and the transformed bacteria were then suspended in 0.9% NaCl at 600 nm (OD₆₀₀ = 0.4). Burkholderia glumae cells carrying pEDV::effector were infiltrated into the right side of N. benthamiana leaves using a 1 ml syringe. As controls, B. glumae cells carrying the pEDV::GFP were infiltrated into the left side of N. benthamiana leaves. Three leaves on 5 plants were used for each biological replicate, and at least three independent experiments were performed. Plants were kept at 25°C. Cell death was observed at 4 days post-inoculation. All primers were synthesized by Sangon and were listed in Supplementary Table S2.

For subcellular localization assays, MiISE2^Δsp, MiISE10^Δsp, MiISE6^Δsp, MiISE6^{Δsp_Δ109_157}, MiISE6^{Δsp_Δ33_108} and MiISE6^{Δsp_MΔ109_118} were amplified by PCR from the M. incognita cDNA, using I-5^TM 2 × High-Fidelity Master Mix (MCLAB, United States). PCR products were digested with Nco I and Spe I and inserted into pCAMBIA1302 digested with the same enzymes. The resulting plasmids harbored the effector gene fused with the GFP gene. The constructs were transformed into Agrobacterium tumefaciens GV3101. After inoculation in LB followed by horizontal shaking for approximately 30 h, A. tumefaciens cultures were diluted to an OD₆₀₀ of 0.8 using infiltration medium (10 mM MES, 10 mM MgCl₂ and 0.1 mM AS). Leaves at the four-leaf stage of N. benthamiana were infiltrated with bacterial inocula using a 1 ml syringe with no needle. After 24–48 h, the infiltrated zones of the leaf were observed under a confocal microscope (Zeiss LSM700, Germany).

RNA samples were extracted from 300 M. incognita at different life stages (pre-J2, parasitic (par)-J2, par-J3, par-J4 and female) for developmental expression analysis using RNAprep pure Tissue Kit (Tiangen, China). The 18S gene (Genbank accession no. U81578) was used as an internal control.

In order to confirm the results of RNA-seq analysis, Arabidopsis plants grow under same condition as those used for the RNA-seq analysis were collected for q-PCR assays. RNA samples were extracted from 100 mg frozen plant tissue using RNAprep pure Plant Kit (Tiangen, China). The Actin2 gene (AT3G18780) was used as an internal control.

Then, 1 μg RNA was used for cDNA synthesis using PrimeScript^TM RT reagent Kit with gDNA Eraser (Perfect Real Time) (Takara, Japan) following the manufacturer’s protocol. Each 25 μl reaction mixture was prepared using SYBR Premix Ex Taq^TM||(TaKaRa, Japan) on a BIO-RAD CFX96 (BIO-RAD, United States) following the manufacturer’s protocol. Cycling conditions for quantitative PCR were: 95°C for 5 min and 40 cycles of 95°C for 30 s and 60°C for 30 s. Gene expression changes were calculated using the 2^-ΔΔCT method. At least three independent experiments were performed, with four technical replicates for each reaction.

Parasitic stage J2s of M. incognita were processed for in situ hybridization to confirm the expression site of MiISE6 as previously described (de Boer et al., 1998). Digoxigenin-labeled sense and antisense cDNA probes were synthesized using the primers DIG-MiISE6-F and DIG-MiISE6-R. Hybridization signals within the nematode were detected with alkaline phosphatase-conjugated anti-digoxigenin antibody and substrate, and observed with an Olympus IX53 microscope.

Fluorescence in situ hybridization was performed using a 28 bp specific cDNA probe (CTTTCTTCAATTCCTT TCATCAACATCC) whose 5′ end was labeled with fluorescein isothiocyanate (FITC, Sangon Biotech) as previously described with some modification (Moter and Gobel, 2000), and the hybridization signals were observed under a confocal microscope (Zeiss LSM700, Germany).

We used the yeast secretion system based on the yeast signal sequence trap vector pSUC2T7M13ORI (Jacobs et al., 1997) to validate the function of MiISE6 signal peptide. PCR products were digested with EcoRI and XhoI and inserted into pSUC2T7M13ORI to generate pSUC2::MiISE6. The primer pairs used for PCR amplification were pSUC2_MiISE6_F and pSUC2_ MiISE6_R. The constructs were transformed into the invertase negative yeast strain YTK12 using the Frozen-EZ Yeast Transformation II kit (Zymo Research, United States). After transformation, yeast was plated on CMD-W (minus Trp) plates for 3 days. Single colonies of YTK12 were replicated onto raffinose-containing YPRAA plates (with raffinose instead of sucrose). Invertase enzymatic activity was detected with TTC as described by Oh et al. (2009).

For MiISE6 overexpression in Arabidopsis, the coding sequence of the MiISE6 gene without signal peptide was cloned into the vector pEGAD, and for host-derived gene silencing, the coding sequence of MiISE6 was inserted into the pFGC5941 vector in both sense and antisense orientations. There was an CHSA intron between sense and antisense fragments. Both pEGAD vector and pFGC5941 vector habor the BASTA selectable marker gene. The transformation of Arabidopsis was performed as previously described (Martínez-Zapater and Salinas, 1998). Transformed T1 seeds were sown on soil and were screened by spraying the herbicide BASTA 2 weeks later to select for transgenic plants. The spraying was repeated three times. Homozygous T3 seeds were used in this research. The expression level of MiISE6 in different transgenic line were confirmed by RT-PCR, RNAi transgenic lines was confirmed by checking the CHSA intron fragment using primer CHSA-F/CHSA-R.

Four-week-old Arabidopsis were inoculated using 400 J2 of M. incognita per plant. 42 days after inoculation, the number of females was counted.

At least 15 plants were used for each treatment, and three independent experiments were performed. Statistically significant differences between treatments were determined by independent samples t-test (P < 0.05) using SPSS.

For M. incognita transcriptome: 50 mg of freshly collected nematodes were used for RNA extraction using TRIzol reagent (Invitrogen) following the manufacturer’s instructions. At least 20 mg of total RNA was sent to the Shenzhen Genomics Institute for Solexa sequencing (commercial service, Shenzhen China), and three independent replicates (Avir-1, Avir-2, and Avir-3) were sequenced. The full genome sequences and gene annotation of M. incognita were downloaded from Wormbase as a reference genome (Abad et al., 2008). The RNA-seq reads were mapped to the reference genome using Tophat and subsequently assembled using cufflinks (Trapnell et al., 2010). The expression abundance of each gene was estimated with fragments per kilobase of transcript per million fragments mapped (FPKM) (Anders et al., 2015).

For A. thaliana transcriptome: The 14-day transgenic Arabidopsis and WT Arabidopsis grown on Murashige and Skoog (MS) medium were used for RNA extraction. RNA samples were extracted from 100 mg frozen plant tissue using RNAprep pure Plant Kit (Tiangen, China). Library construction and RNA-Seq were performed at BerryGenomics (Beijing, China), and four independent replicates were sequenced. The cDNA libraries with average insert sizes of 280 bp were sequenced by Illumina technology. The adaptor sequences and low-quality reads were filtered from the raw data. Tophat2 tools were used to map mRNA reads to the Arabidopsis genome (Langdon, 2015), and htseq-count was then used to calculate the expected FPKM mapped (Anders et al., 2015), and the edgeR package in R was applied to identify DEGs of 2 groups of samples with a threshold criterion of (FDR < 0.01) and |log₂FC|≥ 1. We used AgriGO⁶ and KOBAS (Xie et al., 2011) to perform GO functional enrichment analysis and KEGG enrichment analysis with default parameters.

In total, 8.22 billion clean base pairs were generated from three independent replicates of M. incognita wild type population, with Illumina/Solexa sequencing reads of ∼70x of M. incognita genome. Of these, 13398 genes (95.90% of total) were shared by the three replicates (Supplementary Figure S1); detailed information is shown in Supplementary Datasheet S1. By a series bioinformatic analysis, such as searching out the sequences with signal peptides and discarding the sequences with putative transmembrane-spanning regions, analyzing subcellular location and removing proteins located in mitochondria, NLS motif predicting, totally 110 potential NLS effectors are predicted from M. incognita, and they are taken as candidate effectors with putative immune function for further analysis. Detailed information is shown in Supplementary Datasheets S1, S2.

Previous studies have shown that PPN effectors can suppress the plant hypersensitive response (HR) to facilitate parasitism (Niu et al., 2016; Chen et al., 2017). Here we used a pEDV-based type III secretion system (T3SS) in which B. glumae is used to deliver effectors to plant cells to verify the function of 20 randomly selected novel effector candidates in plant immune suppression. Burkholderia glumae is a bacterial rice pathogen that can cause a fast, localized cell death in non-host N. benthamiana. We also cloned a green fluorescent protein gene (GFP) into the vector as a negative control, and HR symptoms were monitored within 3 days after agroinfiltration. Screening results 3 out of 20 candidate effectors had strong HR suppression ability, including MiISE2 (Minc04520), MiISE10 (Minc08615), and MiISE6 (Minc06775) (Figure 1).

To investigate the subcellular localization of the three candidate effectors in plant cells, the green fluorescent protein (GFP) gene was fused with the effector (without signal peptide and stop codon) at its C terminus with the 35 S promoter, and the effector proteins were then transiently expressed in N. benthamiana. Confocal microscopy revealed distinct subcellular localization of the 3 candidate effectors. MiISE2 shows nuclear-cytoplasmic localization; MiISE10 appears to localize in the cytoplasm; and MiISE6 is enriched in the nucleus (Figure 2A). The subcellular localization of MiISE2 and MiISE6 are consistent with the predictions by NucPred. Moreover, we used q-PCR to quantify the expression levels of three effectors during pre-parasitic J2 (pre-J2), parasitic J2 (par-J2), J3 and Female stages. The candidate effectors MiISE2 and MiISE10 have a high expression level in late parasitism stages (J3 and Female), but they exhibit a very low expression level in early parasitism stages. Conversely, the candidate effector MiISE6 is upregulated in early sedentary stages, reaching a peak in the early par-J2 at 1 dpi, and then undergoing a dramatic reduction in the J3 parasitism stage (Figure 2B). The expression pattern of MiISE6 is similar to other published nematode effectors which play central roles in nematode parasitism (Chronis et al., 2013; Lin et al., 2016; Chen et al., 2017), so we focus on MiISE6 and speculate it may play an important role during the initial establishment of the plant-nematode interaction.

The M. incognita MiISE6 cDNA encodes a 158 amino acid protein with a predicted N-terminal signal peptide (SP, residues 1–32) for secretion; an OGFr_N region (residues 35–62); and two NLS motifs (residues 89–98, 109–118) (Supplementary Figure S2). A search of homologues of MiISE6 from 16 nematodes species (including 4 free living nematodes, 6 animal or insect parasitic nematodes, and 6 plant parasitic nematodes) (Supplementary Table S1) was performed with BlastP (E ≤ 1e-10), and we obtained 20 homologues (Supplementary Table S3). These homologs were not only distributed in plant parasitic nematodes, but also in animal or insect parasitic nematodes and free-living nematodes. Based on the amino acid sequences of MiISE6, a phylogenetic tree was constructed in MEGA5 using the maximum likelihood method with 1000 bootstrap replicates. MiISE6 was closest to the homologue from M. hapla, and they clustered with MiISE6 homologues from Globodera pallida (Figure 3).

Sequence analysis revealed that MiISE6 harbors two putative NLSs including ⁸⁹GMQKRKRALD⁹⁸ and ¹⁰⁹FDRKRRALDM¹¹⁸. To test which NLS play function in nuclear localization, we firstly generate two constructs MiISE6^{Δsp_Δ109_157}:GFP and MiISE6^{Δsp_Δ33_108}:GFP, which contain NLS^89-98 and NLS^109-118, respectively. As shown in Figure 4, the MiISE6^{Δsp_Δ109_157}:GFP localized to the cytoplasm (Figure 4A), whereas MiISE6^{Δsp_Δ33_108}:GFP localized to the nucleus in N. benthamiana leaves (Figure 4B). The results indicate that the predicted NLS^89-98 motif may not responsible for the nuclear localization of MiISE6. To determine whether NLS^109-118 motif play role in nuclear localization, we used site-directed mutagenesis to generate MiISE6^{Δsp_MΔ109_118}:GFP, and the fluorescence accumulation was observed in cytoplasm, and no fluorescence was observed in the nucleus (Figure 4C). These results indicate that the NLS^109-118 lead to the nuclear localization of MiISE6.

The tissue localization of MiISE6 transcripts in the nematode was verified using in situ hybridization. With digoxigenin-labeled antisense cDNA probes of MiISE6, there is a strong signal within the subventral gland cells of the pre-parasitic J2s. No signals are detected using the sense cDNA probes. Consistent with the result of DIG-labeled in situ mRNA hybridization, strong fluorescent signals are observed in the subventral esophageal gland by fluorescence in situ hybridization, and there is no tissue specificity when dyed with DAPI (Figure 5).

We used a yeast system requiring invertase secretion to grow on raffinose media to validate the function of the signal peptide of MiISE6. The predicted signal peptide sequence and the subsequent 24 amino acids of MiISE6 was firstly fused with the SUC2 minus the signal peptide and then inserted in vector pSUC2T7M13ORI (Jacobs et al., 1997). The final construct enables the successful secretion of invertase, and then the invertase mutant yeast strain YTK12 can hydrolyze raffinose and grow on YPRAA media. The invertase secretion was confirmed with 2,3,5-triphenyltetrazolium chloride (TTC), which reacts with monosaccharides to form the insoluble (red colored) triphenylformazan. Compared with the negative control, the fusion of signal peptide of Avrblb2 and MiISE6 to the SUC2 enabled the secretion of invertase, which lead to the yeast cells growing on YPRAA media (Figure 6). These results revealed that the MiISE6 protein habor a functional signal peptide.

To determine whether the expression of MiISE6 in Arabidopsis can change plant phenotype and nematode susceptibility, we generated transgenic Arabidopsis lines expressing ^ΔspMiISE6 driven by the CaMV35S promoter. Three independent transgenic lines (T_1, T_2, T_8) were selected for phenotype observation and the nematode inoculation. We firstly detected the expression of transcripts in independent transgenic lines by RT-PCR. The photo (Figure 7A and Supplementary Figure S3) shows that MiISE6 gene is really expressed in Arabidopsis plants. Then, we counted the number of root-knots in each transgenic plant, and compared the difference between each transgenic line and the wild-type plants by statistical significance test. The results show that there are more root-knots in each transgenic line plants than those in wild-type plants (mean 23–25 vs. 18–19) (Supplementary Figure S4). The difference is significant at a statistic level (all of the three tests with P < 0.001, Supplementary Figure S5). It indicates that all three transgenic lines are more susceptible to M. incognita than control, and there is an increase in females with an average rate of 25.15, 35.26, and 38.80%, respectively (Figure 7B). To further assess the role of MiISE6 in nematode parasitism, we used a host-derived gene silencing method to knock down the expression of MiISE6 in M. incognita during the nematode infection. Three RNAi lines (RNAi1, RNAi5, RNAi11), with at least 15 plants for each treatment, and three independent experiments. By RT-qPCR assays, it is shown that the expression of MiISE6 in M. incognita collected from inoculated on RNAi plants at 3 dpi is significant reduced compared with those from the wild-type (WT) plants (Supplementary Figure S6 and Figure 8A). We also counted the number of root-knots in each RNAi plant and wild-type plant. Statistical analysis showed that there is an obvious reduction on RNAi plants (mean 13-15 knots) (Supplementary Figure S7), compared to those on wild type plants (mean 21). There was a 29.16, 35.26, and 30.86% reduction of females in three RNAi plants, respectively (Figure 8B). The results of t-tests show that the difference is significant (P < 0.001, Supplementary Figure S8). These data indicate that MiISE6 is an important effector in mediating parasitic success of M. incognita.

The difference in expression level of DEGs between WT and MiISE6 transgenic Arabidopsis plants was calculated according to the FPKM method. Under the criteria of |log₂FC|≥ 1 and FDR ≤ 0.01, a total of 852 genes were identified in MiISE6-overexpressing samples compared with WT samples, of which 772 were upregulated and 80 were downregulated. Detailed information is provided in Supplementary Datasheets S3, S4.

To functionally classify the DEGs, gene ontology (GO) terms were assigned to each DEG using AgriGO based on the TAIR 10 database. DEGs are categorized into three groups: Biological Process (BP), Cellular Component (CC), and Molecular Function (MF). The majority of the DEGs are assigned to the BP category, such as cellular process, metabolic process, and response to stimulus. In addition, we find that DEGs are significantly involved in binding, catalytic activity, and transcription regulator activity in the MF category. The results also show that in the CC category, the majority of the DEGs are significantly related with cell part and organelle (Figure 9). The expression of MiISE6 in Arabidopsis leads to changes in multiple plant immune responses, including hormones signaling pathway, ubiquitination proteasome pathway, and calcium signaling pathway (Supplementary Table S4). In “regulation of defense response (GO:0031347)”, most enriched genes work as negative regulator of plant defense responses, such as WRKY40, WRKY18, and NUDT7(Yang et al., 2006; Ge et al., 2007; Pandey et al., 2010). In “protein ubiquitination (GO:0016567),” most enriched genes are plant U-BOX E3 ubiquitin ligases (PUBs). The U-Box domain mediates the binding to the conjugating enzyme (E2) during ubiquitination process. Previous study showed that there was a reduced bacteria infection in pub22/pub24 double mutant and pub22/pub23/pub24 triple mutant compared to control. Although the molecular mechanism is not well clarified, previous studies indicated that those genes may act as negative regulator of PAMP-triggered immunity by altering the activity of 26S proteasome (Stegmann et al., 2012; Cho et al., 2015). In “response to hormone stimulus,” genes involved in modulating multiple hormone signaling pathway, including JA signaling (GO:0009867), SA signaling (GO:0009751) and auxin signaling (GO:0009733). Moreover, we find that a number of genes are enriched in “cell wall modification (GO:0042545),” such as xyloglucan endotransglucosylase-hydrolase (XTHs) and expansin-like gene (EXLA). Those genes involved in the reorganization or loosening of cell walls (Cosgrove, 2000; Olsen and Krause, 2017). Recent studies have suggested that diverse pathogen acquired and exploited these plant derived enzymes for successful parasitism (Balestrini et al., 2005; Nikolaidis et al., 2014; Olsen and Krause, 2017). Our result shows that most of those genes exhibit a significant higher expression level in the MiISE6 transgenic plants compared to control. Based on the present data, we draw a model illustrating the interplay between potential DEGs involved in host defense responses (Figure 10).

Moreover, the results of KEGG enrichment analysis indicate that a total of 228 annotated DEGs are assigned into 25 pathways. As shown in Supplementary Table S5, these pathways are mainly classified into 8 categories (P ≤ 0.05). The most significantly enriched KEGG pathway is “Plant–pathogen interaction,” and most genes enriched in this pathway are calmodulin-related proteins or calcium dependent kinases, such as CML16, CPK29, CPK28. The expression of CPK28 has a 5.3-fold increase (log₂FC = 2.31) in the MiISE6 transgenic plants compared to the WT plants. Previous studies have proven that CPK28 can interact with phosphorylate BIK, and contribute to BIK1 turnover (Monaghan et al., 2014, 2015); the overexpression of CPK28 in Arabidopsis inhibits PTI signaling and immunity (Matschi et al., 2015). Moreover, the NADP oxidase RBOHD is also enriched in the “Plant–pathogen interaction” pathway and presents a higher expression in the transgenic plants, and the expression of SA-responsive genes EDS1, PR1, PR2, and PR5 did not show significant differences in the transgenic plants compared to control. Additionally, in the KEGG pathway “Plant hormone signal transduction,” most enriched genes are JAZ genes, which work as repressors of jasmonate signaling (Chini et al., 2007), and they all showed significantly higher expression in the transgenic plants compared to control. Conversely, the JA-responsive genes PR4 and PDF1.2 showed a significant decrease in the transgenic plants. The expression patterns of these genes were verified by q-PCR, and the results are similar to the transcriptome analysis (Figure 11).