The experiment was designed in the main cotton production fields of Korla (41°45’2.39”N, 85°48’4.69”E), Xinjiang province, Northwest China. The site is located in the temperate continental climate zone, with the same annual mean temperature of 11°C and annual mean precipitation of 58.6 mm. The cotton cultivars used in this study were Zhongzhimian No. 2 (resistant cultivar) and Junmian No. 1 (susceptible cultivar). The cotton plants were sampled in June 2017. Briefly, the 2,500 m² cotton field was selected as the experiment area, and 400 kg of V. dahliae was added when the cotton was planted in 2014 and 2015. The disease symptom was observed in 2014. Cotton plants in the experimental area that displayed no wilt or disease symptoms were classified as healthy plants, while the plants in the experimental area that showed wilt or obvious disease symptoms were classified as diseased plants (Supplementary Figure 1). Three replicates of healthy and diseased cotton plants were collected from three adjacent plots, and each replicate consisted of a composite sample obtained by mixing five individual samples. The leaves of the cotton plant were discarded, and the stem was cut into small parts and stored in a 50 mL tube. All the samples were kept in cold conditions and later transferred to the laboratory and stored at –80°C before further processing. The soil type in the designed field was loam sandy soil. The pH, total nitrogen, total phosphorus, total potassium concentration, and organic matter were 7.98, 0.42 g/kg, 0.67 g/kg, 16.51 g/kg, and 8.94 g/kg, respectively.

For microbial DNA extraction from the endosphere, the plant stem materials were placed in sterile tubes containing phosphate-buffered saline buffer (137 mM NaCl, 2.7 mM KCl, 10 mM Na2HPO4, and 2 mM KH2PO4) and sonicated at 50 kHz for 90 s. This procedure was repeated three times. Then, the plant stem material was rinsed with 75% ethanol for 2 min, then the sample was washed with 1% NaCIO for 1 min, and finally washed again with sterile H₂O for three times, to ensure surface sterilization. Then, the stem material was ground using sterile mortar with liquid nitrogen to extract DNA. Total DNA was extracted from the plant using Plant DNA mini kit (Omega, Beijing, China) according to the manufacturer’s protocol. The quality and quantity of DNA were measured using NanoDrop 2000 spectrophotometer (Thermo Scientific, Waltham, MA, United States). For the 16S rRNA gene library, the V3-V4 region was amplified using the universal primers 338F (ACTCCTACGGGAGGCAGCA) and 806R (GGACTACHVGGGTWTCTAAT) (Huws et al., 2007). The PCR reaction was performed in 10 μL volume, containing 0.3 μL of both primers (10 μM), 50 ng of template DNA, 5 μL of KOD FX Neo buffer, and 0.2 μL of KOD FX Neo (TOYOBO, Shanghai, China). Amplification was carried out with the following PCR program: initial denaturation at 95°C for 5 min, followed by 25 cycles of 95°C for 30 s, 50°C for 30 s, and 72°C for 40 s, and final elongation at 72°C for 7 min. The library quality was assessed using the Qubit 2.0 Fluorometer (Thermo Scientific) and Agilent Bioanalyzer 2100 system. Finally, the library was sequenced on an Illumina Hiseq 2500 platform at Biomarker Co., Ltd (Beijing, China), and 250 bp paired-end reads were generated.

Raw sequences were split according to their unique barcodes, and the adaptors and primer sequences were trimmed using their own script from the company. Paired-end sequences were merged using FLASH software 1.2.11 with the default parameters (Magoc and Salzberg, 2011). Then, the amplicon data were analyzed using a combination of the VSEARCH v2.13.3 and QIIME 1.9.1. software (Caporaso et al., 2010; Zgadzaj et al., 2016). The resultant sequences were quality-filtered, and singletons were removed in VSEARCH. The sequence reads were then clustered into OTUs at 97% similarity level. Chimeric sequences that were identified using the reference-based chimera checking methods were removed from the data. Representative sequences were classified using the blast algorithm with the SILVA (v.13.2) reference database in QIIME (Quast et al., 2013). Mitochondrial and chloroplastic DNA sequences, and OTUs with a total relative abundance of < 0.00001 in all of the samples were discarded. The raw sequencing data have been submitted to the Sequence Read Archive under accession number PRJNA820148.

