Origin of gibberellin-dependent transcriptional regulation by molecular exploitation of a transactivation domain in DELLA proteins

DELLA proteins are land-plant specific transcriptional regulators known to interact through their C-terminal GRAS domain with over 150 transcription factors in Arabidopsis thaliana. Besides, DELLAs from vascular plants can interact through the N-terminal domain with the gibberellin receptor encoded by GID1, through which gibberellins promote DELLA degradation. However, this regulation is absent in non-vascular land plants, which lack active gibberellins or a proper GID1 receptor. Current knowledge indicates that DELLAs are important pieces of the signalling machinery of vascular plants, especially angiosperms, but nothing is known about DELLA function during early land plant evolution or if they exist at all in charophytan algae. We have now elucidated the evolutionary origin of DELLA proteins, showing that algal GRAS proteins are monophyletic and evolved independently from those of land plants, which explains why there are no DELLAs outside land plants. DELLA genes have been maintained throughout land plant evolution with only two major duplication events kept among plants. Furthermore, we show that the features needed for DELLA interaction with the receptor were already present in the ancestor of all land plants, and propose that these DELLA N-terminal motifs have been tightly conserved in non-vascular land plants for their function in transcriptional co-activation, which allowed subsequent exaptation for the interaction with the GID1 receptor when vascular plants developed gibberellin synthesis and the corresponding perception module.


Introduction
DELLA proteins are transcriptional regulators that have been extensively characterized during the past twenty years (Vera-Sirera et al. 2015). They are involved in diverse processes ranging from seed germination to flowering, including legume nodulation, stress responses, or fern sexual reproduction (Peng and Harberd 1997;Floss et al. 2013;Tanaka et al. 2014). In fact, these proteins are responsible for the dwarf phenotype that allowed the development of new crop varieties during the Green Revolution (Peng et al. 1999). DELLAs are one of the main elements that compose the gibberellin (GA) signalling pathway in vascular plants, acting as the negative regulators of the pathway (Dill et al. 2001;Itoh et al. 2002).
As part of the GRAS family of plant specific proteins, they present a highly conserved Cterminal domain, the GRAS domain. Initially, this domain was suggested to be distantly related to the STAT family of metazoan proteins (Richards et al. 2000). More recently, a thorough in silico structural analysis of the domain has evidenced a remarkable similarity to bacterial Rossman-fold SAM-dependent methyltransferases suggesting a bacterial origin of the GRAS domain (Zhang et al. 2012). Even though no chlorophytan alga presents GRAS-like genes, several reported only in vascular plants (Hirano et al. 2007;Yasumura et al. 2007). These early studies indicate that GA-dependent DELLA regulation first appeared in the vascular plants common ancestor, however the analyses were constrained by the limited availability of genomic and transcriptomic resources from only the moss Physcomitrella patens, and the lycophytes Selaginella moellendorffii and S. lepidophylla. In these studies, neither a clear set of late GA synthesis genes nor proper GID1 receptors were detected in non-vascular plants, supporting the idea of GA-mediated regulation of DELLA proteins being vascular plant specific. In fact, no reports are available for the presence of bioactive GAs in mosses, and application of these compounds have no effect on moss growth (Hayashi et al. 2010).
Current knowledge indicates that DELLAs are important pieces of the signalling machinery of late diverging land plants, especially angiosperms, but nothing is known about DELLA function during early land plant evolution or if they exist at all in charophytan algae. The previous lack of data can now be completed with new genomic and transcriptomic sources from earlier diverging land plants and algae to understand the origin and emergence of DELLA proteins. In Table 2), those present in non-vascular plants do not contain the amino-lid sequences necessary for GA perception and DELLA interaction, and resemble those found in Physcomitrella patens (Hirano et al. 2007;Yasumura et al. 2007). Hence, we did not consider them as GID1 receptors.
However, we found two hornwort sequences (an almost complete sequence for Phaeoceros carolinianus and a partial sequence for Paraphymatoceros halli) that represent good candidates for an ancestral state of GID1 receptors, pointing to a possible pre-tracheophytan origin of putative GID1 proteins (Supplementary Fig. 1). Among non-vascular land plants, we were able to find GID2 orthologous sequences only in liverworts (Supplementary Fig. 2, Supplementary Table 3). These data consolidate the presence of GRAS in charophytan algae, not only in Zygnematales but also in Coleochaetales, and is consistent with the idea of GRAS proteins appearing first in charophytes before land colonization. Besides, DELLA proteins most likely appeared in the land plant common ancestor, but the GA signalosome components can only be found simultaneously in tracheophytes (Fig. 1). However, these views may change with further improvement of genomic data quality and availability from Charales and/or hornworts.

