The Fusarium-wilt-suppressive soil was harvested from a field that had remained a fallow for several years in the Châteaurenard region, France (43°53′15″N, 4°51′36″E). The conducive soil was sampled nearby from a temporarily uncultivated field following muskmelon cropping and before lettuce cropping (43°51′53″N, 4°50′04″E). The two soils were geographically very close (3 km apart). Both were collected in October 2011, air-dried and sieved to 4 mm. Both were loamy soils. Physicochemical profiles were respectively 17.1% clay, 47.4% silt, 35.5% sand, 1.98% organic matter (OM), 0.9 N (g/kg), and pH 8.5 for the suppressive soil, and 20.6% clay, 43.2% silt, 36.2% sand, 4.6% OM, 1.5 N (g/kg), pH 8.06 for the conducive soil. Differences in organic matter content, nitrogen content and pH were related to the current uses of the soils. To confirm the biological nature of the soil suppressiveness, the suppressive soil was autoclaved three times on three consecutive days for 20 min at 120°C and stored at room temperature for 1 week before use.

Soil disease suppressiveness was tested under greenhouse conditions with flax as the host plant and the flax pathogen F. oxysporum f. sp. lini MIAE00347 (Collection of Microorganisms of Interest for Agriculture and Environment, INRA Dijon, France). The plants were grown in different types of soil: suppressive soil (S), conducive soil (C), conducive soil amended with 10% (w/w) of suppressive soil (referred to as mixed soil or M), and heat-treated suppressive soil. Soils were inoculated with 10³ conidia of the pathogen per mL of soil, and flax (Linum usitatissimum L. variety Opaline) seeds were sown on the same day. For each modality, three replicates were performed with 20 plants/replicate. Non-inoculated soils were used as controls. Plants were grown in a growth chamber with 70% relative humidity, a 16/8 h daylight/dark photoperiod at 17°C (day) and 15°C (night) the first week, and 25°C (day) and 22°C (night) afterwards. Diseased plants were recorded 25 days after inoculation. Data were subjected to one-way analysis of variance (ANOVA). Pairwise comparisons were performed using Fisher's test (P < 0.05).

Flax seeds were sown in 9/9/9.5 cm pots containing 300 g of soil. Twenty plants per pot were grown under the same conditions of light and humidity as above. Four modalities were prepared: soil S, soil C, soil M, and suppressive soil inoculated with 10³ conidia of F. oxysporum f. sp. lini MIAE00347 per mL of soil (referred to as pathogen-inoculated suppressive soil or IS soil). The inoculum dose used was relatively low (10³ conidia per mL) to ensure that plants do not die too quickly in the conducive soil and prevent to collect the active rhizospheric soil needed to conduct microbial communities comparative analyses. Three replicates per modality were performed. After 25 days, the rhizospheric soil was collected from each replicate. The root systems were isolated and shaken to remove free soil particles. The soil around the roots influenced by root development and plant exudates was considered as rhizospheric soil. It was sieved to 2 mm to remove fine roots and organic debris. Large pieces of roots were removed manually with tweezers. Two grams of soil samples from each modality were placed in 2-ml cryotubes and frozen at −20°C for further DNA extractions.

The protocol used for DNA extractions was described by Plassart et al. (2012). Extraction was performed by mechanical lysis using FastPrep®-24 (MP-Biomedicals, NY, USA) with a lysis buffer and two purification steps with PVPP (polyvinyl polypyrrolydone) minicolumns (BIORAD, Marne-la-Coquette, France) and Geneclean Turbo kit (Q-Biogene, Illkirch, France). Total DNA was extracted from each replicate of the four modalities. Each purified DNA sample was quantified by a fluorometric assay using the PicoGreen dsDNA Assay Kit (Invitrogen Life Technologies) and the StepOnePlus™System (Applied Biosystems®) according to the manufacturer's instructions. DNA extracted from each soil sample served as a template in PCR reactions to amplify a fungal barcode and a bacterial barcode. Three replicates were used for each of the four soil modalities. PCR reactions were performed in quadruplicate for each replicate, with 20 ng of soil DNA.