The sample was prepared as previously reported, with little modification (Wen et al., 2020). For extraction, 50 mg of sample was taken and placed in an EP tube, and then 1,000 μL of extraction solvent containing an internal target was added (V_methanol: V_acetonitrile: V_water = 2:2:1, which was kept at –20°C before extraction). The mixtures were homogenized in a ball mill for 4 min at 45 Hz, and then ultrasound treated for 5 min (incubated in ice water). After homogenization for three times, the sample was incubated for 1 h at –20°C to precipitate the proteins, and then centrifuged at 4°C and 12,000 rpm for 15 min. Next, the supernatant (500 μL) was transferred into a fresh tube, the extracts were dried in a vacuum concentrator without heating, and 100 μL of extraction solvent was added (V_acetonitrile: V_water = 1:1) for reconstitution. The sample was vortexed for 30 s, sonicated for 10 min (4°C water bath), and centrifuged for 15 min at 12,000 rpm and 4°C. Then, the supernatant (60 μL) was transferred to a fresh 2 mL LC/MS glass vial for the UHPLC-QTOF-MS analysis.

The UHPLC-QTOF-MS analysis was performed using a UHPLC system (1290, Agilent Technologies, Palo Alto, CA, United States) with a UPLC BEH Amide column (1.72 μm, 2.1*100 mm, Waters, Milford, CT, United States) coupled to TripleTOF 6600 (QTOF, AB Sciex, Framingham, MA, United States). The mobile phase consisted of 25 mM NH₄OAc and 25 mM NH₄OH in water (pH = 9.75) (A) and acetonitrile (B), and elution was carried out with elution gradient as follows: 0 min, 95% B; 7 min, 65% B; 9 min, 40% B; 9.1 min, 95% B; and 12 min, 95% B, which were delivered at 0.5 mL min^–1. The injection volume was 2 μL. The Triple TOF mass spectrometer was used due to its ability to acquire MS/MS spectra on an information-dependent basis (IDA) during an LC/MS experiment. In this mode, the acquisition software (Analyst TF 1.7, AB Sciex, Framingham, MA, United States) continuously evaluates the full-scan survey MS data as it collects and triggers the acquisition of MS/MS spectra depending on the preselected criteria. In each cycle, 12 precursor ions whose intensity was greater than 100 were chosen for fragmentation at the collision energy (CE) of 30 V (15 MS/MS events with a product ion accumulation time of 50 ms each). The following settings were used: ion source gas 1 at 60 Psi, ion source gas 2 as 60 Psi, curtain gas at 35 Psi, source temperature at 650°C, and ion spray voltage floating (ISVF) at 5,000 V or –4,000 V in positive or negative modes, respectively.

MS raw data (.wiff) files were converted to the mzXML format using ProteoWizard and processed by R package XCMS (version 3.2). The preprocessing results generated a data matrix that consisted of the retention time (RT), mass-to-charge ratio (m/z) values, and peak intensity. R package CAMERA was used for peak annotation after data processing by XCMS. An in-house MS2 database was applied for the identification of metabolites.

The OTU table was rarefied to reads (lowest number of 3,616 for bacterial community) for alpha-diversity and was used as a normalization method for beta-diversity analysis. Bacterial alpha-diversity (observed OTUs and Shannon index) values were calculated using QIIME. The significant differences were evaluated with a one-way analysis of variance (ANOVA). Bacterial beta-diversity was estimated according to the Bray-Curtis distance between samples. Permutational multivariate analysis of variance (PERMANOVA) was performed with 1,000 permutations using the adonis function in vegan R package (Xiong et al., 2021). The Venn diagram was drawn using the OECloud tools¹ to analyze common and specific OTUs for healthy and disease plants of different cultivars. The linear discriminant analysis effect size (LEfSe) was applied (Wilcoxon P-value < 0.05 and logarithmic LDA (linear discriminant analysis) score > 2)² to identify the biomarkers of healthy and diseased plants and different cultivars (Xiong et al., 2021). Welch’s test was used to calculate the significance of differences between different healthy conditions using STAMP software (Parks et al., 2014). Differential abundance analysis was performed using EdgeR’s generalized linear model (GLM) approach, and the differential OTUs with false discovery rate-corrected P-values < 0.05 were identified as indicator OTUs (Robinson et al., 2010).

