Mapping Genetic Variation in Arabidopsis in Response to Plant Growth-Promoting Bacterium Azoarcus olearius DQS-4T

Plant growth-promoting bacteria (PGPB) can enhance plant health by facilitating nutrient uptake, nitrogen fixation, protection from pathogens, stress tolerance and/or boosting plant productivity. The genetic determinants that drive the plant–bacteria association remain understudied. To identify genetic loci highly correlated with traits responsive to PGPB, we performed a genome-wide association study (GWAS) using an Arabidopsis thaliana population treated with Azoarcus olearius DQS-4T. Phenotypically, the 305 Arabidopsis accessions tested responded differently to bacterial treatment by improving, inhibiting, or not affecting root system or shoot traits. GWA mapping analysis identified several predicted loci associated with primary root length or root fresh weight. Two statistical analyses were performed to narrow down potential gene candidates followed by haplotype block analysis, resulting in the identification of 11 loci associated with the responsiveness of Arabidopsis root fresh weight to bacterial inoculation. Our results showed considerable variation in the ability of plants to respond to inoculation by A. olearius DQS-4T while revealing considerable complexity regarding statistically associated loci with the growth traits measured. This investigation is a promising starting point for sustainable breeding strategies for future cropping practices that may employ beneficial microbes and/or modifications of the root microbiome.


Introduction
Plant growth-promoting bacteria (PGPB) benefit plant growth notably by enhancing nutrient uptake, nitrogen fixation, protection from pathogens, stress tolerance and/or boosting plant productivity. Changes in root system architecture and shoot biomass are common plant responses to PGPB, as evidenced on a variety of crops, including maize, rice, wheat, and various bioenergy grasses [1][2][3][4]. Despite an extensive literature documenting the beneficial effect of PGPB on plant growth, we still know relatively little about the molecular details of their mode of action. Genetic variation in the plant host has been reported to modulate the composition of the root-associated microbial population and has been suggested to have important adaptative consequences for plant health [5,6]. The process of domestication has profound consequences on crops, where the domesticate

Plant Growth Conditions and Bacterial Treatment
A collection of 305 natural accessions of Arabidopsis thaliana [30] (Supplementary Material Figure S1) was used to investigate the response to inoculation by Azoarcus olearius DQS-4 T . Arabidopsis seeds were sterilized by vortexing once with 70% ethanol for 1 min, twice with 70% ethanol plus 0.01% Triton X-100 for 2 min followed by once with 100% ethanol for 2 min. After ethanol removal, seeds were allowed to dry in a sterile hood and then 1 mL sterile water was added prior to vernalization at 4 • C for 3 days. Sterilized seeds were sown on square Petri dishes containing agar-solidified 1 2 Murashige and Skoog (MS) medium supplemented with 0.5% sucrose and incubated in a vertical position in a plant growth chamber at 21 • C with 16 h light/8 h dark cycle. After 5 days of germination, seedlings of similar size were transferred to circular Petri dishes containing 1 2 MS agar medium. Plants were inoculated by applying 150 µL of DQS-4 T culture containing 2 × 10 5 cell ml −1 onto the agar medium, approximately 5 cm below the root tip, which allowed the roots to grow into the inoculant. The same procedure was done for control treatments using 150 µL of saline solution (0.9% NaCl) without bacteria. The plates were briefly dried in a sterile hood, sealed with parafilm and placed vertically in a growth chamber until phenotypic analysis.

Phenotypic Response of Arabidopsis Accessions to Azoarcus Olearius DQS-4 T
For each accession, 5 seedlings were grown on a 1 2 MS agar plate. The reference accessions Col-0 and WS were used in each experiment since they represent non-responding and responding ecotypes, respectively. A total of 3 replicate plates (15 seedlings) of control and DQS-4T-treated seedlings were analyzed for each of the 305 accessions tested. Growth parameters were analyzed 7 days upon treatment by counting lateral root number and measuring primary root length (cm), then average primary root length and lateral root number per seedling were determined. Root and shoot fresh weight were analyzed 8 days after treatment. Data were acquired simultaneously for inoculated and control samples from the three replicates. To determine statistical significance between control and inoculated plants, one-way analysis (ANOVA) was used with Tukey's test (p-value < 0.05). Accessions with significant differences between control and inoculated were categorized as: (A) positive-genotypes that showed trait enhancement due inoculation; (B) negative-genotypes where inoculation inhibited growth; (C) non-responsive-genotypes that showed no growth response to bacterial inoculation.