For fungal identification, the variable internal transcribed spacer 1 (ITS1) region was amplified using specific fungal primers ITS1F (5′-CTTGGTCATTTAGAGGAAGTAA-3′) and ITS2 (5′-GCTGCGTTCTTCATCGATGC-3′; White et al., 1990; Gardes and Bruns, 1993; Buée et al., 2009). The primers were tagged with four-base-pair multiplex identifiers (MIDs) at the 5′ and 3′ positions to specifically identify each sample, as recommended by the manufacturer (Beckman Coulter Genomics). PCR conditions were as follows: 95°C for 5 min, 35 cycles of 30 s at 95°C, 40 s at 53°C, and 45 s at 72°C, followed by 7 min at 72°C.

For bacterial identification, a 16S rRNA gene fragment (partial V3, V4 and partial V5) was amplified using primers 530F (5′-ACTCCTACGGGAGGCAGCAG-3′; Acosta-Martínez et al., 2008) and 803R (5′-CTACCNGGGTATCTAAT-3′; Zancarini et al., 2013). Ten-base-pair MIDs at the 5′ and 3′ positions were added to the primers to specifically identify each sample, as recommended by the manufacturer (Beckman Coulter Genomics). PCR conditions were as follows: 95°C for 10 min, 30 cycles of 30 s at 95°C, 30 s at 62°C, and 60 s at 72°C, followed by 10 min at 72°C. All amplifications were performed in a Mastercycler (Eppendorf, Hamburg, Germany). For each sample, amplicons of the four replicated PCRs were pooled and purified using a MinElute PCR Purification Kit (Qiagen, Courtaboeuf, France) following the manufacturer's protocol. Amplicon concentrations were then estimated by fluorometric assay (PicoGreen dsDNA Assay Kit). For each of the two amplified barcodes, an equimolar pooling of all samples was prepared (total DNA amount: ~3.2 μg per library). Adapter sequences were added by ligation as recommended by the manufacturer, and 454-pyrosequencing was carried out by Beckman Coulter Genomics (Danvers, USA) on a Genome Sequencer FLX 454 (Life Sciences/Roche Applied Biosystems).

Concerning the sequenced ITS1 region, sequences were sorted into different files according to their MID using the sfffile program of Roche 454 main software with default parameters. Mismatches were not allowed for MIDs and were removed. Sequences were then cleaned by trimming and denoising using trim.flows in Mothur version 1.20.1 with the default parameters (Schloss et al., 2009), and reverse-complemented if needed. ITS1 was extracted using Fungal ITS Extractor, version 2 (Nilsson et al., 2010) and sequences were filtered based on a minimal length of 100 bp. Operational taxonomic units (OTUs) were generated after two successive clustering steps using Uclust version 3.0 [Usearch version 5.2.32, (Edgar, 2010)] at 97% similarity. The first clustering included all sequences, and the second clustering was conducted with the batch of consensus sequences from the previous clustering.

Concerning the sequenced 16S region, sequences were analyzed using the GnS-PIPE pipeline initially developed by the GenoSol platform (INRA, Dijon, France) and recently optimized (Terrat et al., 2012, 2015). First, all the raw 16S reads were sorted according to their multiplex identifier sequences. Then, a preprocessing step was realized to filter and delete low-quality reads based on (i) their length (<350 bases), (ii) their number of ambiguities (deletion of reads with one N or more), and (iii) their primer sequence(s) (the distal and proximal primer sequences must be complete and without errors). A PERL program was then applied for rigorous dereplication (i.e., clustering of strictly identical sequences with same length). The dereplicated reads were aligned using INFERNAL alignment (Cole et al., 2009; v1.0.2 with chosen parameters: –hbanded, –sub, –dna) to obtain a global alignment. Then, aligned sequences were clustered into OTUs at 93.8% of similarity using the CrunchClust V43 program (Hartmann et al., 2012) that groups rare reads with abundant ones and does not count differences in homopolymer lengths (default parameters were selected). Here, an OTU is defined by the most abundant read, known as the centroid, and every read in the OTU must have similarity above the given similarity threshold with the centroid). The chosen level of similarity was defined by an in silico approach (data not shown). We first extracted all known and reliable microbial sequences (bacteria and archaea) from the SILVA database (release v111) (Quast et al., 2013). We kept from these sequences only the amplified regions of sequences using our primer set, and deleted those with too many mismatches with our primer set (more than 3 mismatches with one primer). Then, all these “artificial reads” with a reliable taxonomy were clustered using our specific program of clustering at various threshold levels (100% to 90%, with a 0.1% step). Regarding the results and the amplified regions, the 93.8% threshold was the best suited to efficiently define the genus level, as it was the closest one to the genus level defined by the SILVA taxonomy of analyzed sequences.

