Salicylic Acid Signals Recruit the Rhizosphere Microbiome Cooperate Resistance to Watermelon Fusarium Wilt

Aims: In order to nd out the mechanism of how salicylic acid signal recruit soil microorganisms to explain their cooperate resistance to Fusarium wilt disease in watermelon plants. Methods: In this experiment, we have examined the soil microbial diversity at three different sampling times after salicylic acid application by metagenomic using Illumina Novaseq6000 platform. Results: The results showed that the incidence of Watermelon Fusarium Wilt was signicantly lower than control group after SA application. Notably, our results indicated that the application of salicylic acid recruit soil microbial communities in watermelon plants. Furthermore, the soil microbial metabolic pathway and microbial function were signicantly different after SA treatment which is benecial to the antagonism of plants against pathogens. Conclusion: Therefore, our results suggested that the induced resistant observed in watermelon against Fusarium wilt might be a case of cooperate resistance dependent on salicylic acid signals to recruit soil microorganisms.


Introduction
Watermelon (Citrullus lanatus) is widely consumed in the summer seasons that is an important horticultural crop around the world (Everts and Himmelstein, 2015). However, watermelon Fusarium wilt disease caused by Fusarium oxysporum f. sp. Niveum (FON), which poses a sever decline to watermelon yield and quality. Firstly, the FON penetrates the root epidermis and then grows intracellularly in the root cortex. Ultimately it usually reaches the roots xylem vessels and spreads to the main root and stem in plants (Gordon, 2017). Since researches into plant microbiota has received substantial attention, since it in uences both plant health and productivity (Berendsen et al., 2018;Levy et al., 2017). Deciphering the various types of interaction between plants and their microbiomes is a hot topic for research in plant sciences and agronomy as well as in ecology (Mendes et al., 2011). Our previous studies have found the important role of soil microbial community structure in controlling the occurrence of watermelon wilt (Zhu et al., 2018;Zhu et al., 2019). More importantly, we have described the altered soil microbial community structure is critical to suppress watermelon Fusarium wilt recently (Zhu et al., 2020). However, there is lack of the molecular mechanism to explain how watermelon plants resistant to Fusarium wilt.
In another sense, plant microbiota is considered as a second genome to a plant (Berendsen et al., 2012).
The establishment of plant and rhizosphere microbiome interaction is a highly coordinated event in uenced by the plant host and soil. For instance, most studies have shown that the plant immune system shapes the microbiome, and these microbiomes can increase the plant immune capacity (Aleklett et al., 2018;Heintz and Mair, 2014;Levy et al., 2017). Reviews of biosynthesis and signal transduction mechanisms of SA have discussed its necessity for activating the hypersensitive response and mediating SAR (Olowe et al., 2020;Pieterse et al., 2012;Yang et al., 2015;Zhang and Li, 2019), as well as its role as a main activator of signal transduction in defense processes to pathogen infection (Berendsen et al., 2018). Studies have also shown that pretreatment with SA enhanced resistant to fungal diseases, such as Fusarium wilt of Arabidopsis , banana (Dou et al., 2019;Wang et al., 2015), chickpea (Mannur et al., 2019) and tomato (Mandal et al., 2009). Lv et al. noticed that FON infection increases the accumulation of salicylic acid in watermelon grown and identi ed some up-regulating disease-and defense-responsive genes in watermelon under wheat intercropping system (Lv et al., 2018). More recently, a genetic approach has shown that SA affects the abundance of speci c bacterial groups in the root via a combination of direct and indirect effects (Trivedi et al., 2020). SA signaling appears to have evolved to regulate plant defense against microbial pathogens after colonization of the land (Berens. et al., 2017;Bowya and Balachandar, 2020;Wang and Sugiyama, 2020). Overall, these results suggested the existence of complex tness trade-offs where the result of the plant-microbe interaction is determined by the speci c combination of plant accession, plant pathogen and microbe strain. Nevertheless, our understanding of the functional role of the plant microbiome remains limited.
At present, it was perceived that the functional versatility and function-based diversity of the microbiome are likely to be dominant factors in niches rather than mere diversity (Delgado-Baquerizo et al., 2018;Liu et al., 2020). As the establishment of plant and rhizosphere microbiome interaction is a highly coordinated, we hypotheses that the salicylic acid signals recruit and remodel rhizosphere microorganisms. Therefore, in this experiment, we have examined the soil microbial diversity at three different sampling times after salicylic acid application by metagenomic sequencing using Illumina Novaseq6000 platform. Notably, our results indicated that the SA signal interplay between microbiota and the plant immune system likely plays a critical role in shaping bene cial plant-microbiota combinations to cooperate resistance to watermelon Fusarium wilt. Moreover, a deeper understanding of the microbial function patterns required for plant colonization will pave the way for the development of highly effective measures to control the disease in agricultural practices.

