Intergrated Metagenomics and Metabolomics Analysis Discovers Nematicidal Microbes, Enzymes and Metabolites From the Plant Rhizosphere Microbiota


 BackgroundRoot-knot nematode Meloidogyne incognita infects root systems of many crops resulting in huge decrease of crop production. Nematicidal microorganisms provides a safe and effective strategy to control M. incognita infection. In order to find more microorganisms with high activity and new nematicidal metabolites, we collected the M. incognita infected tobacco rhizosphere soils (RNI) and non-infected tobacco rhizosphere soils (NS), and investigated their microbial community and network via metagenomics and metabolomics analysis. ResultsMicrobial networks of RNI soils were very different from the NS soils. Many nematicidal microorganisms were enriched in the NS soils, including some isolates such as Aspergillus , Achromobacter , Acinetobacter , Bacillus , Burkholderia , Comamonas , Enterobacter , Lysobacter , Microbacterium , Paenibacillus , Pantoea , Pseudomonas , Streptomyces and Variovorax. Enzymes analysis showed these nematicidal microorganisms can produce proteases, chitinase and lipases. The functions genes belonging to pathways of secondary metabolites biosynthesis and carbohydrate transport and metabolism were overrepresented in the rhizophere microbiota of NS soils comparing with the RNI soils. 102 metabolites contents were significantly different between the RNI and NS rhizosphere microbiota. 35 metabolites were overrepresented in the NS soils comparing the RNI samples, including acetophenone. Acetophenone showed high nematicidal (LC 50 = 0.66 μg/ml) and avoidance activity against M. incognita . A isolate of Bacillus amyloliquefaciens W1 with production of acetophenone can kill 98.8% of M . incognita . ConclusionsIn general, the rhizophere microbiota of NS soils could produce volatile materials, multiple enzymes and secondary metabolites against nematode. Collectively, the microbiota of NS and RNI rhizophere differed significantly in microbial network structure, community composition, function genes and metabolites. Collectively, combination of multi-omics analysis and culture-dependent technology is powerful for finding nematicidal microorganisms and metabolites from soil.


3
The root-knot nematode Meloidogyne incognita can parasitize plant roots and causes dramatic yield losses in many different crops worldwide [1], but the management options available for effectively controlling root-knot nematode are extremely limited. Although traditional control methods such as crop rotation, disease resistance breeding and use of chemical pesticides have reduced nematode diseases indeed, the crop yield losses caused by M. incognita are still enormous. Fortunately, some microbial agents can effectively control root-knot nematodes to reduce M. incognita infection. These microorganisms include rhizosphere bacteria, nematophagous fungi, specialized parasitic bacteria, actinomycetes, etc, [2,3]. Among them, rhizosphere microbiota are especially interested, and currently, there is a growing interest in studying rhizosphere microbiota in the root-knot nematode infected soils [4].
Many rhizosphere microorganisms have been isolated and selected as biological agents for controlling root-knot nematodes [5][6][7]. However, application of microbial agents also raised a novel problem that is the control effect of microbial agents are often unstable in the field [8]. One of the main reasons is that these beneficial microorganisms can't stably colonize in soil and plant rhizosphere for a long time due to their disability to adapt to complex soil environment and climate conditions [9]. In particular, many farmland soils in the Southern China are seriously degraded and the soil ecological environment is deteriorating year by year, resulting in "bad" soil microbial communities and unbalanced networks. Thus, the colonization of beneficial microorganisms will be very difficult in these soils [10,11].
Generally, researchers often pay more attentions to screening of microbes for biocontrol of root-knot nematodes, but the interaction relationships among root-knot nematodes, rhizophere microbiota and biocontrol microbes are rarely understood in the soil ecosystem.
Nowadays, how to control the soil micro-ecological environment and make the beneficial microorganisms colonize more stable and exert their effects furthest in soils is one of the key points to improve their biocontrol effects. Investigation of microbial interactions may help to resolve this bottle-neck problem. Microbial interactions encompass a spectrum ranging from antagonistic and competition (resource competition) to cooperative, syntrophic and mutualistic interactions. Network analysis is favorable for exploring the organisation and dynamics of microbial interactions and niches [12. 13], and help us to find "good friends" (e.g. nematicidal microorganisms) and improve their colonization efficiency in soils. 4 Metagenomics and metabolomics analysis have enabled the study of microbial ecosystem structure to a greater depth and accuracy. Metagenomics sequencing technology is an interesting tool for selection of the desirable biocontrol microbes [14], and metabolomics is favorable for finding active metabolites produced by microorganisms [15]. The main object of this project is to study the important functional genes and key microorganisms related to the occurrence of root-knot nematode disease, and screen novel nematicidal microorganisms and metabolites against root-knot nematodes based on the metagenomic and metabolomic analysis. This study will provide a way for effective control of root-knot nematodes by improving soil micro-ecology, promoting the colonization and efficiency of beneficial microorganisms, and developing new microbial metabolites with anti-nematode activity.

