In this study, we used the OTU6_Psi_1 strain of P. siliginis that has been isolated from the rosette of one individual of A. thaliana collected in spring 2015 in the natural population ESPE-B located in the southwest of France (Bartoli et al., 2018; Ramírez-Sánchez et al., 2022). This strain was demonstrated to have a PGPB effect on A. thaliana under in vitro conditions when inoculated both at seed or seedling stage (Ramírez-Sánchez et al., 2022). In addition, we revealed a high genetic variation in response to this strain among seven A. thaliana accessions located in the southwest of France to be suppressed (Ramírez-Sánchez et al., 2022). Using single-molecule real-time long reads with a PacBio Sequel II system, a de novo genome sequence was obtained for the OTU6_Psi_1 strain, showing a single chromosome containing 5,458 genes (Ramírez-Sánchez et al., 2022).

Fifty-four populations (each represented by three accessions) were selected to represent the ecological and genetic diversity observed among a set of 168 natural populations of A. thaliana from southwest of France (Frachon et al., 2018; Frachon et al., 2019) ( Supplementary Table S1 ). The seeds coming from the maternal plants were harvested in June 2015. To reduce differences in the maternal effects of the 162 seed lots (i.e. 54 populations × three accessions), one plant per accession per generation was grown as followed: (i) several seeds of each accession were sown on October 1^st 2016 in 7 x 7 x 6 cm plastic pots (Soparco^®) filled with damp standard culture soil (PROVEEN MOTTE 20, Soprimex^®); (ii) seeds were stratified at 4°C for four days; (iii) pots were put on November 4^th 2016 to a greenhouse at 22°C with a 16 hours photoperiod; (iv) seedlings were thinned to one on November 25^th 2016; (v) seedlings were transferred to the INRAE campus of Auzeville field station (France) on December 5^th 2016; (vi) when plants started to flower, they were moved to a greenhouse that reproduces outdoor conditions (no extra light or heating) but protects the plants from rain; (vii) aratubes (Arasystem^®) were put on each plant to prevent cross-pollination between accessions; (viii) seeds were collected from late April to early May 2017 and conserved at 4°C.

The sterilization of the seeds was performed with chlorine gas as previously described (Ramírez-Sánchez et al., 2022) and the seeds were then kept at 4°C.

An experiment of 3,888 plants was set up using a split-plot design with two treatments (i.e. inoculation with OTU6_Psi_1 and mock treatment) nested within two blocks. Each ‘block × treatment’ interaction was represented by 24 48-well plates with each well filled with 700 µL of 0.5x MS medium (Murashige and Skoog medium), which contains 2.2 g of MS medium, 0.5 g of 2-(N-Morpholino)-ethanesulfonic acid, 6.0 g of plant tissue culture agar, 1 L of deionized water and a pH adjusted to the range of 5.7-5.8. The six wells of the last column of each plate were sown with seeds from the Col-0 reference accession to control for micro-environmental variation among plates. The 162 accessions were randomly assigned to the remaining columns of the 24 plates, resulting for each accession in a total of 12 replicates (two columns of six wells) in each treatment. The same randomization was done among treatments within a block but was modified between the two blocks.

For the 162 natural accessions and Col-0, one seed was sown in each well. After the 7-day cold treatment, seeds were inoculated and placed in a phytotron (10 hours photoperiod, light intensity ~ 80 µmol m^-2 s^-1, 21°C, 50% hygrometry) with a daily plate randomization.

The OTU6_Psi_1 strain was grown from glycerol stock on solid medium in a 9 cm x 1.5 cm circle polystyrene Petri dish filled with TSA medium for one day. Colonies were diluted in 500 µl of sterile deionized water. 250 µl were then deposited in two new 9 cm x 1.5 cm circle polystyrene Petri dishes filled with TSA and spread with sterile beads. The two plates were incubated at 28°C overnight. Bacterial colonies were then resuspended and diluted in sterile deionized water to an OD_600nm of 0.5 which corresponds to 3.90 x10⁸ CFU/mL. Each seed was inoculated either with 5 µl of OTU6_Psi_1 (P. siliginis treatment) or 5 µl of sterile deionized water (mock treatment). Col-0 seeds that were sown in the last column of each 48-well plate were not inoculated. Plates were sealed with a micropore tape (3M Micropore Surgical x 9.14 m).

