Four soils were collected in October 2012 (Supplementary Table S1): MS8 from Morens, county Fribourg, Switzerland (46° 52′ 04.77″ N and 6° 54′ 15.83″ E; Kyselková et al., 2009), Bmo1 from Béligneux, Ain, France (45° 52′ 22.28″ N and 5° 7′ 53.21″ E), Ysa5 and Ysa8 from Seyssel, Savoie region, France (respectively, 45° 57′ 2.42″ N and 5° 51′ 11.44″ E; Almario et al., 2013), and 45° 58′ 30.42″ N and 5° 51′ 2.43″E). MS8 and Bmo1 have a morainic origin whereas Ysa5 and Ysa8 developed on sandstone material. Bmo1 and Ysa5 were cultivated with maize for at least 3 years, MS8 was an artificial meadow for at least 2 years, and Ysa8 is a natural grassland. Soils were taken from 10 to 30 cm depth at three locations 5–10 m apart in each field, and were sieved at 0.5 mm.

Seeds of hybrid maize cultivar DK315 (Monsanto SAS/Dekalb, USA) and PR37Y15 (Pioneer Semences SAS, France) were surface-sterilized by soaking 1 h in sodium hypochlorite and one wash in 70% ethanol, and were rinsed three times with sterile distilled water. Three seeds of each maize variety were sown each in 2-dm³ pots (2 kg fresh soil/pot; 3 pots/cultivar/soil) and soil water content was maintained at 20% w/w. Unplanted pots (1 pot/soil) served as controls for sampling bulk soils. 21 days after sowing, rhizosphere soils and bulk soils were sampled (i.e., 3 root-adhering soils × 2 cultivars × 4 soils and 4 bulk soils). The root systems were shaken vigorously and rhizosphere extracts were prepared by putting each root system with adhering soil in 20 mL of 0.9% NaCl solution and shaking 1 h at 150 rpm. Four bulk soil extracts were obtained using 10 g soil (in 20 mL).

For isolation of fluorescent Pseudomonas, the 24 rhizosphere and four bulk soil extracts were serially diluted and 20 μL was mixed with 180 μL of King’s B⁺⁺⁺ [i.e., King’s B supplemented with ampicillin (40 mg/mL), chloramphenicol (13 mg/mL); Simon and Ridge, 1974; Landa et al., 2002; Latz et al., 2016] in 96-well microtiter plate following a most probable number (MPN) design with eight wells per dilution. Aliquots from each last positive well were plated on King’s B⁺⁺⁺ agar. At least 55 isolates were randomly selected for each of the 12 conditions [i.e., 4 soils × (2 cultivars and bulk soil)] and all colonies were purified three times successively, giving a total of 698 isolates. Genomic DNA was extracted for all isolates using NucleoSpin^® 96 Tissue kits (Ref - 740454.4; Macherey Nagel, Germany) and identification performed by sequencing the housekeeping gene rpoD (accession numbers LN885567 to LN886065, EMBL-EBI database) using primers rpoDf/rpoDr targeting the rpoD alleles of bacteria from the P. fluorescens group (Frapolli et al., 2007). When rpoD amplification failed, the 16SrRNA encoding rrs gene was amplified with pA/pH (Edwards et al., 1989) and sequenced (accession numbers: LN885368 to LN885566, EMBL-EBI database). rpoD sequences were aligned with MUSCLE (Edgar, 2004). Sequences were manually filtered to discard gaps and aligned regions of low quality. The phylogenetic trees were inferred with PHYML (Guindon et al., 2010) with the GTR model and 500 bootstraps. Isolate redundancy was estimated at 30%, according to both rpoD sequence similarity and functional profiles obtained.

