Soybean cultivar “Williams 82” and “DN50” obtained from the Chinese Academy of Agricultural Sciences were used as WT for all soybean assays. Gmrmi mutants (including Gmrmi1 and Gmrmi2) were mutagenized from the soybean treated with 5% EMS (Wang et al., 2020). All seedlings of soybean were grown in the field from May to October for separate group construction and gene mapping, or in a climate chamber, whose temperature was maintained at 25°C and humidity at 50% with a 14-h photoperiod for nematode infection assays.

We performed the detailed assay of mapping the GmLMM1 gene as the method described by Wang (Wang et al., 2020).

For nematode infection assays toward soybean, 1-mL aliquot of a suspension containing 500 M. incognita pre-J2s was inoculated on the roots of 14-day-old wild-type and mutation soybean at 25°C. Subsequently, the number of nematodes inoculated was recorded at different stages. Fuchsin-stained roots were observed with an Olympus BX53 microscope (Olympus, Japan). For each experiment, 10–20 plants were used per treatment, and each experiment was repeated independently at least three times.

RNA-sequencing was performed by Novogene (Beijing, China). Soybean cultivar “Williams 82” and Gmlmm1-1 were grown for 14 days and then the roots of which were inoculated with 1-mL aliquot of a suspension containing 500 M. incognita pre-J2s. Seven days post inoculation, the roots of soybean were collected and a total amount of 1 μg RNA per sample was used as input material for the RNA-sequencing (RNA-seq) library preparations. Briefly, mRNA was purified from total RNA using poly-T oligo-attached magnetic beads. RNA degradation and contamination were monitored on 1% agarose gels. RNA purity was checked using the NanoPhotometer spectrophotometer (IMPLEN, CA, USA). RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA). Sequencing libraries were generated using NEB Next UltraTM RNA Library Prep Kit for Illumina (NEB, USA) following the manufacturer's recommendations, and index codes were added to attribute sequences to each sample. The clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumia) according to the manufacturer's instructions. After cluster generation, the library preparations were sequenced on an Illumina NovaSeq platform, and 150-bp paired-end reads were generated. Raw data (raw reads) of fastq format were first processed through in-house perl scripts. In this step, clean data (clean reads) were obtained by removing reads containing adapter, reads containing ploy-N, and low-quality reads from raw data. Reference genome and gene model annotation files were downloaded from genome website directly. Index of the reference genome was built using Hisat2 v2.0.5, and paired-end clean reads were aligned to the reference genome using Hisat2 v2.0.5. The fragments per kilobase of exon model per million reads mapped (FPKM) value of each gene was calculated using cuff links (Trapnell et al., 2012), and the read counts of each gene were obtained using htseq-count. DEGs were identified using the DESeq estimating function and nbinom test. Notably, p-value <0.05 and fold changes >2 or <0.5 were used as the criteria for identifying significant DEGs. A hierarchical clustering analysis of the DEGs was performed to explore gene expression patterns. The GO enrichment of the DEGs was estimated using the DESeq R package based on the hypergeometric distribution.

The DNA sequences of the bioactive form of the AtRALF1 peptide (amino acids 72–120) and AtRALF23 peptide (amino acids 89–138) were PCR-amplified from Arabidopsis cDNA, and MiRALF1/3 were amplified from M. incognita cDNAs via PCR (Zhang et al., 2020a). All the RALF/RALF-like genes were constructed into a modified expression vector, pET-28a, with a 6 × His-MBP tag. Then, the constructed vector was transformed into the E. coli strain BL21 Gold (DE3) and cultured in LB medium overnight at 37°C, before transferred to fresh media at a 1:100 dilution. When the OD₆₀₀ value of bacteria solution nearly reached 0.5, a final concentration of 0.5 mM isopropylthio-β-D-1-galactopyranoside (IPTG, Sangon Biotech, Shanghai, China) was used and the cultures were incubated at 28°C for 4 h to induce expression. Subsequently, RALF/RALF-like proteins with the 6 × His-MBP tag were expressed and collected by centrifugation for purification. One volume of the collected cells was resuspended in seven volumes of lysis buffer [50 mM Tris-HCl (pH 8.0), 300 mM NaCl, and 10 mM imidazole] supplemented with 1 mM of the protease inhibitor PMSF (Selleck Chemicals, Houston, USA) and suspended for sonication. The extract was centrifuged at 10,000 × g for 15 min at 4°C, the soluble fraction was obtained, and the desired protein was bound to His-beads (Ni Sepharose HP, GE Healthcare) overnight. Subsequently, the beads were washed sequentially with seven volumes of wash buffer 1 (50 mM Tris-HCl, 300 mM NaCl, and 20 mM imidazole, pH 8.0), wash buffer 2 (50 mM Tris-HCl, 300 mM NaCl, and 50 mM imidazole, pH 8.0), and wash buffer 3 (50 mM Tris-HCl, 300 mM NaCl, and 100 mM imidazole, pH 8.0) for 10 min each prior to eluting the 6 × His-MBP-RALF protein with elution buffer (50 mM Tris-HCl, 300 mM NaCl, and 250 mM imidazole, pH 8.0) for 3 h at 4°C.

