Adaptive growth homeostasis in response to drought in Iberian Arabidopsis

39 40 Plants respond to environmental fluctuations through plastic phenotypic shifts. Whether a plastic 41 response upon environmental variability is adaptive or not has been subject to debate. Using a set of 42 Iberian Arabidopsis accessions, we quantified an interplay between passive plastic reductions in 43 leaf areas that we found typical of accessions from productive environments and homeostatic leaf 44 areas responses to drought typified by accessions originating from unproductive environments. 45 Results from Genome-Wide Association Studies (GWAS) and Transcriptome Wide Association 46 Studies (TWAS) highlight the role of auxin-related processes and, in particular, the possible role of 47 the SMALL AUXIN UP RNA 26 (SAUR26) gene in the regulation of the observed plastic responses. 48 Homeostatic responses in leaf area potential following drought were typical of accessions with 49 lower leaf area potential under well-watered conditions. Transcripts that were negatively associated 50 with leaf area potential and positively associated with homeostatic and positive leaf area plasticity 51 following drought showed functional enrichment in ion transport processes. We hypothesized that 52 the contrasting plastic and homeostatic responses in leaf area potential were associated with 53 differential intrinsic water use efficiency (WUEi). We confirmed this relationship in a metanalysis 54 conducted using previously published δC measurements. Our results highlight the adaptive role of 55 homeostatic leaf area response to water depletion arising from increased WUEi. The concerted 56 utilization of Genome-Wide Association Studies (GWAS), Transcriptome Wide Association 57 Studies (TWAS), and expression Genome-Wide Association Studies (eGWAS) allows integration 58 of phenotype, genotype, and transcript abundance to identify both “plasticity genes” and 59 “homeostasis genes” associated with drought stress responses. 60 61 62