A total of 18 plant-beneficial properties were targeted. Molecular screening for plant-beneficial properties was performed by PCR targeting genes involved in production of 2,4 diacetylphloroglucinol (phlD), pyrrolnitrin (prnD), pyoluteorin (pltC), phenazines (overlapping region of phzC/phzD) and NO (nirS in pseudomonads), as well as ACC deamination (acdS) and nitrogen fixation (nifH). All amplifications were performed with a thermocycler Mastercycler (Eppendorf, Germany). The reaction volumes contained 10x PCR buffer, 50 mM MgCl₂, 2 mM dNTP, 5% DMSO, 10 μM of each primer (Supplementary Table S2), 1 unit of Taq polymerase (Invitrogen, Cergy-Pontoise, France) and 50 ng of DNA. Several amplified fragments were sequenced and data blasted against the NCBI database in order to ascertain that isolates actually harbored the corresponding genes (accession numbers LT607759 to LT607801, EMBL-EBI database). No false-positive PCR results were found.

Screening for isolates with phosphate solubilizing activity was done by measuring the degradation halo on a National Botanical Research Institute’s Phosphate (NBRIP) agar after 6 days at 28°C, according to (Meyer et al., 2011). HCN production (indicated by color orange to red) was assessed after growth (3 days at 28°C) on King’s B⁺⁺⁺ agar, using a Whatman filter paper n°1 previously soaked in 2% sodium carbonate in 0.5% picric acid solution and placed in the lid of the Petri dish (subsequently sealed with parafilm). Production of extracellular protease and chitinase was assessed using, respectively, milk agar and minimum medium supplemented with colloidal chitin (Kim et al., 2003).

An Ultra High Performance Liquid Chromatography (UHPLC) method was developed in order to screen isolates for production of the two auxinic phytohormones indole-3-acetic acid (IAA) and indol-3-butyric acid and the five cytokinin phytohormones trans-zeatin, trans-zeatin riboside, kinetin, 6-benzylaminopurine and isopentenyl adenosine. Briefly, all isolates were grown in 2 mL of King’s B medium supplemented with 250 mM of auxin precursor tryptophan and 0.1 mM of cytokinin precursor adenine (2 days at 28°C, 300 rpm). The cultures were centrifuged at 4500 rpm during 8 min and filtered at 0.2 μm. Supernatants were subjected to UHPLC separation on an Agilent 1290 Series instrument using a 100 mm × 3 mm reverse phase column (Agilent Poroshell 120 EC-C18, 2.7 μm particle size) with a diode array detector. Samples (10 μL) were loaded onto the column equilibrated with water and acetonitrile (98:2). Compounds were eluted by a two-step gradient increasing the acetonitrile concentration to 40% over a 6 min period, then to 100% over 4 min, followed by an isocratic step of 2 min, at a flow rate of 0.5 mL/min. Hormones were detected with an Agilent 6530 Quadrupole Time-of-Flight (Q-TOF) mass spectrometer in positive electrospray ionization, based on comparison with commercial standards on both mass and UV (280 nm) chromatograms, along with accurate mass and UV spectra.

Heatmaps were analyzed using R “pheatmap” package (Kolde, 2012). Clustering analysis was performed using Euclidean distance method or Spearman correlation. For each condition, data were log-transformed for normal distribution and variance homogeneity, and a two-way ANOVA and Tukey’s HSD tests were performed to detect soil or variety impact on fluorescent Pseudomonas population size. Pseudomonas proportions were compared with a χ² test, or Fisher’s exact test when the expected values in any of the cells of a contingency table were below 5. Numbers of Pseudomonas isolates were compared by ANOVA and Fisher’s LSD tests. All analyses were performed at P < 0.05, using R software (Venables and Smith, 2011). Results in text and figures are presented as means ± standard error.

In all conditions, culturable fluorescent Pseudomonas populations were significantly more abundant in rhizosphere soil than in the corresponding bulk soil (Supplementary Figure S1). The mean number of culturable fluorescent Pseudomonas did not differ according to the maize cultivar (7.09 ± 0.26 log cells/g for PR37Y15 and 6.99 ± 0.22 log cells/g for DK315; P = 0.064), but the interaction between soils and cultivars was significant (P = 0.032). Indeed, in Ysa5 soil, pseudomonads were more abundant in PR37Y15 than in DK315 maize rhizosphere, whereas no significant differences were observed between cultivars in the three other soils. In addition, the number of culturable fluorescent pseudomonads in the rhizosphere differed according to the soil for cultivar PR37Y15 but not for DK315 (Supplementary Figure S1).

