MiMIF-2 transgenic A. thaliana lines were previously described (Zhao et al., 2019). Wild-type (WT) A. thaliana (Col-0) was used as a control. Surface-sterilized A. thaliana seeds were placed on Murashige and Skoog Medium (MS) (Sigma, St. Louis, MO, USA) and were incubated at 22°C with a 12-h photoperiod. Seedlings of 15 days old plants were transplanted into soil and were grown in the greenhouse at 22°C and 65% relative humidity with a 16-h light/8-h dark photoperiod.

Tomato plants (Solanum lycopersicum var. “Baiguo”) were used to reproduce M. incognita nematode culture. Infective pre-parasitic second-stage juveniles (pre-J2s) were collected after egg hatching (Zhao et al., 2019).

Homozygous T3 MiMIF-2 transgenic A. thaliana lines and WT were challenged with 200 pre-parasitic J2s of M. incognita, indicated by MiMIF-2+N and WT+N. WT and MiMIF-2 transgenic A. thaliana lines of uninoculated nematodes were used as the controls, indicated by WT and MiMIF-2. All the plants were grown in the laboratory of the Institute of Vegetable and Flower Research of the Chinese Academy of Agricultural Sciences, Beijing, China (116.19.32°E, 39.57.44°N). Rhizosphere soil samples around the A. thaliana roots were collected at 30 days post-inoculation (dpi) with nematodes. After excavating the roots of 30 Arabidopsis plants from the soil, the soil that was loosely attached to the roots was carefully removed. The soil that adhered to the Arabidopsis root system was brushed loose using a delicate brush. This soil was collected in sterile self-sealing bags and mixed thoroughly, and the amalgamated soil was divided into five equal portions. Each treatment was replicated five times with five seedlings per replication (Lundberg et al., 2012).

Microbial community genomic DNA was extracted using 0.5 g of soil per sample employing the E.Z.N.A.^® soil DNA Kit (Omega Bio-tek, Norcross, GA, USA) as per the manufacturer's recommendations. The quality of the DNA extract was confirmed on 1% agarose gel, whereas the DNA concentration and purity were determined using the NanoDrop 2000 UV-vis spectrophotometer (Thermo Scientific, Wilmington, USA). Amplification of the hypervariable region, V3–V4 of the bacterial 16S rRNA gene, was performed using the primer pair 338F (5'-ACTCCTACGGGAGGCAGCAG-3') and 806R (5'-GGACTACHVGGGTWTCTAAT-3'), whereas the primer set SSU0817 forward (5′-TTAGCATGGAATAATRRAATAGGA-3′) and 1196 reverse (5′-TCT GGACCTGGTGAGTTTCC-3′) was used to amplify the V5–V7 region of the 18S rDNA gene through the ABI GeneAmp^® 9700 PCR thermocycler (ABI, CA, USA).

The PCR method used for the amplification comprised of an initial denaturation at 95°C for 3 min, followed by 27 cycles of denaturation at 95°C for 30 s, annealing at 55°C for 30 s, and extension at 72°C for 45 s in each cycle with a final extension at 72°C for 10 min after which the reaction was maintained at 4°C. The PCR mixtures contained 5× TransStart FastPfu buffer 4 μl, 2.5 mM dNTPs 2 μl, forward primer (5 μM) 0.8 μl, reverse primer (5 μM) 0.8 μl, TransStart FastPfu DNA polymerase 0.4 μl, template DNA 10 ng, and ddH₂O to make a final volume of 20 μl. PCR amplifications were performed in triplicate. The PCR products were electrophoresed using 2% agarose gel, and the band of interest was eluted using the AxyPrep DNA Gel Extraction Kit (Axygen Biosciences, Union City, CA, USA) as per the recommendations of the manufacturer, and the eluted DNA was quantified using the Quantus™ Fluorometer (Promega, USA).

