Connection the Rhizomicrobiome and Plant MAPK Gene Expression Response to Pathogenic Fusarium oxysporum in Wild and Cultivated Soybean

Article information

Plant Pathol J. 2019;35(6):623-634
Publication date (electronic) : 2019 December 12
doi : https://doi.org/10.5423/PPJ.OA.04.2019.0111
1Key Laboratory of Mollisols Agroecology, Northeast Institute of Geography and Agroecology, Chinese Academy of Sciences, Changchun 130102, China
2University of Chinese Academy of Sciences, Beijing 100049, China
3Key Laboratory of Vegetation Ecology, Ministry of Education, Institute of Grassland Science, Northeast Normal University, Changchun 130024, China
4College of Life Science, Jilin Agricultural University, Changchun 130118, China
*Corresponding author: Phone) +86-431-85542315, FAX) +86-431-85542298, E-mail) tiancj@neigae.ac.cn
Handling Editor : Kwak, Youn-Sig
Received 2019 April 21; Revised 2019 October 29; Accepted 2019 November 04.

Abstract

Little known the connections between soybeans mitogen-activated protein kinase (MAPK) gene expression and the rhizomicrobiome upon invasion of the root pathogen Fusarium oxysporum. To address this lack of knowledge, we assessed the rhizomicrobiome and root transcriptome sequencing of wild and cultivated soybean during the invasion of F. oxysporum. Results indicated F. oxysporum infection enriched Bradyrhizobium spp. and Glomus spp. and induced the expression of more MAPKs in the wild soybean than cultivated soybean. MAPK gene expression was positively correlated with Pseudomonadaceae but negatively correlated with Sphingomonadaceae and Glomeraceae in both cultivated and wild soybean. Specifically, correlation profiles revealed that Pseudomonadaceae was especially correlated with the induced expression of GmMAKKK13-2 (Glyma.14G195300) and GmMAPK3-2 (Glyma.12G073000) in wild and cultivated soybean during F. oxysporum invasion. Main fungal group Glomeraceae was positively correlated with GmMAPKKK14-1 (Glyma.18G060900) and negatively correlated with GmRaf6-4 (Glyma.02G215300) in the wild soybean response to pathogen infection; while there were positive correlations between Hypocreaceae and GmMAPK3-2 (Glyma.12G073000) and between Glomeraceae and GmRaf49-3 (Glyma.06G055300) in the wild soybean response, these correlations were strongly negative in the response of cultivated soybean to F. oxysporum. Taken together, MAPKs correlated with different rhizomicrobiomes indicating the host plant modulated by the host self-immune systems in response to F. oxysporum.

Most phytopathogenic Fusarium spp. are soil-borne and cause wilting of plants and rotting of roots, stems, leaves and fruits by damaging the vascular system and organs of various plant species (Dean et al., 2012; Garbeva et al., 2004). Therefore, disease caused by Fusarium spp. severely affects both the quality and quantity of different important crops, including soybean (Tu, 1987). Although tremendous progress has been made in recent decades in enhancing plant defence through various classical breeding and modern transgenic approaches, a challenge remains for plant scientists to control this devastating plant disease effectively.

Upon pathogen attack, host plants detect the presence of pathogen- and microbe-associated molecular patterns through plasma membrane-localized pattern recognition receptors, triggering pattern-triggered immunity (PTI), which is the first line of defence in innate immunity (Boller and He, 2009). Activation of PTI results in different immune responses, including mitogen-activated protein kinase (MAPK) cascade activation (Asai et al., 2002; Miya et al., 2007). Studies have shown that plant MAPKs are extremely important because they transduce PTI signals to downstream components (Chisholm et al., 2006; Dodds and Rathjen, 2010; Rodriguez et al., 2010; Segonzac et al., 2011). Therefore, a lot of study has been conducted to illustrate the role of MAPKs in plant-pathogen interactions. Notably, several members of the MAPK family have been revealed to be involved in plant defence responses processes (Zhang and Klessig, 2001).

Over 20 plant pathogenic fungi have been characterized for which MAPKs play a vital role in plant-pathogen interactions (Li et al., 2012; Turrà et al., 2014). For instance, upon pathogen attack, the model plant Arabidopsis thaliana induces the expression of three genes coding for the proteins MPK3, MPK4 and MPK6 (Pitzschke et al., 2009). Several studies have revealed that MAPKs are involved in Rhizobia-legume symbiosis; for example, alfalfa MAPK, called TDY1, is expressed in roots and nodules (Schoenbeck et al., 1999). In addition, MAPK can connect with other beneficial microbes, and a previous study revealed that a MAPK plays an essential role in signal transduction in the Trichoderma asperellum–cucumber system; this protein has been termed Trichoderma-induced MAPK (Shoresh et al., 2006). Although symbiotic and pathogenic interactions appear to be opposite model, common strategies have been found in the early plant response to infection by pathogenic and symbiotic bacteria, such as, MAPKs signal pathway. However, little is known about the interaction between host plant MAPKs and their corresponding root-associated microbiomes.

Soybean (Glycine max L.), a globally important oilseed crop, is a rich source of edible oil, proteins and nutraceutical compounds and originated from the wild soybean (Glycine soja Sieb. and Zucc.) from more 5,000 years ago in ancient China (Li et al., 2008). A previous study revealed that the symbiotic association of arbuscular mycorrhizal fungi with cultivated soybean enhanced plant drought tolerance, an effect that was positively correlated with the molecular dialogue of the plant MAPKs (Liu et al., 2015). Thus, it is reasonable to speculate that diverse microbial taxa might be present in the rhizosphere of wild soybean to connect with MAPKs. However, little information is available highlighting the central role played by MAPK signalling in the molecular dialogue between wild soybean and rhizomicrobiomes during the invasion of fungal pathogens. To obtain such information, in this study, we first aimed to explore the differences in the rhizomicrobiome and transcriptome of cultivated soybean and its wild relative (wild soybean) during Fusarium oxysporum infection. In addition, this study aimed to clarify whether the rhizomicrobiome and GmMAPK gene expression are tightly connected in the response to F. oxysporum invasion.

Materials and Methods

Experimental material