The network analysis based on Pearson correlation scores was performed using the CoNet in CYTOSCAPE v3.8 and visualized in GEPHI (Shannon et al., 2003; Xiong et al., 2021). Briefly, only robust (Pearson’s r > 0.9 and r < –0.9) and statistically significant (p < 0.01) correlations were retained for the next step. The ‘benjaminihochberg’ multiple testing correction (p < 0.01) was performed at bootstrap step restore network step. The topological indices of the networks were calculated in Gephi. Network modular analysis was conducted in Gephi with default parameters. The network complexity was defined according to previous research, and nodes with high degree (>30) and closeness centrality (>0.6) values were identified as hub nodes in the network (Agler et al., 2016; van der Heijden and Hartmann, 2016; Wagg et al., 2019). Analysis was conducted in R 3.0.6 or GraphPad Prism 7.

To investigate the metabolism of different cultivars, if the metabolite feature is detected in <20% of experimental samples or is detected in <50% of quality control samples, then it is deleted from the data analysis. Then the missing values of raw data were filled up by half of the minimum value. In addition, the internal standard normalization method was employed in this data analysis. Finally, features with a relative standard deviation > 30% should be removed from the subsequent analysis. The resulting three-dimensional data involving the peak number, sample names, and normalized peak area were fed into R package metaX for principal component analysis (PCA) (Wen et al., 2017). PCA showed the distribution of origin data. The first principal component of variable importance in the projection (VIP) was obtained. The VIP values summarize the contribution of each variable to the model. The metabolites with VIP > 1 and p < 0.05 (student t-test) were considered as significantly altered metabolites. In addition, commercial databases including KEGG³ and MetaboAnalyst⁴ were utilized to search the pathways of metabolites.

The fresh cotton plants were sampled for the isolation of endophytic bacteria. The method used to isolate the endophytic bacteria was reported previously (Zhang et al., 2017). First, the surface of the cotton plants was sterilized with 75% ethanol for 30 s, immersed in 2% sodium hypochlorite for 3 min, and rinsed four times with sterile distilled water. Then, 3 g of surface-sterilized cotton plants was collected and ground in 5 mL of sterilized water. The resultant suspension was plated on Luria-Bertani (LB, 10 g/L tryptone, 5 g/L yeast extract, and 10 g/L NaCl) and King’s B (KB, 20 g/L protease peptone No. 3, 1.5 g/L K₂HPO₄, 1.5 g/L MgSO₄⋅7H₂O, and 10% glycerol) plate after 10 × dilution. The plates were then incubated at 30°C for 2 days. Colonies of different morphology were transferred onto new LB and KB plates for further characterization. All isolates were purified and stored at –80°C in a nutrient broth containing 15% glycerol. DNA extraction, amplification, and sequencing were performed as described by a previous study (Zhao H. et al., 2019; Zhuang et al., 2021). The antagonistic activity of the isolates against cotton Verticillium wilt was evaluated using the methods reported by Liu et al. (2021b). Briefly, spores of V. dahliae were diluted to 8 × 10⁵ mL^–1 and 200 μL was inoculated on plates containing potato dextrose agar (PDA, 200 g/L potato and 20 g/L glucose). Then, the plates were permitted to dry under a laminar flow hood. Portions of cultures (3 μL) of strains were added to the center of plates, and the zone of inhibition was measured after incubation at 30°C for 4 days.