Soil sampling and properties analysis
Fifteen root-knot nematode M. incognita infected tobacco fields (RNI) and fifteen non-infected tobacco fields (NS) were located in Enshi City, Hubei Province, China (29°22'-29°58'N, 109°20'-109°22'E), with elevations ranging from 990 to 1,127 m and subtropical humid climate. The soil is yellow-brown soil (classified as Alfisols). Rhizosphere soil samples were collected from ten tobacco plants in each field (667 m 2 ) in August 2018 at the maturation stage of tobacco. The tobaccos were removed gently from the field, then the loosely adherent soils were removed by vigorous shaking, and the soils adhering to root were collected with a brush. The collected rhizosphere soil of ten plants in each field was mixed together as a composite soil sample. Soil samples were sieved (2-mm mesh) to remove plant roots, debris and stones. Each soil sample was divided into two parts: one was stored at -80°C for DNA and metabolites extraction within one week, and the other one was stored at 4°C for soil properties and root-knot nematode analysis. Quantification of M. incognita in soil samples was determined via real-time quantitative PCR [16].
Soil organic matter (SOM) content was assayed by the acidified potassium dichromate (K2Cr2O7-H2SO4) heating method [17]. Alkali-hydrolyzable nitrogen (AN) content was determined by alkaline hydrolysis diffusion method. Available phosphorus (AP) and potassium (AK) contents were determined photometrically by a flame spectrophotometer.
After suspending with deionized water (soil : water = 1 : 2.5, w/v), soil pH was determined 5 using a pH meter. Urease activity was detected by colorimetric determination of ammonium.
Catalase activity was determined by colorimetric assay using K2Cr2O7/acetic acid reagent.

Metagenomic sequencing data pretreatment and metagenome assembly
For each soil sample, DNA was extracted from 1.0 g soils using the FastDNA Spin Kit (MP Biomedicals, USA). DNA concentration was measured using Qubit® dsDNA Assay Kit in Qubit® 2.0 Flurometer (Life Technologies, CA, USA). 1 μg DNA per sample was used to construct library. Sequencing libraries were generated using NEBNext® Ultra™ DNA Library Prep Kit for Illumina (NEB, USA) and index codes were added to attribute sequences to each sample. Briefly, the DNA sample was fragmented by sonication to a size of 350 bp, then DNA fragments were end-polished, A-tailed, and ligated with the full-length adaptor for PCR amplification. PCR products were purified by AMPure XP system and libraries were analyzed for size distribution by Agilent2100 Bioanalyzer and quantified using real-time PCR.
Clustering of the index-coded samples was performed on a cBot Cluster Generation System according to the manufacturer's instructions. After cluster generation, the library preparations were sequenced on an Illumina HiSeq platform and paired-end reads were generated.
Preprocessing the raw data using Readfq V8 was conducted to acquire the clean data for subsequent analysis. The reads that are of host origin were filtered using Bowtie2.2.4 software.
The clean data was assembled and analyzed by MEGAHIT software (v1.0.4-beta), then interrupted the assembled Scaftigs from N connection and leave the Scaftigs without N [20].
All samples' clean data were compared to each Scaffolds respectively by Bowtie2.2.4 software to acquire the PE reads not used. All the reads not used in the forward step of all samples were combined and then used the software SOAPdenovo V2.04/MEGAHIT v1.0.4-beta for mixed assembly, broke the mixed assembled Scaffolds from N connection and obtained the Scaftigs, then filtered the fragment shorter than 500 bp in all of Scaftigs for 6 statistical analysis.