For both the sequenced ITS and 16S regions, a filtering step was carried out to check all singletons (OTUs supported by only one sequence likely to be artifacts, such as PCR chimeras) and potentially to delete them from downstream analysis based on the quality of their taxonomic assignments. Then, to efficiently compare datasets and avoid biased community comparisons, high-quality reads were normalized by random selection: each dataset with a large number of reads was randomly cut down to the same number of reads as the dataset with the lowest number of reads (3,451 and 2,505 reads for fungal and bacterial communities, respectively) using a homemade Perl program (similar to the single_rarefaction.py with default parameters from QIIME). The retained high-quality reads were then used for taxonomy-independent analyses, and several diversity and richness indices were determined using the defined OTU composition. More precisely, we used the number of bacterial and fungal OTUs from each sample to determine Shannon and Evenness indices as indicators of soil microbial diversity and structure, respectively. Taxonomy-based analysis was also performed using similarity approaches against: (i) the UNITE database (fungal rDNA ITS sequence database; Kõljalg et al., 2013) using the Basic Local Alignment Search Tool algorithm BLASTn version 2.2.23 (Altschul et al., 1997) with the filter turned off for fungal sequences; (ii) a dedicated reference database from SILVA (Quast et al., 2013), using the USEARCH tool (Edgar, 2010) for bacterial sequences. Finally, global analysis of the soil samples was computed by merging all homogenized high-quality reads from each sample into one global file and defining OTUs as previously described before subsequent analyses.

All raw sequences collected in this study have been deposited in the European Bioinformatics Institute nucleotide sequence database system under the accession number PRJEB24081.

The OTUs defined for fungal and bacterial datasets were used to perform rarefaction analysis using Analytic Rarefaction version 1.3 (Hunt Mountain Software, Department of Geology, University of Georgia, Athens, GA, USA). The differences in fungal and bacterial community compositions (number of genera, OTUs, Shannon and Evenness indices) were assessed by Kruskal–Wallis tests under XLSTAT software (Addinsoft®). The OTU tables obtained from ITS1 and 16S datasets were used to perform an ordination analysis by means of Principal Component Analysis (PCA), using ADE-4 version 2001 (Biometry and Evolutionary Biology Laboratory, University Claude Bernard Lyon 1). For each microbial community, the numbers of unique OTUs and shared OTUs among the different soils were compared in Venn diagrams using VENNY (Oliveros, 2007).

Heat maps were built by applying the heatmap.2 function implemented in “gplots” R package based on the relative abundance values of the dominant fungal OTUs (for which the relative abundance level in one soil at least was higher than 1%) and dominant bacterial OTUs (for which the relative abundance level in one soil at least was higher than 1.5%) across all samples using the global matrices. Datasets were analyzed together, so the mean values of Z-scores were from the total dataset. The relative abundance levels of each row (corresponding to each OTU) were expressed as median-centered Z-scores between all samples, and the colors scaled to standard deviations of the corresponding Z-scores. All these analyses were performed with R free software (version 3.2.3). Moreover, to minimize the probability of incorrect assignments of the dominant OTUs, taxonomic assignments of representative sequences were refined by checking taxonomic assignments of dominant fungal OTUs against the GenBank database. The default algorithm parameters were used, and environmental sequences were excluded. Only sequences associated with a publication and a description of the corresponding strain were considered. Assignments were taken into account at the genus level. Some synonyms were corrected using the online Species Fungorum database (Kirk, 2017). For the few OTUs that had the same percentage of similarity with several fungal genera or families, assignment was made at the family or order level, respectively. Concerning the dominant bacterial OTUs, their taxonomic assignments were checked using SINA online (Pruesse et al., 2012), the SILVA database (version 123), and default parameters. SINA uses the search result to derive a classification with the LCA (lowest common ancestor) method. Each query sequence is assigned the shared part of the classifications of the search results. OTU abundance levels in the different soils were compared based on Kruskal-Wallis tests using XLSTAT software (Addinsoft®).