One-sentence summary: Information on phenotype, genotype, and transcript abundance is 12 integrated to identify both plasticity and homeostasis genes and processes associated with local 13 adaptation to drought stress in Arabidopsis accessions of the Iberian Peninsula. 14 15 16 Author  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  INTRODUCTION   63  64 Earth is undergoing a climate change scenario resulting in extreme weather events that are 65 increasingly frequent and variable (Easterling et al., 2000;Meehl and Tebaldi, 2004;Schär et al., 66 2004). These changing conditions are creating opportunities for some species to spread beyond their 67 native ranges (Chen et al., 2011) while posing a challenge to species that are currently adapted to 68 their local environment (Holt, 1990 as an appropriate response to environmental change during the lifetime of an individual plant. 72 However, simply because a given trait is plastic does not mean that this plasticity is adaptive 73 (Ghalambor et al., 2007). Active plasticity is anticipatory in response to predictable environmental 74 cues. On the other hand, passive plasticity is not anticipatory but a consequence, often proportional, 75 to a resource limitation. Active plastic phenotypic responses to environmental cues are thus more 76 likely to be adaptive than passive plasticity that unavoidably arises from unpredictable resource 77 limitation (Dorn et al., 2000;van Kleunen and Fischer, 2005;Forsman, 2015). However, the 78 division between passive and active plasticity is not straightforward. It not only depends on whether 79 the environmental factor is a reliable cue such as daylength or an unpredictable resource such as 80 water (precipitation) but also on the phenotype under consideration. For example, when plants are 81 subjected to drought, their growth is limited, a passive response. Still, they can also respond 82 actively by accelerating the onset of the reproductive stage (early flowering time). 83 The study of passive plastic responses to changes in resource availability can provide 84 insights into plant responses to climate change-related events that are unpredictable. Water 85 availability is indeed unpredictable, as the frequency and intensity of precipitation vary seasonally, 86 annually, and geographically (Loik et al., 2004). Drought is thus among the strongest factors 87 driving adaptation in plants (Siepielski et al., 2017; Exposito-Alonso et al., 2019). Climate 88 projections for the remainder of the 21 st century predict precipitation extremes to become more 89 severe (Donat et al., 2016) and for this reason, understanding the genetic basis and adaptation 90 potential of phenotypic traits and their plastic or homeostatic responses to drought is paramount. 91 Water availability ranges in a continuum from critically low to abundant and even 92 excessive. The distribution of phenotypic traits and their plasticity in a population is often 93 continuous (Galton, 1894) because they are influenced by variation in a large number of genes. The 94 effect of many of these variants is additive or synergistic and of small effect individually (Fisher,95 1919), and thus difficult to identify (Barton et   responses to drought. We uncover an adaptive homeostatic growth strategy in accessions from less 133 productive environments, dissect its genetic basis, and conclude that, just as has been argued 134 previously for plasticity, homeostasis is under genetic control and of evolutionary importance. 135

136
Natural variation in leaf area and plasticity in response to long-term drought 137 138 To quantify phenotypic variation in Iberian Arabidopsis accessions exposed to sustained drought, 139 we conducted a common garden experiment (Fig. 1A). We used image analysis to derive rosette 140 leaf area for 43 Iberian accessions grown under well-watered or drought conditions (Fig. 1B, C;  141  Table P1). As expected, we observed a reduction in leaf area resulting from drought in most 142 genotypes, as evident from the phenotypic reaction norms, while a small number of accessions 143 exhibited constant or increased leaf area under drought (Fig. 1A,D,E). Based on these data, we 144 calculated leaf area potentials ( Fig.1B,C,S,E) and plasticity indexes. We define leaf area potential 145 as the normalized maximum phenotypic value realized by any given accession under each condition 146 of the experiment (well-watered or drought). We used potential rather than actual values in order to 147 allow comparison of relative leaf area of each accession within the study population between both 148 watering regimes. 149 A leaf area plasticity index was calculated (Bouslama and Schapaugh, 1984) for each 150 accession based on the difference in leaf area potential between drought and well-watered 151 conditions, i.e., based on the reaction norms in Fig. 1F; Table P1. Genetic variants significantly 152 associated with this parameter were identified by GWAS ( Fig. 2A; Table G3). Gene ontology (GO) 153 analysis of functions associated with the genes harboring the variants identified by GWAS revealed 154 an enrichment in auxin-related functions ( Fig. 2A; Table GO3). Among the variants identified from 155 GWAS ( Fig. 2A; Table G3) and TWAS (Table T1,T2) as associated with leaf area plasticity, we  156   highlight two co-varying SNPs in the 5' and 3' UTR of SMALL AUXIN UP RNA 26 (SAUR26,  157 AT3G03850) ( Figure 2B). We also found from eGWAS that the transcript levels of SAUR26 are 158 regulated by this genetic variation (cis-regulation; Fig. 2C; Table G6). 159 We illustrate the haplotypic differences in leaf area plasticity for SAUR26 that we identify in 160 our study ( Figure 2B). Accessions harboring the minor SAUR26 haplotype display negative 161 plasticity in leaf area while the major haplotype is associated with homeostatic and positive 162 plasticity responses in leaf area in response to drought ( Fig. 2D; Wilcoxon P-value < 0.001). These 163 haplotypic differences were also associated with the transcript levels of SAUR26 ( Figure 2C), which 164 we confirmed both in the group of accessions included in our study (Wilcoxon P-value < 0.001) as 165 well as in the 189 Iberian accessions included within the 1001 Genomes Project ( Fig. 2E; Wilcoxon 166 P-value < 0.001). Both the genetic variation defined by this haplotype in SAUR26 and the resulting 167 transcript abundance are associated with the interplay between plasticity and homeostasis that we 168 define in our study (Fig. 2D,F). These results highlight SAUR26 as a candidate regulator of leaf area 169 plasticity in response to drought (Fig. 2G). 170 We conducted a GWA analysis of leaf area potential based on the genome-wide association 171 of individual SNPs for leaf area in response to drought (Fig. S1A, Table G1, G2). In our GWA 172 analysis on leaf area potential (Fig. S1A, Table G1, G2), no SNPs were found that were associated 173 with leaf area potential under both well-watered and drought conditions, i.e., no SNPs that would be 174 depicted as horizontal lines in Fig. S1A. The results from TWA analysis were strikingly different, 175 identifying transcripts associated with both leaf area plasticity (Fig. S1B) and leaf area homeostasis 176 (Fig. 3). A TWA analysis identified five transcripts that are positively correlated with leaf area 177 potential under both well-watered and drought conditions: SOFL2, CPN60B, AT5G02244, CKL8, 178 and IBH1 (Fig. 3A). This means that accessions with higher transcript abundance of these genes 179 display higher leaf area potential under both well-watered and drought conditions. We highlight 180 IBH1 as an example (Fig. 3B). Conversely, we identified four genes negatively correlated with leaf 181 area potential such that accessions with a higher transcript abundance display lower leaf area 182 potentials under both well-watered and drought conditions: RHF1A, AT5G60160, AT4G15260, and 183 AT5G16800 (Fig. 3A). We highlight RHF1A as an example (Fig. 3C). 184 185 Trade-off between growth, plasticity and intrinsic water use efficiency 186 187 We determined that the plastic leaf area responses we observed in our study were typical of 188 accessions from more productive habitats, while homeostatic responses were typical of accessions 189 from less productive habitats (Fig. S2). We next explored the relationship between leaf area 190 potential in the absence of drought and its plastic response when drought is imposed. In Figure 4A, 191 orange points depict accessions that display a reduction in leaf area in response to drought (negative 192 plasticity) but that in the absence of drought display a high leaf area potential. Blue datapoints 193 represent those accessions with homeostatic or positive growth responses to drought and low leaf 194 area potential under well-watered conditions. Red points are accessions with negative plastic 195 responses and low phenotypic potential. Finally, green datapoints are those accessions displaying 196 homeostatic or positive responses to drought that exhibit relatively high leaf areas in the absence of 197 drought. We observe a trade-off between growth, indicated by high leaf area potential in the 198 absence of drought, and leaf area homeostasis in response to drought. Accessions with higher leaf 199 area potential displayed higher plastic reductions in leaf area potential when subjected to drought. 200 On the contrary, accessions with lower leaf area potential under well-watered conditions display 201 homeostatic phenotypic responses to drought (Fig. 4A). 202 We performed a GWA analysis ( Fig. 4B; Table G1,G3) to determine the genetic basis of the 203 trade-off between growth and homeostasis in response to drought illustrated in Fig. 4A Finally, the function of AT5G16990 has not been previously characterized. Gene Ontology analysis 216 on the 17 genes that include the 50 SNPs with the strongest association with leaf area potential 217 reveal an enrichment in responses to osmotic and salt stress (Table GO3). This will become relevant 218 to the tradeoff between growth and plasticity as shown by our analysis of this process using TWAS, 219 described below. 220 Figure 4C depicts the relationship between transcript abundance and the observed trade-off 221 between growth and homeostasis in response to drought (Table T1, Ontology analysis reveals that the transcripts negatively associated with leaf area potential are 227 involved in ion transmembrane processes (Table GO3). This enrichment is also confirmed in the 30 228 transcripts depicted as blue datapoints in Fig. 4C (Table GO5). 229 A relationship between δ 13 C and intrinsic water use efficiency has been confirmed in natural study. The results revealed a trade-off between plasticity and WUE i (Fig. 4D). Accessions that 233 exhibit strong negative plasticity in the presence of drought have a lower intrinsic water use 234 efficiency in the absence of drought (more negative δ 13 C values). Those accessions that maintain 235 leaf area in response to drought display a higher water use efficiency in the absence of drought 236 (more positive δ 13 C values). 237 We found no genetic variants associated with the observed trade-off between intrinsic water 238 use efficiency and leaf area plasticity in response to drought ( Fig. 4E; Table T1,T2), i.e., no 239 datapoints in the upper right quadrant of Fig. 4E. However, TWA analysis (Table S2) provided  240 insight into this trade-off. Red datapoints in Fig. 5F are transcripts with higher greater abundance in 241 those plastic and less water-use efficient accessions depicted in red in Fig. 4D. Conversely, green 242 points in Fig. 4F represent transcripts with greater abundance in the homeostatic and water-use 243 efficient accessions depicted in green in Fig. 4D. GO analysis revealed an enrichment in snRNA 244 and rRNA pseudouridine synthesis ( We determined that in our study, accessions that displayed positive growth or homeostasis in leaf 252 area in response to drought ( Fig 4A) also had higher WUE i in the absence of drought (Fig. 4D). We 253 then conducted a meta-analysis on the 95 Iberian accessions with δ 13 C information from a previous 254 study focused on the global Arabidopsis population (Dittberner et al., 2018). GWA analysis of these 255 accessions rendered a distinctive peak in a single gene, identifying a haplotype consisting of eight 256 co-varying SNPs in ClpX3 (AT1G33360; Fig. 6A; Table G5), an ATP-dependent Clp protease. 257 More informative was the TWA analysis of water use efficiency determined from (δ 13 C) values in 258 the 88 Iberian accessions with both published transcript abundance and δ 13 C information: we 259 uncovered a distinctive enrichment in photosynthesis-related genes ( There is an adaptive trade-off between growth and drought resistance in plants (Zhang et al., 2020). 272 Plants adapted to more productive habitats with higher and more constant precipitation regimes 273 favor decreased water use efficiency and increased carbon capture, producing increased leaf areas 274 to compete for aboveground resources with neighbors (Tilman, 1988). Conversely, plants adapted 275 to unproductive environments typically fix less carbon, have higher water use efficiency, and 276 display low growth rates, even when grown in rich environments (Parsons, 1968;Chapin, 1991), 277 and are typical stress-tolerant (Billings and Mooney, 1968;Grime, 1977). In this study, we observe 278 that accessions with higher leaf area potentials in the absence of drought display strong and 279 negative plastic responses to drought (Fig. 1). On the other hand, accessions with lower leaf area 280 potential in the absence of drought exhibit leaf area homeostasis or even positive leaf area potential 281 plastic responses in response to drought (Fig. 1,4A). 282 Natural variation in the degree of plasticity among different ecotypes (accessions in the case 283 of Arabidopsis) that are subjected to similar depletions in water availability can provide valuable 284 information on the adaptive/neutral/maladaptive role of this plastic response. The variation in leaf 285 area in Arabidopsis accessions in response to drought in this study can be grossly divided into 286 passive and negative phenotypic plasticity in accessions adapted to richer environments (high Net 287 Primary Productivity in spring, NPP spring), vs. adaptive phenotypic homeostasis responses 288 observed in accessions originating from less productive areas of the Iberian Peninsula (low NPP 289 spring) (Fig. S2). We sought to identify the genetic basis of these contrasting ecological strategies. 290 Identifying causal genetic variants related to an abiotic stress response such as drought is 291 challenging given the polygenic nature of these responses (Exposito-Alonso et al., 2018b). 292 Moreover, here, our main interest was not to identify natural variation associated with well-watered 293 or drought phenotypes, but rather natural variation associated with the degree of plasticity evoked 294 in response to drought, i.e., natural variation associated with reaction norms, rather than with 295 phenotypes in one environment or the other (well-watered or drought). Reaction norms in leaf area 296 in response to drought (Fig. 1E, F) illustrate the plasticity index that we calculated based on the 297 difference in leaf area potential between drought and well-watered conditions (Table P1). In 298 exploring the genetic basis of plasticity, we identified the 31 genes corresponding to the top 50 299 candidate SNPs with the highest strength of association with the leaf area plasticity index that we 300 calculated ( Figure 2B; Table G3). Of these 31 genes, 20 have been previously characterized, and 301 many of these are related to auxin regulation and leaf expansion. This enrichment in auxin-related 302 processes was confirmed by GO analysis (Table GO3). The variant with the strongest association 303 with the plastic response to drought was the auxin efflux carrier family protein PIN8 (AT5G15100). 304 The resulting from TWAS analysis on leaf area plasticity index for drought revealed an enrichment of 318 glucosinolate biosynthesis genes (Table S6). 319 We found that among the auxin-related variants associated with leaf area plasticity ( Fig 2019). The two associated variants in SAUR26 were found in the 5' and 3' UTR, respectively ( Fig.  323 2B), regions of the ORF that typically regulate gene expression. Indeed, we found the transcript 324 levels of SAUR26 to be associated with the genetic variation in SAUR26 itself (cis-regulation; Fig.  325 2C), which we confirmed both in the group of accessions included in this study and in the 189 326 Iberian accessions included within the 1001 Genomes Project (Fig. 4E). The minor haplotype is 327 associated with negative plasticity in leaf area in response to drought. Conversely, the major 328 haplotype, associated with higher transcript abundance of SAUR26, is associated with leaf area 329 homeostasis and positive plastic responses (Fig. 2D). 330 Relevant to the non-discrete distribution of the plastic responses in the accessions that we 331 studied, variation in transcript abundance of SAUR26 explained 17% of the observed plastic 332 response in leaf area to drought in our study ( Figure 2F, G). The importance of SAUR26 is 333 supported by previous work that revealed natural genetic variation in SAUR26 associated with 334 growth thermo responsiveness of the rosette compactness index (Wang et al., 2019) and 335 cis elements in the 5' upstream region or the 3'UTR region of SAUR26 that were responsible for 336 this response (Wang et al. 2021). This conclusion is also supported by our meta-analysis of the 337 larger Iberian population with transcriptome data, which also identified SAUR26 transcript 338 abundance as positively correlated with higher water use efficiency (Fig. 6B,C) along with other 339 related auxin-related transcripts such as AUXIN RESISTANT 3 (AXR3, also known as IAA17, 340 AT1G04250) (Leyser et al., 1996)  Genomes collection. We know that Iberian Arabidopsis populations locally adapted to drought 349 conditions in this region out-perform accessions from wetter areas of Northern Europe when both 350 are exposed to drought (Exposito-Alonso et al., 2019). Following this, given that the major 351 haplotype in our study is associated with homeostatic responses to drought (Fig. 2D), and that its 352 frequency became less frequent outside the drought-prone Iberian region, we can speculate that this 353 SAUR26 haplotypic variant is under selection in the Iberian population conferring an adaptive 354 advantage to water-limited conditions in the region. 355 The study of the strength of association of genetic variation with leaf area plasticity provides 356 us with candidates that reveal a binary answer: one allele is associated with plastic reductions in 357 leaf area while the alternative allele will display an association with either homeostatic or positive 358 plasticity (increase in leaf area potential) under drought. We did not identify any SNPs associated 359 with leaf area potential under both well-watered and drought conditions. Therefore, we could not 360 define the genetic basis of growth homeostasis under drought using GWAS (Fig. S1). On the other 361 hand, reaction norms based on the strength of association of transcript abundance with leaf area 362 potential identify five transcripts that are positively correlated with leaf area potential under both 363 well-watered and drought conditions (Fig. 3A). This means that accessions with higher transcript 364 abundance for those genes display higher leaf area potential under both well-watered and drought 365 conditions. Conversely, we identified four transcripts negatively correlated with leaf area potential 366 such that accessions with a higher transcript abundance of those genes display lower leaf area 367 potentials under both well-watered and drought conditions. Therefore, these nine genes are strong 368 functional candidates for homeostatic regulation of leaf area potential. 369 Most of the candidate genes we obtained from this analysis have not been characterized, and 370 for those that have, their function is not necessarily known to be involved in growth processes 371 leading to differential leaf area. Among the homeostatic transcripts positively correlated with leaf 372 area potential, ILI1 binding bHLH 1 (IBH1, AT2G43060) ( Arabidopsis IBH1 in transgenic tobacco plants (Nicotiana tabacum) induces a dwarf phenotype 377 (Nagatoshi et al., 2016). In our study, we observe a direct relationship between the natural variation 378 in the transcript abundance of IBH1 and leaf area potential, which seems counterintuitive based on 379 these latter reports, perhaps reflecting an interaction between genotypic background and IBH1. 380 SOFL2 (Fig. 3C), another candidate positively correlated with leaf area potential and homeostasis. 381 has been previously reported to confer an increase in the levels of specific cytokinins, which 382 positively regulate plant growth (Zhang et al., 2009 To increase leaf area potential and sacrifice WUE i or increase WUE i and 392 sacrifice leaf area potential? 393 We found that accessions with higher leaf area potential displayed greater plastic reductions in leaf 395 area potential (negative leaf area plasticity index) when subjected to drought (Fig. 4A). This 396 response is presumably neither adaptive nor maladaptive but a consequence or passive response to a 397 limited resource rather than an active response. Consistent with this, accessions with plastic 398 responses are typical of productive habitats, while the homeostatic responses that we observed in 399 our study typify accessions collected from less productive habitats (Fig. S2). 400 Gene Ontology analysis on the 17 genes that include the 50 SNPs with the strongest 401 association with leaf area potential (green and yellow datapoints in Fig. 4B) reveals an enrichment 402 in responses to osmotic and salt stress (Table GO1). Gene Ontology analysis also reveals that the 403 transcripts negatively associated with leaf area potential are involved in ion transmembrane 404 processes (Table GO2). This enrichment is also confirmed in the 30 transcripts depicted as blue 405 datapoints in Fig. 4C (Table GO5). We hypothesized that these processes are related to the greater 406 water use efficiency of accessions with reduced growth rates typical of unproductive environments. 407 Indeed, plant species with low relative growth rates are known to exhibit high intrinsic water use 408 efficiency (WUE i ) (Angert et al., 2009;Kimball et al., 2012). WUE i is defined as the ratio between 409 instantaneous carbon gain and transpirational water loss, typically calculated per unit of leaf area 410 (Condon et al., 2002). In support of our hypothesis, the homeostatic response in leaf area potential 411 in response to prolonged drought that we observe in our study correlates positively with a high 412 predicted water use efficiency recorded by Dittberner et al. (2018) for the same accessions (Fig.  413 4D). While we regard the observed plasticity in response to drought as a passive response to 414 resource availability, we propose that the observed leaf area homeostasis is an adaptive response. 415 These area plasticity index and WUE i (Fig. 4E). However, when we compare the candidates obtained from 419 TWAS analysis, we identify transcripts associated with both low WUE i and a negative plastic 420 response in leaf area when subjected to sustained drought. Conversely, we identify a set of 421 transcripts that is significantly associated with both leaf area homeostasis (more positive leaf area 422 plasticity index) and increased WUE i (Fig. 4F). To exemplify these opposed responses to drought, 423 given the enrichment in ion transport genes that we found in the transcripts regulating the tradeoff 424 between leaf area plasticity and potential ( Fig. 4C; Table GO5), we highlight two of the genes 425 identified by TWAS, with opposite correlations between transcript abundance and relative leaf area 426 plasticity (Fig. 6A,B vs. Fig. 6C opening, which results in decreased WUE i . This is supported by the negative correlation between 441 leaf plasticity and AHA4 transcript abundance that we observe, wherein accessions with higher 442 transcript abundance of AHA4 display strong negative plastic responses in leaf area potential in 443 response to drought (Fig. 5C). Accordingly, the increased transcript abundance of AHA4 is also 444 associated with reduced WUE i (Fig. 5D). 445 We have described two opposed processes in response to sustained drought. On the one 446 hand, accessions displaying higher leaf area potential in the absence of drought display a higher 447 plastic reduction in leaf area potential following drought than accessions with lower leaf area 448 potential in the absence of drought. Indeed, species with high relative growth rates are known to 449 exhibit low WUE i ( Table GO8). During photosynthesis, plants discriminate against 13 CO 2 in the 473 reaction center of Ribulose 1, 5-biphosphate carboxylase/oxidase (RuBisCO). As stomata close to minimize water loss in response to drought, the discrimination against 13 C is less intense because 475 the internal CO 2 concentration (C i ) becomes unavoidably reduced, reflecting a trade-off between 476 carbon gain and water loss. As a result of their reduced C i , water-efficient plants experience lower 477 CO 2 concentrations at the sites of carboxylation (C c ) and have more positive δ 13 C values. Among 478 the candidates from our TWAS analysis that are positively related to water use efficiency, we found 479 RuBisCO small subunit 2B (AT5G38420; RBCS2B; Figs. 6B, D). An upregulation in RuBisCO may 480 compensate for suboptimal CO 2 concentrations at the sites of carboxylation that plants from 481 drought-prone environments typically experience. Among these candidates, we also find 482 Sedoheptulose . We can speculate that the large proportion of uncharacterized genes negatively 502 associated with WUE i may be the result of an increased lethality in the loss of function mutations in 503 these genes. (Table GO9). We can further speculate that this enrichment reflects the transport of 504 carbohydrates produced during photosynthesis to promote plant growth. This would be concordant 505 with the tradeoff between growth and WUE i that we observe in our study. 506 507 Limitations 508 509 Our study has revealed genetic variants associated with both plasticity and homeostasis in leaf area 510 following drought. One limitation of our study is that while the number of accessions that we 511 utilized would have been considered ample just a few years ago, with the more recent advent of 512 GWA studies in Arabidopsis such studies typically include more accessions. If the number of 513 accessions utilized is low, low-frequency variants may not be present in the subset of accessions 514 studied. This is less of a concern when we focus on variants that are undergoing selection and thus 515 are at higher frequency in the population. Accordingly, in this study, we did not consider variants in 516 frequencies lower than 10% as we were interested in adaptive genetic variation. In addition, we 517 performed TWAS, which is less susceptible to this concern because it evaluates continuous 518 variation in transcript abundance rather than the presence/absence of SNPs within the study 519 population. Indeed, RNA-seq experiments often are conducted on very few genotypes, frequently 520 comparing transcript abundance between just two genotypes: wild-type and a specific mutant of 521 interest. Another limitation typical of GWAS is the fact that synthetic associations, non-causative 522 markers that are within linkage disequilibrium with causative makers, are difficult to discern from 523 true causative markers ( The validation of candidates from both GWA and TWA studies also suffers from the 539 challenge of restricted genetic backgrounds. Most often, the characterization of knock-out mutants 540 of candidate genes is performed in a single genetic background, typically Col-0. Here, we instead 541 validated our results by meta-analysis of the larger Iberian population. TWAS in this larger sample 542 identified 31 out of the 45 candidate genes that we previously identified in the subset of 20 Iberian 543 accessions in our study for which published δ 13 C data were available (Table T1). Moreover, there 544 are significant correlations between the results from our sample and the larger sample from both 545 GWAS and TWAS (Fig. S4). 546 547

549
In our study, we find a gradient ranging from higher plasticity in response to resource limitation 550 (drought) for generalist accessions originating from productive environments to increased 551 specialization and decreased phenotypic plasticity in accessions originating from more marginal 552 environments. The "passive" plasticity observed here may nevertheless be important over 553 evolutionary time as it gradually allow adaptiveness in unproductive environments as the induced 554 change in resource availability moves the mean phenotype closer to the environment-specific 555 optimum favored by selection (Waddington, 1953;Richards et al., 2006). 556 The adaptive role of plasticity has been at the center of a long-standing debate. Much work 557 and discussion have been focused on the principle that if indeed plasticity is adaptive, it may be 558 under genetic control (Bradshaw, 1965;Via et al., 1995); however, few relevant genes have been 559 identified. Here we propose that not only plasticity but also homeostasis is under genetic control 560 and important from an evolutionary standpoint. Moreover, we dissect the genetic basis of both 561 plasticity and homeostasis in leaf area, identifying specific genes and regulatory variants implicated 562 in each process. Our integrated application of GWAS, TWAS, and eGWAS identifies both 563 candidate "plasticity variants" and candidate "homeostasis variants" that can be further studied. control and drought treatment were well-watered for the first 28 days, by maintaining the soil 589 relative water content over 75%, as measured by monitoring pot weight daily and adding water 590 when necessary. After 28 days, plants within the drought treatment were watered to soil saturation 591 once every 10 days, while the well-watered treatment was maintained near field capacity 592 throughout the duration of the experiment. The experiment continued for the entirety of the life 593 cycle of the plants. Plants were randomized within the greenhouse on a weekly basis throughout the 594 experiment. At the end of the experiment 170 days later, nine accessions within the drought 595 treatment that encountered mortality rates higher than 40% (two of the five replicates; Table S1)  596 were excluded from downstream analysis 597 Photographs of rosettes were taken at the time of flowering and subsequently analyzed using 598 ImageJ (http://imagej.nih.gov/ij/) to determine projected leaf area. For the associations between the 599 phenotypic variables we obtained in this study, and the local environment, we inspected the data 600 using Cook's distance (Cook, 1984;Altman and Krzywinski, 2016) to evaluate the influence of 601 each datapoint on the fit. We considered values greater than 4/n to be overly influential. Based on 602 this, we decided to retain outliers in our analysis except for the phenotypic values we obtained for 603 the accession IP-Coc1 (accession id = 9535), which were consistently unrealistically influential 604 across different phenotypes based on their extreme Cook's distance. We include these removed 605 values in Table S1. 606 607 reflects the leaf area that any given accession displays in the absence of drought stress, divided by 614 the average leaf area of the population in the absence of drought stress. A high leaf area potential 615 reflects a higher vegetative growth of any given accession relative to the set of accessions utilized 616 in this study. We used potential rather than actual values in order to compare the relative leaf area 617 of each accession within the study population between both watering regimes. 618

Determination of phenotypic potential and plasticity
Leaf area plasticity indexes were calculated based on an established definition (Bouslama 619 and Schapaugh, 1984). Leaf area plasticity indexes to drought (LAPI d ) were calculated as: LAPI d = 620 LAP d -LAP ww ; where LAP d = leaf area potential under drought and LAP ww = leaf area potential 621 under well-watered conditions. Positive indexes reflect an increment in leaf area potential value for 622 any given accession after drought treatment, not necessarily an increment in leaf area. For instance, 623 a positive leaf area plasticity index does not necessarily imply leaf areas were higher under drought, 624 but that under drought, a given accession displays a higher leaf area potential relative to the rest of 625 accessions utilized in this study. 626 Reaction norms based on leaf area potential were employed to reflect plastic changes in 627 phenotypic potential between the two treatments (well-watered vs. drought). The difference in leaf 628 area potential between drought and well-watered conditions for the reaction norm of each accession 629 is its leaf area plasticity index. 630 631

633
We used a GWAS approach to identify polymorphisms associated with variation in phenotypic 634 potential and plasticity. The online tool GWAPP (http://gwas.gmi.oeaw.ac.at/) was employed using 635 a linear regression model (Seren et al., 2012). While the linear regression model does not address 636 confounding effects of population stratification, family structure and cryptic relatedness (Price et  637 al., 2010), it was employed here because those issues were minimized by our choice of a regional 638 ( We filtered out low frequency variants (MAF < 10%), to focus on variants less likely to be 640 undergoing negative selection or to be spurious SNPs from sequencing errors. We obtained gene 641 model annotation from the TAIR11 genome release. We defined the 'consensus' transcript variant 642 as the variant producing the longest transcript. We predicted the effect of every SNP using SnpEff 643 (Cingolani et al., 2012) (Table G1-G6). 644 Using the results from GWAS obtained from the association with leaf area potential, we 645 obtained reaction norms that reflect plastic and homeostatic patterns of genome-wide associations 646 of these traits in response to drought. For instance, genetic association plasticity indexes to 647 drought (G d ) were calculated as: GPI d = G d -G ww; where G d = Genome-wide association score (-log  648 p-value) with any given phenotype under drought, and G ww = Genome-wide association potential 649 under well-watered conditions. The association between a phenotype and a SNP was considered 650 more plastic as the genome-wide association score difference in the correlation coefficient between 651 the two conditions reflected in the reaction norm (well-watered vs. drought) became greater. On the 652 contrary, the association between a phenotype and any given SNP was considered to be homeostatic 653 if found significantly associated to that trait under both conditions (well-watered and drought physical distance between the genetic variant. cis-eSNPs were defined as those located within the 690 ORF or the promoter region, which was defined as 5 kb upstream of the most distal transcription 691 start site of the gene. Trans-eSNPs were defined as those located within the ORF/promoter region 692 of a gene different from the one that is transcriptionally regulated by that genetic variant. 693

695
To explore the local environment of this population, we retrieved 204 environmental parameters for 696 these 188 Iberian accessions from AraCLIM (https://github.com/CLIMtools/AraCLIM) ( Analyses 709 All statistical analyses and figures were produced using R (Team, 2013 (Ge et al., 2020) with the following settings: the search species was "Arabidopsis thaliana," the P-715 value cutoff (FDR) was 0.05, and the number of most significant terms to show was 10. To 716 replicate these results, the lists of genes that were input into these analyses based on significance 717 thresholds can be selected from Table G1-G6, T1,T2.  718  719 Supplemental data 720 721 The following supplemental materials are available 722 723 Figure S1. Genotypic and transcriptomic reaction norms based the strength of association of 724 individual SNPs and transcripts highlights the variation of these associations depending on water 725 availability. 726 Figure S2. Accessions from more productive areas display stronger negative plastic reductions in 727 leaf area relative to accessions from less productive environments 728 Figure S3. Illustration of the photosynthesis-related transcript variants that regulate water use 729 efficiency in Iberian Arabidopsis accession. 730 Figure S4.  Table T1. Results from TWAS (r s ). 760 Table T2. Results from TWAS (P-value). 761 Table GO1. GO analysis on the 17 genes corresponding to the top 50 candidate SNPs resulting 762 from GWAS analysis on leaf area potential under well-watered conditions. 763 Table GO2. GO analysis on negatively correlated candidates resulting from TWAS analysis on leaf 764 area potential under well-watered conditions. 765 Table GO3. GO analysis on the 30 genes corresponding to the top 50 candidate SNPs resulting 766 from GWAS analysis on leaf area plasticity index for drought. 767 Table GO4. GO analysis on the positively correlated candidates resulting from TWAS analysis on 768 leaf area plasticity index for drought. 769 Table GO5. GO analysis on the negatively correlated candidates resulting from TWAS analysis on 770 leaf area plasticity index for drought. 771 Table GO5. GO analysis on the transcripts with abundances that overlap negatively with leaf area 772 potential under well-watered conditions, and positively with leaf area plasticity index for drought 773 (blue datapoints in Fig. 4C). 774 Table GO6. GO analysis on the transcripts with abundances that overlap negatively with leaf area 775 potential under well-watered conditions, and positively with leaf area plasticity index for drought 776 (orange datapoints in Fig. 4C). 777 Table GO7. GO analysis on the transcripts with abundances that overlap negatively with leaf area 778 plasticity index for drought and WUEi (red datapoints in Fig. 4F). 779 Table GO8. GO analysis on the positively correlated candidates resulting from TWA analysis of 780 water use efficiency determined from (δ13C) values in 88 Iberian Arabidopsis accessions (Fig. 6A). 781 Table GO9. GO analysis on the negatively correlated candidates resulting from TWA analysis of 782 water use efficiency determined from (δ13C) values in 88 Iberian Arabidopsis accessions (Fig. 6A). 783 784  represent probability densities of leaf area potential among the accessions included in this study. 798 Datapoints depict the leaf area potential of individual accessions, and the reaction norms represent 799 the plastic response in the leaf area of leaf area potential of each accession in response to drought. 800 represent the different chromosomes. The higher the score, the lower the P-value, and the stronger 819 the association between genetic variation and leaf area plasticity or SAUR26 transcript abundance. 820 D. The minor haplotype of SAUR26 is associated with negative plastic responses of leaf area 821 potential in response to drought. Conversely, the major haplotype is associated with homeostatic or 822 positive responses in leaf area potential in response to drought. Violin plots showing significantly 823 different probability densities of leaf area plasticity in response to drought for the major and minor 824 allele haplotypes in SAUR26. E. The haplotype in SAUR26 associated with leaf area 825 plasticity/homeostasis in response to drought is a regulatory haplotype. Violin plots showing 826 significantly different probability densities of SAUR26 transcript abundance for the major and 827 minor allele haplotypes in SAUR26 for both our sample and the entire Iberian population.  (Table S1). A. Detail of the common garden experiment, illustrating the randomization of well-watered and drought treatments, as reflected in soil coloration. B. Leaf areas for accessions grown under optimal conditions. C. Leaf areas for accessions under drought conditions. Accessions are arranged alphabetically. Bars represent means ±SE. D. Images illustrating the plastic and homeostatic responses to drought found among the accessions included in this study. E. Phenotypic reaction norms illustrating an overall reduction of leaf area in response to drought (t-test; P-value < 0.001). Violin distributions represent probability densities of leaf area (cm 2 ) among the accessions included in this study. Datapoints depict the mean phenotypic values of individual accessions, and the reaction norms represent the plastic response in the leaf area of each accession in response to drought. F Phenotypic reaction norms illustrating plastic and homeostatic responses in leaf area potential (t-test; P-value > 0.05). Violin distributions represent probability densities of leaf area potential among the accessions included in this study. Datapoints depict the leaf area potential of individual accessions, and the reaction norms represent the plastic response in the leaf area of leaf area potential of each accession in response to drought. Orange block depicts the single exon, green block the 5' UTR, dark green the 3' UTR. Labels above indicate the position of the two SNPs that co-vary forming a haplotype in our sample. The illustration includes information on the major/minor allele corresponding to each of the SNPs, its frequency in our sample, in 189 accessions from the Iberian within the 1001 Genome Project, and in the whole 1001 Genome Project sample dataset. C. Manhattan plot representing the genome-wide strength of association between the genetic variation present within the 665 accessions from the transcript abundance reference panel used in this study (Kawakatsu et al., 2016) that overlap with the collection of accessions included in the 1001 Genomes project, and the variation in the transcript abundance of SAUR26 (Kawakatsu et al., 2016). The peak in chromosome 3 represents the cis-genetic variants in SAUR26 affecting its transcript abundance. For B and C, the score (yaxis) consists of the negative logarithm of the P-value and provides association strength. Each data point represents an SNP with a color corresponding to the chromosome where the genetic variants are located. The x-axis and the different colors shown in the SNPs depicted in the Manhattan plots represent the different chromosomes. The higher the score, the lower the P-value, and the stronger the association between genetic variation and leaf area plasticity or SAUR26 transcript abundance. D. The minor haplotype of SAUR26 is associated with negative plastic responses of leaf area potential in response to drought. Conversely, the major haplotype is associated with homeostatic or positive responses in leaf area potential in response to drought. Violin plots showing significantly different probability densities of leaf area plasticity in response to drought for the major and minor allele haplotypes in SAUR26. E. The haplotype in SAUR26 associated with leaf area plasticity/homeostasis in response to drought is a regulatory haplotype. Violin plots showing significantly different probability densities of SAUR26 transcript abundance for the major and minor allele haplotypes in SAUR26 for both our sample and the entire Iberian population. F. Scatter plot illustrating the positive relationship between leaf area plasticity index (y-axis) and SAUR26 transcript abundance (x-axis). Grey shadow represents the 95% CI of the fit regression line (dark grey). G. Cartoon illustrating the regulatory role of genetic and transcript variation in SAUR26 and the response of leaf area to drought.