Gene prediction and abundance analysis
The Scaftigs (≥ 500 bp) assembled from both single and mixed assembly were all predicted the ORF by MetaGeneMark V2.10 software. For ORF prediction, the CD-HIT V4. 5.8 software was adopted to redundancy and obtained the unique initial gene catalogue [21]. The clean data of each sample was mapped to initial gene catalogue using Bowtie2.2.4 and get the number of reads to which genes mapped in each sample. We filtered the gene with the number of reads ≤ 2 in each sample and obtained the gene catalogue Unigenes eventually used for subsequently analysis. Based on the number of mapped reads and the length of gene, the abundance information of each gene in each sample was counted.
Taxonomy prediction DIAMOND V0.9.9 software was used to blast the Unigenes to the sequences of bacteria, fungi, archaea and viruses which are all extracted from the NR database of NCBI [22]. The number of genes and the abundance information of each sample in each taxonomy hierarchy (kingdom, phylum, class, order, family, genus, specie) were obtained. The abundance of a specie in one sample equal the sum of the gene abundance annotated for the specie. The gene number of a specie in a sample equal the number of genes whose abundance are nonzero.
Principal component analysis (PCA) and non-metric multi-dimensional scaling (NMDS) analysis were performed based on the abundance of each taxonomic hierarchy. The difference between two groups was tested by Anosim analysis. Metastats and LEfSe analysis were used to find the different species between groups. Permutation test between groups was used in Metastats analysis for each taxonomy and get the P value, then used Benjamini and Hochberg False Discovery Rate to correct P value and acquired q value [23]. LEfSe analysis was conducted by LEfSe software (the default LDA score is 3). Finally, random forest (RandoForest) was used to construct a random forest model, screen out important species by MeanDecreaseAccuracy and MeanDecreaseGin, then cross-validate each model and plot the ROC curve.

Functional database annotations
We adopted DIAMOND V0.9.9 software to blast Unigenes to functional database KEGG, 7 eggNOG, and CAZy database. For each sequence's blast result, the best Blast Hit was used for subsequent analysis. The relative abundance of different functional hierarchy was counted.
The gene number of each sample in each taxonomy hierarchy was obtained. Based on the abundance of each taxonomy hierarchy, not only the counting of annotated gene numbers, the exhibition of the general relative abundance situation, the exhibition of abundance cluster heat map, and the decrease-dimension analysis of PCA and NMDS were conducted, but also the Anosim analysis of the difference between groups based on functional abundance, comparative analysis of metabolic pathways, the Metatat and LEfSe analysis of functional difference between groups were performed.

Network construction and analysis
The relative abundances of species or functional genes were used to construct microbial or functional networks using the Molecular Ecological Network Analysis Pipeline (MENAP), respectively. Briefly, species or gene abundances were normalized to the standardized relative abundances (SRA) [12]. A matrix of species SRA or genes SRA, a matrix of soil variables, and a species or genes annotation file were prepared. The SRA matrix was submitted to MENAP to construct the NS and RNI microbial or functional networks, respectively [24,25].
Modules were detected by the greedy modularity optimization. Three files were generated for network graph visualization by Cytoscape 3.7.2 software. The network graph was represented using different species or genes (nodes) with positive or negative interactions (edges).
Positive interactions indicate the abundances of these species or genes changed following the same trend across different soil samples. Negative interactions indicate the abundances of those species or genes changed following the opposite trend in different soil samples [26,27].
The topological role of each node (specie) was defined by two parameters: within-module connectivity (Zi) and among-module connectivity (Pi). Zi value described how well a node was connected to other nodes within its own module. Pi value described how well a node was connected to different modules [28]. Threshold values of Zi and Pi for categorizing nodes (connectors, module hubs) and network hubs are the keystone microorganisms, which are critical in maintaining the network stability and microbial community structure and function [13].