The plants grown in the non-inoculated soils showed no wilt disease symptoms. Typical symptoms of Fusarium wilt were observed in pathogen-inoculated soils. The percentage of healthy plants was significantly higher (P < 0.05) in soil S (46.6 ± 10.6) and soil M (43.5 ± 19.4) than in soil C (15.3 ± 5.7). The percentage of healthy plants was significantly lower (P < 0.05) in heat-treated suppressive soil (9.2 ± 0.3) than in soils S, M and C, further confirming the biological nature of the soil suppressiveness.

Fungal and bacterial targeted regions were successfully amplified by PCR and sequenced for all soils. The raw sequence libraries were filtered to remove reads originating from sequencing errors or putative chimeric sequences. For the fungal and bacterial datasets, 81,902 reads (88% of raw sequences), and 67,675 reads (89% of raw sequences) passed all quality controls, respectively. The number of high-quality reads ranged from 3,451 to 10,821 for fungal datasets, and from 2,505 to 8,842 for bacterial datasets. After homogenization, a total of 41,412 high-quality fungal ITS sequences were clustered into 1,798 OTUs for all soil samples. Among them, 570 OTUs were considered as singletons, and represented 31.7% of the detected richness. Regarding bacterial communities, 2,280 distinct OTUs, including 1,162 singletons (50.9% of the richness), were observed for 30,060 high-quality bacterial 16S rDNA sequences. Rarefaction curves were drawn for fungal and bacterial datasets based on OTUs (Supplementary Figure 1). Fungal curves showed that: (i) the number of OTUs at 97% similarity increased with the number of sequences, and saturation was not reached for all soils; (ii) based on soil replicates, however, the number of sequences was sufficient to obtain a representative coverage of the major fungal groups; and (iii) differences among soils were recorded in the slope and level of the curves (Supplementary Figure 1A). As regards bacterial richness, rarefaction curves also revealed that the number of OTUs increased with the number of reads, without reaching a plateau, and reads were in sufficient numbers to allow for an accurate description of the major bacterial genera in each of the soil samples (Supplementary Figure 1B). Interestingly, no significant differences were detected in the slope and level of bacterial curves in any of the soils, showing less heterogeneity of bacterial richness among the four modalities.

Fungal and bacterial communities were also evaluated using richness and diversity indices (Table 1). Similar numbers of genera were detected in the different soils. However, significant differences among soils were recorded (P < 0.05) for the number of fungal OTUs, and for Shannon and Evenness indices. Thus, soil M harbored the highest fungal richness (462.0 ± 26 OTUs on average), and soils C and IS the lowest. However, soils C and M had the highest Shannon indices, and soil IS the lowest. Finally, evenness was significantly higher in soil C (P < 0.05) than in soils S and IS. Contrary to fungal communities, bacterial communities did not significantly differ in the number of genera, the number of OTUs, and the Shannon and Evenness indices between the different soils (Table 1).

Fungal and bacterial sequences were assigned at different taxonomic levels (phylum, class, order, family, and genus; Supplementary Table 1). Concerning fungal communities, most assigned sequences belonged to the Ascomycota phylum (90.9%). Basidiomycota and Zygomycota accounted for 2.6 and 1.5% of the sequences, respectively. The proportion of Chytridiomycota and Glomeromycota was <1%. Finally, 4.4% of the sequences were not assigned and were considered as undefined. At the class level, Sordariomycetes were dominant within Ascomycota, followed by Eurotiomycetes, Dothideomycetes, and Pezizomycetes (Figure 1A).

Concerning bacterial communities, the soils had the same overall bacterial composition, with Proteobacteria, Bacteroidetes, Actinobacteria, Firmicutes, and Acidobacteria as major phyla (Figure 1B). Proteobacteria was the dominant bacterial phylum and represented 46–55% of all bacterial DNA sequences. Several minor phyla (Chloroflexi, Gemmatimonadetes, Verrucomicrobia, Planctomycetes, and Nitrospirae) were also identified; they accounted for ca. 7% of bacterial communities. However, a few differences were found among soils, with more Gamma-Proteobacteria in soil M than in soils S and IS (P < 0.05). In contrast, Actinobacteria were more represented in soils S and IS than in soil C (P < 0.05).