Purified amplicons pooled in equimolar volumes were paired-end sequenced (2 × 300) on an Illumina MiSeq platform (Illumina, San Diego, USA) following the standard protocols recommended by Majorbio Bio-Pharm Technology Co. Ltd. (Shanghai, China). The raw reads were deposited into the NCBI Sequence Read Archive (SRA) database (accession number: PRJNA954129).

Splicing and filtering were performed on raw data to secure valid data and to ensure the accuracy and reliability of the results. Operational taxonomic units (OTUs) were clustered using 97% similarity. The Chao and Shannon indices were computed to figure out the alpha diversity and to determine the species richness in each soil sample. A principal coordinates analysis (PCoA) of each sample formulated on the Bray–Curtis distances was accomplished using the R package vegan (version 2.1). Linear discriminant analysis (LDA) effect size measurements (LEfSe) was performed to search for significantly different (P < 0.05) taxa between the two groups. Significantly different taxa were used to generate taxonomic cladograms, which illustrated the differences between sample classes on the website http://huttenhower.sph.harvard.edu/galaxy.

Graphical representations were created using GraphPad Prism 5 (GraphPad Software, Inc., La Jolla, CA, USA). The mean and standard error for each group of data were calculated using the ANOVA and Tukey's honestly significant difference test (P < 0.05) using SPSS 19.0 (SPSS Inc., Chicago, IL, USA).

After Illumina paired-end sequence, quality evaluation, data filtering, and integration, a total of 933,427 bacterial 16S rRNA sequences (Supplementary Table S1) and 1,169,253 effective fungal 18S rDNA sequences (Supplementary Table S2) were generated from 20 samples. A total of 2,299 bacterial OTUs were assigned to 38 different phyla: WT contained 28, MiMIF-2 contained 32, WT+N contained 30, and MiMIF-2+N contained 32 (Table 1). A total of 492 fungal OTUs were assigned to 9 different phyla: WT contained 9, MiMIF-2 contained 9, WT+N contained 9, and MiMIF-2+N contained 8 (Table 2).

By comparing the α-diversities of the fungal species richness (Chao1) index, there was a significant difference between the before and after inoculation nematodes. However, the fungal species richness in MiMIF-2 was significantly reduced compared to WT, when inoculated with M. incognita (P < 0.05, Wilcoxon rank sum test; Figure 1A).

In the case of bacterial species richness index, however, MiMIF-2 lines showed a significant increase compared to WT after inoculation with M. incognita. The bacterial species richness (Chao1) index was also reduced in wild-type Arabidopsis when compared with uninoculated WT (Figure 1B). Interestingly, the species richness index of WT inoculated with M. incognita was significantly lower than that of uninoculated samples (Figure 1B). The increase in bacterial species diversity might be the reason for the increased nematode susceptibility in MiMIF-2.

Statistical analysis of Bray–Curtis distances (β-diversity) was involved in the comparison of the composition of microbial communities. It was performed to determine the community treatment relationship between the samples. Concerning the phylum level, PCoA was performed on each sample, and the distribution of WT+N and other samples was significantly different for bacteria (Figure 1C). At the OTU level, PCoA of Bray–Curtis distances (β-diversity) revealed that WT, MiMIF-2, WT+N, and MiMIF-2+N microbiota in four types of soil exhibited a clear separation (Figure 1D).

OTUs of the four samples were assigned corresponding taxonomies based on the combined search results against the Greengenes and NCBI databases. The relative abundance of different phyla in the four samples was analyzed, of which Proteobacteria, Actinobacteria, and Firmicutes were the dominant bacterial phyla (Figure 2A). Ascomycota and Basidiomycota were the dominant fungal phyla in all soil samples, accounting for more than 90% of the total abundance in each sample (Figure 2B). After inoculation with nematodes, Firmicutes of WT+N were significantly reduced compared to uninoculated samples (Figure 2A). In community composition, Clitopilus, Candida, Meliniomyces, and Pseudeurotium were the dominant fungal genera in all soil samples (Figure 2C) whereas in the bacterium community composition, Komagataeibacter, Lactobacillus, Streptomyces, and Rhodanobacter were the dominant bacterial genera in all the rhizosphere bacteria. Lactobacillus of WT+N were significantly reduced compared to uninoculated samples. On the contrary, Streptomyces, Rhodanobacter, and Actinoallomurus of WT+N species increased compared with other samples (Figure 2D).