Among the 698 isolated Pseudomonas, 209 isolates were isolated from PR37Y15 maize rhizospheres, 255 from DK315 rhizospheres and 234 from bulk soils (Table 1). A total of 18 plant-beneficial properties were targeted among the 698 Pseudomonas isolates obtained. Whereas the nitrogen fixation gene nifH, as well as the ability to synthesize indole-butyric acid, trans-zeatin, trans-zeatin riboside, isopentenyl adenosine or 6-benzylaminopurine, were not found in any isolates from the library, the 12 other plant-beneficial properties studied were found in 2% (extracellular chitinase activity) to 96% (production of IAA) of the 698 pseudomonads. When considering together the 12 displayed plant-beneficial properties, a total of 175 different combinations (i.e., functional profiles) were found among the 698 isolates. When data from all soils were pooled together, no significant differences were observed when comparing the proportion of isolates displaying a given plant-beneficial property, in bulk soil and PR37Y15 and DK315 maize rhizospheres, regardless of the property tested (χ² tests, P > 0.05 for each plant-beneficial property; Figure 1). When soil treatments were analyzed separately (Figure 2), the distribution of these properties in pseudomonads differed according to soil type and maize cultivar, especially for properties involved in the biosynthesis of antimicrobial compounds (i.e., DAPG, pyrrolnitrin, pyoluteorin) found in lower proportion in PR37Y15 and DK315 maize rhizospheres than bulk soil for Bmo1 (P-value of χ² tests from 10^-10 to 0.04). For the same soil, however, the proportion of Pseudomonas harboring acdS gene was higher in both maize rhizospheres than bulk soil (P < 0.01; Figure 2). Conversely, in MS8 soil, proportions of Pseudomonas harboring genes involved in the biosynthesis of antimicrobial compounds (i.e., DAPG, pyrrolnitrin, pyoluteorin and phenazine) were higher in both maize rhizospheres than bulk soil (P-value of χ² tests from 10^-12 to 0.05) whereas acdS was higher in bulk soil than both maize rhizospheres (P < 0.05). In the other two soils, no significant differences were detected for the distribution of each property.

When considering the numbers (rather than relative proportions) of fluorescent Pseudomonas – estimated by multiplying the Pseudomonas proportions with the raw population levels from each bulk and rhizosphere soil, and thereafter pooling together data from all soils – a higher number of pseudomonads harboring each property studied was found in the rhizosphere of one (i.e., PR37Y15 for phlD, prnD, pltC) or both maize cultivars (i.e., phzCD, acdS, nirS, extracellular protease activity, HCN and IAA production, phosphate solubilization) than in the bulk soil, by approximately two orders of magnitude (Supplementary Figure S2). Numbers of kinetin producers and of pseudomonads with chitinase activity were below 2 log cells/g in all three treatments.