Untargeted metabolomics analysis
Metabolites of five NS samples and five RNI samples were measured by untargeted metabolomics analysis. For each soil sample, 100 mg soils were grounded with liquid nitrogen and the homogenate was suspended with chilled methanol and 0.1% formic acid. The samples were incubated on ice for 5 min and then were centrifuged at 15,000 g and 4°C for 5 min. The supernatant was diluted to final concentration containing 60% methanol by LC-MS grade water. The samples were filtrated with 0.22 μm filter and then were centrifuged at 15,000 g and 4°C for 10 min. Finally, the filtrate was injected into the LC-MS/MS system. LC-MS/MS analyses were performed using a Vanquish UHPLC system (Thermo Fisher) coupled with an Orbitrap Q Exactive HF-X mass spectrometer (Thermo Fisher). Samples were injected onto an Hyperil Gold column (100 × 2.1 mm, 1.9 μm) at a flow rate of 0.2 ml/min. The eluents for the positive polarity mode were eluent A1 (0.1% formic acid in water) and eluent B (methanol). The eluents for the negative polarity mode were eluent A2 (5 mM ammonium acetate, pH 9.0) and eluent B (methanol). The solvent gradient was set as follows: Q Exactive HF-X mass spectrometer was operated in positive/negative polarity mode with spray voltage of 3.2 kV, capillary temperature of 320°C, sheath gas flow rate of 35 arb and aux gas flow rate of 10 arb. The raw data generated by UHPLC-MS/MS were processed using the Compound Discoverer 3.0 (CD 3.0, Thermo Fisher) to perform peak alignment, peak picking, and quantitation for each metabolite. Peak intensities were normalized to the total spectral intensity. The normalized data were used to predict the molecular formula based on 9 additive ions, molecular ion peaks and fragment ions, and then peaks were matched with the mzCloud and ChemSpider databases to obtain the accurate qualitative and relative quantitative results. These metabolites were annotated using the KEGG, HMDB, and Lipidmaps databases [29,30].
Principal components analysis (PCA) and partial least squares discriminant analysis (PLS-DA) were performed at metaX. Univariate analysis (t-test) was applied to calculate the statistical significance (P-value). The metabolites with VIP > 1 and P-value < 0.05 and fold change ≥ 2 or FC ≤ 0.5 were considered to be differential metabolites. Volcano plots were used to filter metabolites of interest which based on Log2(FC) and -log10(P-value) of metabolites. For clustering heat maps, the data were normalized using z-scores of the intensity areas of differential metabolites and were plotted by Pheatmap package in R language. The functions of these metabolites and metabolic pathways were studied using the KEGG database.
Bacillus amyloliquefaciens W1 culture was analyzed by untargeted metabolomics. W1 was incubated in 100 ml of M9 medium (MgSO4·7H2O 0.492 g, CaCl2 0.011 g, Na2HPO4·12H2O 17 g, KH2PO4 3 g, NaCl 0.5 g, NH4Cl 1 g, glucose 4 g, 1000 ml ddH2O) at 37°C and 180 rpm for 48 h. The culture was centrifuged at 8,000 g for 15 min, and the supernatant was collected and mixed with equal volume of ethyl acetate, then centrifuged at 8,000 g for 15 min. The extracted liquid was dried by rotary evaporation at 65°C, then dissolved in 2 ml ethanol following with filtration using filter membrane. 100 μl of sample was mixed with 400 μl of 80% methanol and 0.1% formic acid by well vortexing, incubated on ice for 5 min, then centrifuged at 15,000 g and 4°C for 10 min. The supernatant was diluted to final concentration containing 53% methanol by LC-MS grade water. After centrifugation at 15,000 g and 4°C for 10 min, the supernatant was collected for detection by the LC-MS/MS system and the data analysis was conducted as the above mentioned methods.

Analysis of nematicidal activity of metabolites
M. incognita was maintained on tomato roots. A lot of galls formed in the roots after infection.
Nematode eggs were isolated from galls. Before detecting the nematicidal activity of metabolites and microbes, the galls were peeled off from root and placed in water at 20°C until the second-stage juveniles (J2) were hatched. J2 juveniles were used in the following 10 nematicidal bioassays.
The nematicidal activity of metabolites were detected as previously described [31]. The nematicidal activities of different concentrations of acetophenone (0.412, 0.459, 0.515, 0.589, 0.687, 0.824, 1.03, and 2.06 μg/ml) were tested against M. incognita J2 juveniles for 24 h, respectively. Each concentration of acetophenone was detected with three replicates. The LC50 and LC90 values of acetophenone were calculated. Assay for nematode avoidance of acetophenone was conducted as our previous study [32]. Water-agar medium (1.5% agar) was prepared in 10-cm plates. Two or four μL acetophenone were dropped on one side of the plate, and 2 or 4 μL ddH2O were dropped on the other side as a control. Approximately 50 larvae of M. incognita were transferred to the middle of plate. After 2 h, the nematodes in each section were counted and avoidance index was calculated.