The germination date was recorded each day, between three and seven days after sowing. A picture of each 48-well plate was taken on 14, 21, and 28 days after inoculation (dai) using a photo-box designed in the lab and with a mobile camera (Samsung S6 16 Mpx). The vegetative growth was then visually scored for each plant using a scale established in Ramírez-Sánchez et al. (2022) and ranging from one (very small plant) to seven (well-grown plant). A total of 4,460 plants (3,888 plants of the 162 natural accessions and 572 Col-0 plants) were therefore phenotyped after inoculation at 14 dai, 21 dai and 28 dai.

Based on within-population genetic variation previously available for 168 natural populations of A. thaliana (Frachon et al., 2018), a Bayesian hierarchical model (Gautier, 2015) was applied to estimate the standardized allele frequencies corrected for the effect of population structure within each population for 1,638,649 SNPs across the genome (Frachon et al., 2018; Frachon et al., 2019). Standardized population allele frequencies were then retrieved for the 54 populations used in this work. Then, for each of the three PGR traits (i.e. PGR estimated at 14 dai, 21 dai and 28 dai), a genome scan was launched by estimating for each SNP the Spearman’s rho value and associated p-values between standardized allele frequencies and genotypic values obtained at the population level. Manhattan plots and quantile-quantile plots drawn on the p-values associated with Spearman’s rho values indicate an absence of an excess of low p-values. To better describe the genetic architecture associated with PGR, notably the identification of QTLs with small effects, we then implemented a local score approach on the set of p-values (Fariello et al., 2017). The local score allows detection of significant genomic segments by accumulating the statistical signals derived from multiple adjacent SNPs, thereby limiting the number of tests performed while utilizing all the available data (Fariello et al., 2017; Apuli et al., 2021). By following Bonhomme et al. (2019); Aoun et al. (2020); Libourel et al. (2021); Demirjian et al. (2022) and Demirjian et al. (2023), we then implemented a local score approach (with tuning parameter ξ = 2) on these p values. Finally, significant SNP-phenotype associations were found by estimating a chromosome-wide significance threshold for each chromosome (Bonhomme et al., 2019).

