Natural Genetic Variation Shapes Root System Responses to Phytohormones in Arabidopsis

Plants adjust their architecture by modulating organ growth. This ability is largely dependent on phytohormones. While many genes involved in phytohormone pathways have been identified, it remains unclear to which extent and how these pathways are modulated in non-reference strains and whether this is relevant for local adaptation. Here we assess variation of root traits in response to perturbations of the auxin, cytokinin, and abscisic acid pathways in 192 Arabidopsis accessions. We identify common response patterns, uncover the extent of their modulation by specific genotypes, and find that the Col-0 reference accession is not a good representative of the species in this regard. We conduct GWAS and identify 114 significant associations, most of them relating to ABA treatment. The numerous ABA candidate genes are not enriched for known ABA associated genes indicating that we largely uncovered unknown players. We then study two associated regions in detail and identify the CRF3 gene as a modulator of multiple hormone pathways. Finally, we show that natural variation in root traits is significantly associated with climate parameters relevant for local adaption in Arabidopsis thaliana and that, in particular, ABA regulated lateral root traits are likely to be relevant for adaptation to soil moisture. Author Summary The root system is a key component for plant survival and productivity. Apart from anchoring the plant, its architecture determines which parts of the soil are foraged for water and nutrients, and it serves as an interface for interaction with microbes and other soil organisms. Plant hormones play crucial roles in the development of root system architecture and its plasticity. However, while there is substantial natural variation of root architectures, it is not clear to which extent genetic variation in hormone related pathways contribute. Using the model species Arabidopsis thaliana we quantitatively explore the breadth of natural variation in plant hormone responses to three major plant hormones: auxin, cytokinin, and abscisic acid. The Col-0 reference strain can be quite different from a large proportion of the natural accessions of the species, illustrating a severe caveat of relying on a single reference strain in model species and drawing generalizations from it. Using GWAS, we further identify a large number of loci underlying the variation of responses to plant hormones, in particular to abscisic acid, find links between local adaptation and root responses to hormones, and finally using mutants for GWAS candidate genes, identify novel players involved in regulating hormone responses.


INTRODUCTION 50 51
To grow and survive plants need access to sunlight, nutrients and water resources that are 52 not evenly distributed in the environment. As sessile organisms plants adjust their architecture 53 according to the distribution of resources by modulating organ growth and development. 54 response patterns with only some accessions showing strong deviations with regard to specific 127 hormonal perturbations or traits. Surprisingly, Col-0, the reference accession on which most of 128 the previous studies of hormonal effects had been performed, is among the accessions whose 129 root growth responses are frequently quite different from the bulk of Arabidopsis accessions. On 130 average it is in the upper quartile of the accession distribution and belonged to the 1% most 131 extreme phenotypes for some traits in some conditions (Supplemental Figure 3, Supplemental 132 Table 3). Col-0 is therefore not a good representative of the diverse panel of accessions that we 133 have investigated. 134 In conclusion, we have generated the first comprehensive atlas of root responses to perturbation 135 of hormonal pathways. Our data show that while there are clear trends of phenotypic responses 136 to perturbation of specific hormonal pathways, there is significant natural variation in these 137 responses. Overall, this suggests that hormonal signaling pathways within the Arabidopsis 138 thaliana species are largely acting on the same traits, but that the extent of the hormonal control 139 of specific traits is genotype dependent. We therefore conclude that hormonal pathways are 140 subject to natural genetic variation. Importantly, the effects of hormones on root growth in the pathways affected a particular trait correlation, we also calculated the variance of each trait 154 correlation between all conditions. We reasoned that if this variance was high, the observed trait 155 correlation would be strongly dependent on a subset of phytohormone pathways ( Figure 2E). The 156 highest variation of correlations between independent traits was observed for the length of branching zone and the lateral root density of the primary root (σ=0.13). While these traits are 158 highly correlated in control (0.67), CK (0.63) and ABA (0.67) conditions, this correlation 159 reverses upon IAA treatment (-0.31). This strong effect can be traced back to the impact of IAA 160 on root growth rate and its impact on the branching zone, which are key for the second and third 161 ranked correlation variances respectively. Under control conditions, there is a very high positive 162 correlation between primary root growth rate before treatment and growth rate after treatment 163 (r=0.9), indicating that root length after 7 days is generally a very good predictor for root length 164 after 10 days. IAA treatment, however, breaks the correlation of these two traits almost 165 completely (r=0.19). ABA and CK treatments diminish this correlation only slightly, suggesting 166 that the observed effect of IAA on root growth rate is highly specific. Overall, this strongly 167 suggests that the regulation of root growth rate is strongly dominated by the auxin pathway, as 168