Isolating nematicidal microorganisms from rhizospheric soils
Ten grams of the NS rhizospheric soil was mixed with 90 ml of sterile water, then used for spreading TSB agar plates (tryptone 1.5%, soy peptone 0.5%, NaCl 0.5%, and agar 1.5%, pH 7.2) for isolating bacteria, and PDA plates (potato 20%, glucose 2%, agar 1.5%) for isolating fungi after serial dilution with sterile water. After incubation, the bacteria and fungi with different morphology were isolated from plates for screening the ones with nematicidal activity For nematicidal activity assay, the isolated bacteria were incubated in LB medium at 37°C and 200 rpm for 48 h, and the isolated fungi were incubated in PD medium at 28°C and 200 rpm for 5 days. After incubation, the broth was centrifuged at 8,000 g for 5 min for collecting supernatant. 200 μl of supernatant and~50 numbers of M. incognita J2 juveniles were mixed together and then added to each well in 96-well plates. After incubation at 20°C for 24 h, the mortality of nematode was counted. The nematodes without detectable movement were judged as dead. Tests were done with three replicates. LB or PD medium was added in the assay as control [31].
Nematicidal bacterial and fungal isolates were identified based on 16S rRNA or 18S rRNA genes sequence. Briefly, the bacterial or fungal genomic DNA was extracted using the DNA isolation kit (TIANGEN), then the bacterial 16S rRNA gene or fungal 18S rRNA gene 11 was amplified, respectively [33,34]. Sequence alignment was performed with Clustal W, and the phylogenetic tree was constructed using the software MrBayes 3.27, then visualized with iTOL [35,36].
The enzymes activity including chitinase, protease and lipase were detected in the isolated nematicidal bacteria and fungi, respectively. Protease activity was detected using the clearing zones method [37]. Single bacterial or fungal colony was incubated on the surface of selection medium (nonfat dried skimmed milk 5%, agar 2%) plates at 28°C for 48 h, then the diameters of clearing zone (D2) and the bacterial or fungal colony (D1) were measured, respectively. The ratio of D2 to D1 was calculated as an indicator of protease activity.
Chitinase activity was also detected using the clearing zones method on the surface of chitin medium (NaCl 1.5 g, KH2PO4 0.3 g, K2HPO4 0.7 g, FeSO4·H2O 0.02 g, MgSO4 0.5 g, agar 20 g, 250 ml of 2% colloidal chitin, adding H2O to 1000 ml, pH 7.0). Lipase activity was detected on the surface of selection medium (peptone 1%, yeast powder 0.5%, NaCl 1%, glyceryl triacetate 3%, bromcresol purple 0.004%, agar 2%) at 28°C for 48 h. Glycerol triacetate produces free fatty acids under the catalysis of lipase, which reduces the pH of the medium and makes medium change from blue purple to yellow. The diameters of the yellow clearing zone (D2) and the bacterial/fungal colony (D1) were measured, and the ratio of D2 to D1 was calculated as an indicator of lipase activity.

Statistics analysis
Differences in soil properties, compounds, microbial or gene abundances and network indexes between the NS and RNI soils were compared by a least-significant-difference (LSD) test.
Correlations between metabolite and nematode number was analyzed by Pearson correlation analysis with SPSS 20 software [18].