The greenhouse experiments also evaluated the capacity of Pseudomonas sp. to control cotton wilt disease. First, cotton seeds were sterilized in 1.5% NaOCl for 10 min with gentle shaking and washed four times with ddH₂O. Subsequently, the sterilized seeds soaked in the bacterial suspension (28°C, 200 rpm for 12 h) at room temperature for 2 h and seeds soaked in LB or water were used as the control samples. Then, one seed was sown in each pot (the ratio of sterilized nutrient soil to vermiculite is 3:1). Six pots were used for each replicate, and three replicates were used for each isolated bacterial strain. The cotton plants were placed in the growth chamber at 28°C with 12-h light/12-h dark conditions. When the second euphylla appeared, 10 mL spores of V. dahliae (1 × 10⁸mL^–1) were added to the soil of the cotton plant, and this process was repeated every 2 days. In total, 30 mL spores of V. dahliae were added. The disease index of cotton plants was calculated after 40 days. Cotton plants were checked for severity based on the disease index described previously (Gu et al., 2021). Briefly, we rated the disease severity on a scale from 0 to 6 according to the percentage of symptoms. The score of 0 indicated no visible symptoms, and score of 1, 2, 3, 4, 5, or 6 described 0–10%, 11–25%, 26–50%, 51-75%, 76–90%, and >90%, discoloration or defoliation, respectively. The disease index was calculated as follows:

Disease index = [Σ(Rating × Number of diseased plants)/Total number of plants × Highest rating] × 100

The endophytic bacterial community of cotton plants was characterized by Illumina HiSeq sequencing. In total, 803,958 high-quality sequences were obtained, and each sample contained between 65,586 and 68,096 (66,997 ± 643) reads. All the sequences were clustered into 700 operational taxonomic units (OTUs) with > 97% similarity. Meantime, we also checked the disease index of different cultivars, and the results showed that the disease index of susceptible cultivars was significantly higher than that of the resistant cultivars in the experimental area (unpaired t-test, P < 0.001) (Supplementary Figure 2).

Alpha-diversity of healthy and diseased cotton plants was measured by Shannon and observed OTU indices. It showed that the alpha-diversity values of healthy plants were higher than those of the diseased plants (Figure 1A). Meanwhile, we also compared the alpha-diversity of different cultivars and found that the alpha-diversity values of susceptible cultivars were higher than those of the resistant cultivars (Figure 1B). Subsequently, the samples were divided by cultivars and analyzed separately. The results displayed that the values of healthy plants were higher than those of the diseased plants for resistant cultivars. However, the index of diseased plants was higher than that of the healthy plants for susceptible plants, but the difference was not found to be significant (Supplementary Figure 3).

To investigate how cotton Verticillium wilt affects the co-occurrence patterns of the cotton endophytic community, we analyzed the bacterial community network of healthy and diseased plants. Based on the network analysis, we found that the healthy conditions showed an effect on network complexity. The average degree was used to assess the network complexity and showed that the endophytic bacterial community of the diseased plants (8.03 and 4.95) was more complex than that of the healthy plants (5.55 and 2.89, Figure 2). Additionally, the results also displayed that the susceptible cultivars possessed a higher degree than the resistant cultivars in healthy and diseased conditions. Furthermore, the number of “hub nodes” [nodes with high values of degree (>30) and closeness centrality (>0.6) in the network] from the diseased samples was higher than that of the healthy plants (Table 1). The taxonomic composition of the networks differed between healthy and diseased cotton plants, with most nodes annotated to Clostridia for diseased plants and Gammaproteobacteria for healthy plants of susceptible cultivars (Figure 2). Moreover, we recorded a higher ratio of positive to negative correlations in diseased cotton plants for susceptible and resistant cultivars. In addition, a higher average path distance was found in the healthy cotton plants (Table 1).