The structure of the fungal and bacterial communities was analyzed by PCA to obtain an ordination on a factorial map and to compare the soils based on global OTU composition (Figure 2). A clear discrimination of fungal communities based on soil types was observed (Figure 2A). In particular, two groups of soils were separated along the first axis: soils C and M on the one hand, and soils S and IS on the other hand. Moreover, soil C grouped apart from soil M, and soil S apart from soil IS along the second axis. Concerning bacterial communities, sample discrimination was not so obvious (Figure 2B): soils C and M were separated from soils S and IS along the first axis, but no distinction of these two groups was visible along the second axis.

Regarding OTU composition, 189 fungal OTUs and 414 bacterial OTUs were shared between soils C and S (Figure 3). In the fungal community, 205 OTUs were unique in soil C, and 271 in soil S. In the bacterial community, 355 OTUs were unique in soil C, and 293 in soil S. To further focus on the microbial genera associated with disease suppressiveness, we compared the fungal and bacterial community compositions between soils S and M on the one hand, and between soils S and IS on the other hand. 122 fungal OTUs and 100 bacterial OTUs were shared between soils S and M but absent from soil C. In addition, 214 fungal OTUs and 145 bacterial OTUs shared between soils S and IS were potentially promoted by the presence of the pathogen.

To combine taxonomic assignment and soil global OTU composition, we compared microbial community compositions based on the relative abundance of the detected microbial genera using a heat map. Due to the large size of the dataset, only the most dominant fungal OTUs (relative abundance level > 1%) and bacterial OTUs (relative abundance level > 1.5%) were mapped (Figures 4, 5). The 143 dominant fungal OTUs represented in the heat map accounted for 83% of total reads on average, and the 122 dominant bacterial OTUs kept in the heat map accounted for 72% of total reads on average (Figures 4, 5). Firstly, these heat maps showed that biological replicates from a same soil were more similar to each other than replicates from other soils, whatever the soil. Moreover, evidence of clear differences between soils C and M on the one hand, and soils S and IS on the other hand was illustrated by both bacterial and fungal community compositions. Concerning fungal communities, 24 OTUs affiliated to 15 known genera were relatively more abundant in soils S and IS than in soil C (P < 0.05; Figure 4, underlined OTUs). Among them, 17 OTUs were not detected at all in soil C. They were assigned to the Acremonium, Chaetomium, Cladosporium, Clonostachys, Fusarium, Ceratobasidium, Mortierela, Penicillium, Scytalidium, and Verticillium genera. Interestingly, these OTUs were also detected in soil M. All 10 genera mentioned above, exclusively detected in soils S, IS and M, have been described for their antagonistic activity against diverse plant pathogens (Table 2). Moreover, the Acremonium, Chaetomium, Clonostachys, Fusarium, and Penicillium genera include strains used in the biological control of pathogenic F. oxysporum (Table 2).

Concerning bacterial communities, the OTUs detected in soils C and M also differed from those of soils S and IS, with dominant Alpha-Proteobacteria and Gamma-Proteobacteria in soils C and M, and dominant Actinobacteria in soils S and IS (Figure 5). More precisely, several OTUs affiliated to known bacterial genera were more abundant in soils S than in soil C (P < 0.05), e.g., Adhaeribacter, Amycolatopsis, Arthrobacter, Geobacter, Massilia, Microvirga, Paenibacillus, Rhizobacter, Rhizobium, Rubrobacter, and Stenotrophomonas (Figure 5, underlined OTUs). Interestingly, some of these OTUs (i.e., Adhaeribacter, Amycolatopsis, Arthrobacter, Paenibacillus, Rhizobacter, Rubrobacter, and Stenotrophomonas) were also more abundant in soil IS than in soil C (P < 0.05). Finally, some of these OTUs (i.e., Adhaeribacter, Amycolatopsis, Paenibacillus, Rhizobacter, and Rubrobacter) were also more abundant in soil M than in soil C (P < 0.05). Interestingly, four genera mentioned above (Arthrobacter, Paenibacillus, Rhizobium, and Stenotrophomonas) are known to have pathogen suppression potential (Table 2), and two of them (Paenibacillus and Rhizobium) include strains used in the biological control of pathogenic Fusarium (Table 2).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.