Data Analysis
In this study, a natural population of 305 accessions of Arabidopsis thaliana [30] was used to investigate the genetic basis of the growth response to bacterial inoculation. Mapping analysis was performed on the following root parameters: primary root length (∆PRL), lateral root number (∆LRN), root fresh weight (∆RFW), and shoot fresh weight (∆SFW). For all traits, means per seedling (n = 5) per biological replicate (n = 3) were used to calculate the mean per treatment per accession. The mean value of the control treatment was subtracted from the value of the inoculated plants to generate the datasets used for GWAS analysis. Genome-wide association analysis employed a mixed linear model (MLM) using Tassel 5.0 software (http://www.maizegenetics.net/tassel, accessed on 15 November 2022) [36] incorporated with population structure by principal component analysis (PCA) and kinship matrix acquired from 1001 Genomes (https://1001genomes.org/, accessed on 15 November 2022) with minor allele frequency (MAF) = as 0.05, and the data were inferred as a normal distribution by the Kolmogorov-Smirnov test. All accessions were genotyped against the Col-0 reference genome with~214 k single-nucleotide polymorphism (SNPs) markers [37]. Significant SNPs were identified with a strict threshold of significance by Bonferroni correction with p-value = 2.34 × 10 −7 . Annotations of candidate genes were retrieved from TAIR10 (http://www.arabidopsis.org, accessed on 15 November 2022). In order to narrow the list of candidate SNPs to a more focused set of SNPs, the data were also analyzed using the two-stages method implemented in the R package GWAS.BAYES [38]. The two steps are called screening and model selection. The screening step of the GWAS.BAYES method performs a usual GWAS analysis with a linear mixed-effects models with a SNP fixed effect and kinship random effects. The screening step provides the usual list of significant SNPs. The model selection step of the GWAS.BAYES method performs a genetic algorithm search through model space, where candidate models are linear mixed-effects models with kinship random effects and may contain multiple SNPs. The genetic algorithm in the model selection step forces the SNPs to compete to appear in the highest ranked models. As shown in [39], combining a screening step and a model selection step provides a much shorter list of significant SNPs and leads to a much higher true discovery rate.
Manhattan plots and linkage disequilibrium (LD) plots were generated using the R statistical software 4.0 [40].

Validation Using Quantitative Reverse Transcription PCR (qRT-PCR)
Gene expression was evaluated by qRT-PCR. RNA was isolated from roots 7 days after mock or bacterial treatment using Direct-zol RNA kit treated with DNase (Zymo Research, Irvine, CA, USA) following the manufacturer's instructions. cDNA synthesis was performed with M-MLV reverse transcriptase (Promega, Madison, WI, USA). qRT-PCR was carried out as follows: 10 ng cDNA, 3 pM of each primer and PowerUp TM SYBR TM Green master mix (Applied Biosystems, Waltham, MA, USA) were mixed and amplified in Applied Biosystems ABI PRISM TM 7500 detection system (Applied Biosystems). Three biological replicates and three technical replicates for each transcript were analyzed using LinReg PCR 11.1 [41]. Quantitative amplifications were performed for different genes and ubiquitin 10 (At4G05320) was used as an internal reference. Primers used are listed in Supplementary Material Figure S2.

Arabidopsis Thaliana Response to Azoarcus olearius DQS-4 T Inoculation
Previous, published research, including from our own laboratory [23], documented that A. olearius DQS-4 T can produce strong, positive effects on root growth in both rice and Setaria viridis [24]. As a prelude to our larger GWAS study, we initially tested only a few Arabidopsis ecotypes regarding their response to DQS-4 T inoculation. These initial experiments revealed that A. thaliana ecotype Columbia (Col-0) had no measurable response to bacterial inoculation, while ecotype Wassilewskija (Ws) showed a robust and significant increase in all the traits measured (lateral root number, root and shoot fresh weight) (Supplementary Material Figure S3). Although limited, these initial experiments were important in showing phenotypic diversity in the plant response to inoculation, as well as providing both a negative (Col-O) and positive (Ws) control for future experiments.

Natural Variation in the Response of Arabidopsis Accessions to A. olearius Inoculation
In order to investigate the natural variation in the responsiveness of Arabidopsis to A. olearius DQS-4 T , a total of 305 Arabidopsis accessions were analyzed for changes in root and shoot traits upon bacterial treatment. Analysis of correlation to measure the direction and strength between control and DQS-4 T -treated samples showed a low correlation coefficient between control and treated samples for shoot fresh weight, primary root length and root fresh weight (R 2 = 0.452, R 2 = 0.349 and R 2 = 0.106, respectively). However, for lateral root number, the correlation was only slightly stronger (R 2 = 0.501) between control and treated samples ( Figure 1A-H), suggesting that the magnitude of these DQS-4 T -induced growth responses were weakly related to the intrinsic growth capacity for these parameters under the tested conditions. Hence, faster-growing accessions or accessions that form more lateral roots in the experimental setup are not necessarily stronger responders to bacterial treatment.

Natural Variation in the Response of Arabidopsis Accessions to A. olearius Inoculation
In order to investigate the natural variation in the responsiveness of Arabidopsis to A. olearius DQS-4 T , a total of 305 Arabidopsis accessions were analyzed for changes in root and shoot traits upon bacterial treatment. Analysis of correlation to measure the direction and strength between control and DQS-4 T -treated samples showed a low correlation coefficient between control and treated samples for shoot fresh weight, primary root length and root fresh weight (R 2 = 0.452, R 2 = 0.349 and R 2 = 0.106, respectively). However, for lateral root number, the correlation was only slightly stronger (R 2 = 0.501) between control and treated samples ( Figure 1A-H), suggesting that the magnitude of these DQS-4 T -induced growth responses were weakly related to the intrinsic growth capacity for these parameters under the tested conditions. Hence, faster-growing accessions or accessions that form more lateral roots in the experimental setup are not necessarily stronger responders to bacterial treatment.

Response Categories of Arabidopsis thaliana Growth Traits to Treatments
The measured, phenotypic variability, relative to lateral root number (LRN), primary root length (PRL), root and shoot fresh weight (RFW, SFW), among the 305 Arabidopsis accessions classified into 3 categories are: 1. positive responsive: accessions that demonstrate a significant, positive change upon inoculation compared to mock samples; 2. nonresponsive: accessions where bacterial treatment did not affect growth; 3. negative responsive: accessions where treatment inhibited growth ( Figure 2). Within the positive responsive category, eleven genotypes showed growth enhancement in all four parameters analyzed, where most of the changes were statistically significant in lateral root number followed by shoot fresh weight and root fresh weight (Figure 2A,C). The growth of 13 ecotypes clearly benefited from DQS-4 T inoculation ( Figure 2B,C). For instance, the genotype Ws showed an increase in biomass and root growth when inoculated ( Figure 2C).

Response Categories of Arabidopsis thaliana Growth Traits to Treatments
The measured, phenotypic variability, relative to lateral root number (LRN), primary root length (PRL), root and shoot fresh weight (RFW, SFW), among the 305 Arabidopsis accessions classified into 3 categories are: 1. positive responsive: accessions that demonstrate a significant, positive change upon inoculation compared to mock samples; 2. non-responsive: accessions where bacterial treatment did not affect growth; 3. negative responsive: accessions where treatment inhibited growth ( Figure 2). Within the positive responsive category, eleven genotypes showed growth enhancement in all four parameters analyzed, where most of the changes were statistically significant in lateral root number followed by shoot fresh weight and root fresh weight (Figure 2A,C). The growth of 13 ecotypes clearly benefited from DQS-4 T inoculation ( Figure 2B,C). For instance, the genotype Ws showed an increase in biomass and root growth when inoculated ( Figure 2C).
However, the largest number of ecotypes (129) fell within the non-responsive group, which showed no statistically significant positive or negative response to bacterial inoculation for the parameters analyzed ( Figure 2D-F). As mentioned above, Col-0 is a good example of an ecotype within this non-responsive category ( Figure 2F). Interestingly, 82 genotypes showed a negative growth response to bacterial inoculation ( Figure 2G). An example is ecotype PHW-35 ( Figure 2H-I), which showed a negative response in each of the However, the largest number of ecotypes (129) fell within the non-responsive group, which showed no statistically significant positive or negative response to bacterial inoculation for the parameters analyzed ( Figure 2D-F). As mentioned above, Col-0 is a good example of an ecotype within this non-responsive category ( Figure 2F). Interestingly, 82 genotypes showed a negative growth response to bacterial inoculation ( Figure 2G). An example is ecotype PHW-35 ( Figure 2H-I), which showed a negative response in each of the four parameters measured. However, when present, the positive effects on root architecture can be dramatic ( Figure 3); for example, as shown in Figure 3A for ecotype Kr-0 ( Figure 3A). Other ecotypes, such as Ws, CIBC-4, Si-0, Hodvala-2 and wag-3, demonstrated similar increases in root fresh weight upon bacterial treatment ( Figure 3B) and could be easily identified from non-responsive ecotypes, such as Col-0, DUK, N7, Sap-0, Zdrl2-24 and MIB-15 ( Figure 3C,D). However, it is important to note that it was common to find ecotypes that were not consistently responsive for all parameters measured ( Figure 3E,F). For instance, ecotype Aa-0 showed a significant increase in LRN upon inoculation but was non-responsive for the other root and shoot parameters. Another example of trait-related variation is exemplified by ecotype Hod which showed no significant changes to LRN and PRL; however, RFW was significantly reduced by inoculation while increasing shoot biomass. This variability within overall response, responses in individual parameters, and opposing responses (i.e., negative, and positive) suggests significant underlying complexity in the genetics of the plant response, as well as the molecular mechanisms involved. However, unlike Aa-0 and Hod, some ecotypes gave robust responses for all four traits analyzed, including ecotypes Ws, Bla-1, KI-5, Kr-0 or showed a consistent, negative response, such as ecotypes UKNW-06-460, UKSE 06-349, PHW-35 and PHW-37 (Table 1 and Supplementary Material Figure S1 for a complete dataset). Other ecotypes, such as Ws, CIBC-4, Si-0, Hodvala-2 and wag-3, demonstrated similar increases in root fresh weight upon bacterial treatment ( Figure 3B) and could be easily identified from non-responsive ecotypes, such as Col-0, DUK, N7, Sap-0, Zdrl2-24 and MIB-15 ( Figure 3C,D). However, it is important to note that it was common to find ecotypes that were not consistently responsive for all parameters measured ( Figure 3E,F). For instance, ecotype Aa-0 showed a significant increase in LRN upon inoculation but was non-responsive for the other root and shoot parameters. Another example of trait-related variation is exemplified by ecotype Hod which showed no significant changes to LRN and PRL; however, RFW was significantly reduced by inoculation while increasing shoot biomass. This variability within overall response, responses in individual parameters, and opposing responses (i.e., negative, and positive) suggests significant underlying complexity in the genetics of the plant response, as well as the molecular mechanisms involved. However, unlike Aa-0 and Hod, some ecotypes gave robust responses for all four traits analyzed, including ecotypes Ws, Bla-1, KI-5, Kr-0 or showed a consistent, negative response, such as ecotypes UKNW-06-460, UKSE 06-349, PHW-35 and PHW-37 (Table 1 and Supplementary Material S1 for a complete dataset). Statistical analysis was carried out using one-way ANOVA with Tukey's test p-value < 0.05 (* p ≤ 0.05; ** p ≤ 0.01; *** p ≤ 0.001). RFW = root fresh weight. Photographs were taken by scanning the plates using a photo scanner with resolution 640 × 480.

Genome-Wide Association Loci Mapping in the Arabidopsis Population
The dataset collected for DQS-4 T -induced changes in lateral root number (∆LRN), primary root length (∆PRL), root fresh weight (∆RFW) and shoot fresh weight (∆SFW) were averaged and analyzed against the Col-0 reference genome using a mixed linear model (MLM) algorithm in Tassel 5.0 Software. To determine the distribution and quality of the data within the population a quantile-quantile plot (Q-Q plot) was generated for each growth trait (Supplementary Material Figure S4). The GWA mapping results showed highly significant SNPs for two traits ∆PRL and ∆RFW (Figure 4). No significant SNPs were significantly associated to ∆LRN and ∆SFW (Supplementary Material Figure S4). With a threshold of −log 10 (P) > 7 adjusted by Bonferroni correction, a total of 63 loci were detected for root fresh weight and 55 loci correlated to primary root length (Supplementary Material Figure S5). We observed only one SNP associated with both traits ∆RFW and ∆PRL, mapping close to the gene encoding AtFKGP, bifunctional fucokinase/fucose pyrophosphorylase (At1G01220, p-value = 2.19 × 10 −8 and 3.32 × 10 −7 , adjusted by Bonferroni correction, respectively). Given this close association, we measured the expression of AtFKGP by qRT-PCR using mRNA extracted from roots of 7-day-old seedlings of ecotypes representing non-responsive (Col-0, N7), responsive (Ws, Kr-0 and CIBC-4), and negative responding lines (PHW-37 and Lis-1) either mock non-inoculated or inoculated with DQS-4 T . However, this experiment failed to find either down or upregulated gene expression of AtFKGP in response to bacterial inoculation in any of the Arabidopsis accessions tested (Supplementary Material).

Primary Root Length Highly Correlated SNPs and Candidate Genes
We identified 55 SNPs highly correlated to measured changes in PRL (Supplemental Material S5). An example are SNPs mapping close to gene At1G33410, encoding the suppressor of auxin resistance 1 (sar1) gene. SAR1 and SAR3 are proteins similar to vertebrate nucleoporins that are part of the nuclear pore complex (NPC). Plants deficient in either protein exhibit pleiotropic growth defects partially affecting the translocation of proteins involved in hormonal signaling and plant development [42,43]. A large LD resulted in inclusion of many polymorphisms in this candidate region; for example, loci within a gene

Primary Root Length Highly Correlated SNPs and Candidate Genes
We identified 55 SNPs highly correlated to measured changes in PRL (Supplementary Material Figure S5). An example are SNPs mapping close to gene At1G33410, encoding the suppressor of auxin resistance 1 (sar1) gene. SAR1 and SAR3 are proteins similar to vertebrate nucleoporins that are part of the nuclear pore complex (NPC). Plants deficient in either protein exhibit pleiotropic growth defects partially affecting the translocation of proteins involved in hormonal signaling and plant development [42,43]. A large LD resulted in inclusion of many polymorphisms in this candidate region; for example, loci within a gene predicted to encode a tetratricopeptide repeat 9 (TPR9) protein involved in gibberellic acid regulation [44], fascinated stem 4 (Atfas4) protein and a Ring/Ubox super family (At1g01660) protein. Moreover, a SNP highly correlated to primary root length corresponded to the gene encoding a late elongate hypocotyl LHY1, a MYB-related putative transcription factor implicated in circadian regulation of flowering time [45].

Root Fresh Weight Highly Correlated SNPs and Candidate Genes
We identified 65 SNPs highly correlated to changes in RFW Supplementary Material Figure S5). In order to narrow the list of candidate SNPs to a more focused set of SNPs, the data were also analyzed using the two-stages procedure implemented in the R package GWAS.BAYES [38,39]. This analysis reduced the number of SNPs highly associated with ∆RFW to 11 (Table 2); however, no SNPs were highly associated with ∆PRL in this analysis. In most GWA studies, the highly trait-correlated SNPs can present alleles in linkage disequilibrium (LD) at two or more loci in a population. However, because the peak of selection signals is relatively large in the GWA peak region, it is difficult to conclude whether the target of selection is the causative SNP or other alleles are significantly associated with this genome region. To examine the relationship between SNPs and regions, a haplotype block was generated for 10 kb upstream or downstream of the most significantly associated SNP at position 13,459,922 on chromosome 4 observed for ∆RFW ( Figure 5A,B). The selected SNP 13,459,922 (p = 7.55 × 10 −9 ; At4G26690) mapped close to a glycerophosphodiester phosphodiesterase-like 3 GDPDL3 gene, also known as Shaven 3 (SHV3) that is involved in cell wall organization and root hair growth [46]. This SNP was identified as significant in both statistical methods used. SNP 13,459,922 is located at an intronic 5 untranslated region (5 UTR) with a modifier predicted effect [47]. Next, we determined which gene in proximity to this highly associated genomic region underlies the variation. As shown in Figure 5, six haplotype blocks were significantly associated with SNP 13459922, At4G26690. Based on haplotype analysis, seven SNP regions were identified associated, respectively, with genes encoding a member of the vacuolar-type ATPase family (At4G26710) and a protein phosphatase x-1(PP4) (At4G26720), extensin-like protein (Lip5, At4G26750), microtubuleassociated protein 65-2 (MAP65-2, At4G26760), phosphatidate cytidylyltransferase (CDS3, At4G26770), mitochondrial GRPE 2 (At4G26780) and Lipase acyl hydrolase superfamily (GDSL-like, At4G26790) ( Figure 5C). Because the polymorphisms are difficult to identify, we also carried out a matrix analysis, that showed SNPs at position 13,477,249 (gene Lip5) and 13,483,356 (CDS3 gene) highly correlated to SNP 13459922, GDPDL3 corroborating the haplotype block (Supplementary Material). Next, we determined the expression level of the genes correlated to the GWA peak by quantitative RT-PCR (qRT-PCR) using mRNA isolated from root tissue. We assumed that the expression of a candidate gene might be altered in accessions of different responsive categories. Hence, we extracted mRNA from the roots of selected accessions responsive (Ws and Kr-0), non-responsive (Col-0 and N7) and negative response (PHW-37 and Lis-1) regarding the RFW trait. The data show that none of the accessions tested showed a significant change in expression level upon inoculation for the various genes tested (i.e., those encoding GDPDL3, Lip5 and CDS3) (Supplementary Material). Of course, the lack of a transcriptional response to bacterial inoculation does not rule out the possibility that a specific gene could be playing an important role in the response to A. olearius DQS-4 T treatment.