The GmLMM1 extracellular domain (from residues 27 to 452, GmLMM1^ECD) was amplified and inserted into the vector pGEX-4T-1 using the primers GmLMM1-ECD-F and GmLMM1-ECD-R (Supplementary Table S1) to produce a GST-GmLMM1^ECD fusion protein in the E. coli strain BL21 Gold (DE3). Cells containing GST-GmLMM1^ECD were collected by centrifugation. One volume of pelleted cells was resuspended in seven volumes of GST lysis buffer [50 mM Tris-HCI (pH 8.0), 100 mM NaCl, 10 mM MgCl₂, 1 mM dithiothreitol (DTT), and 10% glycerol] supplemented with a 1 × solution of the protease inhibitor PMSF and suspended for sonication. After centrifugation (4000 rpm, 4°C, 15 min), the supernatant was absorbed onto prewashed GST agarose (Glutathione Sepharose 4B, GE Healthcare) overnight. Subsequently, the agarose was washed with GST wash buffer (50 mM Tris-HCI, 150 mM NaCl, 10 mM MgCl₂, 1 mM DTT, 10% glycerol and 0.1% Triton X-100, pH 8.0) three times for 10 min prior to eluting the fusion protein for 2 h in GST elution buffer (25 mM Tris-HCl, 20 mM glutathione, 100 mM NaCl, 10 mM MgCl₂, 1 mM DTT, and 10% glycerol, pH 7.5). For GST pull-down assays, the agarose beads were prewashed three times with 1 mL of GST-binding buffer (20 mM HEPES, 40 mM KCl, and 1 mM EDTA, pH 7.5). The GST-GmLMM1^ECD fusion protein and GST protein were separately mixed with AtRALF1, AtRALF23, and MiRALF1/3 and then incubated with the prewashed agarose beads (20 μL) in 400 μL of GST-binding buffer at 4°C overnight. Subsequently, the mixture was separated by centrifugation at 500 × g for 1 min, and the precipitate was washed three times with GST wash buffer for 10 min. The washed precipitate was then eluted with 50 μL of GST elution buffer for 2 h. After low-speed centrifugation, the eluate was mixed with SDS loading buffer, boiled, and analyzed using SDS-PAGE and immunoblotting with a GST-Tag (12G8) Mouse mAb (Abmart, Shanghai, China) and a His-Tag (2A8) Mouse mAb (Abmart, Shanghai, China). GST pull-down assays for other RALF proteins were performed using the same procedure.

Proteins were extracted from 7-day GmLMM1-HA transgenic plants (Wang et al., 2020) using NEBT buffer, incubated with AtRALF23 or MiRALF1 at a concentration of 1 μM. After immunoprecipitation with HA beads (Bimake, Houston, USA) overnight at 4°C, the beads were washed three times with NEBT buffer supplemented with 1% TritonX-100 for 15 min at 4°C and then boiled with protein loading buffer. GmLMM1-HA was tested using an anti-HA antibody [HA-Tag (C29F4) Rabbit mAb], and RALFs-His was tested by His-Tag (2A8) Mouse mAb.

