Responses of Soil Microbial Communities to the Invasion of Fusarium Oxysporum and Screening for Fusarium Antagonistic Microorganisms From Pathogenic Soil

2.1 Fugal strains and plant material

The fungal pathogens Fusarium oxysporum f. sp. momordicae (Fom) strain SG-15 (Agricultural Culture Collection of China with accession ACCC 39204) was isolated from diseased plant of bitter gourd wilt in Shouguang, Shandong Province in P. R. China. Other fungal pathogens besides Fom, such as Fusarium graminearum, F. moniliforme, F. oxysporum f. sp. cucumarinum, F. oxysporum f. sp. niveum, Pyriculariaoryzae, Trichothecium roseum, Cryphonectria parasitica, Cytospora chrysosperma, Phytophthora parasitica, Rhizoctonia solani, used in this study were preserved in our laboratory. The bitter gourd cultivar is Ruyu-41.

2.2 Study sites and soil sample collection

We collected soil samples (the standing plant was bitter gourd) at a long-term experimental station of the Shandong Agricultural University, in Taian of Shandong Province, China, August to November 2019. The artificial disease nursery of Fusarium wilt was established in the experimental station in 2015 to investigate the disastrous mechanisms and control strategies of exogenetic invasive F. oxysporum. The incidence of Fusarium wilt of bitter gourd is 100% from 2017 to 2019. Three treatments with five replicates were set up, including healthy soil without Fusarium inoculation (Control), pathogenic soil with Fusarium inoculation (S-FOM), and the pathogenic soil with Fusarium inoculation was treated with antagonistic fungus (TP). We collected 1 kg soil samples from the rhizosphere of bitter gourd in each plot. The collected soil samples were sieved (2 mm). A part of the soil for isolation and screening of beneficial microorganisms, a part of the soil for DNA extraction was stored at -80 ℃, and the rest soil for soil fertility was stored at 4 ℃. The soil physicochemical characterizations were measured based on the norm system, the soil enzyme activities (i.e. soil fluorescein diacetate hydrolase, soil urease, soil polyphenol oxidase and soil acid phosphatase) were measured according to instructions of the kit (Beijing Solarbio Science & Technology Co., Ltd.).

2.3 Isolation of biocontrol agent from the pathogenic soil

Approximately 10 g soil was incubated in 100 mL of sterile distilled water and constantly shaked at 200 rpm for 30 min. 0.2 mL diluted soil extract was plated on PDA medium Petri plates. Plates were then incubated at 28℃ for 3–5 days and the growing colonies were transferred to new Petri dishes with PDA medium in order to establish bacterial or fungal pure cultures. More than 100 colonies were collected and used in subsequent screenings. All microbial collections isolated from the pathogenic soil were tested for antagonism against Fom. Antagonistic examination of microbial isolates against Fom was performed in vitro using PDA plates. Agar disks of Fom was applied on the Petri dish 6 cm apart from each microbial culture spot. The Petri dishes were then incubated at 28℃ for 7 days. The percentage of mycelia growth inhibition I (%) by isolated strain was calculated according to the following formula: I (%) = (1 − Dn/Do) × 100, where I (%) represents the mean inhibition of mycelia growth, Dn represent the mean diameter of fungal pathogens growth in the presence of isolated strain, and Do represents the mean diameter of fungal pathogens growth in the absence of isolated strain (Control). 14 microbial isolates showing high Fom inhibition were named SK1-SK6, KG-1, Q1-Q2, and M1-M5 (Supplementary Table S1).

2.4 Morphological and molecular identification of Talaromyces sp. strain Q2

Talaromyces sp. strain Q2 was inoculated and grown on Petri dishes (90 mm) with potato dextrose agar (PDA) medium at 25°C in the dark, which were performed in triplicate. 2–5 days later, the morphology observation of conidiophores and conidia was performed under Nikon Eclipse 90i microscope (Nikon, Tokyo, Japan). The pure strains were incubated on different media including Czapek agar (CA), Czapek yeast extract agar (CYA), Sabouraud maltose agar (SMA), Corn meat agar (CMA), Martin media (MM), Sabouraud sucrose agar (SSA), Sabouraud dextrose agar (SDA), Oatmeal agar (OA) at 25°C for 7 days, and the colony morphology and microscopic characters were examined, measured and photographed. Morphological identification of Talaromyces sp. strain Q2 was conducted according to the references of Visagie and Su [19, 20]. β-1, 3-glucanase and chitinase produced by strains Q2 were determined qualitatively by the clear zone method on β-D-glucan and colloidal chitins-containing agar plate. At the same time, we also tested the changes of the enzyme-producing ability of strain Q2 induced by the different substrate. The β-1, 3-glucanase and chitinase were assessed quantitatively by chitinase enzyme Activity Determination Kit and β-1, 3-glucanase enzyme Activity Determination Kit (Beijing Solarbio Biotechnology Co., Ltd.), respectively.

Fungal mycelium was harvested from the surface of PDA medium and frozen by liquid nitrogen. Then the DNA samples of strain Q2 were extracted using the Fungal Genomic DNA Kit (BioTeKe Biotech Co. Ltd, Beijing, China) according to the manufacturer’s instructions and stores at -20˚C. Primers of ITS1/ITS4, Bt2a/Bt2b, and cmdAD1/cmdQ1 were used to identify the strain Q2 [21, 22, 23]. The primers, PCR program, and the reference sequence used in the current study have been described in the Supplementary Table S2. The genomes of strain Q2 was sequenced using the whole genome de novo sequencing at the Illumina NovSeq platform and Pacbio Sequel platform.

2.5 The solid fermentation of strain Q2

On the basis of our previous research, we selected wheat grain, corn straw, and humic acid as the main material of solid fermentation medium. Orthogonal test table L₉(3⁴) was used for studying the optimal solid-state fermentation media formula of strain Q2. The optimal solid-state fermentation media formula of strain Q2 was 66.67% wheat grain, 14.29% corn straw, 19.05% humic acid, and 65.0% moisture and the number of spores can reach 6 ×10⁹ cfu·g^− 1 (Supplementary Table S3). After drying, culture material was crushed and adjusted spores concentration containing 2×10⁷ cfu·g^− 1.