The principal coordinate analysis (PCoA) plot of microbial communities revealed a clear separation between the cotton plants of different cultivars. The two main coordinates extracted explained 41.62% of the variation, of which PC1 explained 23.63% of the variation and PC2 explained 17.99% of the variation (Figure 3A). Moreover, the cultivars had a significant influence on the endophytic bacterial community (R² = 14.98%, P < 0.05). Furthermore, hierarchical clustering analysis revealed a separate clustering between healthy and diseased plants (Supplementary Figure 4). The results also showed that the healthy conditions had a strong effect on the bacterial endophytic community.

The bacterial composition at the phylum level was analyzed in the cotton samples. We found that the most abundant bacterial phyla were Alphaproteobacteria (45.4%), Gammaproteobacteria (22.3%), Actinobacteria (21.8%), Bacteroidetes (5.5%), and Firmicutes (3.7%) (Figure 3C). In addition, taxonomic classification showed that Rhizobiaceae was the most abundant bacteria in healthy and diseased cotton plants at the family level. Furthermore, we found that the population of Rhizobiaceae was significantly increased in the healthy cotton plants at the family level, whereas Streptococcaceae was significantly decreased in the healthy cotton plants compared with diseased cotton plants (Figure 3D). Meanwhile, the LEfSe results showed that Rhizobiaceae and Patescibacteria were the biomarker taxa for healthy and diseased cotton plants (Figure 3E).

To determine which OTUs are responsible for the healthy conditions and cultivar effect, we performed OTU analysis with OTUs identified from healthy and diseased plants, and different cultivars (Figure 3B). In the resistant cultivars, we identified 8 and 25 OTUs that were specifically observed in healthy and diseased plants, respectively. The OTUs that were specifically noticed in healthy plants primarily correspond to Clostridia, Bacteroidia, Actinobacteria, and Verrucomicrobiae. In contrast, the diseased plants specifically showed the presence of Clostridia, Alphaproteobacteria, Gammaproteobacteria, Bacteroidia, Deinococci, and Negativicutes at the class level. For the susceptible cultivars, nine OTUs consisting of Clostridia, Bacteroidia, Gammaproteobacteria, Actinobacteria, and Verrucomicrobiae were specifically selected by healthy plants, and 35 OTUs belonging to Clostridia, Bacteroidia, Alphaproteobacteria, Deltaproteobacteria, Gammaproteobacteria, Acidimicrobiia, Nitrospira, Verrucomicrobiae, Nitriliruptoria, and Melainabacteria were specifically selected by diseased plants (Figure 3B and Supplementary Table 1). Furthermore, the Venn diagram analysis of all samples indicated that the diseased plants (84) had higher specific OTUs than the healthy plants (24). The OTUs selected by healthy plants annotated as Ruminococcaceae, Opitutaceae, Sphingobacteriaceae, Weeksellaceae, Microscillaceae, and Nocardioidaceae, while the OTUs selected by diseased plants belonged to Ruminococcaceae, Lachnospiraceae, Muribaculaceae, Burkholderiaceae, and Chitinophagaceae (Supplementary Figure 5A). In addition, we found that the susceptible cultivars (81) possessed more specific OTUs than the resistant cultivars (49). The 81 OTUs mainly corresponded to Ruminococcaceae, Muribaculaceae, Desulfovibrionaceae, Flavobacteriaceae, Burkholderiaceae, Xanthomonadaceae, Sphingobacteriaceae, and Chitinophagaceae, and the 49 OTUs mainly corresponded to Ruminococcaceae, Lachnospiraceae, Sphingobacteriaceae, and Weeksellaceae at the family level (Supplementary Figure 5B).

Remarkably, the different abundance analysis indicated that two OTUs which corresponded to Xanthomonadaceae and Saccharimonadaceae at the family level were significantly enriched in the diseased plants. In contrast, one OTU, belonging to Pseudomonas, was significantly depleted in the diseased plants. However, there were no enriched or depleted OTUs in the susceptible cultivars and all the samples (Supplementary Figures 5, 6). For the different cultivars, we found that three and six OTUs were significantly enriched in resistant and susceptible cultivars, respectively. In the resistant cultivars, the OTUs belonged to Burkholderiaceae, Moraxellaceae, and Opitutaceae at the family level, and in the susceptible cultivars, the OTUs corresponded to Pseudomonadaceae, Xanthomonadaceae, Cellvibrionaceae, and Ruminococcaceae (Supplementary Figure 5D).