For the current study, we selected healthy seeds of wild soybean (Glycine soja, ZYD 01-289) and cultivated soybean (Glycine max, Williams 82) varieties. Briefly, sterilized the soybeans seeds for 5 min with 75% ethanol and then treated by 0.5% sodium hypochlorite for 1 min. Next, all soybean seeds were rinsed three times with sterilized water. Sterilized seeds were then incubated in the dark at 28°C in an incubator until germination (about 3 days). Germinated seeds of wild and/or cultivated soybean were planted in pots (18 cm dia. × 15 cm high), containing 1 kg autoclaved soil. Soybean seedlings were then grown in a growth chamber with 75% relative humidity with a 16 h light period at 28°C and an 8 h dark period at 25°C. Three replicate pots of each genotype were used. The wild and cultivated soybeans were watered with sterilized water every other day.

The soil used in the pots was collected at the Changchun Agricultural Station of the Northeast Institute of Geography and Agroecology in Jilin Province, China (125°23′44″E, 43°59′58″N) in July 2016. The soil was sieved and homogenized, and 1,500 g of the soil was stored in a plastic bag and autoclaved for 60 min twice. The soil physicochemical characteristics, including soil organic matter, total N, available N, P, and K, were 32.8 g/kg, 651.9 mg/kg, 109.2 mg/kg, 7.5 mg/kg, and 88.7 mg/kg, respectively, and the pH value was 7.0.

Fungal strain, culture conditions, and infection processes

Fusarium oxysporum Schltdl. is one of the phytopathogenic fungal species most frequently associated with soybean root rot (Jimenez, 2017). The fungal pathogen F. oxysporum was originally isolated from soybean roots. Molecular identification was performed based on fungal internal transcribed spacer (ITS) region fragment amplification by using the universal ITS primers with ITS1 and ITS4 (White et al., 1990). Upon the appearance of trifoliate leaves, the seedlings of both soybean genotypes were inoculated with an F. oxysporum spore suspension (1.0 × 107 spores/ml) according to the method previously described (Scandiani et al., 2015). We then tested and observed the phenotype of wild and cultivated soybean without or with pathogen invasion after inoculated F. oxysporum 7 days ago (Chang et al. 2018).

Collection and storage of rhizosphere soils

Rhizosphere soils were collected from wild and cultivated soybean roots 7 days after inoculation with F. oxysporum. The soil loosely adhering to the roots was collected by gently shaking the roots and was stored at −80°C until DNA extraction. In addition, wild and cultivated soybean roots were harvested, immediately frozen in liquid nitrogen after washing, and then stored at −80°C for RNA extraction. For each treatment, the roots of five soybean seedlings were collected and pooled together, with 3 biological samples per treatment (n = 3).

Metabarcoding of the bacterial V3–V4 region and the fungal ITS2 region

Total soil DNA was extracted with a FastDNA SPIN Kit (MPBio, Santa Ana, CA, USA) according to the manufacturer’s instructions and quantified with a NanoDrop 2000 spectrophotometer (NanoDrop Technologies, Inc., Wilmington, DE, USA). The V3–V4 region of the bacterial 16S rRNA gene was amplified using the primers 341F and 806R (Luo et al., 2018). The fungal ITS2 region was amplified by the primers ITS3F and ITS4R (Zhou et al., 2017). Both the bacterial and fungal primer pairs were ligated with Illumina adapter sequences and barcodes (Caporaso et al., 2010). PCR was performed in a 30 μl mixture containing 3 μl each of primer (2 μM), 10 μl of template DNA (1 ng/μl), 15 μl of Phusion High-Fidelity PCR Master Mix (BioLabs, Inc., LA, USA) and 2 μl of water. The following thermal programme was used for amplification: 95°C for 1 min followed by 30 cycles of 98°C for 10 s, 50°C for 30 s, and 72°C for 30 s and a final extension step at 72°C for 5 min. PCR was repeated four times for each sample, and all amplicons were pooled and purified using a Qiagen Gel Extraction Kit (Qiagen, Hilden, Germany). The sequencing libraries were generated using a TruSeq DNA PCR-Free Sample Preparation Kit (Illumina, San Diego, CA, USA) according to the manufacturer’s instructions and pooled at an equimolar ratio. The Illumina HiSeq 2500 platform at ORI-GENE Biotech Company was used for 250 bp paired-end (PE) sequencing.