2.6 Antagonism assay in vitro and in vivo

The ability of strain Q2 to inhibit the growth of pathogenic fungi listed in Sect. 2.1 has been investigated. Antagonism of strain Q2 against pathogenic fungi was conducted in vitro using PDA plates as described in Sect. 2.3.

For the study on the control efficacy of strain Q2 in greenhouse and field, healthy bitter gourd seedlings with 2–4 leaves were selected and transplanted to health soil and pathogenic soil containing fungal pathogen inocula of Fusarium wilt, the preparation of soil infested with the pathogenic fungus was conducted according to Wang [24]. 15 g of microbial agents of strain Q2 were poured for each plant in health soil (H-TP) and pathogenic soil (TP). And the equal amounts of solid fermentation medium without strain Q2 was poured for each plant health soil (Control) and pathogenic soil (S-FOM), all the seedlings were cultured at the same condition. After transplantation for 30 days, diseased incidence and disease index of bitter gourd seedlings were investigated. Control efficacy of Fusarium wilt of bitter gourd was calculated. Six classes were divided according to Chen [25]. The density of Fusarium spp. in the bitter gourd rhizosphere was estimated by metagenomic sequencing and dilution plating (using Fusarium selective medium Komada’s medium[26].

2.7 The mechanism of strain Q2 inhibiting the growth of FOM

In order to verify the mechanism of strain Q2 inhibiting the growth of Fom. Two treatments with three replicates were set up, including the sample singly-cultured of Fom (Control treatment) and the sample co-cultured with strain Q2 and Fom (T). Spores of strain Q2 and Fom were washed with sterile water, the final concentration of 4×10⁷cfu·ml^− 1 and 2.5×10⁷ cfu·ml^− 1, respectively. Subsequently, 0.8 mL of conidial suspensions of Fom combined with 0 ml (Control), 0.05 ml (T1), 0.5 ml (T2), 1.0 ml (T3), 2.5 ml (T4), and 5.0 ml (T5) of conidial suspensions of Q2 in multiple proportions were inoculated into 100 mL PDB, co-cultured at 180 rpm·min^− 1 and 28°C for 3–5 days. The growth of Fom in different treatments was recorded by measuring the spore number of Fom. The axenic culture was filtered through three layers of gauze, and then filtrates were centrifuged at 8000 rpm·min^− 1 for 10 min, and supernatant was discarded. Spores were precipitated and then were used to analyze the effects of strain Q2 on the metabolic activity of Fom (i.e. N-acetylglucosamine content of fermentation supernatant liquid, intracellular glycerol content and malondialdehyde content, they were measured according to instructions of the ELISA kit (Jiangsu Jingmei Bioechnology Co., Ltd.)). Meanwhile, we had chosen the sample singly-cultured of Fom (Control treatment) and the sample co-cultured with strain Q2 and Fom (T3 treatment group) for transcriptome sequencing, there were named FOM and Q2-FOM, respectively.

2.8 Identification of the role of a CFEM domain-containing protein orthologue gene from Fom

Generation of FomCFEM deletion mutants FomCFEM deletion mutants were generated by the split-marker approach. Briefly, a 760-bp upstream sequence and a 1027-bp downstream sequence fanking the FomCFEM gene were amplified with primer pairs CFEM-AF/R and CFEM-BF/R, respectively. The hygromycin phosphotransferase (hph) gene was amplified from pCB1003 plasmid with primer pairs HYG-F/HY-R and YG-F/HYG-R. Two equal quantities of the fragments (a total of 20 µg) were transformed into protoplasts of Fom, and the resulting transformants selected on TB3 medium containing 100 µg/mL hygromycin B were identified by polymerase chain reaction (PCR) with primers (CFEM-F/R, CFEM-K1F/R, and CFEM-K2F/R) and southern blot analysis. For southern blot assay, genomic DNA was extracted from mycelia of each strain as described previouslyand digested with Pvu I [27, 28]. The DNA probe was amplified by PCR using the primer pairs tCFEM-F/R. Southern blot analysis was performed with a DIG High Prime DNA Labeling and Detection Starter kit I (Roche Diagnostics, Mannheim, Germany) according to the manufacturer’s protocol. According to the size difference between the wild-type and ΔFomCFEM strain hybrid fragments, the correctness of the converter was judged. The primer sequences for each gene product are listed in Supplementary Table S4.

Complementation of FomCFEM For complementation assays, a 2.6-kb fragment bearing 1.1 kb of the full-length FomCFEM gene and 1.5 kb of its promoter was amplified with the primer pair Ech1-CF/Rand then cloned into pFL2 vector by gap repair. This fusion plasmid was reintroduced into protoplasts of the ΔFomCFEM mutant to obtain complemented strains (ΔFomCFEM-C). G418-resistant transformants were screened and confirmed by PCR.

Characterization of Fom transformants To assay growth rate, the WT strain, FomCFEM deletion mutants, and FomCFEM-C transformants were cultured on PDA at 28℃ for 5 days. Then, the colony morphology and conidia of these strains were examined after being grown on PDA plate for 7 days. For stress tolerance assay, each indicated strain was inoculated onto CM plates containing 0.2 g·L^− 1 Congo red (CR), 1 mol·L^− 1 sorbitol, 0.05% SDS, 0.7 mol·L^− 1 NaCl, 0.7 mol·L^− 1 KCl, and 5 mmol·L^− 1 H₂O₂ as previously described [27, 29]. To analyze the conidial germination, 50 µL conidial suspension (1 × 10⁵ conidia·mL^− 1) of each strain was spotted on hydrophobic coverslips and incubated at 28℃ for 12 h. Then, the conidial germination was examined by the Nikon Eclipse 90 i microscope (Nikon, Tokyo, Japan). Virulence assay of each strain on bitter gourd was performed as described in Sect. 2.6.

2.9 DNA and RNA extraction and sequencing

Metagenomic sequencing Total microbial genomic DNA from each soil sample was extracted from 0.5 g of soil using a FastDNA Spin kit (MP Biomedical, Santa Ana, California, USA) or DNeasyPowerSoil Kit (QIAGEN, Inc., Netherlands) according to the instruction provided by the manufacturer. The DNA quality was assessed on a NanoDrop spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). The extracted microbial DNA was processed to construct metagenome shotgun sequencing libraries with insert sizes of 400 bp by using the Illumina TruSeq Nano DNA LT Library Preparation Kit. Each library was sequenced by the Illumina HiSeq X-ten platform (Illumina, USA) with PE150 strategy. Raw sequencing reads were processed to obtain quality-filtered reads for further analysis. Co-occurrence networks were constructed according to Zhang [30].