Soil properties were different between the NS and RNI soils
It was found that the number of M. incognita was very different between the non-infected tobacco rhizosphere soil (NS) and the root-knot nematode infected tobacco rhizosphere soil (RNI). The number of M. incognita in NS was 0.08 to 3.56, with an average value of 0.98; whereas the number of M. incognita in RNI was 325.21 to 3070.75, with an average value of 1091.76 (table S1). The number of M. incognita in RNI was significantly (P < 0.05) higher than that in NS.
Soil properties were also different between NS and RNI. The activities of invertase and catalase, pH, CEC, Fe, Ca, SWC, STP, SCMC and MWD were all significantly (P < 0.05) higher in the NS soils compared to the RNI soils. Conversely, the phosphatase activity, AK and BD were all significantly (P < 0.05) higher in the RNI soils than the NS soils. Catalase activity, CEC, Fe, Ca, STP and MWD were significantly (P < 0.05) negatively correlated with the number of M. incognita, while AP and BD were significantly positively correlated with the number of M. incognita (table S2).
Phanerochaete chrysosporium is a decomposer in carbon cycle (such as lignin degradation) [40]. These results suggested that the composition of microbial networks had varied considerably in the RNI soils when compared with the NS soils. It seems that the microbial network topology structure and composition were substantially shaped in the root-knot nematode infection soils in this mountain ecosystem.
We further analyzed the keystone microorganisms in bacterial networks ( fig. 2A

Functional diversity of the rhizosphere microbiota
Of the genes retrieved from the metagenome data, 58 to 61% were assigned to a known function (fig. S1B, C). Among 979,916 genes with associated functions, 3,641 genes and 1611 functions were significantly (P < 0.05) enriched in the rhizophere microbiota of the non-infected soils. These genes belonged mainly to the pathways classified as "Energy production and conversion", "Amino acid transport and metabolism", "Signal transduction mechanisms", "Secondary metabolites biosynthesis, transport and catabolism" and "Carbohydrate transport and metabolism".
Co-occurrence function network analysis revealed the average path distance (GD) of the NS function network was smaller than the RNI function network, which means that all the nodes in the NS network are closer than those in the RNI network. Average clustering coefficient (avgCC) and modularity of the NS network were higher than the RNI network, indicating the NS function network has more stability and resilience than the RNI network.
AvgK of the RNI function network was higher than the NS network (table S7) For more detailed resolution of the specific functions associated with cluster of orthologous groups (COG) Category Q and Category G (carbohydrate transport and metabolism), we searched for carbohydrate-active enzymes (CAZymes) and secondary metabolite bio-synthetic gene clusters within the metagenome sequences using dbCAN [41] and antiSMASH [42], respectively. Among Category Q, 41 genes were significantly   S10). APE Vf and bicornutin were significantly overrepresented (P < 0.05) in the NS microbiota comparing with the RNI microbiota.

Untargeted Metabolomics analysis of the rhizosphere microbiota
PCA analysis showed that metabolites composition of the NS rhizosphere microbiota were different from the RNI rhizosphere microbiota ( fig. S11A, B). Partial least squares discrimination analysis (PLS-DA) showed the RNI samples were well separated from the NS samples ( fig. 4A, B). Different metabolites analysis showed 102 metabolites were significantly different between the RNI and NS rhizosphere microbiota ( fig. 4 C, D; fig. S12A fig. S13D).