The raw sequencing data were de-multiplexed, quality-filtered and processed by FLASH (Magoč and Salzberg, 2011) as the following steps: (1) chimeric sequences were identified and removed using USEARCH 6.1 in QIIME (Caporaso et al., 2010), and all sequences shorter than 200 bp in length, with mismatching barcodes, with more than two base pair mismatches, with ambiguous bases, or with an average quality score < 25 in the raw reads were discarded; and (2) UPARSE (version 7.0.1001, http://drive5.com/uparse) was used to cluster quality-trimmed sequences at a 97% sequence similarity cut-off and to form operational taxonomic units (OTUs). (3) The representative sequence of each OTU was taxonomically classified by using the Ribosomal Database Project (RDP) Classifier software (version 2.2, http://sourceforge.net/projects/rdp-classifier/). The reference databases at this step were used Silva for bacteria communities (DeSantis et al., 2006) and the UNITE databases for fungal communities, which precluster at 97% similarity cutoff (Abarenkov et al., 2010). OTUs belonging to a functional group were identified using the Functional Annotation of Prokaryotic Taxa database (FAPROTAX database) (Louca et al., 2016).

RNA-sequencing (RNA-seq) and differential genes expression profiling

Root samples (100 mg each) were ground into a fine powder in liquid N2, and total RNA was isolated by using a Promega RNA extraction kit (LS1040, Promega, China) following the manufacturer’s instructions. The extracted RNA was treated with RNase-Free DNase (Qiagen) to digest and remove any contaminating genomic DNA. RNA quantity and quality were assessed by using a NanoDrop spectrometer (NanoDrop 2000, NanoDrop Technologies, Inc.) and electrophoresis on a 1% (wt/vol) agarose gel. The RNA samples were then stored at −80°C until subsequent use. cDNA library construction and RNA-seq were performed by the ORI-GENE Biotech Company (Beijing, China) using the Illumina PE sequencing method following the Illumina protocol. PE sequencing (2 × 100 bp) was carried out using the Illumina HiSeq 2000 platform. Reads that were not aligned with the soybean reference genome (e.g., fungal reads) were filtered out and were not included in further analyses. Gene expression (FPKM, Fragments Per Kilobase of transcript Per Million fragments mapped) levels were estimated using Cufflinks software (version 2.1.1) (Trapnell et al., 2012), and differentially expressed genes (DEGs) were selected by using the criteria q < 0.05 and |log2(fold change)| ≥ 2.

Statistical analyses

The OTU abundance matrix was rarefied before performing downstream analysis. The functional microbial communities involved in carbon and nitrogen metabolism were annotated by the FAPROTAX database (Louca et al., 2016). The correlation network was visualized using Cytoscape (Smoot et al., 2010). Statistically significant differences in the expression of each gene among the treatments were assessed by SAS version 9.1 (SAS Institute, Cary, NC, USA) with Tukey’s honestly significant difference test at the P < 0.05 level.

Availability of data and materials

The raw sequencing data for the microbe diversity and transcriptome were submitted to the Sequence Read Archive (SRA) at the NCBI of the accession numbers SRP079341 and SRP111527, respectively.

Results

In total, 4 treatments, the non-inoculated wild (W) and cultivated (C) soybeans and the inoculated wild (W + F) and cultivated (C + F) soybeans, were enrolled and included in the metagenomics sequencing analyses (Fig. 1). For the metagenomics analysis, a total of 245,263 qualified bacterial reads were obtained after filtering the original 295,190 raw reads (ranging from 12,060–35,975 bp), resulting in 2,584 total OTUs (ranging from 584–1,313 bp) (Supplementary Table 1). A total of 149,929 qualified reads out of 156,955 raw sequences were obtained (ranging from 6,472–19,814 bp), resulting in 1,279 total OTUs (ranging from 72–150 bp) during the ITS region sequencing (Supplementary Table 1). Rarefaction and Shannon curves indicated that the sequencing coverage was valid to quantify most of the species and OTU numbers produced by Illumina sequencing (Supplementary Fig. 1).

Fig. 1

Experimental design for examining the plant-microbe interactions of wild and cultivated soybean in response to Fusarium oxysporum. Soybean cultivars were selected from wild (01–289) and cultivated (Williams 82) soybean. MAPK, mitogen-activated protein kinase.

Taxonomic characteristics of bacterial and fungal communities

Overall, the 16S rRNA metabarcodes were classified into 28 bacterial phyla, and 2.28 ± 0.65% of the bacterial sequences were unclassified at the phylum level. The dominant bacterial phyla (relative abundance > 1%) accounted for more than 90% of the abundance and included Proteobacteria (mean ± SD, 48.06 ± 2.20% of the total sequences), Cyanobacteria (34.28 ± 2.72%), Actinobacteria (4.38 ± 1.02%), Bacteroidetes (3.93 ± 1.32%), Planctomyces (2.25 ± 0.59%), and Verrucomicrobia (1.30 ± 0.31%). The groups Firmicute, Chloroflexi, Acidobacteria, Armatimonadetes, Candidate_division_TM7, Candidate_ division_OD1, Chlamydiae, and Gemmatimonadetes were less abundant (relative abundance > 0.5%) but were still found in all soil samples (Fig. 2A). In addition, the relative abundance of Bradyrhizobium spp. was significantly higher in the W + F group than in the C + F group, but there were no significant (P < 0.05) differences between the W and W + F groups (Fig. 2B). Interestingly, upon analysis with the FAPROTAX database, bacteria involved in nitrogen fixation were most abundant in the W + F group (Supplementary Fig. 2).

Fig. 2

The relative abundance of bacteria (A, B) and fungi (C, D) in wild and cultivated soybean in response to Fusarium oxysporum. Values (mean ± standard error) with different letters are significantly different at P < 0.05 (Tukey’s honestly significant difference test). The four treatments were non-inoculated wild soybean (W), non-inoculated cultivated soybean (C), wild soybean inoculated with F. oxysporum (W + F), and cultivated soybean inoculated with F. oxysporum (C + F). Panels B and D represent the different relative abundance levels of Proteobacteria and Glomeromycota, respectively.

The most dominant fungal phyla were Ascomycota (ranging from 31.67–82.45%), Basidiomycota (0.23–8.31%) and Glomeromycota (0.25–8.15%), accounting for more than 70% of the fungal community (Fig. 2C). Chytridiomycota (0.53–0%) and Zygomycota (0.03–0%) were less abundant. The relative abundance of Zygomycota was recovered only in the C and C + F groups. More interestingly, Glomeromycota was more highly enriched in the W + F group than in the other three groups. Furthermore, taxonomical classification at the genus level detected 138 genera, among which the genus Fusarium was dominant across all the soil samples (Fig. 2D).

Host plant transcriptome profile analysis of wild and cultivated soybean responses to F. oxysporum

In total, transcriptional profiling obtained 49.28, 37.36, 49.28, and 45.57 million clean reads from the W, C, W + F, and C + F groups, respectively. The average error rate of the sequences ranged from 0.0128–0.0132%, and more than 96% of the bases had error rates < 0.1% (Supplementary Table 2) (Chang et al., 2018). The short read sequences were blasted against the soybean reference genome (Wm82.a2.v1), and the total mapping rates ranged from 54.96% to 74.76%. Significant DEGs were determined by using a |fold change| ≥ 2 and a q-value < 0.05. Upon using these criteria, a total of 1,060, 1,372, 5,300, and 6,505 DEGs were identified in the C + F vs. C, W + F vs. W, W + F vs. C + F and W vs. C comparisons, of which 211, 614, 2,344, and 2,810 DEGs were upregulated and 849, 758, 2,956, and 3,695 were downregulated, respectively (Chang et al., 2018). Seven DEGs were selected for validation of the RNA-seq data using quantitative reverse transcription polymerase chain reaction (RT-qPCR) (Supplementary Table 3). The RT-qPCR data for these genes were consistent with the RNA-seq results from the four group samples, which indicated a high degree of reproducibility between the DEG fold changes assayed using RNA-seq and the expression profiles revealed by RT-qPCR (Chang et al., 2018).

Expression profiles of the MAPK signalling pathway in soybean

Thirty-eight GmMAPK, 11 GmMAPKK, and 150 GmMAPKKK orthologous genes were identified across the published data sources from Arabidopsis, rice and poplar (Neupane et al., 2013). In the present study, the transcriptome data showed that 35 MAPK genes were differentially expressed in cultivated and wild soybean in response to F. oxysporum (Fig. 3). More specifically, among the genes of the MAPK cascade pathway 6 MAPK-related DEGs (4 up- and 2 downregulated), were observed in the W + F vs. W comparison, while 4 downregulated genes were observed in the C+F vs. C comparison (Fig. 3, Supplementary Table 4). In addition, Glyma.14G099300/MAPK3-2, a Glyma.14G099300 gene encoding an ACT-like protein tyrosine kinase family protein, was downregulated in the C + F vs. C comparison (Fig. 3, Supplementary Table 4). Furthermore, 23 MAPK-related DEGs (11 up- and 12 downregulated) and 21 DEGs (14 up- and 7 downregulated) were observed in the W vs. C and W + F vs. C + F comparisons, respectively (Fig. 3, Supplementary Table 4). The protein kinase-related gene MAPK 4 (MPK4; Glyma.01G222000) was strongly upregulated in the W + F vs. C + F comparison, while there was no significant difference in the W vs. C comparison (Fig. 3, Supplementary Table 4). We assumed that these transcriptional changes stimulated protein changes in both wild and cultivated soybeans in response to F. oxysporum infection.

Fig. 3

Differentially expressed genes related to the mitogenactivated protein kinase (MAPK) family in wild and cultivated soybean in response to Fusarium oxysporum infection. The four treatments were non-inoculated wild soybean (W), non-inoculated cultivated soybean (C), wild soybean inoculated with Fusarium oxysporum (W + F) and cultivated soybean inoculated with F. oxysporum (C + F). The red, blue and white colours represent the upregulated (fold changes ≥ 4 with a q-value < 0.05), downregulated (fold changes ≤ −4 with a q-value < 0.05) and unchanged genes, respectively. Additional information pertaining to MAPK signalling pathway-related genes is presented in Supplementary Table 2.

Correlation analysis between the rhizomicrobiome and plant MAPK gene expression

In Spearman correlation analysis, between abundance of rhizosphere microbial communities (bacterial and fungal communities) and the level of plant MAPK-related DEGs with correlation coefficients of R > 0.75 or R < −0.75 with P-values < 0.05 were considered to be correlated with each MAPK (Supplementary Table 5). Transcriptional co-expression networks revealed correlations between 35 MAPK-related DEGs and 66 bacterial microbes (consisting of 246 OTUs) (Fig. 4, Supplementary Table 5). More specifically, 44 positive correlations and 75 negative correlations were observed in the W + F vs. W comparison, whereas 21 positive correlations and 72 negative correlations were observed in the C + F vs. C comparison. For instance, Pseudomonadaceae and Streptomycetaceae were positively correlated with GmMAPKKK 13-2 (Glyma.14G195300) in the W + F vs. W comparison, while Planctomycetaceae was strongly negatively correlated with GmRaf6-4 (Glyma.02G215300), GmRaf35-2 (Glyma.17G065700) and GmRaf53-2 (Glyma.20G165800) in the W + F vs. W comparison (Fig. 4, Supplementary Table 5). In the C + F vs. C comparison, Pseudomonadaceae and Planctomycetaceae were both positively correlated with GmMAPK3-2 (Glyma.12G073000), and Streptomycetaceae was strongly negatively correlated with GmRaf49-3 (Glyma.06G055300) (Fig. 4, Supplementary Table 5). In addition, 228 positive correlations and 233 negative correlations were found in the W vs. C comparison, whereas 184 positive correlations and 306 negative correlations were found in the W + F vs. C + F comparison (Supplementary Table 5).

Fig. 4

Correlations between bacterial taxa and the mitogen-activated protein kinase (MAPK) family in wild and cultivated soybean in response to Fusarium oxysporum. (A) Wild soybean inoculated with F. oxysporum (W + F) vs. non-inoculated wild soybean (W). (B) Cultivated soybean inoculated with F. oxysporum (C + F) vs. non-inoculated cultivated soybean (C). (C) Wild soybean (W) vs. cultivated soybean (C). (D) Wild soybean inoculated with F. oxysporum (W + F) vs. cultivated soybean inoculated with F. oxysporum (C + F). The nodes represent the operational taxonomic units, and the colourful nodes indicate the microbial taxa at the phylum level. The connecting lines indicate a positive (red line) or a negative (blue line) correlation between the bacterial taxa and the MAPK genes. Additional information pertaining to the correlations is presented in Supplementary Table 4.

Regarding the fungal community and MAPKs, correlations existed between 34 MAPK-related DEGs and 8 fungal microbes at the family level (consisting of 9 OTUs), with 18 positive and 16 negative correlations in the Spearman analysis (Fig. 5, Supplementary Table 5). 15 positive and 39 negative correlations were found in the W + F vs. W comparison, while 6 positive and 11 negative correlations were found in the C + F vs. C comparison. For instance, positive correlations between Glomeraceae and GmMAPKKK14-1 (Glyma.18G060900) and between Chaetomiaceae and GmRaf53-2 (Glyma.20G165800) were found in the W + F vs. W comparison, and negative correlations between Glomeraceae and GmRaf6-4 (Glyma.02G215300), between Hypocreales_family_ Incertae_sedis and GmRaf35-2 (Glyma.17G065700), and between Hypocreales_family_Incertae_sedis and GmRaf53-2 (Glyma.20G165800) were also observed. In addition, in the C + F vs. C comparison, Hypocreaceae and Ceratobasidiaceae were positively correlated with GmMAPK3-2 (Glyma.12G073000), while Glomeraceae was strongly negatively correlated with GmRaf49-3 (Glyma.06G055300), GmRaf49-1 (Glyma.14G099300), and GmRaf49-2 (Glyma.17G225400). Furthermore, 39 positive and 40 negative correlations were found in the W vs. C comparison, while 34 positive and 41 negative correlations were found in the W + F vs. C + F comparison.

Fig. 5

Correlations between fungal taxa and the mitogen-activated protein kinase (MAPK) family in wild and cultivated soybean in response to Fusarium oxysporum. (A) Wild soybean inoculated with F. oxysporum (W + F) vs. non-inoculated wild soybean (W). (B) Cultivated soybean inoculated with F. oxysporum (C + F) vs. non-inoculated cultivated soybean (C). (C) Wild soybean (W) vs. cultivated soybean (C). (D) Wild soybean inoculated with F. oxysporum (W + F) vs. cultivated soybean inoculated with F. oxysporum (C + F). The nodes represent operational taxonomic units, and the colourful nodes indicate the microbial taxa at the phylum level. Connecting lines indicate a positive (red line) or a negative (blue line) correlation between the bacterial taxa and the MAPK genes. Additional information pertaining to the correlations is presented in Supplementary Table 5.

Discussion

In this study, the functional associations between the rhizomicrobiome and the host plant immune system were determined by analysing both the host root-associated microbiome and the transcriptome in an investigation of the root microbes of soybeans (G. max and G. soja) and their potential impact on the MAPK signalling pathway in response to F. oxysporum infection (Fig. 1). MAPK signalling has been demonstrated to regulate cellular processes involved in activating host defence responses against attacking pathogens (Andreasson and Ellis, 2010). Little is known thus far about whether MAPKs are active in the plant-symbiont interaction; however, certain rhizobia have the capacity to interfere with and even activate the host MAPK pathway (Bartsev et al., 2004; Fernandez-Pascual et al., 2006). During the arbuscular mycorrhiza (AM)-plant interaction, the only relevant data in the literature indicate that there is increased transcription of a single MAPK gene during the pre-contact and the appressorium stages between the barrel medic Medicago truncatula and Glomus mosseae (Weidmann et al., 2004). In our study, correlation analyses between diverse microbes and the MAPK signalling pathway revealed the answer to a key question about the response of wild and cultivated soybeans to F. oxysporum by indicating the possible interactions between the microbiome and MAPKs (Figs. 4 and 5). However, sequencing analyses explore only the taxonomical composition of the plant microbiome without taking into account the whole bacterial genomes or addressing the functions these microbes are performing.

Previous studies have fully demonstrated that different genotypes or varieties of the same species can develop distinct microbial communities, as has been observed for Arabidopsis, rice and maize (Edwards et al., 2015; Micallef et al., 2009; Peiffer et al., 2013). In addition, our study revealed that greater abundance of rhizosphere bacterial and fungal microbe groups occurred in wild soybean than in cultivated soybean, especially in response to F. oxysporum (Fig. 2). Interestingly, we found that wild soybean possessed more Bradyrhizobium spp. and Glomus spp. than cultivated soybean (Fig. 2B and D), and specific microbiota involved in carbon and nitrogen metabolism, such as Bradyrhizobium spp., were enriched in wild soybean (Supplementary Fig. 2). Wild soybean is a source of genetic diversity to identify novel genes underlying responses to biotic and abiotic stresses (Kim et al., 2011; Zhang et al., 2016). In our study, the RNA-seq results revealed more DEGs in the W + F vs. W comparison than in the C + F vs. C comparison, which implied that wild soybean exhibited a stronger response to F. oxysporum infection (Fig. 3). Furthermore, our study revealed more MAPK-related DEGs in the W + F vs. W comparison than in the C + F vs. C comparison (Fig. 3). Previous studies on crops such as rice have revealed that wild relative crops have a stronger ability to resist pathogen infection than cultivated varieties (Tian et al., 2018).

A better understanding of the various interactions between microbes and host plants may lay a valuable foundation of knowledge. To our knowledge, this study is the first to reveal the correlations between the rhizomicrobiome (bacterial and fungal) and soybean transcriptional expression; most previous studies have focused on specific plant microbiomes and have placed more emphasis on microbial diversity than on gene function (Bulgarelli et al., 2015; Ofek-Lalzar et al., 2013). In our study, the correlations revealed that root-associated microbiota mediated MAPK transcriptional expression; however, the microbiota-dependent MAPKs were different in wild and cultivated soybeans under F. oxysporum infection (Figs. 4 and 5). For example, Pseudomonadaceae was positively correlated with GmMAPKKK 13-2 (Glyma.14G195300) in the W + F vs. W comparison but was strongly negatively correlated with GmMAPK3-2 (Glyma.12G073000) in the C + F vs. C comparison. Pseudomonas spp. have been identified as plant growth-promoting rhizosphere bacteria and can elicit induced systemic resistance (Pieterse et al., 2014). However, the majority of OTUs detected in our investigation have unknown functions, making it difficult to correlate the microbial community with MAPK profiles. Future analysis of knock-out mutants and ectopic overexpression of the candidate genes will be performed to examine their role in the activity of the root-associated microbiota.

As it is known, the seed endophytic are often defined as bacteria or fungi that lived in plant tissues without expressing symptoms or visible signs (Jimtha et al., 2014). Lundberg et al. (2012) studied on the core root microbiome reveal that the seed endophytic microbial diversity is less than the rhizomicrobiome (Lundberg et al., 2012), however, some strains can be spread out into the rhizosphere and soil, which can be induced the host plant’s defence mechanisms. Plant-endophyte interactions were studied in different defence reactions that indicated the seed endophytes played a vital role in host plant defence mechanisms (Rosenblueth and Martínez-Romero, 2006). In our study the soybean seeds surface were sterilized, but some fewer microbiome were still remained on the seed surfaces. Rhizosphere microbiome tightly attach to plant roots, thus the rhizomicrobiome extracted from plant tissue possibly contained surface bacteria (Edwards et al., 2015; Lundberg et al., 2012). Future study should be stressed on the possibility of microorganisms from seed endophytic, and make a fully elucidation on the relationship between the seed endophytic and rhizomicrobiome with MAPK signalling pathway.

Conclusion

In summary, our study provides a preliminary analysis of the associations between the rhizomicrobiome and MAPK expression in wild and cultivated soybean responses to F. oxysporum. In this study, wild and cultivated soybean presented differential correlations with MAPK and rhizomicrobiome; hence, to better utilize soybean’s transcriptional information is important to develop multiple disease resistant cultivars for soybean. This specificity of the effect of the rhizomicrobiome on soybean transcriptional networks may result from host plant-rhizomicrobiome coevolution.

Acknowledgements

This research was financially supported by the National Key Research and Development Program of China (2016YFC0501202), Science Foundation of Chinese Academy of Sciences (XDB15030103, XDA23070501), the Key Research Program of CAS (KFZD-SW-112), the National Natural Science Foundation of China (41571255), Cooperative Project between CAS and Jilin Province of China (2019SYHZ0039), the Science and Technology Development Project of Changchun City of China (18DY019) and the Science and Technology Development Project of Jilin Province of China (20180519002JH, 20190303070SF). We thank the China Scholarship Council for support to CCL.

Notes

Electronic Supplementary Material

Supplementary materials are available at The Plant Pathology ournal website (http://www.ppjonline.org/).

References

Abarenkov K., Henrik Nilsson R., Larsson K.H., Alexander I.J., Eberhardt U., Erland S., Høiland K., Kjøller R., Larsson E., Pennanen T., Sen R., Taylor A.F., Tedersoo L., Ursing B.M., Vrålstad T., Liimatainen K., Peintner U., Kõljalg U.. 2010;The UNITE database for molecular identification of fungi: recent updates and future perspectives. New Phytol 186:281–285. 10.1111/j.1469-8137.2009.03160.x. 20409185.
Andreasson E., Ellis B.. 2010;Convergence and specificity in the Arabidopsis MAPK nexus. Trends Plant Sci 15:106–113. 10.1016/j.tplants.2009.12.001. 20047850.
Asai T., Tena G., Plotnikova J., Willmann M.R., Chiu W.-L., Gomez-Gomez L., Boller T., Ausubel F.M., Sheen J.. 2002;MAP kinase signalling cascade in Arabidopsis innate immunity. Nature 415:977–983. 10.1038/415977a. 11875555.
Bartsev A.V., Deakin W.J., Boukli N.M., McAlvin C.B., Stacey G., Malnoë P., Broughton W.J., Staehelin C.. 2004;NopL, an effector protein of Rhizobium sp. NGR234, thwarts activation of plant defense reactions. Plant Physiol 134:871–879. 10.1104/pp.103.031740. 14966249. 344561.
Boller T., He S.Y.. 2009;Innate immunity in plants: an arms race between pattern recognition receptors in plants and effectors in microbial pathogens. Science 324:742–744. 10.1126/science.1171647. 19423812. 2729760.
Bulgarelli D., Garrido-Oter R., Münch P.C., Weiman A., Dröge J., Pan Y., McHardy A.C., Schulze-Lefert P.. 2015;Structure and function of the bacterial root microbiota in wild and domesticated barley. Cell Host Microbe 17:392–403. 10.1016/j.chom.2015.01.011. 25732064. 4362959.
Caporaso J.G., Kuczynski J., Stombaugh J., Bittinger K., Bushman F.D., Costello E.K., Fierer N., Peña A.G., Goodrich J.K., Gordon J.I., Huttley G.A., Kelley S.T., Knights D., Koenig J.E., Ley R.E., Lozupone C.A., Mc-Donald D., Muegge B.D., Pirrung M., Reeder J., Sevinsky J.R., Turnbaugh P.J., Walters W.A., Widmann J., Yatsunenko T., Zaneveld J., Knight R.. 2010;QIIME allows analysis of high-throughput community sequencing data. Nat Methods 7:335–336. 10.1038/nmeth.f.303. 20383131. 3156573.
Chang C., Tian L., Ma L.N., Li W., Nasir F., Li X., Tran L.-S.P., Tian C.. 2018;Differential responses of molecular mechanims and physiochemical characters in wild and cultivated soybeans against invasion by the pathogenic Fusarium oxysporum Schltdl. Physiol Plant 166:1008–1025. 10.1111/ppl.12870. 30430602.
Chisholm S.T., Coaker G., Day B., Staskawicz B.J.. 2006;Host-microbe interactions: shaping the evolution of the plant immune response. Cell 124:803–814. 10.1016/j.cell.2006.02.008. 16497589.
Dean R., Van Kan J.A.L., Pretorius Z.A., Hammond-Kosack K.E., Di Pietro A., Spanu P.D., Rudd J.J., Dickman M., Kahmann R., Ellis J., Foster G.D.. 2012;The top 10 fungal pathogens in molecular plant pathology. Mol Plant Pathol 13:414–430. 10.1111/j.1364-3703.2011.00783.x. 22471698. 6638784.
DeSantis T.Z., Hugenholtz P., Larsen N., Rojas M., Brodie E.L., Keller K., Huber T., Dalevi D., Hu P., Andersen G.L.. 2006;Greengenes, a chimera-checked 16S rRNA gene database and workbench compatible with ARB. Appl Environ Microbiol 72:5069–5072. 10.1128/AEM.03006-05. 16820507. 1489311.
Dodds P.N., Rathjen J.. 2010;Plant immunity: towards an integrated view of plant-pathogen interactions. Nat Rev Genet 11:539–548. 10.1038/nrg2812. 20585331.
Edwards J., Johnson C., Santos-Medellín C., Lurie E., Podishetty N.K., Bhatnagar S., Eisen J.A., Sundaresan V.. 2015;Structure, variation, and assembly of the root-associated microbiomes of rice. Proc Natl Acad Sci USA 112:E911–E920. 10.1073/pnas.1414592112. 25605935. 4345613.
Fernandez-Pascual M., Lucas M.M., de Felipe M.R., Boscá L., Hirt H., Golvano M.. 2006;Involvement of mitogen-activated protein kinases in the symbiosis Bradyrhizobium–Lupinus. J Exp Bot 57:2735–2742. 10.1093/jxb/erl038. 16868044.
Garbeva P., van Veen J.A., van Elsas J.D.. 2004;Microbail diversity in soil: selection of microbial populations by plant and soil type and implications for disease suppressiveness. Annu Rev Phytopathol 42:243–270. 10.1146/annurev.phyto.42.012604.135455. 15283667.
Jimenez D.R.C.. 2017. Soybean root rot caused by Fusarium oxysporum and Fusarium graminearum: interactions with biotic and abiotic factors. PhD thesis Iowa State University; Ames, IA, USA:
Jimtha J.C., Smitha P.V., Anisha C., Deepthi T., Meekha G., Radhakrishnan E.K., Gayatri G.P., Remakanthan A.. 2014;Isolation of endophytic bacteria from embryogenic suspension culture of banana and assessment of their plant growth promoting properties. Plant Cell Tissue Organ Cult 118:57–66. 10.1007/s11240-014-0461-0.
Kim M., Hyten D.L., Niblack T.L., Diers B.W.. 2011;Stacking resistance alleles from wild and domestic soybean sources improves soybean cyst nematode resistance. Crop Sci 51:934–943. 10.2135/cropsci2010.08.0459.
Li G., Zhou X., Xu J.-R.. 2012;Genetic control of infection-related development in Magnaporthe oryzae. Curr Opin Microbiol 15:678–684. 10.1016/j.mib.2012.09.004. 23085322.
Li Y., Guan R., Liu Z., Ma Y., Wang L., Li L., Lin F., Luan W., Chen P., Yan Z., Guan Y., Zhu L., Ning X., Smulders M.J., Li W., Piao R., Cui Y., Yu Z., Guan M., Chang R., Hou A., Shi A., Zhang B., Zhu S., Qiu L.. 2008;Genetic structure and diversity of cultivated soybean (Glycine max (L.) Merr.) landraces in China. Theor Appl Genet 117:857–871. 10.1007/s00122-008-0825-0. 18587557.
Liu Z., Li Y., Ma L., Wei H., Zhang J., He X., Tian C.. 2015;Coordinated regulation of arbuscular mycorrhizal fungi and soybean MAPK pathway genes improved mycorrhizal soybean drought tolerance. Mol Plant-Microbe Interact 28:408–419. 10.1094/MPMI-09-14-0251-R. 25390189.
Louca S., Parfrey L.W., Doebeli M.. 2016;Decoupling function and taxonomy in the global ocean microbiome. Science 353:1272–1277. 10.1126/science.aaf4507. 27634532.
Lundberg D.S., Lebeis S.L., Paredes S.H., Yourstone S., Gehring J., Malfatti S., Tremblay J., Engelbrektson A., Kunin V., del Rio T.G., Edgar R.C., Eickhorst T., Ley R.E., Hugenholtz P., Tringe S.G., Dangl J.L.. 2012;Defining the core Arabidopsis thaliana root microbiome. Nature 488:86–90. 10.1038/nature11237. 22859206. 4074413.
Luo S., Tian L., Chang C., Wang S., Zhang J., Zhou X., Li X., Tran L.-S.P., Tian C.. 2018;Grass and maize vegetation systems restore saline-sodic soils in the Songnen Plain of northeast China. Land Degrad Dev 29:1107–1119. 10.1002/ldr.2895.
Magoč T., Salzberg S.L.. 2011;FLASH: fast length adjustment of short reads to improve genome assemblies. Bioinformatics 27:2957–2963. 10.1093/bioinformatics/btr507. 21903629. 3198573.
Micallef S.A., Shiaris M.P., Colón-Carmona A.. 2009;Influence of Arabidopsis thaliana accessions on rhizobacterial communities and natural variation in root exudates. J Exp Bot 60:1729–1742. 10.1093/jxb/erp053. 19342429. 2671628.
Miya A., Albert P., Shinya T., Desaki Y., Ichimura K., Shirasu K., Narusaka Y., Kawakami N., Kaku H., Shibuya N.. 2007;CERK1, a LysM receptor kinase, is essential for chitin elicitor signaling in Arabidopsis. Proc Natl Acad Sci USA 104:19613–19618. 10.1073/pnas.0705147104. 18042724. 2148337.
Neupane A., Nepal M.P., Piya S., Subramanian S., Rohila J.S., Reese R.N., Benson B.V.. 2013;Identification, nomenclature, and evolutionary relationships of mitogen-activated protein kinase (MAPK) genes in soybean. Evol Bioinform Online 9:363–386. 10.4137/EBO.S12526. 24137047. 3785387.
Ofek-Lalzar M., Sela N., Goldman-Voronov M., Green S.J., Hadar Y., Minz D.. 2013;Niche and host-associated functional signatures of the root surface microbiome. Nat Commun 5:4950. 10.1038/ncomms5950. 25232638.
Peiffer J.A., Spor A., Koren O., Jin Z., Tringe S.G., Dangl J.L., Buckler E.S., Ley R.E.. 2013;Diversity and heritability of the maize rhizosphere microbiome under field conditions. Proc Natl Acad Sci USA 110:6548–6553. 10.1073/pnas.1302837110. 23576752. 3631645.
Pieterse C.M.J., Zamioudis C., Berendsen R.L., Weller D.M., Van Wees S.C.M., Bakker P.A.H.M.. 2014;Induced systemic resistance by beneficial microbes. Annu Rev Phytopathol 52:347–375. 10.1146/annurev-phyto-082712-102340. 24906124.
Pitzschke A., Schikora A., Hirt H.. 2009;MAPK cascade signalling networks in plant defence. Curr Opin Plant Biol 12:421–426. 10.1016/j.pbi.2009.06.008. 19608449.
Rodriguez M.C., Petersen M., Mundy J.. 2010;Mitogen-activated protein kinase signaling in plants. Annu Rev Plant Biol 61:621–649. 10.1146/annurev-arplant-042809-112252. 20441529.
Rosenblueth M., Martínez-Romero E.. 2006;Bacterial endophytes and their interactions with hosts. Mol Plant-Microbe Interact 19:827–837. 10.1094/MPMI-19-0827. 16903349.
Scandiani M.M., Luque A.G., Razori M.V., Ciancio Casalini L., Aoki T., O’Donnell K., Cervigni G.D.L., Spampinato C.. 2015;Metabolic profiles of soybean roots during early stages of Fusarium tucumaniae infection. J Exp Bot 66:391–402. 10.1093/jxb/eru432. 25336687.
Schoenbeck M.A., Samac D.A., Fedorova M., Gregerson R.G., Gantt J.S., Vance C.. 1999;The alfalfa (Medicago sativa) TDY1 gene encodes a mitogen-activated protein kinase homolog. Mol Plant-Microbe Interact 12:882–893. 10.1094/MPMI.1999.12.10.882. 10517028.
Segonzac C., Feike D., Gimenez-Ibanez S., Hann D.R., Zipfel C., Rathjen J.P.. 2011;Hierarchy and roles of pathogen-associated molecular pattern-induced responses in Nicotiana benthamiana. Plant Physiol 156:687–699. 10.1104/pp.110.171249. 21478366. 3177268.
Shoresh M., Gal-On A., Leibman D., Chet I.. 2006;Characterization of a mitogen-activated protein kinase gene from cucumber required for trichoderma-conferred plant resistance. Plant Physiol 142:1169–1179. 10.1104/pp.106.082107. 16950863. 1630744.
Smoot M.E., Ono K., Ruscheinski J., Wang P.-L., Ideker T.. 2010;Cytoscape 2.8: new features for data integration and network visualization. Bioinformatics 27:431–432. 10.1093/bioinformatics/btq675. 21149340. 3031041.
Tian L., Shi S., Nasir F., Chang C, Li W., Tran L.-S.P., Tian C.. 2018;Comparative analysis of the root transcriptomes of cultivated and wild rice varieties in response to Magnaporthe oryzae infection revealed both common and species-specific pathogen responses. Rice 11:26. 10.1186/s12284-018-0211-8. 29679239. 5910329.
Trapnell C., Roberts A., Goff L., Pertea G., Kim D., Kelley D.R., Pimentel H., Salzberg S.L., Rinn J.L., Pachter L.. 2012;Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nat Protoc 7:562–578. 10.1038/nprot.2012.016. 22383036. 3334321.
Tu J.C.. 1987;Integrated control of the pea root rot disease complex in Ontario. Plant Dis 71:9–13. 10.1094/PD-71-0009.
Turrà D., Segorbe D., Di Piertro A.. 2014;Protein kinases in plant-pathogenic fungi: conserved regulators of infection. Annu Rev Phytopathol 52:267–288. 10.1146/annurev-phyto-102313-050143. 25090477.
Weidmann S., Sanchez L., Descombin J., Chatagnier O., Gianinazzi S., Gianinazzi-Pearson V.. 2004;Fungal elicitation of signal transduction-related plant genes precedes mycorrhiza establishment and requires the dmi3 gene in Medicago truncatula. Mol Plant-Microbe Interact 17:1385–1393. 10.1094/MPMI.2004.17.12.1385. 15597744.
White T.J., Bruns T., Lee S., Taylor J.. 1990. Amplification and direct sequencing of fungal ribosomal RNA genes for phylogenetics. PCR protocols: a guide to methods and applications In : Innis M.A., Gelfand D.H., Sninsky J.J., eds. p. 315–322. Academic Press. New York, USA:
Zhang J., Wang J., Jiang W., Liu J., Yang S., Gai J., Li Y.. 2016;Identification and analysis of NaHCO3 stress responsive genes in wild soybean (Glycine soja) roots by RNA-seq. Front Plant Sci 71842;10.3389/fpls.2016.01842. 28018382. 5161042.
Zhang S., Klessig D.F.. 2001;MAPK cascades in plant defense signaling. Trends Plant Sci 6:520–527. 10.1016/S1360-1385(01)02103-3. 11701380.
Zhou X., Tian L., Zhang J., Ma L., Li X., Tian C.. 2017;Rhizospheric fungi and their link with the nitrogen-fixing Frankia harbored in host plant Hippophae rhamnoides L. J Basic Microbiol 57:1055–1064. 10.1002/jobm.201700312. 28902963.

Article information Continued

Fig. 1

Experimental design for examining the plant-microbe interactions of wild and cultivated soybean in response to Fusarium oxysporum. Soybean cultivars were selected from wild (01–289) and cultivated (Williams 82) soybean. MAPK, mitogen-activated protein kinase.

Fig. 2

The relative abundance of bacteria (A, B) and fungi (C, D) in wild and cultivated soybean in response to Fusarium oxysporum. Values (mean ± standard error) with different letters are significantly different at P < 0.05 (Tukey’s honestly significant difference test). The four treatments were non-inoculated wild soybean (W), non-inoculated cultivated soybean (C), wild soybean inoculated with F. oxysporum (W + F), and cultivated soybean inoculated with F. oxysporum (C + F). Panels B and D represent the different relative abundance levels of Proteobacteria and Glomeromycota, respectively.

Fig. 3

Differentially expressed genes related to the mitogenactivated protein kinase (MAPK) family in wild and cultivated soybean in response to Fusarium oxysporum infection. The four treatments were non-inoculated wild soybean (W), non-inoculated cultivated soybean (C), wild soybean inoculated with Fusarium oxysporum (W + F) and cultivated soybean inoculated with F. oxysporum (C + F). The red, blue and white colours represent the upregulated (fold changes ≥ 4 with a q-value < 0.05), downregulated (fold changes ≤ −4 with a q-value < 0.05) and unchanged genes, respectively. Additional information pertaining to MAPK signalling pathway-related genes is presented in Supplementary Table 2.

Fig. 4

Correlations between bacterial taxa and the mitogen-activated protein kinase (MAPK) family in wild and cultivated soybean in response to Fusarium oxysporum. (A) Wild soybean inoculated with F. oxysporum (W + F) vs. non-inoculated wild soybean (W). (B) Cultivated soybean inoculated with F. oxysporum (C + F) vs. non-inoculated cultivated soybean (C). (C) Wild soybean (W) vs. cultivated soybean (C). (D) Wild soybean inoculated with F. oxysporum (W + F) vs. cultivated soybean inoculated with F. oxysporum (C + F). The nodes represent the operational taxonomic units, and the colourful nodes indicate the microbial taxa at the phylum level. The connecting lines indicate a positive (red line) or a negative (blue line) correlation between the bacterial taxa and the MAPK genes. Additional information pertaining to the correlations is presented in Supplementary Table 4.

Fig. 5

Correlations between fungal taxa and the mitogen-activated protein kinase (MAPK) family in wild and cultivated soybean in response to Fusarium oxysporum. (A) Wild soybean inoculated with F. oxysporum (W + F) vs. non-inoculated wild soybean (W). (B) Cultivated soybean inoculated with F. oxysporum (C + F) vs. non-inoculated cultivated soybean (C). (C) Wild soybean (W) vs. cultivated soybean (C). (D) Wild soybean inoculated with F. oxysporum (W + F) vs. cultivated soybean inoculated with F. oxysporum (C + F). The nodes represent operational taxonomic units, and the colourful nodes indicate the microbial taxa at the phylum level. Connecting lines indicate a positive (red line) or a negative (blue line) correlation between the bacterial taxa and the MAPK genes. Additional information pertaining to the correlations is presented in Supplementary Table 5.