Experiment description and Sampling
This study was conducted at Hunan Academy of Agricultural Sciences in the city of Changsha, Hunan Province, China (112 °5 8'42''E,28 °1 1'49''N). The soil was sandy loam with same background, sterilization before use (LDZM-80KCS-3 Vertical pressure steam sterilizer, ZHONGAN, Shanghai, China). The trial crop was watermelon strain 8424, which was provided by Xinjiang Farmer Seed Technology Co., Ltd. China. The watermelon seedling nutrition bowl was cultivated and grown in biochemical incubator at Tm 25°C, light 16h/ Tm 18°C, dark 8h. After 30 days, we transplanted each plant into the pots separately (LRH-300, ZHUJIANG, Taihong, Shaoguan, Guangdong, China). 100uM exogenous Salicylic acid (BBI Life Sciences, Shanghai, China) were prepared. FON1 cultured on PDA at 28°C in the dark for 7days. A bam plug from a PDA plate was placed into 300ml of potato dextrose broth in a ask and placed on a rotary shaker at 200rpm. SA, 5ml aliquots of 100uM exogenous SA were irrigated into the root zone of each plant and added 5ml again 24h later. The result of screening optimum concentration of exogenous SA was list in Supplementary material 1. After 2 days, 5ml aliquots of 10 6 conidia/ml FON1 were irrigated into the root zone of each plant. SF, 5ml aliquots of 10 6 conidia/ml FON1 were irrigated into the root zone of each plant. S, 5ml sterile deionized H 2 O were irrigated into the root zone of each plant.
The rhizosphere soils of 15 watermelon plants were selected as one repetition set three replicates for each sample group. Samples named as S0 (Control group), S3 (Control group, 3 dpi: days post infection), S5 (Control group, 7dpi), SA3 (SA + FON1 treatment, 3dpi), SA5 (SA + FON1 treatment, 7dpi), SF3 (FON1 treatment, 3dpi), SF5 (FON1 treatment, 7dpi). At each sampling time, every treatment of rhizosphere soil sample was pooled from 15 plant roots as one repetition, three independent replicates. Samples were placed in 5ml sterile centrifuge tubes and then we divided each sample into three parts. Two of them were placed in sterile centrifuge tubes, stored at -80°C for sequencing analysis. While the other was used for measuring the soil properties, stored at room temperature. We have collected total of 21 samples in three different sampling times.