Phylogenetic analysis of GRAS proteins
To elucidate the origin of the DELLA subfamily, we analysed the phylogenetic relationship between GRAS sequences in algae and land plants. For this, we used the GRAS domain of the sequences found, and added previously described eubacterial GRAS sequences to use as outgroup in a phylogenetic tree (Zhang et al. 2012). Interestingly, algal and land plant sequences formed two independent and statistically supported clades ( Fig. 2A). This suggests that all land plant GRAS genes arose from a single gene present in an algal and land plant common ancestor. Further expansion and loss of GRAS subfamilies has occurred independently several times during plant evolution, and no clear correlation between the number of GRAS sequences and factors such as biological complexity seems to exist (Fig. 2B). Among land plants, we found sequences from twelve known GRAS subfamilies in Arabidopsis thaliana (SCL9, SHR, PAT1, SCL16/32, SCL29, SCL3, DELLA, SCL28, SCR, LAS, SCL4/7 and HAM), and sequences with no clear A.
thaliana match in at least two early diverging land plants that resemble RAM1 sequences (Fig.   2C). We conducted phylogenetic analysis of these GRAS domain sequences and obtained highly supported clades for these groups in all land plant lineages (Fig. 2D). In fact, these groups greatly coincide with those recently published in the Marchantia polymorpha genome (Bowman et al. 2017). Altogether, these analyses indicate that previously known GRAS subfamilies are land plant specific and appeared early in a land plant common ancestor. Consequently, we consider that only land plant genomes may contain DELLA genes.