Transcriptome sequencing Total RNA of fungal mycelium (filtered through a filter paper) was isolated using the Trizol Reagent (Invitrogen Life Technologies), after which the concentration, quality and integrity were determined using a NanoDrop spectrophotometer (Thermo Scientific). Sequencing libraries were generated using the TruSeq RNA Sample Preparation Kit (Illumina, San Diego, CA, USA). First strand cDNA was synthesized using random oligonucleotides and SuperScript II. Second strand cDNA synthesis was subsequently performed using DNA Polymerase I and RNase H. After adenylation of the 3′ ends of the DNA fragments, Illumina PE adapter oligonucleotides were ligated to prepare for hybridization. To select cDNA fragments of the preferred 200 bp in length, the library fragments were purified using the AMPure XP system (Beckman Coulter, Beverly, CA, USA). The sequencing library was then sequenced on a Hiseq platform (Illumina) by Shanghai Personal Biotechnology Cp. Ltd. To quantify the identified gene expression in Fom, reference genome index was built by Bowtie2 (2.2.6) and the filtered reads were mapped to the Fusarium oxysporum f. sp. cubense race 1 (https://www.ncbi.nlm.nih.gov/genome/707? genome_assembly_id = 202328) reference genome using Tophat 2 (2.0.14), the default mismatch was no more than 2. We use HTSeq (0.9.1) statistics to compare the Read Count values on each gene as the original expression of the gene, and then used FPKM (fragments per kilobase of exon model per million mapped reads) to standardize the expression. Then we used DESeq (1.30.0) to analyze the genes of difference expression with screened conditions as follows: expression difference multiple | log₂FoldChange| > 1, significant P-value < 0.05.

2.10 Validation of DEGs by RT-PCR

To confirm the transcriptome data, the genes were randomly selected and verified and their primers were designed using the Primer Premier 5.0 software listed in Supplementary Table S5. Quantitative reverse-transcription PCR (RT-PCR) was performed on a Light Cycler® 96 instrument (Roche) using SYBR Green qPCR kits (TaKaRa) according to the manufacturer’s protocol. Reactions were performed at 95°C for 2 min, 45 cycles of 95°C for 10 s, 58°C for 30 s and 72°C for 30 s. EF1α was used as the reference genes. The relative expression levels of the target genes were calculated using the 2^−△△CT approach.

2.11 Statistical Analysis

All these experiments were carried out using three replicates and were repeated at least twice, with reproducible results. For multiple comparisons, two-way ANOVA with post hoc LSD and Duncan were performed using the SPSS 22.0.P < 0.05 was considered as significant.

3.1. Responses of soil physical-chemical properties and soil enzyme to the invasion of Fom

In this study, modifications in physical-chemical properties and enzyme of the soil in which it was infected by Fom were described in Supplementary Table S6. The soil organic matter content (SOC), pH, electrical conductivity (EC), available nitrogen content (AN), available phosphorus content (AP) and available iron content (AFe) increased in S-FOM treatment. The significant increase in the value of SOC (57.7%), pH (2.0%), and AFe (27.7%) were obtained in S-FOM soil. The control group had the higher values of soil fluorescein diacetate (FDA) hydrolase and polyphenol oxidase (PPO). Although the inhibitory effects of different enzymes were different, the invasion of Fom significantly inhibited these enzymes. For example, the obvious decreases in the activities of FDA hydrolase (47.5%) and PPO (17.4%) were obtained in S-FOM soil. Soil FDA hydrolase and PPO activity were negatively correlated with the relative abundance of F. oxysporum, pH, organic matter content, and available Fe. This result indicated that the inhibition of F. oxysporum on enzyme activities of the soil might play critical roles in the conversion of organic matter and the change of pH.

3.2. Responses of community structure and diversity of soil to the invasion of Fom

The data set including 10 samples comprised a total of 22,492,589 contigs with 2,005,895 to 2,511,389 sequences per sample (mean 2,161,810.6 (Control) and 2,336,707.2 (S-FOM)) and 21,972,887 scaffolds with 1,931,781 to 2,469,478 sequences per sample (mean 2,106,715 (Control) and 2,287,862.4 (S-FOM)). A total of 6,923 microbes were assembled. The higher observed species and Chao1 index were detected for S-FOM, whereas the difference was not significant. The Fom infected samples showed a significantly higher Simpson index and Shannon index of bacteria than Control, however, the Fom infected samples showed lower Simpson index and Shannon index of fungi than Control (Fig. 1A). When the beta diversity was assessed, the Control and S-FOM was significantly different (R² = 0.43,P = 0.013) (Supplementary Fig. S1A). Redundancy analysis (RDA) revealed that the environmental factor variable possesses significant correlation with community structure. The variable of pH (R² = 0.65, P = 0.017) produced significant positive correlations with S-FOM microbes communities, the variable of FDA hydrolase (R² = 0.84, P = 0.002) and PPO (R² = 0.58, P = 0.043) produced significant negative correlations with S-FOM microbes communities (Supplementary Fig. S1B).

The soil was predominantly colonized by Proteobacteria, Actinobacteria, Acidobacteria, Gemmatimonadetes, Planctomycetes, Firmicutes, Nitrospirae, and Bacteroidetes, accounting for 63.0% of the total scaffolds, when the microbe was assessed at the phylum level (Fig. 1B). Ascomycota and Basidiomycota were the largest fungal taxonomic group, but accounted for a substantially lower proportion of microbe-assigned reads (Fig. 1D). The higher abundance of increased Proteobacteria, Firmicutes, and Nitrospirae were detected in the S-FOM samples that were subjected to invasion of F. oxysporum. Meanwhile, the obvious decreases in the relative abundance of Acidobacteria, Gemmatimonadetes, Planctomycetes, and Basidiomycota were obtained in S-FOM soil (Fig. 1B and D). The modifications in the composition of the core microbe at the genus level of the S-FOM samples were described in Fig. 1C, Fig. 1E, and Supplementary Fig. S1C. The invasion of Fom significantly influenced these microbes, although the influence differed in different microbes. For example, the obvious decreases in the relative abundance of Aspergillus, Pseudomonas, Sphingobium, Lysobacter, Variovorax, Sporisorium, Methylobacterium, Mycobacterium, Micromonospora, Actinoplanes, Mycolicibacterium, Amycolatopsis, Frankia, Pseudonocardia, and Ilumatobacter were obtained in S-FOM samples. Simultaneously, the relative abundance of Fusarium, Agromyces, Ramlibacter, Streptomyces, Conexibacter and Kribbella were significantly increased in S-FOM samples. Sphingomonas, Trichoderma, and Bacillus had no difference between the control sample and the S-FOM samples. At the species level, the relative abundance of Trichoderma asperellum, T. atroviride, Bacillus anthracis, B. kochii, Talaromyces stipitatus and so on showed an increased tendency in response to the invasion of F. oxysporum (Supplementary Fig. S2). Meanwhile, Fusarium sp. was composed of F. incarnatum, F. culmorum, F. delphinoides, F. dimerum, F. equiseti, F. fujikuroi, F. graminearum, F. langsethiae, F. oxysporum, F. pseudograminearum, F. solani, F. venenatum, F. verticillioides, and Unclassified Fusarium. The F. oxysporum, F. solani, and F. fujikuroi showed significantly higher proportions than others. A significant increase of the relative abundance of F. oxysporum (51.7%) was observed in S-FOM samples, while there was no difference in other species.

3.3. Co-occurrence networks analysis of microbial communities

To study the effect of the invasion of Fom on microbial interactions, the networks were constructed (Supplementary Fig. S3). The co-occurrence network was constructed from 10 samples (Control and S-FOM) based on the strong and significant (P < 0.05) correlations of microbial species. In total, 5035 nodes (5025 microbial nodes and 10 environmental nodes) and 17912 edges were shared among co-occurrence networks, and the co-occurrence networks contains 89.54%-92.13% microbial biomass (Supplementary Table S7). The degrees of the networks followed power-law distributions, the average path length was significantly larger than the corresponding random networks, which indicated that co-occurrence networks was a small world (Supplementary Fig. S2 and Table S7). Meanwhile, modularity values were higher than the random networks and greater than 0.34. This indicated that the co-occurrence networks had modular structures. For the networks, approximately 82.75% of the nodes belonged to bacteria, 7.98% of the nodes belonged to eukaryota, 7.36% of the nodes belonged to archaea, 1.91% of the nodes belonged to viruses. Proteobacteria (33.29%), Actinobacteria (18.23%), and Firmicutes (7.02%) accounted for more links in the network. There are 353 modules among the networks, of which 15 modules had more than 50 nodes, such as module 1141, module 676, module 1458, module 1465, module 168, module 1460, module 1229, module 1163, module 1124, module 1281, module 1250, module 1238, module 1225, module 1286 and module 1283, the nodes account for 37.87% of the total number of nodes and the relative abundance of microbe the nodes account for 70.39–75.23% of the total number of microbe, including Proteobacteria (35.26%), Actinobacteria (20.86%), Firmicutes (6.25%), Euryarchaeota (2.94%), Thaumarchaeota (2.73%), Ascomycota (2.47%), and Bacteroidetes (2.15%) (Supplementary Fig. S3 C and D). We investigated the distribution of soil microbial species by Zi-Pi plot (Supplementary Fig. S4). In the network, 94.22% of the nodes belong to the peripherals, 4.39% of the nodes belong to the connectors, 0.71% of the nodes belong to the module hubs, 0.08% of the nodes belong to the network hubs, among them, the connectors, module hubs and the network hubs were considered as key microorganisms. Among the 261 key microbial species, 159 microbial species belong to Proteobacteria, 27 microbial species belong to the Actinobacteria, and the relative abundance of 13 microbial species were increased in S-FOM samples (Table 1). The relative abundance of sp1571 (Kribbella flavida), sp1615 (Microlunatus soli), sp2002 (Streptomyces sp. 4F), sp4492 (Devosia sp.), sp5385 (Burkholderia sp.), sp5454 (Paraburkholderia sprentiae), sp5496 (Acidovorax sp.), sp5501 (Alicycliphilus denitrificans) and sp5675 (Chromobacterium vaccinii) were significantly positively correlated with the relative abundance of F. oxysporum. The relative abundance of sp5951 (Haliangium ochraceum), sp744 (Acidimicrobium ferrooxidans) and sp4354 (Bradyrhizobium erythrophlei) were negatively correlated with the relative abundance of F. oxysporum.

Table 1 The topological features of key microorganisms in co-occurrence networks and the correlation coefficient with F. oxysporum and Talaromyces purpureogenus

^# means connectorssp; ^& means module hubs. * indicate a significant difference or correlation (P<0.05); ** indicate a significant difference or correlation (P<0.01). 1227: Arsenicicoccus sp.; sp1571: Kribbella flavida; sp1615: Microlunatus soli; sp2002: Streptomyces sp. 4F; sp4492: Devosia sp. A16; sp4982: Azospirillum thiophilum; sp5385: Burkholderia sp.; sp5454: Paraburkholderia sprentiae; sp5496: Acidovorax sp. RAC01; sp5501: Alicycliphilus denitrificans; sp5675: Chromobacterium vaccinii; sp5951: Haliangium ochraceum; sp744: Acidimicrobium ferrooxidans; sp1175: Modestobacter marinus; sp4354: Bradyrhizobium erythrophlei; sp7444: Unclassified Bacteria. the same below

The soil invaded by F. oxysporum was used for the discovery of beneficial microbial species antagonism against Fom. More than 100 microbial isolates were screened and successfully revealed 14 biocontrol agents with strong antagonistic activity to Fom on the plate, including 7 bacterial strains, 2 Penicillium sp. (1 Talaromyces sp.), and 5 Trichoderma sp.. Fom inhibition rate was 52%-83% (Supplementary Table S1). The secondary screening was conducted in pot trials and successfully revealed 4 biocontrol agents that have strong control effects on Fusarium wilt of bitter gourd in the greenhouse, including bacterial strain SK2 and SK6, Talaromyces sp. strain Q2, and Trichoderma sp. strain M2, their control efficacy on Fusarium wilt were 68%-79% (Supplementary Table S8). Meanwhile, the relative abundance of sp1571 (Kribbella flavida), sp2002 (Streptomyces sp. 4F), sp4492 (Devosia sp.), sp5496 (Acidovorax sp.), sp5501 (Alicycliphilus denitrificans), sp5675 (Chromobacterium vaccinii) and F. oxysporum were significantly positively correlated with the relative abundance of strain Q2, and the correlation coefficient between strain Q2 and F. oxysporum reaches 0.79 (Table 1).

3.5. Micromorphology, biocontrol potential, and phylogenetic analysis of strain Q2

Strain Q2 had strong antagonistic activity to Fom on the plate (Fig. 2A ). Colonies of strain Q2 growing moderately fast, reaching 70.46 mm on PDA, 49.44 mm on CA, 76.68 mm on CYA, 76.22 mm on OA, 70.66 mm on CMA, 75.32 mm on MM, 67.06 mm on SSA, 61.00 mm on SDA and 66.56 mm on SMA at 25 ℃ for 7 days (Fig. 2D and Supplementary Table S9 ). Colonies of strain Q2 on PDA moderately deep, plane; margins low, entire; marginal mycelia yellow and white; conidia pile greyish green to dull green; sometime exudates minute reddish droplets; reverse dark ruby (Fig. 2D). Under light microscope, hyphae of strain Q2 were septate, smooth to coarse, and had branched. Conidiphore of strain Q2 produces from hyphal string. Conidiophores biverticillate; stipes smooth walled, top not swollen (Fig. 2B and C), metulae three to five, divergent, 8.0–13 µm × 2.3-3.0 µm; phialides acerose, three to six per metulae, 10–20 µm×1.8-3.0 µm; conidia oval or round, smooth walled, 3–3.5 × 2–2.5 µm. Conidia chain loose, similar cylinder shape. Ascomata not observed. The optimal growth temperature was 30 ℃, maximum 55 ℃; optimal growth pH was 7; optimal carbon source was glucose (Supplementary Fig. S5).

Strain Q2 had great potential for biological control. For example, it was able to produce β-1, 3-glucanase and chitinase (Fig. 2E and F). The β-1, 3-glucanase activity of the crude culture filtrate was 2.51 U·mL^− 1 (Q2) and that of the culture filtrate of substrate induced was 3.13 U·mL^− 1 (Q2-CW), 4.20 U·mL^− 1 (Q2-Lam) and 4.38 U·mL^− 1(Q2-CW-Lam), respectively (Fig. 2G). The chitinase activity of the crude culture filtrate was 1.79 U·mL^− 1 (Q2) and that of the culture filtrate of substrate induced was 2.22 U·mL^− 1 (Q2-CW), 1.92 U·mL^− 1 (Q2-cChi) and 2.40 U·mL^− 1 (Q2-CW-cChi), respectively (Fig. 2H). Strain Q2 was able to inhibit mycelia growth of eleven pathogenic fungi such as F. oxysporum f. sp. momordicae (Fom), F. graminearum, F. moniliforme, F. oxysporum f. sp. cucumarinum, F. oxysporum f. sp. niveum, Pyricularia oryzae, Trichothecium roseum, Cryphonectria parasitica, Cytospora chrysosperma, Phytophthora parasitica, Rhizoctonia solani (Supplementary Table S10). After co-culture for 5 days, obvious inhibitory zone was occurred between the colony of strain Q2 and Fom, which inhibited the mycelia growth of Fom by 65.86% (Fig. 2A and Supplementary Table S10).

The DNA sequencing was used to identify the strain Q2. The lengths of nucleotide sequence for ITS, beta-tubulin gene (BenA), and calmodulin gene (CaM) was 585 bp, 446 bp and 705 bp, respectively and be submitted to GenBank with the accession number KX432212, KY047419 and KX781300. The maximum likelihood (ML) phylogenetic tree showed that strain Q2 was the same as the type strain Talaromyces purpureogenus CBS 286.36, which indicated that strain Q2 is T. purpureogenus strains (Supplementary Fig. S6A). Together with the morphological identification, we determined the strain Q2 to be T. purpureogenus, anamorphic type Penicillium purpurogenum. The genome of strain Q2 was assembled with PacBio technology, has 57 scaffolds with a length of 27.73 Mb, including 9,600 genes with the average length of 1,633 bp (Supplementary Table S11). Although, 8968 genes of strain Q2 could be annotated in the KOG database, the function of 41.1% was unknown (Supplementary Fig. S6B). Meanwhile, a total of 463 CAZymes (Carbohydrate active enzymes) have been identified, encoding 200 glycoside hydrolases (GHs), 83 glycosyl transferases (GTs), 84 carbohydrate esterases (CEs), 79 auxiliary activities (AAs), 15 carbohydrate-binding modules (CBMs), and 2 polysaccharide lyases (PLs) (Supplementary Fig. S6C). The number of CE1 (Acetyl xylan esterase), CE10 (arylesterase), GH18 (chitinase), GH3 (β-glucosidase), GH16 (xylanase), and GH109 (α-N-acetylgalactosidase) were significantly higher than other CAZymes (Supplementary Fig. S6D).

3.6. The control efficacy of strain Q2 on Fusarium wilt of bitter gourd

The control efficacy of strain Q2 on Fusarium wilt of bitter gourd was 63.4% (greenhouse) and 60.2% (field) (Table 2). The disease incidence and disease index of bitter gourd seedlings were significantly decreased after treatment with strain Q2. Besides, the time of bitter gourd seedling infected by Fom was obviously postponed, which demonstrated that the resistance to Fusarium wilt of bitter gourd seedlings was enhanced after inoculation with strain Q2. Meanwhile, the obvious modifications the composition of the fungal communities at the genus level of the Fom pathogenic soil treated with strain Q2 were discovered in this study (Supplementary Fig. S7). The lower observed species and Chao1 index were detected in TP samples (The pathogenic soil with Fusarium inoculation was treated with antagonistic fungus T. purpureogenus strain Q2), whereas the difference was not significant, the TP samples showed a significantly lower Simpson index and Shannon index of fungi than S-FOM samples, but, the TP samples showed a significantly higher Simpson index and Shannon index of bacteria than S-FOM samples (Supplementary Fig. S7A). The relative abundance of Actinobacteria, Firmicutes, Ascomycota, Chytridiomycota, Streptomyces, Conexibacter, Sphingomonas, Rhodoplanes, Gemmatirosa, Sphingobium, Nocardioides, Mesorhizobium, Lysobacter, Mycobacterium, Micromonospora, Talaromyces, Aspergillus, Thermothelomyces, and Chaetomium were increased in TP samples. The decreases in the relative abundance of Fusarium, Cordyceps, Colletotrichum, Anthracocystis, Trichoderma, Thielavia, Cladophialophora, Gaeumannomyces, Sporisorium, Verticillium, Fonsecaea, Exophiala, Trametes, and Diplodia were observed in TP samples (Supplementary Fig. S7 B-E). Talaromyces became the largest fungal taxonomic group in the TP samples compared with S-FOM samples, the relative abundance of strain Q2 that obtained by comparing the genome of the strain Q2 to the metagenomic sequencing data in the TP samples was 2828.01 reads compared to S-FOM samples with 51.87 reads and control samples with 27.09 reads (Supplementary Fig. S7E and I). A random forest analysis was then applied to identify the major microbes that contribute to the variation in soil microbiota. Ascomycota, Thaumarchaeota, Euglenida, Firmicutes, Chytridiomycota, Actinobacteria, Nitrospirae, Bacteroidetes, Microsporidia and so on were the primary factor that affected features of soil of microbial community (Supplementary Fig. S7 F). Spearman’s correlation network analyses confirmed that Actinobacteria, Ascomycota, Nitrospirae, Tenericutes, Bacteroidetes, Bathyarchaeota, Beckwithbacteria, Campbellbacteria, Saccharibacteria, Tectomicrobia, Woesebacteria, Wolfebacteria, Elusimicrobia, Thaumarchaeota and T. purpureogenus obviously influenced the Simpson index of soil microbiota (Supplementary Table S12). Meanwhile, T. purpureogenus seems to have a strong positive correlation with Actinobacteria, Firmicutes, Thaumarchaeota, Ascomycota, Bacillariophyta, Bacteroidetes, Beckwithbacteria, Campbellbacteria, Saccharibacteria, Woesebacteria, Wolfebacteria, Elusimicrobia, and Tenericutes (Supplementary Table S12). The relative abundance of Fusarium decreased in the TP sample (Supplementary Fig. S7G), but the number of Fusarium in the rhizosphere of bitter gourd had decreased significantly on the Fusarium selective medium (Komada’s medium) (Supplementary Fig. S7H). The density of Fusarium under TP treatment was lower, with only 10,700 cfu·g^− 1 soil compared to 17, 600 cfu·g^− 1 soil in S-FOM treatment. Meanwhile, strain Q2 could colonize in soil, but its increasing trend in soil was gradually decreasing (Supplementary Fig. S7J).

Table 2 The biocontrol efficacy of strain Q2 on Fusarium wilt of bitter gourd in pot and field

The co-occurring network was built based on S-FOM and TP, in order to more clearly observe the influence of T. purpureogenus on the network interactions among microbes. The co-occurrence network included 4954 nodes (4946 microbial nodes and 8 environmental nodes), 16275 edges, 1752 modules and 87.79%-89.00% microbial biomass (Supplementary Fig. S9A). Proteobacteria, Actinobacteria, Acidobacteria Planctomycetes and Firmicutes were dominant phyla (Supplementary Fig. S9B). 13 modules had more than 50 nodes, sunch as module 897, module 1725, module 1491, module 1441, module 1448, module 1468, module 1368, module 543, module 1529, module 1450, module 1409, module 1492 and module 1505 (Supplementary Fig. S9A). There are 1512 nodes, 33 phyla and 62.77%- 64.63% microbial biomass in 13 larger modules (Supplementary Fig. S9 C). T. purpureogenus participated in the construction of co-occurrence network as environmental factor. It showed significant negative correlations with the Acidobacteria, Firmicutes, Chloroflexi Armatimonadetes, Chromerida, Chytridiomycota and Chlamydiae, obvious positive correlations with Actinobacteria and Dictyoglomi (Supplementary Fig. S9 E). In contrast, F. oxysporum showed positive correlations with the Acidobacteria, Firmicutes, Chloroflexi Armatimonadetes, Chromerida, Chytridiomycota and Chlamydiae, negative correlations with Actinobacteria and Dictyoglomi. Meanwhile, T. purpureogenus tended to co-occur with Actinobacteria and Proteobacteria of module 897. Module 897 was a module with Actinobacteria as the main component (Supplementary Fig. S9 D). The microbial species, such as Kribbella flavida, Frankia casuarinae, Kitasatospora setae, Nocardiopsis alba, Streptomyces niveus and Xylanimonas cellulosilytica, in module 897 were significantly positively correlated with the relative abundance of T. purpureogenus. The results indicated that T. purpureogenus strain Q2 involved in the resistance of soil microbial communities to the invasion of F. oxysporum, and the growth and development of F. oxysporum was inhibited by T. purpureogenus and its “friends”.

3.7 Inhibition of the growth and development of Fom by strain Q2 in co-culture

To understand the mechanisms of inhibition of growth and development of Fom by strain Q2 under natural condition. In this study, the liquid co-cultivation of strain Q2 and Fom was firstly attempted, and the growth of strain Q2 and FOM in the axenic and co-culture conditions in PDB were significantly different at 48 h. It was observed that the culture liquor of FOM appeared beige, thicker, hazy mixture, while the color of culture liquor of strain Q2 was relatively lighter and appeared red white after they were singly-cultured for two days. However, when the inoculation ratio of strain Q2 and FOM was greater than 2:1, after co-cultivation for two days, the culture broth color was red and deepened with time, and the mycelia of strain Q2 were dominant (Supplementary Fig. S10). Fom hyphae from the control group was not changed under the light microscope and hyphae grown and developed normally (Fig. 4A). Strain Q2 hyphae of all the groups were normally grown and developed, and the color of strain Q2 was not changed (Fig. 4A). There were morphological changes that hypha swelling appears round, the appearance of septa digestion, cytoplasma concentration and more than globular material in hypha cells of Fom compared with control group (Fig. 4A). Propidium iodide (PI) staining showed that the Fom conidia with red pigment accumulation lost vitality in T1 to T4 samples (Fig. 4B and C). The number of conidia of Fom in culture liquor was lower in T1, T2, T3, T4 and T5 treatment, with only 8.20 × 10⁷ cfu·mL^− 1, 6.73 × 10⁷ cfu·mL^− 1, 9.83 × 10⁶ cfu·mL^− 1, 4.33 × 10⁶ cfu·mL^− 1 and 1.77 × 10⁶ cfu·mL^− 1, respectively, compared to 1.20 × 10⁸ cfu·mL^− 1 in control treatment at 5th days of cultivation (Fig. 4D). The results showed that the hyphae growth and spores production of Fom was inhibited significantly in the co-culture by strain Q2, and the red pigment produced by Q2 penetrated cell wall and cell membrane of the dead cell of Fom and deposited.

The presence of propidium iodide and the red pigment in Fom cells indicated that the integrity of the cell wall and cell membrane of Fom were damaged. We investigated the effect of strain Q2 on malondialdehyde (MDA) concentration and triglyceride content of Fom cell. The results showed that the MDA content and triglyceride content of Fom cell increased in T1 to T4 samples compared to control (Fig. 4F and G). The MDA content and triglyceride content in T2 treatment were significantly higher than those in the control. At the same time, it is worth noting that the small fragments increased in the culture filtrate of T3 samples (Fig. 4B), and the N-acetylglucosamine (N-acetylglucosamine is a hydrolysate of chitin) content in the culture filtrate of T1 to T3 samples was significantly higher than that in the control (Fig. 4E), and the content of N-acetylglucosamine gradually increases with inoculation concentration of strain Q2 increases. The results indicated that strain Q2 maybe inhibit the growth and development of Fom by destroying the integrity of the cell wall and cell membrane of Fom.

3.8 The influence of strain Q2 on the metabolic function of Fom at transcriptional level

A total of more than 336.47 million high quality clean reads were obtained after removing duplicate reads, trimming adapters and low quality sequences from six libraries. Q30 percentage was over 92% for all the samples. These data showed that the RNA-Seq quality was applicable for further analysis. In FOM samples, 89.20 % to 89.73 % of the total reads were mapped to the genome of F. oxysporum f. sp. cubense race 1; and more than 95 % of the reads in each sample were uniquely mapped to the genome. In Q2-FOM samples, 28.25 % to 47.17 % of the total reads were mapped to the genome of F. oxysporum f. sp. cubense race 1; and 96.86 % to 97.15% of the reads in each sample were uniquely mapped to the genome (Supplementary Table S13).

Differentially expressed genes (DEGs) were identified between FOM and Q2-FOM groups through comparing the expected read counts normalized using the DEseq package. To compare the expression levels of the DEGs between two groups, sequences count based differential expression analysis was performed. It was observed that 4,893 genes expression levels in FOM were identified as differentially expressed in FOM and Q2-FOM groups, including 2,435 up-regulated and 2,458 down-regulated (Supplementary Fig. S11A). The Gene Ontology analysis (GO) and KEGG analysis of DEGs was further conducted (Supplementary Fig. S11B and C). Following GO annotation, the DEGs were classified into different groups belonging to biological process (BP), molecular functions (MF) and cellular component (CC) categories. In the all category, the highly represented groups were ‘intrinsic component of membrane’ (410 up-regulated genes and 748 down-regulated genes), ‘integral component of membrane’ (410 up-regulated genes and 744 down-regulated genes), ‘membrane part’ (416 up-regulated genes and 748 down-regulated genes), ‘membrane’ (428 up-regulated genes and 762 down-regulated genes), ‘oxidoreductase activity’ (355 up-regulated genes and 207 down-regulated genes), ‘transmembrane transport’ (149 up-regulated genes and 256 down-regulated genes) and ‘transmembrane transporter activity’ (65 up-regulated genes and 146 down-regulated genes) groups. In the Q2-FOM groups, the genes associated with membrane and transport appeared to be the most significantly enriched. To further investigate the biochemical pathways of these DEGs in the process of T. purpureogenus inhibiting the growth of Fom, we mapped all DEGs identified in the RNA sequencing to terms in the KEGG database. The results showed the DEGs could be classified into four pathways, including metabolism, genetic information processing, environmental information processing and cellular processes. Top 20 enriched KEGG pathways associated with DEGs affected by T. purpureogenus were showed in Supplementary Fig. S11B, including glyoxylate and dicarboxylate metabolism, glutathione metabolism, valine, leucine and isoleucine degradation, glycerophospholipid metabolism, phenylalanine metabolism, nitrogen metabolism, propanoate metabolism, pyruvate metabolism, ABC transporters, galactose metabolism, tyrosine metabolism and glycerolipid metabolism. These results further indicated that stress of strain Q2 may greatly influence carbohydrate metabolism, amino acid metabolism, lipid metabolism, energy metabolism, component of membrane and membrane transport of F. oxysporum.

Based on the STRING output, the interaction of 371 coding genes was described (Fig. 5A), including cell cycle relating genes, starch and sucrose metabolism relating genes, amino sugar and nucleotide sugar metabolism relating genes, glycan biosynthesis and metabolism relating genes and biosynthesis of secondary metabolites relating genes (Supplementary Table S14). Real-time PCR was used to validate the expression levels of thirteen genes (Fig. 5C), including cell cycle relating genes, starch and sucrose metabolism relating genes, amino sugar and nucleotide sugar metabolism relating genes and glycan biosynthesis and metabolism relating genes. Among these genes, FOC1 g10015395 (CDC5, Cell cycle serine/threonine-protein kinase), FOC1 g10015375 (CG21, G2/mitotic-specific cyclin-B), FOC1 g10014984 (MPG1, Mannose-1-phosphate guanyltransferase), FOC1 g10014983 (HXKG, Glucokinase), FOC1 g10014811(GNPI1, Glucosamine-6-phosphate isomerase 1), FOC1 g10014809 (NAGA, N-acetylglucosamine-6- phosphate deacetylase), FOC1 g10014793 (INV, Invertase), FOC1 g10014440 (GPI1, N-acetylglucosaminyl- phosphatidylinositol biosynthetic protein), FOC1 g10013913 (ARK1, Serine/threonine-protein kinase), FOC1 g10013617 (BGL1B, Beta-glucosidase 1B), FOC1 g10012541 (TREA, Trehalase), FOC1 g10012205 (MALT, Alpha-glucosidase), FOC1 g10012008 (GPI18, GPI mannosyltransferase 2), FOC1 g10010259 (BIR1, Protein), FOC1 g10009667 (GPI15, Phosphatidylinositol N-acetylglucosaminyltransferase subunit), FOC1 g10009617 (WSC1, cell wall integrity and stress response component 1), FOC1 g10009158 (TPSY, Putative alpha, alpha-trehalose-phosphate synthase), FOC1 g10008820 (AGDC, Probable alpha/beta-glucosidase), FOC1 g10008111 (ASE1, Anaphase spindle elongation protein 1), FOC1 g10006897 (GPI2, Phosphatidylinositol N-acetylglucosaminyltransferase subunit), FOC1 g10005945 (TPS2, Trehalose-phosphatase), FOC1 g10005528 (CHS7, Chitin synthase 7), FOC1 g10003572 (AMYG, Glucoamylase) and FOC1 g10003114 (BUB1, Checkpoint serine/threonine-protein kinase) were down-regulated in Q2-FOM. FOC1 g10015825 (GUN, Endoglucanase), FOC1 g10015543 (AGM1, Phosphoacetylglucosamine mutase), FOC1 g10012466 (GAL10, Bifunctional protein) and FOC1 g10011050 (YBB2, Probable mannose-1-phosphate guanyltransferase) were up-regulated in Q2-FOM. Four GPI biosynthesis and one chitin synthases genes were all down-regulated in Q2-FOM. The gene expression results from the real-time PCR were consistent with the transcriptome data (Fig. 5B).

We found a gene with CFEM domains (the characteristic eight cysteine residues at N terminus) and putative GPI-anchored site (FOC1 g10014283) in further analysis of differentially expressed genes. The expression level of FOC1 g10014283 was significantly decreased in Q2-FOM samples compared to FOM samples, although his function is unknown in F. oxysporum (Fig. 6A and Supplementary Table S14). We experimentally isolated and characterized an FOC1 g10014283 orthologue, FomCFEM, in F. oxysporum f. sp. momordicae.

3.9 FomCFEM is important for full virulence and stress tolerance

The SMART-PFAM revealed that FomCFEM consists of a CFEM domain, a signal peptide, and the potential GPI-modification site (the amino acids N²⁹¹ and A²⁹²) (Supplementary Fig. S12A and Supplementary Fig. S13). A phylogenetic analysis of the putative FomCFEM, homologous sequences of Fusarium and other fungal genes with CFEM domain indicated that homologous sequences of FomCFEM were common in Fusarium and the FomCFEM is closely homologous to the gene with CFEM domain from Trichoderma asperellum (GFP54046.1) (Supplementary Fig. S12B). To assess functions of FomCFEM, in F. oxysporum f. sp. momordicae by homologous recombination, FomCFEM deletion transformants harboring a hygromycin-resistance gene instead of the entire FomCFEM were generated. After PCR and Southern blot assays, the predicted FomCFEM separate mutants were obtained (Supplementary Fig. S12C and Fig. S14B). One of the mutants (designated as ∆FomCFEM) was selected for complementation by reintroducing the entire FomCFEM gene under the control of its native promoter into the ∆FomCFEM mutant (Supplementary Fig. S12D).

To determine the role of FomCFEM in pathogenicity, the rhizosphere soil of bitter gourd seedlings were inoculated with conidial suspensions of WT, ∆FomCFEM mutant, and ∆FomCFEM-C strains (the inoculation amount was 1×10⁶ spores/g·soil). The disease incidence and disease index of bitter gourd seedlings were significantly decreased after treatment with ∆FomCFEM mutant (Fig. 6B and C). The ∆FomCFEM mutant resulted in a 58.8% reduction in virulence compared to the WT strain. These results suggest that FomCFEM is essential for full virulence in F. oxysporum f. sp. momordicae. Although conidial production of ∆FomCFEM was decreased compared with WT and ∆FomCFEM-C strains, there was no difference in the conidium morphology and conidial germination of ∆FomCFEM, WT, and ∆FomCFEM-C (Fig. 6F-G). In addition, to investigate whether FomCFEM is involved in external stress tolerance, growth assays of indicated strains on CM media supplemented with salt stress (0.7 mol·L^− 1 NaCl and 0.7 mol·L^− 1 KCl), osmotic stress (1 mol·L^− 1 sorbitol), H₂O₂ (5 mmol·L^− 1) and cell wall stress (0.05% SDS, and 0.2 g·mL^− 1 CR) were analyzed, respectively. The results showed that the ∆FomCFEM strain was more hypersensitive to the cell wall stress compared with WT and ∆FomCFEM-C strain (Fig. 6D and E). These results indicated that ∆FomCFEM may play an important role in the integrity and function of the cell wall of F. oxysporum f. sp. momordicae.