Discussion
The benefits of PGPB in promoting plant growth, improving nutrient uptake and plant resilience to biotic and abiotic stress, and boosting crop production are documented by an expansive literature [3,[48][49][50]. While the molecular mechanisms and specific pathways that underlie the growth-promoting responses in plants by PGPB have been investigated to some extent regarding the bacterial functions, left largely unexplored are the plant functions involved. Better defining these functions is important since they may help address the problems of consistency and efficiency that are found commonly when PGPB are used under field conditions to enhance crop yield and sustainability [51][52][53]. GWAS is now a popular method to harness natural genetic variation in a population to identify genetic loci critical for specific agronomic traits and in support of breeding improvement programs [13,[54][55][56]. Although used with great success in many studies, classical GWAS relying on SNPs has its limitations due to 'missing heritability' [57]. Failure to capture rare variants, allelic heterogeneity, epistasis, and/or epigenetic variation often decreases the detecting capacity of GWAS [58][59][60][61][62][63]. To test the feasibility of this approach to investigate PGPG-plant interactions, we applied GWAS to map loci within the model plant Arabidopsis crucial for the beneficial response to the PGPB Azoarcus olearius DQS-4 T . Such an

Discussion
The benefits of PGPB in promoting plant growth, improving nutrient uptake and plant resilience to biotic and abiotic stress, and boosting crop production are documented by an expansive literature [3,[48][49][50]. While the molecular mechanisms and specific pathways that underlie the growth-promoting responses in plants by PGPB have been investigated to some extent regarding the bacterial functions, left largely unexplored are the plant functions involved. Better defining these functions is important since they may help address the problems of consistency and efficiency that are found commonly when PGPB are used under field conditions to enhance crop yield and sustainability [51][52][53]. GWAS is now a popular method to harness natural genetic variation in a population to identify genetic loci critical for specific agronomic traits and in support of breeding improvement programs [13,[54][55][56]. Although used with great success in many studies, classical GWAS relying on SNPs has its limitations due to 'missing heritability' [57]. Failure to capture rare variants, allelic heterogeneity, epistasis, and/or epigenetic variation often decreases the detecting capacity of GWAS [58][59][60][61][62][63]. To test the feasibility of this approach to investigate PGPG-plant interactions, we applied GWAS to map loci within the model plant Arabidopsis crucial for the beneficial response to the PGPB Azoarcus olearius DQS-4 T . Such an approach has been used previously; for example, to examine the response of an Arabidopsis natural population to inoculation with the PGPB Pseudomonas fluorescens. This study identified 10 potential genes candidates involved in changes root architecture and shoot biomass but found none strongly correlated to growth responses to bacterial inoculation with no common gene that could be correlated to a PGPB mediated effect [21].
A few general conclusions can be made from our study. Consistent with published reports, including some from our own studies [2,23,64] plant genotype largely determines whether a given PGPB strain will or will not enhance or inhibit plant growth. This could be a major factor in field-to-field variation in published PGPB studies [51][52][53]. However, perhaps most impactful, is that the data point to considerable complexity in the mechanisms that underlie a beneficial plant response to PGPB inoculation. For example, within the Arabidopsis population, 27% of ecotypes showed no response to inoculation, while others showed either a negative (12% PRL, 4% LRN, 6% RFW and 6%SFW) or positive (8% PRL, 13% LRN, 10% RFW and 11% SFW) response to a specific trait. Considering the four growth parameters tested, 11 ecotypes showed consistently positive response whereas 6 accessions responded negatively. Indeed, some ecotypes showed different responses regarding a specific phenotypic parameter. For instance, ecotype Hod showed increased shoot fresh weight while root biomass was negatively affected. This complexity correlates well with our recent, mutational analysis of two PGPB strains that suggested that the gene functions necessary for plant root colonization are unique to a given strain, with only a few genes appearing essential for both strains tested. This large variation, coupled with normal issues found when applying biological inoculation to cropping systems, could, in large part, explain why it is not uncommon to find very variable, inconsistent results when PGPB are used under field conditions [51][52][53].
The power of the GWAS method is the ability to statistically correlate specific genetic loci with the phenotypic variation measured across the entire population [65]. Hence, we used two, orthogonal methods for data analysis with increased statistical stringency. For example, our initial analysis using a mixed-effects model identified several candidate SNPs associated with changes in primary root length and root fresh weight. However, in root length none of those SNPs were significant when submitted to a stricter statistical analysis. We encountered two major challenges to identify statistically significant SNP associations with the phenotypes measured. First, although several statistically robust models have been developed, false positives can still arise from population structure. Second is the large extended, linkage disequilibrium (LD) which results in the inclusion of many candidate genes within a single LD block, making it if very hard to ultimately identify the underlying, causative gene. For example, this might explain why none of the candidate genes associated with ∆PRL were found to respond transcriptionally to PGPB inoculation.
Among the candidate genes identified were a phosphate transporter essential for arbuscular mycorrhizal symbiosis [66] and glutamine synthetase (GS) known for its role in assimilation of nitrogen and biosynthesis of glutamine in plants. Transcriptional levels of the GS2 gene showed no changes in sugarcane leaves inoculated with several PGPB bacteria strains. However, these authors did detect biochemical differences relative to concentrations of glutamine and glutamate [67]. The ankyrin repeat family protein, was previously found in a study that characterized genes that contribute to bacterial adaptation isolated from Arabidopsis, maize, and poplar trees [68]. In rice, ankyrin-repeat protein was shown to be upregulated in response to rice leaf bright pathogen Xanthomonas oryzae (xoo) suggesting its role as a positive regulator in basal defense pathways [69].
Unfortunately, ecotype Col-0 did not respond to inoculation for any of the growth parameters measured. Hence, all the genetic resources available for this ecotype were of little use for follow-up genetic studies. For example, the gene encoding Lyst-interacting protein 5, (LIP5) was within the haplotype associated with ∆RFW. Arabidopsis LIP5 is part of the endosomal sorting complexes required for transporters (ESCORTs) for sustained protein trafficking. Lyst-interacting protein 5 (LIP5) interacts with MAP kinase 3 (MPK3) and MPK6 in response to pathogen infection, playing a critical role in plant basal resistance [70,71]. Mutational disruption of LIP5 expression had little effect on pathogen associated molecular pattern (i.e., flagellin) or salicylic acid-induced defense responses but compromised basal resistance in response to the bacterial pathogen Pseudomonas syringae [70]. We obtained the available T-DNA insertion line of AtLIP5 within the Col-O background and found, not surprisingly, that this mutation had no effect on the response of seedlings to A. olearius DQS-4 T inoculation (data not shown). Hence, confirmation of those candidate genes implicated by our GWAS analysis in the PGPB response awaits the ability to use gene-editing to create mutations in those specific ecotype backgrounds that do respond to inoculation.
In summary, our study of the natural genetic variation within an Arabidopsis population showed considerable variation in the ability of plants to respond to inoculation by A. olearius DQS-4 T while revealing considerable complexity regarding statistically associated loci with the growth traits measured and in the patterns (positive, neutral, negative) of those responses. Considering all candidate genes with SNP-trait associations in the GWA analysis, several have known or predicted functions that hold promise for being functional in mediating PGPB growth effects.
Supplementary Materials: The following supporting information can be downloaded at: https://www. mdpi.com/article/10.3390/microorganisms11020331/s1, Figure S1: Differential growth responses of Arabidopsis accessions to the beneficial bacterium A. olearius DQS-4 T . A. Bacterial cells recovered from roots of Col-0 and WS five days after treatment with DQS-4 T . No bacterial cells were recovered from the mock treatment. Results are shown in log of colony forming units (CFU) per gram of fresh root tissue. Bars are the average of root samples per treatment (n = 15). B. Growth parameters analyzed 10 days after DQS-4T treatment in Col-0 and WS accessions. Bars are average of samples (n = 32). Statistical significance is represented by ** p > 0.01, *** p > 0.001 calculated by one-way Anova; Figure S2: Manhattan plot representing SNPs upon GWA mapping. A. Lateral root number (LRN) and B. Shoot fresh weight (SFW), no significant SNPs were found at a threshold -log 10 (P) = 10 −7 for both traits. C. Q-Q plot representation of data distribution within 305 Arabidopsis accessions. Each color represents a different set of data in -log 10 (p-value), red = PRL, yellow = RFW, blue = LRN and green = SFW. According to Q-Q plot PRL and RFW showed high variability in distribution within the data, unlike LRN and SFW; Figure S3: Relative gene expression of At1G01220, L-fucose gene. RT-qPCR analysis of mRNA expression from root tissue of non-responsive genotypes (Col-0 and N7) and responsive genotypes (Ws, Kr-0 and CIBC4) 7 days after inoculation. The values are expressed as the relative quantity with respect to the control. The data represent the average of three independent experiments; Figure S4: Pearson's correlation matrix 10 kb region surrounds a highly significant SNP 13459922 at chromosome 4. SNPs in LD with 13459922 are shown in the table with SNPs physical location; Figure S5: Relative gene expression of At4G26690, At4G26750 and At4G26770. RT-qPCR analysis of mRNA expression from root tissue of non-responsive genotypes (Col-0 and N7), responsive (Ws, Kr-0 and CIBC4) and negative response (PHW-37 and Lis-1) 7 days after inoculation. Lis-1 showed no expression detected in At4G26690 and AtG426770 genes. The values are expressed as the relative quantity with respect to the control. The data represent the average of three independent experiments.