Early evolution of DELLA proteins
To elucidate DELLA evolution, we generated a new phylogenetic tree adding previously undetected DELLA proteins from species belonging to different clades across land plant phylogeny (Fig. 3A). The N-terminal domain was also excluded from this analysis because the high level of divergence in this region yielded trees that were in conflict with known taxonomic relationships ( Supplementary Fig. 3). We confirmed the previously reported major DELLA clades corresponding to the eudicot clades RGA and RGL (DELLA1 and DELLA2, respectively), which are fused into a single DELLA1/2 clade in non-eudicot tracheophytes, and the DELLA3 clade, also named DGLLA/SLRL, that is present in all vascular plant lineages. These three clades are found in every major clade analysed, with the sole exception of DELLA1/2 being absent in ferns. This data suggests the occurrence of two main duplication events in DELLA genes coinciding with the appearance of tracheophytes and the emergence of eudicots. However, multiple duplications and loses have occurred in specific groups and species, such as the lack of DELLA1 and DELLA3 in Solanum lycopersicum.
Due to the scarce knowledge of DELLA function in non-vascular land plants, we expanded our search for DELLA sequences in liverworts, mosses and hornworts and analyzed their 6 phylogenetic relationship (Fig. 3B). As suggested by the previous tree, no ancient major duplications in DELLAs seem to have been maintained in land plants prior to vasculature emergence. We also detected a type of DELLA proteins sequences that we named DELLA-like proteins (DLP), whose phylogenetic position is unclear and most likely represents a liverwortspecific duplication with no resemblance in their N-terminal domains to those of DELLA proteins ( Supplementary Fig. 3 To determine the precise motifs and residues that provide high significance value for the identity of the N-terminal domain, we aligned the corresponding regions of representative sequences from each clade (three liverworts, ten mosses, three hornworts, and three tracheophytes) (Fig. 4). The three important motifs for the interaction with the GID1 receptor (and, therefore, for GA signalling) in tracheophytes were differentially conserved among non- This distribution of motifs suggests that the N-terminal domain of DELLA proteins was established early in land plant ancestors and maintained during evolution. To confirm this hypothesis, we performed ancestral protein reconstruction by maximum likelihood methods, excluding late divergent moss sequences to avoid bias towards DELLA or VHYNP motif loss.
To avoid the lack of consensus in bryophyte relationships and early land plant evolutionary history, we performed the analysis using different phylogenetic relationships among the three early diverging land groups (hornworts, liverworts and mosses), with almost identical results. The predicted ancestral sequences for the land plant common ancestor DELLA N-terminal domain harbours the canonical motifs and is strictly conserved compared to those known in tracheophytes and in hornworts (Fig. 4). Three of the four alpha helices harboured in the N-terminal domain show an ancestral state highly conserved in late diverging DELLA proteins, coinciding with the alpha helices that form the surface interacting with GID1.
Two more pieces of evidence support that the putative function of the ancestral N-terminal domain of DELLA proteins needed to be maintained during evolution: the Ka/Ks ratio is particularly low (around 0.2) precisely in the region that interacts with GID1 in higher plants

DELLA interaction with GA-signalling components
The solid conservation of the N-terminal domain of DELLA proteins even in land plants that lack the necessary GA signaling elements suggests that the ancestral DELLA was preadapted for subsequent interaction with the GA receptors. To gather additional experimental evidence, we first modelled the structure of putative complexes between AtGID1a and DELLAs from species of different land plant taxa based on the previously described GA4-AtGID1a-AtGAI structure.
For comparison, the DELLAs were selected from A. thaliana (Angiosperma), S. moellendorffii (Lycophyta), Nothoceros vincentianus (Anthocerotophyta), P. patens and T. lepidozioides (IEELD). This mutation which should decrease the interaction with K28 (Fig. 5) does not seem important because the degradation of SmDELLA1 by GAs has already been described (Hirano et al. 2007;Yasumura et al. 2007).This result can be easily understood by the mobility of the side chain of Lysine (K) 28 which could likely interact with the aspartate residue of SmDELLA1. On the contrary, despite a similar fold, T. lepidozioides TlDELLA displayed critical changes in the LEQLE motif (absent in the alignment, substituted by the subsequent LGGAQ sequence in the model) of the αB helix which are predicted to prevent interaction with GID1. In addition, M. polymorpha and P. patens DELLAs present not only modifications in the αB helix, but also in the DELLA and VHYNP motifs preventing interaction with AtGID1a. In summary, the models indicate that, unlike lycophytes and hornworts, moss and liverwort DELLAs should be unable to interact with AtGID1a.
To experimentally test the models, we performed yeast-two-hybrid assays between the six full-length DELLAs and the Arabidopsis GID1 proteins in the presence or absence of GA3 ( Fig. 5B; Supplementary Fig. 5B). As expected, both tracheophytan DELLAs interacted with AtGID1s when GA was present. Interestingly, N. vincentianus DELLA was also able to interact with the receptors in a GA-dependent manner. However, moss and liverwort sequences did not show interaction with the receptors either in a GA-dependent or independent way. These results indicate that N-terminal domain conservation is necessary but not sufficient for GID1 interaction.
Moreover, the ability of hornwort DELLAs to recognize the Arabidopsis GA receptors in a GAdependent manner (Fig. 5B), and the presence of putative GA receptor sequences in some Anthocerotophyta genomes ( Fig. 1; Supplementary Table 2) suggests that the origin of GA signaling might predate the separation between hornworts and vascular plants. included those from non-vascular land plants ( Fig. 6B; Supplementary Fig. S6A). Previously, it has been suggested that both DELLA and VHYNP motifs are involved in this activity (Hirano et al. 2012). However, the N-terminal domain of PpDELLAa is strongly capable of transcriptional activation, but lacks both of these motifs. The most conserved region among all land plant DELLA N-terminal ends is the αD helix (Fig. 4), and we found that deletion of this region (ΔαD) in PpDELLAa prevented the induction of reporter expression in yeast, while the αD helix alone was also capable of promoting transactivation in a yeast two-hybrid assay ( Fig. 6C and Supplementary Figure 6B). To study if this also happens in planta, we confirmed these results by performing transient expression assays of a dual luciferase reporter in Nicotiana benthamiana (Fig. 6D, Supplementary Fig. 6C). The activity of certain transactivation domains from different origins (viral VP16, yeast GAL4 or PHO4, mammalian p53 or NFAT, and others) have been proposed to harbor a particular 9-amino acid transactivation domain (9aaTAD) (Piskacek et al.

2007) which directly interacts with the KIX domain of general transcriptional coactivators like
Mediator's MED15 subunit (Piskacek et al. 2016). Interestingly, the αD helix of the DELLA Nterminal domain displayed a high score in a 9aaTAD evaluation (Supplementary Fig. 7).
Our results suggest that the main driving force for the conservation of the N-terminal regions of DELLA proteins has been its role in transcriptional activation, and its eventual co-option by the GID1 receptor allowed hormonal regulation of DELLA stability. Although molecular exploitation has been described as an evolutionary strategy to expand hormone-receptor complexity (Bridgham et al. 2006;Baker et al. 2013), the function of the evolving receptor was equivalent to the ancestral one (i.e., interaction with the ligand). However, the origin of the GA signaling pathway illustrates how molecular exploitation can occur upon domains with completely unrelated functions. Curiously, the co-existence of degron and transactivation motifs in the same stretches of residues has been described in several mammalian transcriptional activators (Salghetti 2000;Salghetti et al. 2001). In contrast, this degron-TAD overlap has not been studied in plants but some examples can be identified, such as MYC2, in which a degron is found within the MID domain, the transactivation domain (Zhai et al. 2013). Therefore, DELLAs would have these two functions encoded in a single region, and interaction with GID1 has been reported not only to promote DELLA degradation, but also prevent transactivation by DELLAs (Hirano et al. 2012). In summary, the coincidence of transactivation and protein stability regulation in a single protein domain is a widespread property and has independently emerged several times during evolution and through different molecular mechanisms. Original sequences, raw and trimmed alignments, and trees are available at https://data.mendeley.com/datasets/bjcp6ggjk9/draft?a=6c3474c7-6b37-47d0-ad3d-19ad9c062627 Ancestral sequence reconstruction 11 DELLA N-terminal domain ancestral state for each codon position was determined using MEGA7 (Kumar et al. 2015). Translated nucleotide coding sequences for the previously mentioned DELLA were first aligned as mentioned above, and then untranslated to codon-aligned nucleotide sequences. Ancestral sequence inference was then performed using Maximum Likelihood, including different predefined tree topologies comprising monophyletic Bryophyta, and paraphyletic "bryophyta", both with a Setaphyta clade (liverworts and mosses), and either hornworts as sister clade to the rest of land plants or to vascular plants. A fourth topology included the step-wise paraphyly of "bryophyta" with liverworts as sister clade to the rest of land plants, followed by mosses, and hornworts. Finally, Tamura-Nei model of nucleotide substitution was used for ancestral state inference. Gap residues were trimmed if absent in more than 90% of the sequences.
Yeast-two-hybrid assay Arabidopsis GID1 and GID2 were fused to the Gal4-DNA binding domain (DBD) in the pGBKT7-GW vector as bait, and DELLA full length ORFs from the different species were fused to the Gal4-Activation domain (AD) in the pGADT7-GW vector. DELLAs, GID1 and GID2 ORFs were either amplified by PCR using sequence-specific primers (Supplementary Table 9) or synthesized as gBlocks ® (I.D.T) and transferred to pDONR221 or pDONR207 via BP Clonase II (Invitrogen), or to pCR8 via TOPO ® TA cloning ® (Invitrogen) to create entry vectors. Final constructs were made by recombining entry clones to GATEWAY destination vectors via LR Clonase II (Invitrogen). Direct interaction assays in yeast were performed following Clontech small-scale yeast transformation procedure. Strain Y187 was transformed with pGADT7 derived expression vectors, while strain Y2HGold was transformed with pGBKT7 vectors, and selected in SD medium without Leu or Trp, resp. Subsequently, diploid cells were obtained by mating and selected in SD medium lacking both Leu and Trp. Interaction tests were done in SD medium lacking Leu, Trp and His, in the presence of different concentrations of 3-aminotriazol (3-AT) (Sigma-Aldrich). GID1 interaction tests were supplemented with 100 μM GA3.
Yeast transactivation assay DELLA ORFs were obtained as described above. DELLA N-end clones and C-end clones were obtained by PCR amplification using sequence-specific primers (Supplementary Table 9) and transferred to pDONR207 via BP Clonase II (Invitrogen). All the entry clones were used to create Gal4-DBD fusions in the pGBKT7-GW vector via LR Clonase II (Invitrogen), transformed into yeast strain Y2HGold and selected in SD medium lacking Trp. Transactivation tests were done as yeast two hybrid assays, but with Leu supplementation.

In planta transient transactivation assay
The reporter construct containing 2xGal4 UAS upstream of a 35S minimal promoter previously and DELLA were counted either when Pfam (PF12041.7) was positive or BLASTP e-value was lower than 1E-20 when using either AtRGA or PpDELLAa. GID1 were counted by the presence described N-lid residue presence and α/β-hydrolase active site conservation. GID2 were selected by thresholding BLAST results with 1E-15 based on similarity to S. moellendorffii GID2 proteins.

Conserved to tracheophytes
Conserved to predicted ancestral residues  (A) Predicted structural model for DELLA-AtGID1a interaction using AtRGA, SmDELLA1, NvDELLA, PpDELLAa, TlDELLA, and MpDELLA. DELLA structure is shown in yellow and AtGID1a structure in light blue. AtGID1a residues involved in DELLA interaction are written in dark blue. Residues different to that of AtRGA/GAI in main motifs are presented in red. Possible residue to residue interactions affected are pointed with a red circle (B) Yeast-two-hybrid assay results between DELLA proteins and the three Arabidopsis GID1 receptors with or without GA 3 . Positive interactions are accounted when yeast growth occurs in 5 mM 3-aminotriazol. Constitutively expressed Renilla luciferase (REN) for normalization. Data shown are normalized to NE value and represent the average of three biological replicates. Error bars represent standard deviation. Letters indicate significant differences between groups (P < 0.01, one-way ANOVA, Tukey's HSD post hoc test).