Pseudomonas isolates possessed from 0 (11 isolates) to a maximum of 9 (1 isolate) co-occurring plant-beneficial properties (i.e., functions harbored by a same isolate) among the 18 plant-beneficial properties targeted, with a mean of 3.6 properties per isolate. To assess whether maize preferentially selects fluorescent Pseudomonas harboring a high number of co-occurring plant-beneficial properties, the number of co-occurring plant-beneficial properties per isolate was analyzed according to their origin of isolation (rhizospheres or bulk soils). When data from all soils were pooled together, most rhizosphere isolates (433 of 464 rhizosphere isolates, i.e., 93%) displayed one to five co-occurring plant-beneficial properties per isolate, whereas 81 of 234 bulk soil isolates (i.e., 35%) displayed six to nine of them (Figure 3A). When considering each soil separately, the same proportion distribution of Pseudomonas sharing 0 to 9 co-occurring plant-beneficial properties per isolate was found in the maize PR37Y15 and DK315 rhizospheres (from P = 0.32 for Ysa8 soil to P = 0.56 for Bmo1 soil; Supplementary Figure S3). In addition, the distribution was similar for bulk soil and rhizosphere in MS8 and Ysa5 soils (Supplementary Figures S3A,B) whereas it differed in soils Ysa8 and Bmo1 (Supplementary Figures S3C,D). When population levels were estimated for fluorescent Pseudomonas harboring various numbers of co-occurring plant-beneficial properties and data from all soils were pooled together, it appeared that the numbers of culturable fluorescent Pseudomonas harboring up to five plant-beneficial properties reached about 6 log cells/g in the maize rhizosphere of both PR37Y15 and DK315, but only 4 log cells/g in bulk soil (P < 0.001; Figure 3B). In contrast, pseudomonads with six to nine plant-beneficial properties were found in the maize rhizospheres in similar substantial numbers, with 2.25 ± 0.47 log cells/g for PR37Y15 and 2.77 ± 1.42 log cells/g for DK315 compared with bulk soil (2.70 ± 0.88 log cells/g; Figure 3).

Among the 698 isolates, the ratio of the number of functional profiles (i.e., different associations of plant-beneficial properties in isolates) to the number of isolates was significantly higher (P = 0.030) in PR37Y15 (0.56 ± 0.02) compared to DK315 (0.45 ± 0.03) maize rhizosphere and bulk soil (0.45 ± 0.05). However, correspondence analysis indicated that bulk soil and rhizospheres overlapped completely when assessing Pseudomonas isolate profiles (Supplementary Figure S4). Results thus point to a lack of selection of pseudomonads with specific assortments of plant-beneficial properties in the rhizosphere.

Some plant-beneficial properties preferentially co-occurred in Pseudomonas isolates. Genes for synthesis of organic antimicrobial compounds (i.e., DAPG, pyrrolnitrin, pyoluteorin and phenazines) grouped together following hierarchical clustering analysis (Figure 2) and significantly correlated with one another (Supplementary Figure S5). Likewise, properties involved in modulation of plant hormonal balance (i.e., ACC deaminase activity gene acdS and NO production gene nirS, which correlated positively; r = 0.70, P < 0.011, Supplementary Figure S5), phosphate solubilization, extracellular protease activity and HCN production grouped together (Figure 2).

In order to characterize taxonomically the Pseudomonas isolates, rpoD amplification and sequencing were performed using specific primers for the P. fluorescens group, which was successful for 498 isolates (60% of PR37Y15 isolates, 82% for DK315 isolates and 81% of bulk soil isolates; i.e., 70% overall Pseudomonas isolates). The 16S rRNA rrs gene sequencing on the 200 remaining isolates affiliated all of them to the P. putida group. rpoD phylogeny (Supplementary Figure S6) enabled affiliation of the rpoD-sequenced isolates to 7 of the 9 subgroups of the P. fluorescens group proposed by Garrido-Sanz et al. (2016), i.e., the subgroups P. corrugata (3.6% of isolates), P. chlororaphis (5.9%), P. protegens (9.3%), P. fluorescens (20.2%), P. mandelii (6.9%), P. jessenii (11.5%), or P. koreensis (14.0%; Figure 4A; Supplementary Figure S6; Table 1). According to the phylogenies of Gomila et al. (2015) and Garrido-Sanz et al. (2016), two clusters, hereafter termed clusters FMJK (for the P. fluorescens, P. mandelii, P. jessenii, P. koreensis subgroups) and CPC (for the P. chlororaphis, P. protegens, P. corrugata subgroups), were defined in this study.

The distribution of plant-beneficial properties was assessed according to their affiliation to the 7 P. fluorescens subgroups recovered in this study (Figure 4; Table 2). The P. chlororaphis subgroup shared the highest functional diversity with 12 distinct plant-beneficial properties recovered among isolates from this subgroup (regardless of their co-occurrence) whereas strains affiliated to P. koreensis subgroup shared the lowest functional diversity (i.e., 8 distinct plant-beneficial properties; Table 2). The number of plant-beneficial properties found per isolate (i.e., co-occurring properties) was significantly lower (P < 0.001) for the FMJK cluster (3.36 ± 0.07 properties per isolate) than the CPC cluster (5.76 ± 0.13 properties; Figure 4B).