The root exudates of resistant and susceptible cotton cultivars were analyzed by UHPLC-QE-MS after collecting from cotton root samples. A total of 3,777 peaks were detected, and 212 compounds were identified across all the samples, including lipids and lipid-like molecules, benzenoids, organic acids and derivatives, organic oxygen compounds, organoheterocyclic compounds, phenylpropanoids and polyketides, and nucleosides, nucleotides, and analogs at the superclass level. The overall exudation patterns of the resistant and susceptible cotton plants were found to be distinct. The first two components of the PCA explained 93.9% of the total variance (79.4 and 14.5% for PC1 and PC2, respectively), and this analysis clearly separated the susceptible and resistant samples (Figure 4A). Furthermore, the data revealed the presence of multiple metabolites that showed significant differences in abundance between susceptible and resistant cultivars according to the VIP thresholds (VIP > 1) and the student’s t-test (P < 0.05). A total of 76 compounds were found to be significantly different between these two cultivars. Specifically, 18 compounds showed a higher relative abundance in the resistant cultivars, and 58 compounds showed a higher relative abundance in the susceptible cultivars (Figure 4B). To further reveal the relative content of all the annotated and significantly differentially expressed metabolites in the two cultivars and determine their relationships, cluster analysis was conducted, and the results were presented as a heatmap (Figure 4C). The susceptible cultivars had more number of differentially expressed metabolites. Then, the pathway analysis based on the significantly different metabolites was carried out using the KEGG enrichment analysis. The results displayed that pathways related to plant hormone signal transduction (salicylic acid), biosynthesis of various secondary metabolites, biosynthesis and metabolism of amino acids, including tryptophan, alanine, aspartate, glutamate, cysteine, and methionine, and zeatin biosynthesis (Figure 4D) were significantly altered.

To further characterize the antagonist effect of Pseudomonas sp. on the cotton Verticillium wilt, we isolated 30 candidate bacteria from the healthy cotton plants. In total, two bacterial isolates were obtained and phylogenetically characterized based on the distinct 16S rDNA gene sequence, allowing us to recover strains for further studies. The results indicated that these two strains belonged to the genus Pseudomonas and inhibited the growth of cotton Verticillium wilt on the PDA plate, and the diameters of the inhibition zone were 29.67 ± 2.08 and 31.33 ± 1.53 for strains N1 and L4, respectively (Figures 5A,B). Furthermore, we were curious to know if the Pseudomonas sp. played a role in the control of the cotton wilt under greenhouse conditions. The results displayed that the disease index of cotton plants treated with the strains L4 and N1 was significantly lower than that of the control plants. There were no significant differences in the disease index between control and LB-treated plants (Figures 5C,D). Importantly, all inoculation treatments could promote plant growth. The height of plants treated with Pseudomonas sp. strains N1 and L4 was significantly higher than that of the control plants (Figure 5E). The results displayed that Pseudomonas sp. promoted the growth of cotton plants and protected the plants by inhibiting the growth of V. dahlia.

The microbiome of the plant plays a crucial role in plant health. Microbial diversity and distinct changes in the microbial community were previously shown to be related to the disease incidence (Berg et al., 2017). In the current study, we found that the alpha-diversity of the diseased plants was higher than that of the healthy plants for the endophytic bacterial microbiome (Figure 1). The previous studies found that the bacterial diversities of infected root and stem samples were significantly higher than those of the healthy plants. The pathogen disrupted the plant defense system, which allowed entry to more number of organisms into the plant (Hu et al., 2020). Moreover, the results also showed that the infected samples had lower bacterial diversity compared to the healthy samples based on the healthy and decaying sugar beets (Kusstatscher et al., 2019; Liu et al., 2019). The decrease in the microbial diversity of plants was probably related to the invasion by the pathogenic species into the communities. The plant microbiome was found to be affected by many factors, such as soil properties, plant genotype, plant developmental stage, and background microbial composition of the soil (Qiao et al., 2017). We also found that the different cotton cultivars affected the composition of the bacterial microbiome (Figure 3A). Furthermore, the plant developmental stage had an impact on the microbial community. The richness indices of the bacterial species in the shoot and roots increased with the growth of the banana plant (Liu et al., 2019). The endosphere microbiome of cotton during plant growth and disease development should be further studied.

