GAI MoRFs Regulate Cleft and Channel Binding Pathways for Gibberellin in GID1A

Abstract The hormone gibberellin (GA) promotes arabidopsis growth by enhancing binding between GA Insensitive DELLA transcriptional repressors and GA Insensitive Dwarf 1 (GID1) receptors to regulate DELLA degradation. The binding mechanism for GA was elucidated by employing a computational study of dissociations of the N-terminus of the DELLA family member GAI (GA Insensitive transcriptional repressor) from the GID1A receptor in the presence and absence of bound GA, and of GA from GID1A in the presence and absence of GAI. The tRAMD method was employed to deduce egression pathways for a diverse set of GA molecules (GA(x)). Two pathways in the form of a newly identified cleft and a previously identified channel are prevalent. The cleft pathway is open in the absence of GAI. Upon GAI binding, the cleft route is blocked, resulting in a slower process for GA(x) to exit and enter the binding pocket through the channel. Several binding pocket residues are identified as gate-keepers to the channel. Molecular recognition features (MoRFs) found in the disordered signaling protein GAI affect GA(x) binding and GID1A dynamics. A three-step synergistic binding cycle is proposed where GAI MoRFs regulate the process. Rapid binding takes place through the cleft where little to no distinctions are made between major and less active forms of GA(x). After GAI is bound to the GA(x) · GID1A complex, the channel supports a rectification process that increases the retention of major active forms of GA within the binding pocket. Both the cleft and channel contact residues to GA(x) are markedly conserved in a GID1 phylogeny, suggesting this binding process in the GID1 · DELLA GA-receptor complex represents a general paradigm for GA binding. Non-specific GA binding assists binding of GAI, which then helps to select the major active forms of the hormone and induce a downstream signalling cascade in response to bioactive GA. Non-expert Summary Statement Gibberellins are plant hormones essential for growth and development. The DELLA proteins are a disordered family of repressors that transcriptionally repress GA responsive genes. Degradation of DELLA proteins in response to GA results in GA-responsive genes being upregulated. Binding of GA to the GA-Insensitive Dwarf 1 receptor (GID1) facilitates binding of DELLA to the GA · GID1 complex. Through computational modelling and phylogenetic analyses, we identified a new GA binding cleft that is blocked by DELLA binding and a three-step mechanism for the GA · DELLA · GID1 complex that also involves the known GA binding channel. We propose a dual (cleft/channel) pathway that allows access to the binding pocket as a paradigm for selection of specific GA forms among a mixture of major active and inactive forms. The cleft is less selective, but preference for active GA in the binding pocket of GID1A is amplified by expunging inactive GA forms, followed by recruiting active forms through the more selective channel. This mechanism allows plants to sense concentration changes of GA with high specificity to enable certain GA variants to trigger specific signalling events. These novel insights into the receptor mechanism in part may explain the large number of different GA forms that exist in nature. Graphical Abstract

147 the DELLA partners, are identified with respect to these pathways. Conceptualized in Figure 1, a 148 complete kinetic cycle is proposed that provides a mechanism for preferential binding of major 149 bioactive GA in the receptor complex. This model provides significant new insights that explain 150 how specific GA regulation can occur within the presence of a large number of GA variants that 151 exist in nature. These insights will also inform strategies to manipulate the protein complex to 159 To cover the spectrum of activities that may be inherent in this binary classification scheme, we 160 selected a diverse set of GA (x) that includes GA variants widely referred to as bioactive {GA (a) } 161 and as inactive {GA (i) }. Additional types of GA are also considered to increase the diversity in 162 molecular moieties. Known GA (a) GA 1 , GA 4 , GA 3 , GA 7 , and recently reported GA 12-16ox (33) 163 were chosen to assess active GA containing systems. We also analysed a set of putative inactive 164 forms (GA (i) ) that represent variants with no evidence of biological activity. Most of these were 165 chosen as a progressive step from GA 12 , the main precursor to all GA, to the various GA (a) 166 examined. In particular, GA 4MeO , GA 4-16/17ox , and GA 34 are all oxidation products of GA (a) and 167 are of interest to understand the signalling control in this system. 9 168 GA binding pocket dynamics during simulation 169 All production runs for each gDG system retained the bound GA (x) ligand within the binding 170 pocket. Counts of the recurrent hydrogen bonds (H-bonds) and hydrophobic contacts between 171 GA variants and the GID1A residues were monitored. Table 1 summarizes the binding pocking 172 interactions to GA (x) that are largely in common across all the GA variants investigated here.
173 Highlighting these pocket residues outlines the proximity to key parts of GA 3 . TYR247 (see 174 methods for residue numbering) interacts with GA (x) C3-OH, SER116 and GLY115 wrap around 175 the C7 carboxylate. TYR31 has interactions both with OH groups in those GA (x) that have this 176 moiety around C13 and interactions with all GA vicinal C16-C17 bonds, appearing to be aided 177 by pi-pi interactions. This interaction aids in aligning GA inside the GID1A pocket and is a 178 shared moiety across most GA. 179 In the GA (x) • GID1 (gG) systems a subset of common interactions was found. With no additional 180 interactions identified, the gG system maintains the same set of interactions except for 181 interactions with TYR31 and TYR247. The loss of the TYR31 and TYR247 interactions enable 182 there to be a less compact environment among the surrounding residues in gG systems compared 183 to gDG systems. Interestingly, the contact residues shown here are highly conserved across the 184 phylogeny of GID1 across diverse plant families (Tables S3-5)