Discussion
Soil properties of the NS soils were very different from the RNI soils. Invertase involved in soil carbon cycle is a characterization of soil fertility level. Catalase is associated with soil respiration and microbial activity. Here, these two enzymes activities in the NS rhizosphere soils were higher than those in the RNI soils, suggesting that microbial activity in the non-infected soils is higher than that in the root-knot nematode infected soils. Soil acidification was favorable for the hatching of root-knot nematode eggs and could increase the density of second instar larvae in soil, so aggravated the occurrence of root-knot nematode disease [43]. The pH value of the RNI soils was much lower than the NS soils, thereby soil acidification in this mountain region possibly aggravated the root-knot nematode disease. Soil total porosity and mean weight diameter of the NS soils were higher than the RNI soils, implying the non-infected soils were more loose to be favorable for improving microbial activity. Instead, the root-knot nematode infected soils showed compaction phenomenon (e.g. increasing soil bulk density, decreasing porosity and water content), thus lowered crop performance via stunting the aboveground growth coupled with reduced root growth [44].
Alterations in soil porosity after compaction strongly limited the air and water conductivity, so soil compaction could reduce microbial abundance and persistently alter the structure of microbiota [45]. Soil calcium plays roles in plant tolerance to the stresses like acidosis, toxic metals, osmosis and pathogen infection [46,47]. Thereby, lower calcium content in the RNI soils may aggravate nematodes infection in the tobacco roots. Iron (Fe) is an essential element, which acts as an important co-factor for the enzymes involved in regulating various metabolic processes in plants (e.g. chlorophyll synthesis), for plant growth and development [48] Alcaligenes faecalis showed toxicity against M. incognita by producing an extracellular serine protease, which can damage the intestines of nematode [37]. Enterobacter asburiae displayed promising nematicidal activity against M. incognita, and its nematicidal virulence factor was also determined as the proteolytic enzymes [58]. Rhizobacteria such as Burkholderia gladioli, Delftia acidovorans, Rhizobium etli, Pseudomonas spp. and Bacillus spp. could significantly reduce the hatching of eggs, the numbers of second stage juvenile of Meloidogyne and root galls [6,[59][60][61]. Escherichia coli secreted small molecules including indole, indole-3-carboxaldehyde and indole-3-acetic acid, as virulence factors to kill nematode [62].
incognita [63]. Many nematocidal antibiotics and metabolites (e.g. jietacin, polyketides, lactones) were produced by Streptomyces spp. to reduce root galls and inhibit egg hatching [64][65][66]. Metagenomics analysis showed these nematicidal microorganisms were enriched in the NS soils comparing with the RNI soils, and from the NS soils we got 62 bacterial isolates and two fungal isolates belonging to these genera with high nematicial activity. Bacillus spp.
and Pseudomonas spp. were more frequently isolated, indicating they may play key roles in protecting plants against nematode infection. The bacterial community richness and diversity in root-knot nematode diseased soils was significantly different from that in healthy soils [67].
This was also observed in our findings that a significant difference in the microbial community composition was found between the NS and RNI rhizospheric soils, suggesting the rhizospheric microbial community structure was shifted in the RNI soils compared to the 21 NS soils.
Higher average connectivity (i.e. average links per node in the network) means a more complex network. A smaller average geodesic distance means all the nodes in the network are closer [12,25]. We found the NS fungal and bacterial networks were more complex and their nodes were closer than the RNI networks. The increased complexity of the NS networks was reflected by the increased average connectivity, as well as the shorter geodesic distances, indicating the NS rhizosphere microbiota has a greater potential for interactions and can alter ecological networks [68,69]. Resource and food availability are important drivers of network structures [70], so the increased complexity of microbial networks was likely due to the increased amount of carbon input in soil [12]. Large number of nematode in the RNI soils (1091 nematodes per gram of soil) consumed large amounts of food and energy which may lead to consequent shortage of resource and food for microbial growth, thus decrease the complexity of microbial network. Interestingly, more positive edges were found in the NS networks than the RNI networks, indicating the potential for extensive mutualistic interactions among microorganisms assembles in the NS rhizosphere [71]. Negative co-occurrence patterns (co-exclusion) predominated in the RNI networks, suggesting that root-knot nematode infection triggered microbial competition and antagonistic interactions in the RNI soils [72].
Ten keystone microorganisms were found in the NS bacterial network. These keystone species are critical in maintaining the NS network structure and function, and ecosystem stability [27,73]. Terriglobus saanensis and Granulicella mallensis both belong to Acidobacteria and can hydrolyze diverse carbohydrates including complex polysaccharides [74]. Xylanosome produced by Oerskovia turbata exhibits powerful degradative ability against xylan and hemicellulose [75]. Streptomyces thermoautotrophicus is able to fix N2 [76]. 22 These keystone microorganisms participate in carbon and nitrogen cycles, possibly supply carbon and nitrogen resources for other microorganisms in the NS network. Occurrence of these keystone species may be indicative of the influence of soil nutrients and root exudates on microbial co-occurrence relationships in soil. The actinomycetes such as Rubrobacter radiotolerans produces two dimeric indole derivatives with acetylcholinesterase inhibitory activity [77]. The acetylcholinesterase inhibitors, e.g. oxamyl and avermectin, cause initial hyperactivity of nematode juveniles followed by a progressive decline in movement and then a final loss of activity [78]. Thereby, as keystone specie, R. radiotolerans may play important roles in control of root-knot nematode disease. The keystone taxa shifted as conditions (land uses, pathogen infection) changed [27,79]. Here, the keystone species of the RNI networks were different from those of the NS networks, indicating the keystone species of microbial network were altered in the RNI soil.
Microorganisms can kill nematodes by producing extracellular enzymes, including serine protease, metalloproteinase, collagenase, lipase and chitinase, which can degrade the nematode cuticle and eggshell [80,81]. We found the metabolites of the NS microbiota were different from those of the RNI microbiota based on the untargeted metabolomics analysis. Alternariol, a kind of mycotoxin produced by the phytopathogenic fungi Alternaria, acts as a virulence and colonization factor during their infection of plants [83]. Interestingly, alternariol was enriched in the RNI microbiota, indicating the tobacco plants were also infected by Alternaria in the RNI soils that was consistent with our investigations (data unshown). The co-infection of Alternaria and root-knot nematode may aggravate the death of plants.
Acetophenone was overrepresented in the NS samples. Two derivatives of acetophenone, 4-nitroacetophenone and 4-iodoacetophenone, were previously reported with nematicidal activity against M. incognita, with EC50 values of 12 and 15 mg/ L respectively [84]. In this study, we found that acetophenone was also very effective for killing nematodes. Not only that, its nematicidal activity (LC50 = 0.66 mg/L) was even much higher than its derivatives of 4-nitroacetophenone and 4-iodoacetophenone. Indeed, acetophenone showed a strong nematicidal activity. After treatment for 24 h, 0.66 μg/ml and 1.01 μg/ml of acetophenone could kill 50% and 90% of M. incognita J2 juveniles, respectively, but how acetophenone kill nematode is needed to be further investigated. The nematicidal activity of acetophenone was even higher than the popularly used commercial nematicide avermectin (LC50 = 2.0 mg/L) [85]. Moreover, as nematodes resistance to avermectins is on the rise, alternatives or new products with enhanced potency and high activity are urgently needed [86,87]. Considering that acetophenone has the advantages of simple structure, high efficiency and safety [88], it is a very potential nematicide candidate in the future. Besides for directly killing nematode, we found that acetophenone could strongly repel M. incognita as a repellent of root-knot nematode. Interestingly, acetophenone was also reported with strong repellency activity to some pests such as the kissing bug Rhodnius prolixus at a low dose (1%) [89], and as an anti-attractant for the pest western pine beetle, Dendroctonus brevicomis LeConte [90]. 24 Addition to acetophenone, indole-3-acetic acid (IAA) was also overrepresented in the NS samples. It is known that IAA produced by microbes can promote plant growth. Thereby, overrepresent of IAA in the NS rhizosphere soils was favorable for improving tobacco growth.
Co-occurrence network analysis found Bacillus was positively interacted with acetophenone. We isolated nematicidal Bacillus spp. from the NS rhizosphere soils. Among these isolates, B. amyloliquefaciens W1 culture could kill 98.8% of M. incognita J2 juveniles and could produce acetophenone. This study firstly found the nematicidal activity of acetophenone produced by B. amyloliquefaciens. This new discovery has broken new ground in exploration of the nematicidal metabolites from microbes through combining multi-omics analysis as well as culture-dependent technology. Interestingly, several intermediary metabolites in phenylalanine metabolism pathway were also found in the W1 culture, including 2-phenylacetamide, 2-hydroxyphenylacetic acid, phenylpyruvic acid and phenylacetaldehyde. We supposed that acetophenone may be one of the metabolites in phenylalanine metabolism, and the pathway of acetophenone production was needed to be investigated asap. Pantoea and Pseudomonas could induce systemic resistance against root-knot nematodes [54]. Here, we found that Pantoea and Pseudomonas positively interacted with Bacillus, indicating these nematicidal bacteria could cooperate together to kill nematode by multiple nematicidal mechanisms. Previous study has demonstrated that a native consortium of five native bacterial isolates protected their host plant from a fungal sudden-wilt disease [91]. Here, a consortium of Bacillus, Pantoea and Pseudomonas isolates may effectively control M. incognita via complementary traits.

Conclusion
We compared the non-infected and root-knot nematode infected tobacco rhizosphere microbiota by metagenomics and metabolomics analysis. Rhizosphere microbiota was altered accompanying with root-knot nematode infection. The RNI microbial community, microbial network, functional network and metabolites were quite different from the NS microbiota.
The NS microbial networks were more complex than the RNI networks. Nematicidal microbes were enriched in the NS rhizosphere microbiota. Potential nematicidal metabolites, antibiotics and extracellular enzymes, which are the nematicidal virulence factors, were 25 overrepresented in the NS microbiota. Consistently, nematicidal bacterial and fungal isolates were got from the NS rhizosphere soils, which could produce extracellular enzymes.
Importantly, the nematicidal B. amyloliquefaciens W1 could produce acetophenone to kill and repel root-knot nematode. Collectively, combination of multi-omics analysis and culture-dependent technology is powerful for finding nematicidal microorganisms and metabolites from soil.

Ethics approval and consent to participate
Not applicable.

Consent for publication
Not applicable.

Availability of data and material
All sequencing data has been submitted to the NCBI repository and can be accessed via the following accession numbers: whole-genome shotgun sequencing PRJNA657468, PRJNA657384, PRJNA657444.

Competing interests
The authors declare that they have no competing interests.

Funding
This work was supported by the National Natural Science Foundation of China (31870030).