The relative abundances of major fungal phyla and genus, including Basidiomycota, Clitopilus, and Candida, were different in MiMIF-2 transgenic A. thaliana (Figure 3A). Phylum Basidiomycota abundance in MiMIF-2 transgenic A. thaliana decreased compared with WT; however, the relative abundances increased after inoculation with M. incognita. Similar results were obtained during the analysis at the genus level. Clitopilus behaved similarly as Basidiomycota. Candida was abundant in MiMIF-2 transgenic A. thaliana; however, no differences were observed after inoculation with M. incognita (Figure 3B).

The analysis of bacterial phylum-level relative abundances revealed significant differences among Firmicutes, Actinobacteria, Acidobacteria, Epsilonbacteraeota, Spirochetes, and Chloroflexi. Firmicutes, Epsilonbacteraeota, Spirochetes, and Fusobacteria were significantly higher in MiMIF-2+N compared to those in WT+N. However, the richness of the phylum of WT was significantly reduced after inoculation with nematodes. Actinobacteria, Acidobacteria, and Chloroflexi were reduced in MiMIF-2+N compared to WT+N (Figure 4A). A similar phenomenon occurs at the genus-level relative abundances. More than 15 genera showed significant differences. For instance, Komagataeibacter, Lactobacillus, Prevotella, Moritella, Vibrio, Escherichia-Shigella, and Pseudomonas were reduced in MiMIF-2 lines compared to those in WT. However, these genera were significantly higher in MiMIF-2+N compared to those in WT+N (Figure 4B).

The LEfSe was used to identify the dominant phylotypes responsible for the differences among the samples. The groups were shown in cladograms, and LDA scores of 4 (Supplementary Tables S3, S4) or greater were confirmed using LEfSe (Figure 5). In MiMIF-2+N, 11 groups of bacteria were significantly enriched, while no fungi were detected at a significant level. The bacteria with the highest LDA value in MiMIF-2+N was Firmicutes (logarithmic LDA score = 4.73). In WT+N, 4 groups of bacteria and 3 groups of fungi were significantly enriched. The bacterium and fungus with the highest LDA value in WT+N were Gammaproteobacteria (logarithmic LDA score = 4.57) and Pleurotaceae (logarithmic LDA score = 4.87), respectively. In WT, 5 groups of bacteria and 5 groups of fungi were detected at a significant level. The bacterium and fungus with the highest LDA value in WT were Proteobacteria (logarithmic LDA score = 4.98) and Entolomataceae (logarithmic LDA score = 4.96), respectively. In MiMIF-2, 4 groups of bacteria and 6 groups of fungi were detected at a significant level. The bacterium and fungus with the highest LDA value in MiMIF-2 were Actinobacteriota (logarithmic LDA score = 4.79) and Saccharomycetes (logarithmic LDA score = 4.84), respectively. Community structure and predominant taxa of different samples may account for differences in sensitivity to nematodes.

The relative abundance of PICRUSt inferred function is shown in Figure 6A. All soil rhizosphere microorganisms consist mainly of energy production and conversion, amino acid transport and metabolism, carbohydrate transport and metabolism, cell wall/membrane/envelope biogenesis, and general function prediction. Whether or not the nematode is inoculated, the relative abundance of all the COG function classifications in WT was slightly larger than that in MiMIF-2 lines (Supplementary Table S5). For example, defense mechanisms, secondary metabolites biosynthesis, transport, and catabolism may be related to nematode infection.

We used a manually curated set of designations based on FunGuild, a recently established fungal classification tool with rigorously defined and referenced trophic group assignments. The relative abundance of animal endosymbiont-undefined saprotroph was higher after inoculation with nematodes; however, it was almost none after transferring to the MiMIF-2 effector (Figure 6B).