220
GAx tRAMD analysis: Channel and cleft pathways 221 After 20 ns production runs (computational simulations) were completed, this information was 222 used in tRAMD to identify the GA (x) egression pathways during the dissociation process. We 223 observed that there are essentially two egression pathways, as shown in Figure 3. The longer 224 pathway uncovers the previously reported channel pathway by Hao et al. (28 245 For each GA variant, the set of GID1A residues that it interacts with when leaving through the 246 channel are summarized in Table 2. In most systems the previously identified hydrogen bonding 247 residues from Hao et al. are recurrent. It is notable that most of the hydrophobic contact residues 248 are simply adjacent to these residues, but are highly conserved across the phylogeny of GID1. 249 The only recurrent residue with large discrepancy between the GA (a) and GA (i) subsets is ASN32.
250 Unlike the pocket residue interactions, the GA (i) and GA (a) interacted similarly with channel 251 residues along the egressions. In simulations without the GAI N-terminus (gG systems), GID1A 252 was found to be flexible enough to allow egression from this same pathway without a large 253 portion of the interactions reported in Hao et al., Table S1. Residue interactions we identified in 254 the channel pathway in both the gG and the gDG systems were THR240, ASN32, and MET220.
267 A second pathway is also uncovered, as shown in 285 When the DELLA protein is bound to GID1A, the cleft opening is blocked. Nevertheless, as 286 Figure 3 indicates, GA (x) can egress from the cleft pocket opening even in the presence of the 287 DELLA protein, albeit in a heavily perturbed pathway. The binding pocket residence times for 288 the set of GA variants considered here for the gDG system are summarized in Table 3. On 289 average the egression pathways through the channel is of higher propensity than the cleft 290 pathway. However, among the 12 GA variants considered here, GA 9 , GA 4 and GA 34 are 291 exceptions that exhibit a higher propensity to leave the binding pocket through the cleft pathway.    311 Even within the known bioactive set of GA ligands there is variation in average residence times 312 particularly between GA 7 and GA 12-16ox . To glean a mechanistic explanation, each egression 313 pathway was analysed separately where the counting of hydrogen bonds and hydrophobic 314 contacts during production runs found that GA 7 had significantly more hydrogen bonding to 315 residues TYR247, GLY115, and VAL319. During these progression runs the GA 7 ligand 316 exhibited greater mobility than any of the other putative active GA (x) ligands, in both the gG and 317 gDG systems. A simple mechanistic explanation for markedly short egression times in GA 7 318 could be differences in mobility within the pocket of the DELLA • GID1A complex due to the 319 lack of an alcohol group at C13. However, no correlation was found between RMSD values and 320 egression times across the GA (x) to support such an explanation, Figure S5. 321 322 Table 5: Hydrogen bonding and contact recurrent residues during Cleft pathway egressions gG. Recurrent hydrogen bonds, in 323 blue, and hydrophobic contact residues, in black, during production runs of the gDG system. Recurrent set generated from GA (a) 324 subset and any pocket residue reported in Table 1 is removed for clarity, and otherwise, X denotes absence of interaction in 325 given system, while checks denote presence of interaction.

GA Res
327 Analysis of the contact residues in the cleft pathway are shown in Table 5. In the complexed 328 gDG system the cleft pathway is more complicated as the GAI MoRF 1 largely covers the cleft, 329 but the space for this pathway is intrinsically present in the binding pocket of GID1A. As a GA (x) 330 leaves GID1A in the gDG through the cleft, it must directly pass the DELLA residues in MoRF 1 331 near GAI-LEU32, as highlighted in Table S2, and the nearby GID1A residue LEU217.

332
DELLA dynamics and tRAMD uncapping 333 As a complementary analysis to removing the GA ligand from the complex pocket, tRAMD 334 egressions on the GAI N-terminus were also performed to explore the coordination of DELLA 335 onto GID1A. This was performed for all gDG systems previously discussed including the apo 336 DELLA • GID1A system (aDG). Important residue subsets described as MoRFs are known to be 337 conserved features of the GAI N-terminus and considered to be responsible for recognizing and 338 binding to GID1A [22]. The residues for MoRF 1 and 2 were included in the Hao modelling but 339 the residues linking MoRFs 1 and 2 and beyond were missing from their model. Nevertheless, 340 Hao et al. concluded there are differences between holo, with GA 4 , and the apo DELLA • GID1A 341 complex at the incomplete protein interface. DELLA protein residues, most of which are within 342 or near the later part of MoRF 1, were found in that study to change interaction by conformation 343 changes induced when GA 4 was present in GID1A. These residues interacted with GID1A N-344 terminal residues LEU9 and ASN32 on the alpha helix lid of GID1A. Although it was concluded 345 that GA 4 recognition leads to a conformationally stabilizing response, the previous work was 346 mainly informed by the static protein structure from X-ray crystallography. In light of these 347 critical differences we analysed both residence time (Table 6) and key contact residues for the 348 MoRFs in tRAMD uncapping modelling experiments. 349  364 The residues previously identified as important in GID1A LEU9 and ASN32 were found to be 365 entirely conserved when analysing the hydrogen bond and contact interactions with MoRF 1 366 during production runs of the complex containing the GA (a) ligands 1, 3, 4, 7, and 12-16ox.
367 Hydrogen bonding to GID1A residues 9 and 32 were not statistically different between GA (a) • 368 DG and the apo aDG systems. Comparing the interaction behaviour of the residues of MoRF 1 369 between GA (a) and GA (i) , only LEU9 was found to have conserved contacts.
370 Analysing the uncapping of DELLA (Table 6) Table 1 were recurrent in these uncapping 377 simulations as well.
378 MoRF 1 also showed long lasting hydrogen bonds to ASN19, LYS28, ARG51, and ASN326 379 residues of GID1A. Interaction with the N-terminal of GID1A were found to be largely 380 stabilized by hydrophobic interactions. While both of these interaction sets persist for a large 381 number of the uncapping events, the interactions often break before the end of the DELLA 382 egression. It is notable that the DELLA-GLU26 from the eponymous DELLA motif of MoRF 1 383 directly interacts with ASN19 of GID1A and was observed to be long-lasting during uncapping.
384 Analysis of MoRF 2 shows recurrent interactions of the active GA systems are largely intra 385 chain interactions, suggesting MoRF 2 controls the tertiary structure of the DELLA subunit.
386 MoRF 2 also maintains long-lasting hydrogen bonds to TYR48 and ARG51 of GID1A.
387 Comparisons between the GA (x) systems and the APO system hydrogen bonding of MoRF 2 388 showed inter chain residue interactions to GID1A residue ARG51 to interact significantly more 389 in the APO system. GID1A's ARG51 likely swings between MoRF 1 and 2 during normal 390 dynamics. 391 MoRF 3 was expected to only interact with intra-chain residues within DELLA due to its 392 position above the entire complex (Figures 1 and 5). Interactions of the main bioactive GA 393 systems showed the expected intra-chain interactions to consist of highly localized residues, 394 within plus or minus five residues of MoRF 3. However, between the gDG and aDG systems the 443 Furthermore, any GA ligand present in the DELLA • GID1 complex is able to stabilize the GAI 444 subunit in its attachment to GID1A. This is suggested by results from comparing uncapping 445 residence times, finding that they are, on average, higher for all tested gDG systems than for the 446 aDG systems. This result applies to each GA (x) case studied separately. Taking this result alone, 447 we cannot discriminate any critical differences between GA (x) that would classify them as 448 bioactive or inactive. This is not to say there are no differences. Consistent with our production 449 run results, it is likely that the different GA (x) create differences in conformational dynamics 450 within the gDG complexes (not studied here) that affect the degradation of the attached DELLA 451 protein.
452 Relatively lower RMSD was found in the MoRFs compared to the linkers between the MoRFs in 453 both gDG and aDG systems. Overall mobility in the DELLA protein is less in gDG compared to 454 aDG for each GA (x) . As Figure 2 shows, the region between MoRF 2 and MoRF 3 has greater 455 mobility with known bioactive GA compared to known inactive GA variants. The RMSF for the 456 carbon alpha atoms in the key interacting pocket residues of GID1A are found to be less mobile 458 has greater mobility in aDG than found within the gDG systems, see Figure 2 and Figures S8 and 459 S9 for additional RMSF data. These comparisons were made across our GA (x) set. Shifts in 460 mobility from one region to another is easily understood by an enthalpy/entropy compensation 461 mechanism driven by rigidity changes (2). Conformational entropy increases in certain regions 462 that become more flexible, while other regions become more rigid with more favourable 463 enthalpy due to conformational shifts in a protein. The greater RMSF found in the DELLA 464 protein when complexed in the aDG system is consistent with experimental data that without GA 465 the DELLA protein does not bind to GID1A. In essence, any of the GA (x) forms can act as a 466 cofactor. Furthermore, if GA leaves gDG, it is expected that aDG becomes an unstable structure 467 with a relatively short lifetime relative to the gDG lifetime. However, our data show it has a 468 relatively long lifetime relative to GA binding to aDG through the channel pathway.  Table 5 when compared to channel egressions in Table S1 in the GID1A 491 system. Along with the easier arrangement of GA carboxylate to organize into the pocket and the 492 higher likelihood of observation of cleft egression (see Table 4), this suggests an easier route in 493 the absence of DELLA motif. However, when GAI is bound to GID1A, by MoRF 1 and 2, the 494 contact residues to GA at the surface of the cleft pathway and the side chain flexibility needed to 495 open the cleft pathway are hindered. Egressions along the cleft pathway were still observed with 496 GAI present but drastically altered. In the complexed state GA was forced to take paths directly 497 under the DELLA MoRF 1 segment (see GAI-LEU32 contacts in Table S2). This suggests that 498 GA (x) egressions from the closed cleft will destabilize the gDG structure, as hydrophobic 499 interactions at this interface are disturbed. This in turn could help control the necessary GA (a) 500 recognition, as the more bioactive GA ligands are less likely to leave or to interact with this 501 aliphatic region. MoRF 2 interactions are found to be similar between the GA (a) and GA (i) 502 systems, removing it as a possible source for discrimination. We also note that a number of the 503 channel contacts evident in GA (x) egression pathways are missing when the ligand is GA (i) . Thus, 504 there is a higher propensity for GA (i) to be expunged from the closed cleft than a GA (a) . 505 Dynamical differences in the intra-chain interactions within DELLA are observed when 506 bioactive gDG and aDG are compared.
507 The opening of the cleft pathway in GID1 is a hole that can fit a GA (x) shaped molecule, and it 508 forms when the unstructured loop that LEU323 is attached to becomes part of a partial helix.
509 This releases helix B of GID1 to tilt outward as it lacks interaction with that loop (see Fig S2).
510 This motion appears as a natural breathing motion within GID1 when the DELLA protein is 511 absent. It is worth mentioning that based on structural similarities, this cleft might be the remnant 512 of the clade IV carboxylesterases (CXE) binding pocket in ancient plants (36). Assisting 513 interactions with HIS119 allow the GA (a) ligands to twist, having the alkyl backbone face 514 LEU323 and the aliphatic residues adjacent. These events were observed as essential dynamics 515 of this localized region based on a principal component analysis. In Figures S3, S4, S6, and S7 516 we compare the influence of DELLA binding on these essential motions (37). From analysing 517 the pocket residues that contact the ligand when fully docked, it is likely that key interaction to 518 TYR322 prevent the LEU323 loop from opening up and allowing a facile egression. This 519 tyrosine may also be important for bioactivity because it directly interacts with the key hydroxyl 520 moiety on C3 of GA. 521 Based on the observation that the cleft opening is relatively large with a shorter passage to the 522 binding pocket compared to the narrower and longer channel passage, we suggest GA (x) binding 523 to GID1A through the cleft in absence of DELLA protein is a significantly more rapid process 524 than through the channel. The GAI is assumed to bind to gG with little distinction in GA (x) . 525 Although it has been shown that mutations in GID1B can be made to facilitate GAI binding 527 is to act as a cofactor to facilitate the gDG complex to form. In other words, it is not the case that 528 an aDG complex cannot be stable, but rather, it is important for biological function to have the 529 presence of GA as a facilitating condition in order to control growth regulation in a plant.
530 Furthermore, GA (x) dissociates from gDG on a much faster time scale than GAI can dissociate 531 from GID1, implying there will be some degree of binding of DELLA to GID1 in the absence of 532 GA.

533
Proposed rectification mechanism that can select for bioactivity 534 We propose a novel kinetic model based on the newly found cleft pathway that also takes into 535 account the role played by the DELLA N-terminus that includes three MoRFs. Our kinetic 536 model elucidates how this system recognizes bioactive GA variants under competitive binding. 537 Furthermore, we performed a phylogenetic analysis that shows there is conservation of key 538 residues within the binding pocket, cleft and channel pathways -suggesting that the same 539 general mechanism in the GAI-GID1A interaction will apply to all DELLA • GID1 complexes. 540 There is uncertainty in whether a particular GA is bioactive for a given receptor. Therefore, we 541 introduce a more appropriate language to describe GA variants in terms of major and less active 542 forms of GA (x) relative to a given receptor. Subject to a wide spectrum of activities, GA (a) 543 denotes a subset of major active forms, while GA (i) denotes a subset of less active forms. For 544 simplicity of this discussion, the two subsets represent extremes in the activity spectrum for the 545 GID1A receptor (active versus inactive). The rank order of GA variants from most bioactive to 546 inactive would likely be different with respect to a different receptor. If the sequence of another 547 receptor is highly conserved with respect to GID1A in the binding pocket as well as the two 548 binding pathways, then the rank ordering of GA activity will presumably be similar. To 549 emphasise generality of the kinetic model, we discuss its components in terms of the GAI and 550 GID1 family of proteins. 551 Our results suggest that a molecule resembling the columbic and geometrical properties of GA (x) 552 can bind into GID1, particularly via the cleft pathway. Notably this includes GA variants that are 553 often only described as "inactive" or "intermediates" in the literature. Therefore, a degree of 554 caution should be exercised when classifying GA (x) as bioactive or not. With this in mind it is 555 axiomatic that recognition of bioactivity must involve GID1 in some capacity. If GID1 can bind 556 an array of GA (x) that likely includes both bioactive and inactive examples, then, this suggests a 557 mechanism of action removes inactive (or less active) forms of GA from the receptor or affect 558 the DELLA protein binding in other ways. Differentiation between GA (a) and GA (i) can occur 559 before and/or after DELLA is bound to the gG system. As shown in Figure 6, we construct a 560 generic 1st order kinetic model for a proposed rectification process. 561 562 Figure 6: Synergetic three-stage kinetic model of DELLA protein assisted rectification of GA in GID1A receptor. 563 The GA (x) in Figure 6 refers to any GA variant, but once a GA variant interacts with GID1, it 564 then becomes possible to classify the GA as a major or less active form. As a simplification, only 565 the extremes are tracked to show the possibility for competition between active and inactive 566 forms of GA, denoted as GA (a) and GA (i) respectively. This simplification demonstrates salient 567 features of the model by considering kinetic rates for only two types of GA (active and inactive), 568 but in reality, different kinetic rates will likely apply to each GA variant comprising a broad 569 spectrum of activities.
570 In state A, we have a free GID1 in an environment containing a mixture of GA variants, where 571 GID1 has opportunity to bind to any GA variant, represented by GA (x) . Since the open cleft does 572 not have a high level of specificity for the type of GA (x) to bind, the cumulative concentration of 573 all GA variants drives GA association to state B. The GA • GID1 complex is now primed to bind 574 with GAI. Since GAI association is a relatively slow process, it is expected that many cycles of 575 association and dissociation of different GA variants through the open cleft in GID1 occurs in 576 the interim. The key role of the cleft opening is to multiplex the kinetics of all GA (x) associations 577 that can take place to populate state B while depleting state A as much as possible depending on 578 a variety of details. While details depend on environmental conditions such as concentration of 579 repressor, receptor, protease, and GA variants, the qualitative features remain robust. Eventually, 580 GAI binds to GID1 to form the GAI • GA • GID1 complex that defines state C. 581 The details of which GA variant is bound affects the stability and conformational dynamics of 582 the GAI • GA • GID1 complex, which in turn affects the propensity for the E3 ubiquitin ligase to 583 exercise proteolysis to degrade GAI. The GAI • GA (a) • GID1 complex is primed for GAI to 584 undergo degradation, but the GAI • GA (i) • GID1 complex with less active GA varaints is