Genotypes can Alter Phenotypic Responses to Hormones 183
While our correlation analysis revealed common trends over a large number of genotypes, our 184 atlas of root responses to perturbation of hormonal pathways also allowed us to study which 185 genotypes modulate the responses to perturbations of hormone pathways. To do so, we 186 performed a hierarchical-clustering of our ten traits separated by condition. Different clusters 187 contain groups of accessions with genotypes that generate a similar RSA under the respective conditions ( Figure 3, and Supplemental Table 4). We observed that perturbation of specific 189 hormone pathways shift RSA to a distinct state (for instance, IAA causes roots to become shorter 190 and more branched) but that the genotype determines the extent to which this state is reached. 191 For instance, IAA treatment cluster 5 is shifted to the IAA RSA state to a lesser degree than IAA 192 cluster 6. Overall, our analysis partitions the phenotypic space that we explored using the 193 systematic perturbation of hormone pathways and illustrates the tremendous impact of the 194 genotype in determining how root development responds to hormones. However, while this 195 analysis is highly useful for identifying classes of accessions 196 and visually illustrates the effects of hormones and how these effects differ, it was based on 197 capturing RSA using individual traits, some of which are dependent on each other and some of 198 them not. Moreover, hierarchical clustering gives insight into classes, rather than continuous 199 similarities. To gain a more integrative insight into the action of hormones, we conducted a accessions and four treatments). 99% of RSA variation (given by our 10 traits) could be captured 202 by five PCs, while 9 PCs are needed to explain the complete variation (Supplemental Table 5). 203 The first PC captures 62% of the total variation. Several traits contribute to it, including root 204 number (LR.No), lateral root density (LRD_P and LRD_R), total lateral root length (TLRL), and 205 length ratio (LRR). The second PC accounts for 26% of the variation and represents mainly 206 primary root length (P) and total root length (TRL) (Supplemental  Table 6). For each 217 accession, we then calculated the average Euclidian distance of all conditions (Supplemental 218 Table 6). This average represents a measure of how profoundly a genotype differs from the norm 219 in altering its RSA in response to all hormonal perturbations that we had performed. Eight 220 accessions (Sh-0, Lu-1, Col-0, Sorbo, UKSE06-373, Kz-1, Or-0, Uk-1) were more than 2 221 standard deviations from the mean of all accessions in this regard ( Figure 4A, Supplemental 222 Table 6), indicating that they show significant alterations of their hormonal response from the 223 average Arabidopsis thaliana accession. Again, the Col-0 reference accessions was among these 224 accessions, again highlighting the conclusion that much of what has been learned in studies on 225 Col-0 may not represent the species as a whole. Taken together, these data demonstrate that 226 hormone perturbations generally lead to specific patterns of developmental responses that shift 227 RSA from one configuration to another (i.e. IAA shifts RSA to shorter and more branched 228 configurations). However, these response patterns can be significantly altered by genotypic 229 variation. This genotypic variation must therefore be in the pathways that perceive or respond to 230 phytohormone cues. Overall, this demonstrates that natural genetic variation in phytohormone 231 pathways can be a major contributor in shaping the Arabidopsis root system. 232 233