In addition, we found that the alpha-diversity in the susceptible cultivars was higher than that of the resistant cultivars regardless of the healthy conditions in the bacterial community, which is consistent with the findings in a previous study, which showed that cucumbers susceptible to Fusarium oxysporum f. sp. cucumerinum exhibited higher bacterial community evenness than the resistant variety of cucumber (Supplementary Figure 3, Wen et al., 2020). High diversity was observed in the susceptible cultivars, which could be attributed to a relatively high variety of exudates supporting more microbial niches (Wen et al., 2020). We also found that the susceptible cultivars had more different abundant compounds compared with the resistant cultivars (Figure 4). Furthermore, we displayed that the network complexity of diseased plants was more complex than that of healthy plants, and the results were consistent with the previous report (Hu et al., 2020, Figure 2). The pathogenic produce positive feedback that enhances their competitiveness and subsequent interactions with the neighbors by changing the composition. Importantly, the highly connected and modular microbiota could prime the plant immune system for accelerated activation of defense against the pathogen (Hu et al., 2020).

Plants can alter their microbiome composition upon pathogen infection and recruit some disease-resistance and growth-promoting beneficial microbes for their survival. For instance, the relative abundances of Pseudomonas, Bacillus, and Falsibacillus, which are often considered to be plant-beneficial bacteria, showed a significant decrease in the diseased plants when compared to the healthy plants. These genera colonize different plant compartments and play a vital role in modulating host performance, particularly in the suppression of plant pathogens (Berendsen et al., 2018; Hu et al., 2020). In the present study, the comparative analysis of the composition of the healthy and diseased plants showed that the healthy plants were enriched in Rhizobiaceae, which is also the biomarker taxa of healthy plants, and diseased plants were enriched in Streptococcaceae (Figures 3D,E). Moreover, we found that the healthy plants were enriched in one OTU, annotated as Pseudomonas, and the diseased plants were enriched in two OTUs, classified as Xanthomonadaceae and Saccharimonadaceae, of resistant cultivars (Supplementary Figure 6). The enriched OTUs of healthy plants probably play an important role in protecting the plants from cotton Verticillium wilt. This finding is consistent with the previous report (Yuan et al., 2020; Zhang et al., 2021). Furthermore, the results displayed that the healthy plants and resistant cultivars contained specific OTUs corresponding to Sphingobacteriaceae, which included the biocontrol genus of Sphingobacterium (Xu et al., 2020). We also need to confirm the capable of suppression cotton Verticillium wilt. However, recent studies have reported that several potential beneficial bacteria, such as Pseudomonas, Streptomyces, and Bacillus, were enriched in diseased plants (Gao M. et al., 2021). The opposite results observed in different plant species need further exploration.