GmLMM1^ECD and GmLLG2 were docked onto the 3D structure of the RALF23-LLG2-FER^ECD complex (PDB code: 6A5E) by mutating the corresponding residues of FER^ECD and LLG2 in COOT (Emsley and Cowtan, 2004), guided by sequence alignment of GmLMM1^ECD with FER^ECD and GmLLG2 with LLG2. In the final models, no residues fell in the disallowed regions in the Ramachandran plot, and the complex was perfectly assembled with no obvious clashes in individual interfaces. Structural images were prepared in PyMOL (The PyMOL Molecular Graphics System, version 1.7 Schrödinger, LLC). The structure prediction was performed via AlphaFold (https://www.alphafold.ebi.ac.uk).

Amino acid sequences of GmLMM1 and GmRALFs homologs were obtained by a BLASTP search against Phytozome (version 12.1; http://www.phytozome.net). The phylogenetic analysis was performed as described previously with some modifications (Wang et al., 2020).

Statistical significance was determined by conducting Student's t-tests or one-way ANOVA using SPSS 23.0 software (SPSS).

To identify genes that negatively regulate soybean resistance against M. incognita, we inoculated 600 M2 seedlings that were mutagenized with EMS (Wang et al., 2020) and screened for M. incognita resistant mutants. We found that two mutants were more resistant against M. incognita compared with wild-type soybean (Figures 1A,B), which were named resistant to M. incognita 1 (rmi1) and rmi2, respectively. Precisely, there were significant differences in the nematode number and nematode development between rmi1/2 and wild-type plants. At 3 days post infection (dpi), the number of nematodes was significantly lower in the rmi1 mutant compared with that in the wild-type plants, but no differences in rmi2 (Figures 1A,B). At 7 dpi, most of the nematodes parasitizing wild-type plants turned to the parasitic third stage (par-J3, morphologically thick and short), whereas the nematodes in both rmi1/2 still remained in the parasitic second stage (par-J2 stage, morphologically thin and long) (Figure 1A). Subsequently, at 30 dpi, most of the nematodes reached the reproduction stage, the features of which is the formation of galls. Notably, the galls number in rmi1/2 was significantly lower compared with wild-type plants (Figure 1C). We also calculated female RKNs (morphologically round basal knobs) in rmi1/2 and wild-type plants and found a significant decrease in rmi1/2 compared with wild-type plants (Figure 1D). Since the obvious differences in developmental stages ranging from par-J2 to the parasitic fourth stage (par-J4), statistical analysis of the developmental stage ratio was conducted at 14 dpi. The ratio of nematodes in par-J2 was <70% in wild-type plants compared to nearly 80% in rmi1/2 (Figure 1E). Moreover, the ratio of nematodes in par-J3 almost reached 30% in wild-type plants, while it only reached 15% in rmi1/2 (Figure 1F). Similarly, the ratio of nematodes in par-J4 was nearly 5% in wild-type plants and only 2% in rmi1/2 (Figure 1G). Collectively, there was a significant developmental delay ranging from par-J2 to par-J4 of M. incognita in rmi1/2 mutants compared with wild-type plants, suggesting that rmi1/2 mutants were more resistant against M. incognita parasitism.

Interestingly, rmi1 and rmi2 mutants have not only nematode resistance phenotype, but also spontaneous lesions on their leaves, namely Gmlmm1-1 and Gmlmm1-2 in the recent article Wang et al. (2020) we published. Thus, rmi1 and rmi2 carried two different mutation types in the same predicted gene, Glyma.13G054400 (GmLMM1). rmi1 (Gmlmm1-1) had a C to T substitution within the coding region of Glyma.13G054400, leading to an early stop codon and thus an incomplete protein (Supplementary Figure S1A). Whereas rmi2 (Gmlmm1-2) harbors a T to A substitution in Glyma.13G054400, this substitution only changed the amino acid from histidine to leucine at position 407 (Supplementary Figure S1A) (Wang et al., 2020). For consistency, we named the gene-controlling rmi1 and rmi2 nematode-resistant phenotypes GmLMM1 as previously described (Wang et al., 2020).

GmLMM1 encodes a RLK that contains a malectin-like domain in its extracellular domain, a transmembrane domain, and an intracellular kinase domain in its C-terminal region (Supplementary Figure S1A). The closest homology of GmLMM1 in Arabidopsis is the CrRLK1L FER which shares a high degree of amino acid identity, particularly in the extracellular domain (Supplementary Figures S1B, S2). To confirm the role of GmLMM1 in regulating nematode resistance, we generated heritable mutation in the coding sequence of GmLMM1 by CRISPR-Cas9 mutagenesis (Wang et al., 2020). Among them, line #9 harbored a 4-bp deletion in GmLMM1 (Figure 2A) and showed an autoimmunity-related root phenotype (Figure 2B). Thus, line #9 was selected to further confirm the role of GmLMM1 in the nematode parasitism. Indeed, there were significant differences in the nematode infection number and nematode development between line #9 and wild-type plants (DN50). At 7 dpi, the number of nematodes was significantly lower in line #9 compared with that in the wild-type plants (Figure 2C). Moreover, at 30 dpi, the formation of galls in line #9 was repressed (Figure 2D), and the female number in line #9 was significantly reduced compared with wild-type plants (Figure 2E). Altogether, these results demonstrate that GmLMM1 negatively regulates nematode resistance.

To understand how GmLMM1 negatively regulates soybean resistance at the transcriptional level, we performed RNA-seq on roots of wild-type plants and Gmlmm1-1 with or without M. incognita infection. We observed profound changes in the plant gene expression profile after infection of M. incognita in wild-type plants or after GmLMM1 mutation. A total of 4, 091 differentially expressed genes (DEGs) (p < 0.05, fold change >2 or fold change < 0.5) were detected after infection with M. incognita in wild-type plants (Figure 3A). Intriguingly, a proportion of the DEGs that detected in wild-type plants under the infection of M. incognita showed similar differential expression pattern after GmLMM1 mutation (p < 0.05, Figure 3A). Clustering analysis showed that the expression changes of most GmLMM1 and nematode infection-regulated genes were correlated (r = 0.79, p < 0.001, Figure 3B). Notably, most DEGs observed seem to be unchanged in Gmlmm1-1 after M. incognita infection (Figure 3A), indicating the key role of GmLMM1 in regulation of gene expression in response to nematode infection. Furthermore, Venn diagrams were conducted to explore the intersection of GmLMM1 mutation and nematode infection-regulated genes. Compared with in wild-type soybean, 4,315 DEGs were upregulated and 3,184 DEGs were downregulated in Gmlmm1 mutants, which are marked with yellow or purple, respectively (Figure 3C). Meanwhile, 2869 DEGs were upregulated and 1,222 DEGs were downregulated in wild-type soybean after M. incognita inoculation, which were marked with green or pink, respectively (Figure 3C). Nearly half of DEGs detected in wild-type plants after M. incognita inoculation were observed in the DEGs detected after GmLMM1 mutation, among of which 1,337 DEGs were intersected upon the groups of Gmlmm1 Up and M.i. Up and 571 DEGs were intersected upon the groups of Gmlmm1 Down and M.i. Down (Figure 3C). Gene ontology (GO) analysis revealed that the top enriched GO processes of both upregulated and downregulated intersection DEGs, such as defense response-associated processes, cell wall organization, and response to biotic stimulus (Figure 3D). Notably, previous studies have shown that these processes play an important role in defensing nematode infection (Wang et al., 2019; Zhang et al., 2020a), suggesting that GmLMM1 is indeed involved in transcriptional regulation of soybean resistance against M. incognita. To further get insight into GmLMM1-regulated genes, we selected the intersected DEGs above and performed DNA-binding motifs prediction in their promoters (Figure 3E). Based on the predicted binding motifs, we searched for transcription factors that recognize such DNA-binding motifs and further predicted their potential downstream targets of these transcription factors. GO enrichment analysis showed that the predicted target genes are mostly involved in JA response, salicylic acid response, and ethylene-activated signaling pathway (Figure 3E; Supplementary Figure S3), all of which are key immune signaling pathways against nematodes (Nahar et al., 2011; Martinez-Medina et al., 2017; Wang et al., 2019). Taken together, RNA-seq analysis revealed that mutation of GmLMM1 transcriptionally upregulates defense-associated genes to enhance nematode resistance.

Having ascertained that GmLMM1 regulates nematode resistance, we sought to understand the molecular mechanism behind it. Our recent work identified several RALF-likes from M. incognita (MiRALFs) and found that MiRALFs can mimic the function of plant RALF by interacting with Arabidopsis FER, inhibiting plant immunity and facilitating nematode parasitism (Zhang et al., 2020a). Given that GmLMM1 is a homology of FER, we wondered whether MiRALFs can also interact with GmLMM1. It has been shown that the N-terminal region of Arabidopsis RALF23 is perceived by LLG2 to nucleate the assembly of RALF23-LLG2-FER complex (Xiao et al., 2019). To test whether GmLMM1 shares a similar RALF23-LLG2-FER-binding mode, the extracellular domain (ECD) of GmLMM1 and GmLLG2 was docked onto the 3D structure of the RALF23-GmLLG2-GmLMM1 complex. The result showed that GmLMM1^ECD finely matched the electron density of FER^ECD without causing significant steric clashes in their models (Figure 4A). Next, we expressed four RALFs (AtRALF1/23 and MiRALF1/3) and GmLMM1^ECD with N-terminal His-Tag or GST-Tag, respectively, in E. coli and purified the proteins. Pull-down assay indicated that only AtRALF23 and MiRALF1 were pulled down with the ECD of GmLMM1 (Figure 4B). Notably, AtRALF1 did not interact with GmLMM1. To validate the interaction between RALFs and GmLMM1 in planta, we generated HA-labeled GmLMM1 transgenic Arabidopsis thaliana lines and coimmunoprecipitation (Co-IP) assay was conducted. The result showed that AtRALF23 and MiRALF1 indeed can interact with GmLMM1 in planta (Figure 4C). Meanwhile, there were about 24 RALFs in soybean, and 2 RALFs clustered together with AtRALF23 (Supplementary Figure S4). GST pull-down assay showed that Glyma.07G247500 and Glyma.17G026400 clustered together with AtRALF23 could bind to GmLMM1^ECD (Supplementary Figure S5). Furthermore, we have conducted the infection assays after exogenous MiRALF1 and AtRALF23 peptides treatment in Williams 82 and Gmlmm1-1. The result showed that exogenous MiRALF1 treatment enhanced RKN parasitism 3 dpi, whereas this promotion was abolished in Gmlmm1-1 mutant, suggesting that this process depends on GmLMM1 (Supplementary Figure S6). However, exogenous AtRALF23 treatment had no significant impact on nematode infection (Supplementary Figure S6). Taken together, these results suggest that AtRALF23 and MiRALF1 potentially work as the ligands of the RLK GmLMM1, which also echoes the sight of our recent work (Zhang et al., 2020a). We assume that targeting of GmLMM1 by MiRALF1 suppresses soybean immunity and promote nematode infection.

To investigate the relationship between GmLMM1 and soybean resistance against M. incognita, we examined the variation in the coding region of GmLMM1. Four SNPs were identified in the coding sequence of GmLMM1, and three of them directing amino acid substitutions (p. T9P, p. T40N, and p. N359H) (Figure 5A, Supplementary Table S2). Based on the three nonsynonymous mutation, five haplotypes (H1-H5) were defined among 232 soybean varieties (Supplementary Table S2). To further investigate which GmLMM1 haplotype was related to the resistance against M. incognita, 27 cultivars were selected as representatives of the indicated five haplotypes for further nematode infection assays (Figures 5A–D). A total of 11 varieties, including Williams 82, belong to H1, which shows no difference in the three SNPs of GmLMM1 compared with reference sequence (Glycine max W82.a2.v1) (Figure 5A, Supplementary Table S2). After inoculation with M. incognita, the varieties belong to H3 or H5 showed enhanced resistance against M. incognita. In additional, the numbers of nematode parasitizing on H3 or H5 haplotypes were significantly lower compared with the other three haplotypes (Figure 5B). Furthermore, the representative varieties of H1 to H5 were selected to calculate the galls and female number. At 30 dpi, the galls in H3 and H5 plants were significantly lower than other varieties (Figure 5C), and the numbers of female RKNs in H3 and H5 lines were also significantly reduced at 30 dpi (Figure 5D), suggesting that the cultivars carried haplotypes 3 or 5 would improve the resistance against RKN. In addition, we test the polymorphisms of the three SNPs in GmLMM1 contribution to the resistance to RKN. Of these, the p. N359H substitution in particular affects obviously to the resistance to RKN (Supplementary Figure S7), implying that the SNP played an important role for soybean resistance against to RKNs. Overall, these results suggest that GmLMM1 is essential for the parasitism of M. incognita in soybean.