For each of the three PGR traits, the candidate genes underlying the QTLs were retrieved using a custom script written under the R environment (Libourel et al., 2021). The lists of the candidate genes were then submitted to the classification SuperViewer tool on the University of Toronto website (http://bar.utoronto.ca/ntools/cgibin/ntools_classification_superviewer.cgi) using the MapMan classification, to allow the identification of biological pathways significantly over-represented (P< 0.05).

To test whether the SNPs underlying the QTLs identified by BHM-LS (hereafter named top SNPs) have suggestive signatures of local adaptation, we followed the method previously described in Brachi et al. (2015) for each of the three PGR traits. We looked for an over-representation of the top SNPs in the extreme upper tail of the XtX distribution obtained for the set of 168 natural populations of A. thaliana (Frachon et al., 2018). For a given SNP, XtX is a measure of the variance of the standardized population allele frequencies, which results from a rescaling based on the covariance matrix of population allele frequencies (Gautier, 2015). The formula used to calculate the fold enrichment in suggestive signatures of local adaptation was:

Here n is the number of SNPs in the upper tail of the XtX distribution. Here, we considered XtX statistic values as suggestive of local adaptation if they were among the top 1% of genome-wide XtX statistic values (i.e. 16,386 SNPs). n_a is the number of top SNPs that were also in the upper tail of the XtX distribution. N is the total number of SNPs tested genome-wide and Na is the total number of top SNPs. Following the methodology described in Hancock et al. (2011), the statistical significance of enrichment was assessed by running 10,000 null circular permutations across the five chromosomes of A. thaliana.

Despite the phyllosphere representing 60% of the total biomass on Earth and concentrating 10²⁶ bacteria (Vorholt, 2012), most GWAS carried out with non-pathogenic bacteria have focused on symbiotic bacteria or non-symbiotic PGPB isolated from the belowground compartment of plants (Stanton-Geddes et al., 2013; Wintermans et al., 2016; Curtin et al., 2017; Vidotti et al., 2019; Cotta et al., 2020; Torkamaneh et al., 2020). In this study, extensive genetic variation was observed among 162 natural accessions of A. thaliana in response to one strain of the native PGPB P. siliginis, which is an abundant and prevalent bacterial species in the leaf and root compartments of natural populations of A. thaliana located in the southwest of France (Ramírez-Sánchez et al., 2022). Because other strains of P. siliginis have been isolated from the leaf compartment of A. thaliana and characterized at the genomic level (Ramírez-Sánchez et al., 2022), it would be interesting to test whether the level of genetic variation of PGR and the underlying genetic architecture are similar among P. siliginis strains when inoculated at the seed stage. In addition, because P. siliginis and other Pseudomonas species with a PGPB effect, such as Pseudomonas moraviensis (Ramírez-Sánchez et al., 2022), belong to the subgroup Pseudomonas koreensis (Girard et al., 2021), it would be informative to check whether the genetics of PGR to P. siliginis extends to other phylogenetically close Pseudomonas species. Finally, in agreement with seed coating as an efficient way of introducing PGPB to seedlings (Ma, 2019; de Souza et al., 2020), we set up our GWAS by inoculating P. siliginis on seeds. Because the strength of the PGPB effect of P. siliginis on A. thaliana can depend on the developmental stage of the plants (Ramírez-Sánchez et al., 2022), it would be complementary to set up a GWAS by inoculating P. siliginis at the seedling stage.

Importantly, we observed a strong negative trade-off between plant growth in absence of P. siliginis and PGR to the OTU6_Psi_1 strain. To our knowledge, such a negative trade-off has not been reported in the literature. Identifying the mechanisms underlying this negative trade-off deserves further investigation. For instance, these mechanisms might rely on differences in seed size and physiology among the 162 accessions tested in this study. Beyond identifying the mechanisms, such a negative trade-off should promote the maintenance of genetic diversity at the underlying candidate genes, with the selection of growth-inductor responsive A. thaliana genotypes in presence of P. siliginis and the selection of growth-inhibitor responsive A. thaliana genotypes in absence of P. siliginis. Because the trade-off was observed both at the among-population and within-population levels, it suggests that the dynamics of A. thaliana - PGPB interactions should be studied at the metapopulation level rather than at the population level, as previously evidenced by studies on natural plant pathosystems such as Plantago lanceolata -Podosphaera plantaginis (Laine, 2005; Thrall et al., 2012) and A. thaliana - Pseudomonas syringae (Karasov et al., 2014).

Theoretical predictions suggest that the temporal regulation of QTLs often drives phenotypic changes in ontogenetic time, typically time-to-event or time-to-failure traits such as flowering time or death time (Johannes, 2007). Accordingly, previous GWAS performed on plant response to pathogens revealed temporal patterns in the detection of QTLs along the infection stages, with association peaks being detected only either at the earlier or at the later stages of infection (Aoun et al., 2017; Bartoli and Roux, 2017; Aoun et al., 2020; Demirjian et al., 2022; Demirjian et al., 2023). For instance, the atypical meiotic cyclin SOLO DANCERS gene was functionally validated in A. thaliana as conferring susceptibility to the bacterial pathogen Ralstonia solanacearum but only at the early stages of the infection (Aoun et al., 2020). Another A. thaliana gene, BWS1 (bacterial wilt susceptibility 1), was also revealed by GWAS as a susceptibility factor with a temporal dynamic in response to R. solanacearum (Demirjian et al., 2023). In this study, a similar dynamic in the detection of QTLs was observed for PGR to the OTU6_Psi_1 strain, suggesting that the PGPB effects conferred by P. siliginis depend on the time specificity of the genetic effects of A. thaliana. Importantly, in line with the mean PGPB effect of the strain OTU6_Psi_1 that increases over time, the enrichment in suggestive signatures of local adaptation of the candidate genes also increases over time, thereby highlighting the eco-evolutionary relevance of this native A. thaliana - P. siliginis interactions, similarly to the native interactions between A. thaliana and the bacterial pathogen P. syringae (Karasov et al., 2014; Roux and Bergelson, 2016). Our population genomics approach for identifying suggestive signatures of local adaptation across the genome of A. thaliana allows taking into account both the effect of selective processes at all life stages of A. thaliana while controlling for the effect of local demographic history (Frachon et al., 2017; Roux et al., 2023). This indirect approach is based on the calculation of the XtX statistics, analogous to F _ST but explicitly corrected for the covariance matrix of allele frequencies among populations (Gautier, 2015), and has been used to identify candidate genes associated with suggestive signatures of local adaptation in diverse species such as the European white oak (Leroy et al., 2020), the Aleppo pine Pinus halepensis (Ruiz Daniels et al., 2019), Caenorhabditis elegans (Crombie et al., 2022), the white-footed mice Peromyscus leucopus (Harris and Munshi-South, 2017) and the fungal wheat pathogen Zymoseptoria tritici (Hartmann et al., 2018). The PGPB effect observed on vegetative growth should therefore translate to fitness proxies such as total seed production in natural conditions, but remains to be tested. However, although a direct approach for testing local adaptation is based on setting up field experiments, in particular reciprocal field experiments, some fitness components such as germination rate can be hard to estimate under natural conditions (Savolainen et al., 2013). In addition, we must caution that estimating the effect of an adaptive allele on fitness proxies such as total seed production does not always predict the fate of the evolution of the frequency of this allele, as previously demonstrated for alleles conferring herbicide resistance (Roux et al., 2004; Roux et al., 2005; Roux et al., 2006; Vila-Aiub et al., 2011).

Of the 21 candidate genes highlighted, i.e. 14 genes in common between two or three scoring time points and seven genes of the three enriched biological processes, more than half of them are potentially involved in cell wall proliferation during seed and root development. Among these candidates, we observed both primary (i.e. expansins) and secondary (i.e. endoglycosylase/hydrolases) wall-loosening factors that are key players in cell wall structuring. For instance, expansins are cell wall-loosening proteins that directly induce cell wall extension by breaking non-covalent bonds between cellulose micro-fibrils and associated matrix polysaccharides in the cell wall (Lee and Kim, 2013). This study highlights EXP17 and xyloglucan endotransglucosylases/hydrolases, which facilitate LR emergence (Maris et al., 2009; Lee and Kim, 2013; Xu et al., 2020). XTH isoenzymes also strengthen the side-walls and cell walls of root hairs in the root differentiation zone after the completion of cell expansion (Maris et al., 2009). In addition, we identified GTs that catalyze protein glycosylation, a major post-translational modification of proteins, which significantly affects protein folding, conformation, distribution, stability and activity (Poulsen et al., 2015). Specifically, we identified a galactosyltransferase (Geshi et al., 2013) and three UGTs (UDP-glycosyltransferases), which are described to glycosylate various phytohormones and metabolites in response to biotic and abiotic stress in plants (Li et al., 2015; Rehman et al., 2018). Finally, studying PDX1 revealed that vitamin B₆ is essential for root development and stress tolerance (Chen and Xiong, 2005; Boycheva et al., 2015). Other candidates are involved in both growth development and plant defense. For instance, the candidate gene SBT3.5 may have a direct or indirect role in root and/or root hair development, particularly via the processing of PME7 in planta (Sénéchal et al., 2014), as PME are ubiquitous cell wall enzymes involved in important developmental processes (Micheli, 2001). Beyond the role of SBT3.5 in root development, another subtilase, SBT3.3, plays a role in immune priming during plant-pathogen interactions (Ramírez et al., 2013). In addition, the candidate gene JAZ11 inhibits A. thaliana hypersensitivity to the key phytohormone JA and represses susceptibility to Pst DC3000 (Liu et al., 2021).

The next step to understand the genetic and molecular mechanisms underlying the adaptive negative trade-off in response to P. siliginis observed in this study would be to phenotype (i) the mutant lines of the candidate genes for PGR of both leaves and roots, and (ii) the ability of the OTU6_Psi_1 strain to multiply in the leaf and root compartment of seedlings. In addition, it would be of particular interest to study the expression profiles with the spatial and subcellular localization of the candidate genes after inoculation with the OTU6_Psi_1 strain. Finally, exploiting the haplotypic diversity of the candidate genes among the 162 natural accessions of A. thaliana used in this study may help to identify the polymorphisms that have been selected in nature to respond to the PGPB P. siliginis.