Genetic Variation leads to Changes in the Molecular Regulation of Hormone Pathways 234
To test whether the phenotypic differences we observed are reflected at the molecular level, we 235    Figure 5). 265

266
Despite the large number of significant associations for the ABA pathway, the genes close to 267 these associations didn't contain obvious core ABA signaling related genes. To test whether this 268 absence of bona fide ABA candidate genes was statistically significant, we identified all genes 269 for which a role in the ABA pathway had been assigned based on reported experimental 270 evidence or mutant phenotypes (Supplemental table 8). There was only one gene that overlapped 271 between GWAS-based candidate gene list and the annotation-based ABA bona fide candidate 272 gene list. We then tested whether this overlap is expected by chance, we found that our observed 273 overlap was at the lower limit (towards our GWAS list being depleted of bona fide candidate 274 genes) of what was to be expected by chance (p<0.16; Supplemental Figure 6). This clearly 275 showed that bona fide candidate genes are not overrepresented in our list of GWAS candidates 276 and indicates that, for complex traits such as root responses to hormonal treatment, genes 277 targeted by natural variation are frequently not the genes found by traditional experimental To identify high confidence candidate genes whose variants underlie the associated traits, we 281 focused on two cases. The first case involved the trait LRR (ratio of total lateral root length and 282 primary root length) where a significant association on chromosome 3 was found under both  Figure 7A). When analyzing linkage disequilibrium (LD) in this region, we found 290 that SNPs in only two genes (AT3G17890 and AT3G17900) were in LD (R 2 >0.4) with the SNP 291 that we had identified through GWAS. We tested whether the expression of any gene in this 292 region was correlated with the LRR trait. Only the three genes in LD showed expression levels 293 correlated with the LRR phenotype. Of these three genes, AT3G17910 showed the highest 294 correlation ( Figure 5C, Supplemental Figure 7B). JAI3, the most likely candidate gene based on 295 its known role, did not show any correlation. To test whether one or more of these genes are 296 involved in regulating lateral root related traits, we obtained available T-DNA insertion mutant 297 lines for AT3G17890 (SALK_050958), AT3G17910 (SALK_045920) and JAI3 (jai3-1) (Chini 298 et al., 2007). No insertion line in the genic region of AT3G17900 was available. The expression 299 levels of the respective genes were significantly reduced in the SALK lines that we obtained 300 (Supplemental Figure 7C). We examined the phenotypes of these three mutants under control 301 and cytokinin treatment conditions. While the jai3 line showed an overall reduced root growth in 302 both conditions, the response to CK was altered for multiple traits in T-DNA mutants of the 303 other two genes ( Figure 5D and 5E). In particular, alterations were observed in both mutants for 304 the branching zone (R) and lateral root number (LR.No) traits, as well as for primary root length 305 (P) and total root length (TRL) in one mutant (SALK_045920). While compared to Col-0, the 306 SALK_045920 (AT3G17910) line showed a decreased branching zone under control and an perturbation of gene expression in the SALK_045920 line. Overall, this strongly suggests that 312 the yet uncharacterized genes AT3G17890 and AT3G17910 are involved in regulating aspects of 313 CK hormonal signaling that are relevant for shaping the same root traits in a complex manner. 314 We note that we can't evaluate the role of the third gene, AT3G17900, which is also a potential 315 candidate based on the LD pattern and the expression level correlation. 316 In contrast to the first case, a well-known candidate gene was identified for another trait. Here, a 317 CRF3 has been shown to be a key link between cytokinin signaling and auxin transport in root 321 development, acting downstream of cytokinin perception and transcriptionally regulating auxin 322 transport (PINs) (Simaskova et al., 2015). As this region was devoid of other genes and a SNP in 323 the coding region of CRF3 was in LD (R 2 >0.4), this was the only and best candidate gene. To 324 test whether we would obtain phenotypes that consistent with the GWAS phenotypes in our 325 conditions, we obtained a loss of function T-DNA line (SAIL_240_H09) for CRF3 that has been 326 described previously (Simaskova et al., 2015). We found that this line had an overall reduced 327 root growth, affecting both root length and lateral root number, when grown under control and 328 abscisic acid (ABA) conditions ( Figure 6C and 6D). While the most pronounced LR effect upon 329 ABA in the wild type is a strong reduction of the lateral root length (LRL and TLRL), these traits 330 were not affected by ABA in the crf3 mutant line, leading to a highly significant interaction of 331 the crf3 mutation and ABA treatment ( Figure 6C and 6D). Moreover, expression levels of genes 332 related to auxin, cytokinin and ABA hormonal pathways were significantly altered in the crf3 333 line (Supplemental figure 8B). Taken together, this shows that CRF3 not only links auxin and 334 cytokinin responses, but also impinges on ABA dependent root growth regulation. 335 336