The disease investigation
We calculated disease incidence which listed in Fig. 1B. The disease incidence (%) = (number of infected plants / total number of surveys) x 100%. DNA extraction, library preparation, and Metagenomic sequencing A total amount of 1ug DNA per sample was used input material for the DNA sample preparations. Brie y, a size of 350bp DNA sample was fragmented by sonication, then DNA fragments were end-polished, Atailed, and ligated with the full-length adaptor for Illumina sequencing with further PCR ampli cation. Then PCR products were puri ed (AMPure XP system) and libraries were analyzed for size distribution by Agilent2100 Bioanalyzer. 1% agarose gel was monitored on DNA degradation degree and potential contamination, the OD value is between 1.8-2.0. Preprocessing the Raw Data obtained from the Illumina Novaseq6000 sequencing platform using Readfq V8, https://github.com/cj elds/readfq was conducted to acquire the Clean Data for subsequent analysis. Clean Data need to be blast to the host database which default using Bowtie2.2.4 software (Bowtie2.2.4, http://bowtiebio.sourceforge.net/bowtie2/index.shtml) to lter the reads that are of host origin. Then all the reads not used in the forward step of all samples are combined and then use the software of SOAPdenovo (V2.04) for mixed assembly with the parameters same as single assembly; Break the mixed assembled Scaffolds from N connection and obtained the Scaftigs. Filter the fragment shorter than 500 bp in all of Scaftigs for statistical analysis both generated from single or mixed assembly.
The Scaftigs (≥ 500 bp) assembled from both single and mixed are all predicted the ORF by MetaGeneMark (V2.10, http://topaz.gatech.edu/GeneMark/) software, and ltered the length information shorter than 100 nt (Nielsen et al., 2014) from the predicted result with default parameters. For ORF predicted, CD-HIT (Fu et al., 2012) software (http://www.bioinformatics.org/cd-hit ) is adopted to redundancy and obtain the unique initial gene catalogue (Sunagawa et al., 2015;Zeller et al., 2014). The Clean Data of each sample is mapped to initial gene catalogue using Bowtie2.2.4 and get the number of reads to which genes mapped in each sample with the parameter setting (Li and Jia, 2014) are --end-toend, --sensitive, -I 200, -X 400. Filter the gene which the number of reads ≤ 2 (Qin et al., 2012) (Cantarel et al., 2009) database (http://www.cazy.org/). For each sequence's blast result, the best Blast Hit is used for subsequent analysis (Bäckhed et al., 2015). Based on the abundance table 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 are conducted, but also the Anosim analysis of the difference between groups (inside) based on functional abundance, comparative analysis of metabolic pathways, the Metastat and LEfSe analysis of functional difference between groups are performed.
Use Resistance Gene Identi er (Pedraza et al.) software to align the Unigenes to CARD database (https://card.mcmaster.ca/) (Jia et al., 2017;Martínez et al., 2015) with the parameter setting are blastp, evalue ≤ 1e-30. Based on the aligned result, count the relative abundance of ARO. Based on the abundance of ARO, the abundance bar charts, the abundance cluster heatmap and the resistance genes' number difference between groups are displayed. In the same way, the resistance genes' abundance distribution in each sample, the species attribution analysis of resistance genes and the resistance mechanism of resistance genes analysis are also conducted.

Statistical analysis
The exhibition of abundance cluster heat map, PCA (R ade4 package, Version 2.15.3) and NMDS (R vegan package, Version 2.15.3) decrease-dimension analysis are based on the abundance table of each taxonomic hierarchy. The difference between groups is tested by Anosim analysis (R vegan package, Version 2.15.3). In addition, functional-enrichment analysis including GO and KEGG were performed to identify which DEGs were signi cantly enriched in GO terms and metabolic pathways at Bonferroni-corrected p-value ≤ 0.05 compared with the whole-transcriptome background. Statistical analysis was performed using SPSS version 20.0 (SPSS Inc., Chicago, IL, USA) and Microsoft O ce 2010 (Microsoft Corporation, Redmond, WA, USA). Differences between two treatments were tested by independent sample t-test at P ≤ 0.05. ANOVA, followed by Tukey's test (P ≤ 0.05) was used to analyze the differences among treatment means. All the values were expressed as mean ± standard error. The gures of the microbial diversity indices and relative abundance of functional pro les were prepared using Adobe Illustrator CS5 (Adobe Systems Incorporated, San Jose, CA, USA).

Results
Effect of salicylic acid application on watermelon Fusarium wilt Compared to SA-treated plants, the plants without SA application showed disease symptoms of rotted, discolored root at 7 days post infection of FON1 (Fig. 1A) and the incidence of watermelon Fusarium wilt after SA application was signi cantly lower than control group (Fig. 1B).

Sequencing and metagenome assembly
The depth analysis of the 15 soil samples sequencing gene length distribution data indicated that the ORF length of gene length was mostly around 200-600 nt ( Fig. 2A). The Fig. 2B of QC analysis index re ects the evenness of the number of observed ORFs in the samples tends to be consistent and the GC quality all above 60%. The Venn diagram showed that there were 624730 overlapped genes in all groups of samples and there were unique genes distribution among each groups as well (Fig. 2C). And we found gene number of different groups at all sampling times had signi cant differences showed on the box plot (Fig. 2D). In addition, we found that the structures of soil microbial communities have dynamics through all sampling times based on beta diversity analysis performed by NMDS (Fig. 3A on phylum level and Comparative analysis of soil microbial community structure at different sampling times Therefore, to nd the dynamic changes of dominant soil microbial communities during all three sampling times with different treatments, we used community bar plot analysis to identify the top10 abundance soil microbial communities both on phylum and genus levels ( Fig. 3B and Fig. 3D). The study demonstrated that they phylum of Actinobacteria, Proteobacteria, Bacteroidetes, Gemmatimonadetes, Firmicutes, Acidobacteria, Verrucomicrobia, Deinococcus-Thermus, Elusimicrobia (Meheust et al., 2020) and Chloro exi were consistently present throughout plant development. Additionally, the heatmap showed different genus have dynamic changes of different samples (Fig. 4A) as well. Furthermore, we used student's t-test to detect most abundance soil microbial communities' signi cant differences at 3days post inoculation ( Fig. 4B and Fig. 4D) and 7 days post inoculation ( Fig. 4C and Fig. 4E) respectively after SA application both on phylum and genus levels. This suggests that while the structure of the overall soil microbial community does not change, speci c microbes may be changing through development.

Common functional database annotations and taxonomy prediction
In order to nd out the molecular mechanism to explain how the rhizosphere microbial community cooperate resistance to watermelon Fusarium wilt, we blasted the unique genes to functional database (KEGG, eggNOG, CAZy, CARD) to annotation their abundances, visual display and statistical analysis.
By analyzing the enrichment of KEGG pathways cluster on level 1, four major pathways were signi cantly activated after SA application at 3 days post inoculation, including environmental information processing, cellular processes, metabolism and human diseases (Fig. 5A). The 43 pathways detailed information on cluster level 2 have displayed in Fig. 5B. And the heatmap analysis of KEGG ortholog group (Fig. 5C) showed that the K06042 (precorrin-8X/cobalt-precorrin-8 methyl mutase), K02055 (putative spermidine/putrescine transport system substrate-binding protein), K16164 (acyl pyruvate hydrolase) and K18930 (D-lactate dehydrogenase) had highly enrichment on SF3 compared with SA3.
Furthermore, the distribution of LDA value of KEGG ortholog group between SA3 and SF3 indicated that the K02014 (iron complex outer membrane receptor protein) had increased in SA3 while the K03088 (RNA polymerase sigma-70 factor, ECF subfamily) had increased in SF3 (Fig. 5D). Surprisingly, our KEGG metabolic pathways results have con rmed the key role of rhizosphere microbial in resistance of watermelon to Fusarium wilt through several signal transduction pathways to activate a series of defensive feedback under environment stress. At the same time, the eggNOG classi cation has merged those signi cantly expressed genes to 12 groups which showed in Fig. 6A. The histogram of LDA value distribution shows the resistance genes with statistical difference between groups. The Fig. 6B showed that genes in groups of O: Post translational; U: Intracellular tra cking, secretion and N: Cell motility and M: Cell wall: membrane envelope biogenesis had increased signi cantly in SA3, while the K: Transcription; Q: Secondary metabolites biosynthesis had increased signi cantly in SF3. The Fig. 6C showed that genes in groups of J: Translation, ribosomal structure and biogenesis D: Cell cycle control, cell division; G: Carbohydrate transport and metabolism; U: Intracellular tra cking, secretion; C: Energy production and conversion; V: Defense mechanisms; L: Replication, recombination and repair had increased signi cantly in SF5. Heatmap analysis of 6 CAZy classes relative abundance in samples. The NMDS results showed that the distribution of resistance genes has dynamics through all sampling times (Fig. 7A). We noti ed that the heatmap analysis indicated that the class of GH (Glycoside Hydrolases) and PL (Polysaccharide Lyases) have highly accumulated in SA3 (Fig. 7B). The detailed unigenes signi cantly differenced relative abundance listed in Supplementary material 2. Furthermore, to identify the relationship between resistance genes and species we used two circle graphs showed in Fig. 8. According to the species annotation results of all samples in the group, the relative abundance of Proteobacteria corresponding to resistance gene is the highest in SA3 which means more resistance genes may came from these kinds of species. Notably, our results indicated that the relative abundance of Chloro exi, Acidobacteria and Cyanobacteria have increased in SA5, and the Ascomycota only in SF groups, besides the Firmicutes only observed in SF5.
Therefore, our results explained that the SA treatment which is bene cial to the antagonism of plants against pathogens,which help con rmed our hypothesis of that SA signal recruit soil microorganisms to cooperate resistance to Fusarium wilt disease in watermelon plants.

Discussion
Effect of salicylic acid on rhizosphere microbial community The rhizosphere plays fundamental role in microbe-microbe and plant-soil-microbe interactions. The microbes and their interactions can extend the plant capacity for disease resistance as well as improving their nutrient use e ciency (Berendsen et al., 2018). Emerging evidence indicates that some microbial symbionts communicate with plant immune system based on mutiple feedback mechanisms, such as pathogen resistant and maintain plant growth and development (Liu et al., 2020). In agree with the theory, our Venn graph, box plot of gene number and NMDS plot results indicate that the structure of the 10 dominant soil microbial communities have dynamics through all sampling times after the application of SA. Likewise, we observed that the p__Proteobacteria had accumulated signi cantly in SA3 but the p__Actinobacteria had signi cantly increased in SF3. Moreover, the p__Blastocladiomycota and p__Elusimicrobia had signi cantly accumulated in SA5, while the p__Dictyoglomi, p__Gemmatimonadetes and p__Ascomycota had accumulated signi cantly in SF5. Indeed, we identi ed that the community abundance of Rhodanobacter had increased signi cantly however the Actinoplanes, Burkholderia, Chryseolinea, Luteimonas, Micromonospora, Nitrolancea, Ohtaekwangia, Sorangium, Sphingomonas and Xanthomonas had decreased signi cantly in SA3 compared with SF3 showed in Fig. 4D. And in Fig. 4E, the results showed that the Stenotrophomonas, Sphingomonas, Ramlibacter, Herbaspirillum, Polaromonas and Azospirillum had increased signi cantly, however the Dokdonella, Gemmatirosa, Castellaniella, Gemmatimonas, Chitinophaga, Myxococcus, Mizugakiibacter, Trichoderma and Moritella had decreased signi cantly in SA5 compared with SF5.
Therefore, these signi cantly altered microbial community con rmed our hypothesis that the salicylic acid signal recruit and remodel rhizosphere microorganisms. Another scienti c question came out is that whether the change of rhizosphere micro ora under the in uence of salicylic acid bene cial to the antagonism of plants against pathogens? Is salicylic acid directly responsible for the increased resistance of plants to pathogens? Or is it mediated by rhizosphere microbiota? Or do they work together?
Salicylic acid remolding rhizosphere micro ora of watermelon is bene cial to plant antagonism to pathogens There were considerable researches have emphasized the role of soil microbial communities enhancing plant growth and health (Aleklett et al., 2018;Toju et al., 2018). For example, Xin's study showed that plant immune system was needed to maintain the normal growth of commensal bacteria in Arabidopsis (Xin et al., 2016). Some recent work reported that assemblages of host-speci c microbiomes in the rhizosphere. For instance, some disease-resistant crop varieties enrich speci c sets of bacterial species in the rhizosphere functions contribute to suppress pathogens (Badri et al., 2009;Delgado-Baquerizo et al., 2018). It is not known how plant roots normally select and maintain a healthy rhizosphere microbiota. SA plays important roles in regulating plant immunity which is necessary for systemic acquired resistance (SAR) (Trivedi et al., 2020). However, the mobile signals are still up in the air. Therefore, we assumed that SA may activate the immune signal in watermelon plants by mediating rhizosphere microbial community to cooperate enhancing resistance to disease.
Likewise, in our experiment, the Proteobacteria which is the second largest phylum of hydrogenogenic CO oxidizers (Badger and Bek, 2008;Wang and Sugiyama, 2020) had accumulated signi cantly at 3dpi after SA application. Furthermore, our results identi ed that there were signi cant increasing community abundance of Rhodanobacter at 3dpi and Azospirillum, Herbaspirillum, Stenotrophomonas and Sphingomonas at 7dpi after SA treatment on the genus level. Surprisedly, these genus all belong to Proteobacteria. Rhodanobacter for ammonia oxidation and denitri cation (De Clercq et al., 2006;Huo et al., 2018). Some members of the genus Azospirillum have biocontrol activity of phytopathogens which recognized as biofertilizers owing to their plant growth-promoting activities, such as the biologic nitrogen xation, hormones production, phosphate solubilization, and siderophore production (Mendez-Gomez et al., 2020). Sphingomonas is a kind of abundant new microbial resources, which can be used for biodegradation of aromatic compounds. It has great potential in environmental protection and industrial production because of its high metabolic capacity and multifunctional physiological characteristics (Vorholt et al., 2017). Stenotrophomonas maltophilia widely exists in water, soil and animals. It is an opportunistic pathogen. With the wide and large dose application of clinical antibiotics and immunosuppression (Hassan and Bano, 2016;Mendes et al., 2013;Messiha et al., 2007). Herbaspirillum species have been described in close association with plants, both endophytically and epiphytically because its nitrogenase activity to promote plant growth (Wang;. For instance, histochemical analysis of seedlings of maize, sorghum, wheat and rice grown in vermiculite showed that H. seropedicae colonizes root surfaces and inner tissues (Ramos et al., 2020). Collectively, our ndings suggest that the application of SA may recruit bene cial rhizosphere microbial community to active watermelon plant defense activity.
In order to nd out the molecular mechanism to explain how the rhizosphere microbial community cooperate resistance to watermelon Fusarium wilt, we blasted the unique genes to annotation their function abundances, visual display and statistical analysis. Surprisingly, after SA application at 3 days, the enrichment of KEGG pathways of environmental information processing, cellular processes and metabolism pathways were signi cantly activated. Furthermore, the K02014 (iron complex outer membrane receptor protein) had signi cantly increased in SA3 compared with SF3. At the same time, the genes in groups of Post translational; Intracellular tra cking, secretion and Cell motility and Cell wall: membrane envelope biogenesis had increased signi cantly in SA3. The results con rmed our hypos of rhizosphere microbiome assemblage is affected by SA signal (Chaparro et al., 2014;Delgado-Baquerizo et al., 2018;Mendes et al., 2011;Vorholt et al., 2017) and the 3dpi is a key time point for watermelon plants to activate disease resistance Ren et al., 2016;Xu et al., 2015). Notably, the resistance genes from Glycoside Hydrolases and Polysaccharide Lyases have highly accumulated in SA3, especially the result indicated that AAC6_IIC (drug class of aminoglycoside antibiotic, resistance mechanism of antibiotic inactivation) has a signi cant accumulation constantly after SA application compared with control group. Furthermore, we noticed that the relative abundance of Proteobacteria corresponding to resistance gene is the highest post inoculation which means more resistance genes may came from these kinds of species. In the same way, our results con rmed studies about that that plant genetics and agricultural practices can potentially impose selective pressures on speci c microbes and microbial communities. For instance, Xin's assumed that a plant genetic network guarding phyllosphere bacteria microbiota. Roeland identi ed that SA or JA systemically accumulates in response to infection by biotrophs or necrotrophy, respectively, and can affect rhizosphere microbiome assembly after foliar infection by pathogens in Arabidopsis (Chen et al., 2020). Thus, our genomic basis of plantmicrobe interactions results may inform the key role of rhizosphere microbial in resistance of watermelon to Fusarium wilt through several signal transduction pathways to activate a series of defensive feedback under environment stress.

Conclusion
In conclusion, this experiment demonstrates that the mechanism of how salicylic acid signal recruit soil microorganisms to explain their cooperate resistance to Fusarium wilt disease in watermelon plants. Collectively, our ndings suggest that the incidence of watermelon Fusarium Wilt was signi cantly lower than control group after SA application. Notably, our results indicated that the application of salicylic acid recruit soil microbial communities in watermelon plants. In addition, the soil microbial metabolic pathway and microbial function were signi cantly different after SA treatment which is bene cial to the antagonism of plants against pathogens. This work advances the current understanding of plant-microbe interaction research and supports the coevolution theory of mutualism between the plant and microbes. It is vital that we identify the microbial functions required for plant colonization in order to effectively implements them potentially in management and agricultural practices of plant healthy.