Deciphering the composition of the healthy and diseased plants is critical for harnessing the plant microbiome to enhance plant growth and health. The previous studies showed that Acidobacteria, Actinobacteria, Bacteroidetes, Planctomycetes, Proteobacteria, and Verrucomicrobia were the major taxa in the cotton rhizosphere (Qiao et al., 2017), while we found that Proteobacteria, Actinobacteria, Bacteroidetes, and Firmicutes were the main taxa in the cotton endophytic bacterial community (Figure 3C). There are some differences between rhizosphere and endophytic communities, and the crop microbiome is mainly derived from soil and is gradually enriched and filtered at different plant compartment niches (Xiong et al., 2021). The soil microorganisms likely affect the pathogen migration, and the interaction between plant endophytes and soil microorganisms facilitates the pathogen invasion process. It was found that pathogen invasion may begin in the bulk soils, then transfer to the plant roots, and further infect the plant stems, causing bacterial wilt by Ralstonia (Hu et al., 2020). Under stress conditions, a plant could attract distant beneficial microbes by actively releasing root exudates, such as amino acids, nucleotides, and long-chain organic acids. Based on the metabolism results, we found that the susceptible cultivars had more differentially expressed metabolic pathways (Figure 4B). Pathway enrichment analysis showed that differentially expressed metabolic pathways were related to plant hormone signal transduction, biosynthesis of various secondary metabolites, and biosynthesis and metabolism of amino acids (Figure 4D). Importantly, the amino acids could induce plant growth through rhizosphere bacteria (Zhao et al., 2019). The relationship between rhizosphere and root exudates needs further study.

Pseudomonas sp. are widely used as biocontrol agents to combat plant diseases. In the previous study, Pseudomonas spp. could directly inhibit the growth of cotton Verticillium wilt under laboratory and greenhouse conditions. In the meantime, the abundances of Pseudomonas in the soil microbial community with less incidence had higher abundances change are higher (Erdogan and Benlioglu, 2010; Zhang et al., 2022). The Pseudomonas was the key plant growth-promoting rhizobacteria in the rhizosphere of garlic plants and had the ability to promote plant growth (Zhuang et al., 2021). The genera of Bacillus and Pseudomonas are the most dominant biocontrol strains in the market, and these two genera could coexist and cooperate with each other. The bio-organic fertilizers stimulate the increase of the Pseudomonas spp. in the soil microbiome to enhance plant disease suppression, and B. velezensis SQR9 induced the enrichment of Pseudomonas sp. in the rhizosphere and promoted the growth of cucumber plants synergistically (Tao C. et al., 2020; Yin et al., 2021; Sun et al., 2022). Moreover, the cotton rhizosphere soil was temporarily altered after applying B. axarquiensis strain. Among them, Streptococcus and Pontibacter were decreased, and Lysobacter and Bacillus were significantly increased (Gao C. et al., 2021).

In the present study, we found that two Pseudomonas strains could inhibit the cotton Verticillium wilt under the laboratory and greenhouse conditions (Figure 5). However, it has been reported that complex microbial community assemblages rather than single microbial strains protected plants from pathogen infection based on multiomics sequencing technologies. More number of studies revealed that the synthetic communities could be applied to achieve sustainable agriculture. The seven microbiome strains together showed better protection of maize roots against Fusarium verticillioides infection than the individual strains (Niu et al., 2017). The members of Chitinophagaceae and Flavobacteriaceae were enriched in the plant endosphere following pathogen invasion, and reconstruction of a synthetic community of Flavobacterium and Chitinophaga consistently suppresses fungal root disease against sugar beet Rhizoctonia damping-off disease (Carrion et al., 2019). Importantly, the enriched bacteria acted as an early warning agent by activating the salicylic acid and jasmonic acid pathways in the presence of the pathogen (Liu et al., 2021a). Also, the results revealed that beneficial bacteria showing high abundance and low abundance performed different functions. The high-abundance bacteria mainly protected the plants by promoting plant growth and inhibiting pathogen invasion, while the low-abundance bacteria controlled the diseases by inducing resistant systems (Ritpitakphong et al., 2016; Jack and Nelson, 2018; Li et al., 2021). Pseudomonas spp. could colonize the roots of plants and produce antibiotics to control the disease (Tao X. et al., 2020). Moreover, the endophyte P. fluorescens PICF7 could trigger defense responses and compete in the niche with the pathogen in roots (Prieto et al., 2009; Gómez-Lama Cabanás et al., 2014). In this study, it displayed strong inhibition of cotton Verticillium wilt when the cotton plant was treated with isolated Pseudomonas strains, and it also promoted the plant growth (Figure 5). We further need to confirm the mechanism by which isolated Pseudomonas strains prevent the cotton Verticillium wilt.