Natural Variation in Hormone Signaling Pathways and Local Adaptation 337
As natural variation in pathways that affect root growth responses to phytohormones has a 338 tremendous impact on root system architecture related traits, it was tempting to speculate that the 339 showed a significant (FDR < 5%) overlap (Supplemental Table 9). We then tested whether the 353 root traits and the climate traits were significantly correlated when accounting for population 354 structure and correcting for multiple testing. We found that out of these 15 trait combinations 355 that had a significant overlap of associated SNPs, 9 displayed a significant population structure 356 corrected trait correlation (FDR < 5%) and another 2 were marginally significant (FDR < 10%) 357 (Supplemental Table 9). Overall this suggests that the root traits that we measured in our study 358 are involved in local adaptation. 359 The most significant overlap between SNPs was found between ABA induced lateral root 360 density and "Mean moisture index of coldest quarter" (FDR < 0.2%), consistent with the well-361 known role of ABA in response to water related responses. The trait correlation of this most 362 overlapping combination was negative (r= -0.189610096), suggesting that less lateral roots are 363 produced in an ABA dependent manner in accessions that are derived from areas in which the 364 soil contains high moisture in the coldest quarter of the year. As this trait correlation is 365 marginally significant when correcting for population structure (FDR < 7.5%), this suggests that 366 variation in genes that determine the impact of ABA on lateral root development is adaptive in 367 environments with contrasting soil moistures. Interestingly, none of the genes in proximity to the 368 19 overlapping SNPs is a bona fide ABA gene (Supplemental Table 10). 369 While the most significant overlap at the SNP level clearly suggested an important role of ABA 370 for local adaption for soil moisture, auxin and cytokinin treatments were also among the list of 371 root traits that were significantly overlapping with environmental parameters at the level of SNP 372 overlap and at least marginally significant correlated to the same climate trait (FDR < 8%). However, 7 out of the 11 combinations displaying significant overlaps at the level of SNPs and 374 traits did not associate with hormonal treatment, but with RSA traits under non-treated 375 conditions. Finally, only one of these overlapping traits is related to primary root growth. 376 However, as only 2 out of 10 traits that we had scored related specifically to the primary root, 377 this is still in the realm of expectation. Overall, these data clearly show that natural variation in 378 root traits is significantly associated with parameters that are highly relevant for local adaption in 379 Arabidopsis thaliana and that, in particular, ABA regulated lateral root traits seem to be relevant 380 for adaptation to soil moisture.  Table 3, Figure 3 and Figure 4). This is not simply a 453 peculiarity of our investigation; for instance, in a study dissecting 107 diverse phenotypes 454 (Atwell et al., 2010), Col-0 is located in the 1% lower tail of the distribution of all accessions for 455 14% of the traits (Supplemental Figure 11). While these findings don't affect the validity of the 456 fundamental mechanisms for hormone responses that have mainly been discovered in Col-0, they 457 suggest that, at a certain level of detail, studies will uncover relations and mechanisms that are 458 only found in specific genotypes within a species. We note that we don't suggest that Col-0 is 459 always an outlier. In fact, for a small number of traits under some conditions in this study 460 (Supplemental table 3), and for root traits in other studies, it has been described to be rather variation of responses to hormones. Using our data, we were able to search for common genetic 469 variants that associate with altered hormone responses using GWAS. Using the large number of 470 associations with ABA responses, we found that genes already implicated in the ABA pathway 471 were not enriched, and possibly even underrepresented, in our set of GWAS candidate genes 472 (Supplemental Figure 6). This supports the notion that genes that are targeted by natural 473 variation are frequently not the genes that are found by traditional means (such as forward 474 genetic mutant screens). 475 476 As more than 100 significant associations were detected, there were many good candidate genes 477 and we needed to focus on the most interesting ones. We chose to study two cases with known 478 candidate genes in regions associated with root traits upon hormone treatment. These turned out 479 to be a interesting cases with regard to candidate gene approaches. The first case was an 480 association for LRR that was detected under both control and CK treatment. We chose to 481 investigate this region because a gene, JAI3 that had been described to be involved JA signaling 482 (Chini et al., 2007) was located within 20KB of the associated SNP ( Figure 5B). As many 483 hormonal signaling pathways are tightly intertwined (Vanstraelen & Benkova, 2012), it was 484 tempting to assume JAI3 was also involved in CK signaling. However, based on LD patterns, 485 other genes were even better candidates, leading us to treat JAI3 as one of 5 candidate genes in 486 this region. Indeed, much like the LD pattern, gene expression analysis suggested that three other 487 genes, and not JAI3, were the prime candidates in this region, based on the correlation of their 488 expression pattern with the LRR trait ( Figure 5C, Supplemental Figure 8B). The root traits 489 observed in the available T-DNA insertion lines suggest that while JAI3 has an effect on root 490 length, this is not impacted by CK treatment and it is therefore not a likely candidate gene. In 491 contrast, at least two other genes, At3g17890 encoding an unknown protein and At3g17910 492 encoding a protein similar to a cytochrome c oxidase assembly factor in human, are within the 493 associated region and regulate root traits in response of CK ( Figure 5D and 5E). These genes 494 (and a third gene for which we could not obtain insertion lines) were the prime candidates based 495 on LD and expression correlation with the trait and amongst each other (Supplemental Figure 8B). Consistently, both available lines showed an altered expression level of multiple hormone 497 related genes. While cytokinin related genes were among those (most notably ARR12), in 498 particular PIN1 expression was altered strongly (Supplemental Figure 8D)