The ecological functions conferred by plant-beneficial properties harbored by isolates also differed according to the taxonomic affiliation of Pseudomonas isolates. When regrouping plant-beneficial properties according to their type of action, differences between Pseudomonas subgroups were observed regarding the proportion of (i) biocontrol properties (i.e., genes for antimicrobial compounds DAPG, pyrrolnitrin, pyoluteorin and phenazine, and production of extracellular protease, chitinase, and HCN) and (ii) properties related to plant hormonal balance modulation and plant nutrition (i.e., ACC deaminase and NO genes, production of kinetin and IAA, and phosphate solubilization; Figure 4A). Overall, higher proportions of biocontrol properties were found in isolates belonging to the CPC cluster (68% for P. chlororaphis, 67% for P. protegens and 50% for P. corrugata subgroups) than among most members of the FMJK cluster (29% for P. fluorescens, 20% for P. mandelii, 30% for P. jessenii and 56% for P. koreensis subgroups; Figure 4A). Conversely, phytostimulatory properties (i.e., plant hormone modulation and nutrition properties) were found, for the most part, in higher proportion in FMJK isolates (71% for P. fluorescens, 80% for P. mandelii, 70% for P. jessenii and 44% for P. koreensis subgroups) than CPC isolates (32% for P. chlororaphis, 33% for P. protegens and 50% for P. corrugata subgroups; Figure 4A).

The proportions of all individual biocontrol properties in the pseudomonads were also higher (P < 0.001 each) in cluster CPC than FMJK (Table 2). Contrariwise, all properties related to plant hormonal balance modulation and plant nutrition were found in similar proportions in CPC and FMJK clusters.

Isolates belonging to the P. putida group displayed more plant-beneficial properties related to plant hormonal balance modulation and plant nutrition than to biocontrol (Figure 4A). Moreover, these isolates harbored a low number of co-occurring plant-beneficial properties, i.e., 2.5 ± 0.1 (Figure 4B).

The FMJK and CPC clusters were present in the same relative proportion (40%) in bulk soil, but in the rhizosphere, FMJK isolates were retrieved in a significantly higher proportion than CPC isolates (up to 60% for FMJK versus below 10% for CPC; Figure 5). The P. putida group was found in the same proportion in bulk soil and in the rhizosphere. The proportion of this group was significantly more important than that of the CPC cluster in PR37Y15 maize rhizosphere (Figure 5). The same trends were observed when considering the estimated number of Pseudomonas from FMJK and CPC clusters, in bulk soils or rhizospheres (Supplementary Figure S7).

When considering the amounts of co-occurring plant-beneficial properties per isolate, FMJK isolates with one to five plant-beneficial properties were more abundant than FMJK isolates with six to nine plant-beneficial properties in all soils (P < 0.001; Figure 6A). This was also the case when considering estimated Pseudomonas population levels (Supplementary Figure S7B).

In parallel, the proportions of isolates belonging to CPC cluster with six to nine plant-beneficial properties were smaller in the rhizospheres of both cultivars (P = 0.002) than in bulk soil, whereas those with one to five plant-beneficial properties were found in similar proportions in bulk soils and maize rhizospheres (P = 0.28; Figure 6B). When considering estimated Pseudomonas population levels, the numbers of CPC isolates with one to five or six to nine plant-beneficial properties were similar between bulk soil and maize rhizospheres, except in the PR37Y15 maize rhizosphere where a higher number of CPC isolates with one to five plant-beneficial properties was found than in the other conditions (Supplementary Figure S7C).

In conclusion, the selection of Pseudomonas strains in the maize rhizosphere depends more on the number of plant-beneficial properties they harbor than the Pseudomonas cluster they belong to.