in both lines and 499
several other genes involved in the SHY2 network were altered in the second mutant line 500 (At3g17910). However, without further studies, it is difficult to clearly attribute the causality for 501 these expression changes to a perturbation of hormone signaling or altered root morphology in 502 these lines. Nevertheless, our data demonstrate that the yet uncharacterized genes At3g17890 and 503 At3g17910 are involved in cytokinin dependent root growth regulation. This shows that an 504 annotation based candidate gene approach can be misleading and that data-driven approaches 505 (e.g. LD, expression) seem to be more efficient in identifying genes with a function in the 506 GWAS trait of interest. 507 508

CRF3 is an integrator of multiple hormonal pathways 509
While further analysis of JAI3 revealed that this is almost certainly not the causal gene, our 510 second example showed that known genes should not be disfavored for consideration as 511  Figure 8B). An opposite effect on PIN1 and PIN7 had been previously observed 519 in root tips in this same line (Simaskova et al., 2015), suggesting that the effect of CRF3 strongly 520 depends on the tissue type (we used whole roots for our transcript measurements). Overall, our 521 results suggest the possible involvement of CRF3 as an integral part of 3 major hormone 522 pathways, auxin and cytokinin, which were already known (Simaskova et al., 2015), and abscisic 523 Lynch, 1995). Hormonal pathways are major factors in the regulation of root traits and root 530 system architecture (RSA). Here we showed that perturbation of hormonal pathways can shift 531 RSA through the phenotypic space but the extent of this regulation can strongly vary between 532 natural accessions of Arabidopsis thaliana (Fig. 4). As these accessions had been collected from 533 different parts of the world, ranging from northern Sweden to southern Spain (Supplemental 534 figure 1), we tested whether the variation in the root growth responses to hormones is associated 535 with the climate at the places from which the accessions were derived (Supplemental table 9). thaliana, ABA-responsive cis-regulatory elements (ABRE) showed significantly higher 545 frequencies of polymorphisms in genes that responded to drought in a genotype dependent 546 fashion than in drought responsive genes that responded similarly in all genotypes, suggesting 547 that natural genetic variation in ABA responsive transcriptional control contributes to genotype 548 specific drought responses (Lasky et al., 2014). 549 Despite 75% of all traits in our data were measured upon hormonal treatments, the majority of 550 significant overlaps at the level of SNPs that were supported by significant trait correlations did 551 not associate with hormonal treatments but with RSA traits under non-treated conditions. What 552 does this mean for the importance of hormonal pathways for local adaptation? It is difficult to 553 come up with a simple answer for this, as, also under control conditions, hormonal signaling 554 pathways are involved in shaping root traits. Certainly, in the case of the link between ABA and 555 Bio35, the data support an outstanding role for ABA for adapting to different soil moisture 556 regimes in the coldest quarter, as no control condition for the same root trait is significantly of wettest quarter [W/m 2 ]), as well as for P2 upon cytokinin treatment and Bio02 (Mean diurnal 559 temperature range (mean(period max-min)) [°C]). Overall, we can conclude that natural variation 560 in root traits is significantly associated with parameters that are highly relevant for local adaption 561 in Arabidopsis thaliana and that, in particular, ABA regulated lateral root traits seem to be 562 relevant for adaptation to soil moisture. Finally, we note that both SNP overlaps and trait 563 correlations clearly indicate that we are far from explaining a major proportion of adaptation to 564 climate variables, and that we expect numerous other processes and pathways to contribute to 565 local adaptation. We would hope that similar systematic studies in the future will lead to a more 566 complete picture that can then be tested in the wild and which might yield highly valuable 567 insights into which traits and which underlying pathways contribute to local adaptation. inhibition of root growth in the tested concentrations (Supplemental figure 9). Thus, we decided 592 to use the lowest concentration tested. Next, we tested the effect of placement in the upper vs. 593 lower row in the plate, because we intended to fit 4 accessions on one plate. Thus, it was 594 necessary to test if there is plate-position effect leading to differences in root growth between the 595 two rows. Our results show that there was no significant difference between the rows in the plate 596 in the tested accessions/concentrations (Supplemental Figure 10). To determine a list of ABA bona fide candidate genes, we used the GO SLIM annotation of the 621 TAIR10 release (ftp://ftp.arabidopsis.org/Ontologies/Gene_Ontology/ATH_GO_GOSLIM.txt; 622 creation date: 4/9/15). We then filtered for GO annotations that contained the key word "abscisic 623 acid" and were based on either IDA ("inferred by direct assay") or IMP ("inferred by mutant 624 phenotype"). The resulting list (Supplemental table 8

Computing the overlap of signals between GWASs 634
To determine whether the top associations between hormone GWAS and GWAS for 635 environmental variables was higher than expected by chance, we first obtained the 0.5% most 636 significantly associated SNPs in each GWAS. To avoid biases caused by outlier accessions, only 637 SNPs with a minor allele count equal to or greater than 6 were included in the list of SNPs from 638 which the top SNPs were derived. For each possible combination of root trait and climate 639 variable GWAS, we computed the overlap of these top SNPs. We then used the hypergeometric 640 distribution to compute a P-value for the observed overlap. This was done in python using the 641 scipy.stats.hypergeom function. The list of all computed P-values for all possible combinations 642 was corrected for multiple testing using the p.adjust function in R applying the "fdr" correction 643 (Benjamini & Hochberg, 1995