A gene regulatory network critical for axillary bud dormancy directly controlled by Arabidopsis BRANCHED1

A gene regulatory network critical for axillary bud dormancy 1 directly controlled by Arabidopsis BRANCHED1 2 Sam W. van Es1,2,7#, Aitor Muñoz-Gasca3,#, Francisco J. Romero-Campero4, Eduardo 3 González-Grandío3, Pedro de los Reyes4, Carlos Tarancón3, Aalt D.J. van Dijk5,6, Wilma 4 van Esse1,2, Gerco C. Angenent1,2, Richard Immink1,2,*, Pilar Cubas3,* 5 1. Bioscience, Wageningen Plant Research, Wageningen University & Research, 6708 PB 6 Wageningen, The Netherlands 7 2. Laboratory of Molecular Biology, Wageningen University & Research, 6708 PB 8 Wageningen, The Netherlands 9 3. Department of Plant Molecular Genetics, Centro Nacional de Biotecnología/Consejo 10 Superior de Investigaciones Científicas, Campus Universidad Autónoma de Madrid, 28049 11 Madrid, Spain 12 4. Plant Development Unit, Institute for Plant Biochemistry and Photosynthesis, 13 Universidad de Sevilla - Consejo Superior de Investigaciones Científicas, E-41092 14 Department of Computer Science and Artificial Intelligence, Universidad de Sevilla, E- 15 41012 Seville, Spain 16 5. Bioinformatics Group, Department of Plant Sciences, Wageningen University & 17 Research, Wageningen, 6708 PB, The Netherlands 18 6. Biometris, Department of Plant Sciences, Wageningen University & Research, 19 Wageningen, 6708 PB, The Netherlands 20 7. Current address: Umeå Plant Science Centre, Department of Plant Physiology, Umeå 21 University, SE-901 87 Umea, Sweden 22 23 *For correspondence (email: [email protected]; [email protected]) 24 #These authors have contributed equally to this work 25 preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted December 15, 2020. ; https://doi.org/10.1101/2020.12.14.394403doi: bioRxiv preprint 26 Abstract 27 The control of branch outgrowth is critical for plant fitness, stress resilience and 28 crop yield. The Arabidopsis thaliana transcription factor BRANCHED1 (BRC1) plays 29 a pivotal role in this process as it integrates signals that inhibit axillary bud growth 30 to control shoot branching. Despite the remarkable activity of BRC1 as a potent 31 growth inhibitor, the mechanisms by which it promotes and maintains bud 32 dormancy are still largely unknown. 33 Here we combine ChIP-seq, transcriptomic and systems biology approaches to 34 characterise the BRC1-regulated gene network. We identify a group of BRC1 direct 35 target genes encoding transcription factors (BTFs) that orchestrate, together with 36 BRC1, an intricate transcriptional network enriched in abscisic acid signalling 37 components. The BRC1 network is enriched in feed-forward loops and feed-back 38 loops, robust against noise and mutation, reversible in response to stimuli, and 39 stable once established. This knowledge is fundamental to adapt plant architecture 40 and crop production to ever-changing environmental conditions. 41 42 Introduction 43 In multicellular organisms, growth and development occur by the sequential 44 activation of gene regulatory networks (GRNs) that influence patterns of cell 45 division expansion and differentiation in developing tissues, which results in 46 species-specific morphologies. Sessile organisms, such as plants, need to adapt 47 these developmental programs to resource availability and environmental 48 conditions for fitness and survival. Thus, the activity of these GRNs, controlled by 49 preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted December 15, 2020. ; https://doi.org/10.1101/2020.12.14.394403doi: bioRxiv preprint transcriptional master regulators, is modulated by external signals, such as light 50 quality, photoperiod and temperature, and endogenous cues that signal nutrient 51 and energy limitations or stress. 52 In higher plants, one key developmental program that determines aboveground 53 plant architecture is the formation of branches, which are lateral shoots derived 54 from axillary meristems. Axillary meristems are groups of undifferentiated cells 55 formed in the leaf axils, which generate axillary buds or branch primordia. Under 56 favourable conditions, these axillary buds continue to grow and elongate to form 57 branches. In contrast, under actual or potentially energy-limiting situations, the 58 buds remain dormant1. The status of each bud, either active or quiescent, is thus 59 influenced by the environment and coordinated at the plant level by sugar 60 availability, strigolactone (SL) hormone signalling, and auxin transport2. However, 61 the molecular operations required for cells in axillary buds to integrate those signals 62 locally, and implement the pathways of dormancy or growth are still poorly 63 understood. 64 In Arabidopsis thaliana, BRANCHED1 (BRC1) is a central local regulator of bud 65 activity that encodes a TEOSINTE BRANCHED1, CYCLOIDEA, PCF (TCP) 66 3transcription factor (TF), whose function is widely conserved in 67 angiosperms. BRC1 and its orthologs are expressed in axillary buds and promote 68 bud growth arrest: brc1 mutants have increased axillary bud activity and excessive 69 branching in a wide variety of species2. Due to the spatially restricted expression 70 patterns of BRC1, limited to axillary buds, this gene only prevents lateral shoot 71 outgrowth. However, when ectopically expressed in Arabidopsis seedlings, it can 72 also cause growth arrest of shoot and root apical meristems and leaf primordia4. In 73 preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted December 15, 2020. ; https://doi.org/10.1101/2020.12.14.394403doi: bioRxiv preprint hybrid aspen, BRC1-like genes expressed in the shoot apex in short days also 74 promote growth cessation5. Despite these remarkable growth-inhibiting effects of 75 BRC1, and its major role in the control of plant architecture, the gene targets and 76 GRNs controlled directly by BRC1 remain largely unknown in Arabidopsis. Three 77 closely related HD-zip-I genes (HB21, HB40, and HB53) are direct BRC1 targets 78 and mediate BRC1-induced abscisic acid (ABA) synthesis and response in buds. 79 ABA signalling in turn promotes inhibition of bud growth6. Furthermore, two 80 clusters of co-expressed genes -one enriched in genes related to cell division and 81 DNA replication, another in chloroplast and ribosomal genesare downregulated 82 in response to BRC1, although it is unclear whether this control is direct4. In maize, 83 the closely related TCP factor Teosinte Branched1 (TB1) directly influences 84 hormone signalling (ABA, gibberellic acid (GA), and jasmonate (JA)) and sugar 85 signalling7 but it is unclear whether the TB1 targets are shared with Arabidopsis 86 BRC1, as there is abundant evidence of the divergent regulatory pathways 87 associated with TB1 and BRC1. 88 In this work, we have characterised the transcriptional network that operates 89 directly downstream of Arabidopsis BRC1 to promote axillary bud growth arrest. 90 For this, we have compared the transcriptome of seedlings after BRC1 induction, 91 the genome-wide binding sites of BRC1 by Chromatin Immunoprecipitation 92 sequencing (ChIP-seq), and the transcriptional profiling of active vs dormant 93 axillary buds. We have identified a group of direct targets of BRC1 encoding TFs 94 that we have termed BTFs, which play a critical role in the regulation of the 95 network. We have integrated this BRC1 ChIP-seq data with ChIP-/DAP-seq data of 96 the BRC1-dependent TFs into a regulatory network, BRC1NET. To better explore 97 preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted December 15, 2020. ; https://doi.org/10.1101/2020.12.14.394403doi: bioRxiv preprint and characterise the network properties we have developed an open interactive 98 web application. Our findings reveal a pivotal role of BRC1 as a master regulator 99 of an intricate transcriptional network enriched in ABA-related components. This 100 transcriptional network, robust against noise and mutation, is reversible in the initial 101 steps but becomes stable once established. These properties help fine-tune the 102 transcriptional responses promoting bud dormancy in response to limiting 103 conditions, an adaptation critical for plant fitness and survival. 104 105 Results 106 BRC1 regulates nine coexpressed gene clusters 107 To identify the genes regulated by BRC1, our first approach was to perform a 108 thorough in silico analysis of existing expression data. González-Grandío et al. 109 (2013) carried out transcriptional profiling of axillary buds of wild-type and brc1 110 mutant plants exposed for eight hours to either white light (WL), or white light 111 enriched in far-red light (WL+FR). The latter treatment simulates a canopy shade, 112 which leads to BRC1 mRNA accumulation and promotes bud dormancy. We used 113 the collection of genes differentially expressed (DE, FDR<0.05) in response to BRC1 114 (BRC1-dependent genes, Supplemental Figure 1) as a starting point for this study. 115 To investigate their potential functional relationships, we studied their degree of co- 116 expression in > 15,200 publicly available microarray experiments and clustered 117 them by their global similarity of transcriptional responses. This analysis revealed 118 six co-expression clusters of up-regulated genes and three co-expression clusters of 119 down-regulated genes (Fig. 1a). We then searched in silico for additional genes co- 120 expressed with members of each cluster and tested whether they were also DE in 121 preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted December 15, 2020. ; https://doi.org/10.1101/2020.12.14.394403doi: bioRxiv preprint dormant (WL+FR-treated wild-type) buds, but not in active (WL-treated wild-type, 122 or brc1 mutant) buds (Fig. 1b). Genes that met this criterion were added to 123 the BRC1-dependent category (Supplemental Fig. 1 and 2). We termed these 124 clusters of co-expressed genes, UP_C1 to C6 and DOWN_C1 to C3, for the up- and 125 down-regulated genes in dormant buds, respectively (Fig. 1b, Supplemental 126 Dataset 1). 127 Each co-expression cluster was enriched in a different set of gene ontology (GO) 128 terms9 (Fig. 1a, Supplemental Dataset 1): UP_C1 was enriched in genes involved 129 in ethylene signalling, senescence, and catabolism; UP_C2, in genes related to 130 water transport; UP_C3, brassinosteroid and auxin signalling; UP_C4, response to 131 light; UP_C5 did not show an statistically significant enrichment in GO terms, but 132 contained genes related to cell wall remodelling and GA signalling, among 133 others; UP_C6, was enriched in genes related to ABA signalling and response to 134 abiotic stress; DOWN_C1, in genes related to chloroplast activity and 135 organization; DOWN_C2, in DNA replication and cell division; 136 and DOWN_C3 in genes related to biosynthesis of secondary metabolites such as 137 flavonoids, glucosinolates and pectins. 138 This collection of 947 genes constitutes a core transcriptional response associated 139 with BRC1 activity and to the growth-to-dormancy transition in Arabidopsis axillary 140 buds. 141 142 Identification of BRC1 direct targets 143 We next investigated how BRC1 orchestrates the expression of these BRC1- 144 dependent gene clusters in buds. BRC1 could directly regulate the expression of a 145 preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted December 15, 2020. ; https://doi.org/10.1101/2020.12.14.394403doi: bioRxiv preprint majority of genes from each cluster, or it could control a limited number of key 146 regulators, which in turn would mediate a wider transcriptional response. To 147 discern between these possibilities, we looked for direct BRC1 targets in the 148 Arabidopsis genome using ChIP-seq9. BRC1 expression levels are low and limited 149 to axillary buds10, which complicate the study of protein-DNA binding events11. We 150 therefore performed ChIP in ten-day-old seedlings of the estradiol-inducible 151 line LEXA:GFP:BRC1ind; brc1-2 (BRC1ind)6. We collected material five hours after the 152 beginning of estradiol induction, when the GFP:BRC1 protein was detectable6. 153 ChIP-seq was performed for three independent biological replicates (Supplemental 154 Fig. 3, Supplemental Dataset 2). We found that 4339 loci contained BRC1 peaks 155 within 3 kb upstream of the TSS. A large proportion (~60%) of the BRC1 peaks were 156 located in the proximal upstream region of protein-coding genes, on average within 157 300 bp of the transcription start site (TSS) (Fig. 2a). Around 4000 loci (~60%) 158 contained peaks within 3 kb upstream of their TSS in the two replicates with the 159 highest sequencing and mapping depth. The other peaks were positioned in introns, 160 exons, or UTRs (totally ~ 9%), and the remaining ~10 and ~20% in the downstream 161 and intergenic regions, respectively (Fig. 2b). This distribution was comparable to 162 that found for ChIP-seq performed on other Arabidopsis TFs (e.g.12,13), but divergent 163 from that described for cucumber TEN, a closely related TCP factor, for which a 164 large proportion of binding sites are localised in gene bodies14. 165 BRC1-bound peaks in the two most enriched ChIP-seq samples showed a 166 significant overrepresentation of several motifs (Supplemental Table 1, DREME15 167 most notably GDCCCA and BACGTGKC (D=AGT; B=CGT; K=GT)(Fig. 2c). Around 168 60% (57.5 and 65.1% for samples 1 and 2, respectively) of all peaks contained the 169 preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted December 15, 2020. ; https://doi.org/10.1101/2020.12.14.394403doi: bioRxiv preprint GDCCCA motif, which closely resembles a TCP binding site16,17 and the BRC1 170 binding motif GGgcCCmc, inferred by Protein-binding microarray assays6. 171 Remarkably, 50% of all peaks (48.3 and 51.9% for samples 1 and 2, respectively) 172 contained the BACGTGKC motif, which closely resembles a G-box (CACGTG) 173 recognised by bHLH and bZIP TFs18. Interestingly, in approximately 30% of the 174 peaks (28.5 and 34.5% for samples 1 and 2, respectively), the TCP binding site was 175 accompanied by a G-box-like motif. In these cases, the TCP motif was centrally 176 enriched, whereas the G-box was positioned either up- or downstream of the TCP 177 binding site (Fig. 2d). This suggests that BRC1 frequently acts in concert with G-box 178 binding proteins to control gene expression. The overrepresentation of a G-box-like 179 motif was not reported in the recently published ChIP-seq experiment for maize 180 TB17 and this prompted us to reanalyse that dataset. Our analysis confirmed that 181 the G-box is also overrepresented in the maize TB1 targets (Supplemental Results), 182 revealing that this feature is indeed conserved in maize, a monocot distantly related 183 to Arabidopsis. 184 In parallel to the ChIP study, we did a gene expression analysis (RNA-seq) in 185 Arabidopsis of similar estradiol-treated material, to identify genes directly 186 responding to BRC1 activity (differentially expressed genes, DEG, Supplemental 187 Dataset 3). We found 5718 DEGs of which 1845 were upregulated and 3873 were 188 downregulated after induction (false discovery rate, FDR <0.05). We considered 189 BRC1 direct targets those genes whose regulatory regions were bound by BRC1 (i.e. 190 BRC1-bound genes) and were DE after BRC1 induction in seedlings (Supplemental 191 Fig. 1). We identified 1438 BRC1 direct targets according to this criterion. 192 193 preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted December 15, 2020. ; https://doi.org/10.1101/2020.12.14.394403doi: bioRxiv preprint BRC1-dependent direct gene targets 194 To further refine the BRC1 target gene list, we compared the BRC1 direct targets 195 identified by ChIP-seq and RNA-seq with the BRC1-dependent genes identified in 196 the active-vs-dormant bud transcriptomics. Genes present in the three lists were 197 termed bona fide BRC1 direct targets (Supplementary Fig. 1, Supplemental Dataset 198 4). We found 97 genes; 90 induced and 7 repressed in response to BRC1, both in 199 seedlings and in buds entering dormancy. This strongly indicates that BRC1 acts 200 mainly as a transcriptional activator. The list was enriched in genes containing ABA- 201 related terms (ABA signalling, FDR= 0.045; ABA response, FDR= 7.4 10-4; ABA 202 transport, FDR=0.037), and genes responding to water deprivation and salt/osmotic 203 stress (FDR=0.01) (Supplemental Dataset 4). Another significantly enriched 204 category (FDR=0.037) was that of regulation of transcription: 19 bona fide BRC1 205 direct targets encoded TFs, many of which had ABA-related annotations (Table 1). 206 The list of bona fide BRC1 targets also contained induced genes encoding enzymes 207 controlling hormone metabolism and homeostasis such as MAX1, NCED3 and 208 ST2A (that regulate SL, ABA, and JA synthesis respectively,19–21) Beta-carotene 3- 209 hydroxylase 2 (BCH2, which catalyses the synthesis of zeaxanthin, a precursor of 210 SL and ABA); CYP707A2 and GA2OX2 (involved in ABA and GA degradation 211 respectively); and WES1/GH3.5 (that control auxin amino acid conjugation). 212 Among the 6 downregulated genes was UGT76C2 (that promotes cytokinin 213 glycosylation and inactivation). Direct targets of BRC1 are also genes encoding 214 ubiquitination-related proteins, as well as water, ion, amino acid and hormone 215 membrane transporters. 216 preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted December 15, 2020. ; https://doi.org/10.1101/2020.12.14.394403doi: bioRxiv preprint physically interact. High-throughput CrYTH-seq assays revealed a potential 358 interaction between BRC1 and ABI531 but this has not yet been confirmed in vivo. 359 Notably, although BRC1 expression is sufficient to activate BTFs and their 360 downstream targets6 (this work), these genes also belong to a well-characterised 361 ABA network, induced in tissues where BRC1 is not usually expressed (e.g. ABA- 362 treated seedlings24,32. This strongly suggests that BRC1 has co-opted a pre-existing 363 ABA signalling network to be activated (under BRC1 control) in axillary buds. This 364 could have been achieved by acquisition of BRC1 binding sites at promoters of top- 365 tier ABA master regulators to become BTFs (i.e. ABI5, ABF3, GBF2, GBF3, 366 ANAC032, HB21, HB40, HB53), which could in turn activate additional ABA- 367 related targets, including other TFs. BRC1 binding sites also evolved in these 368 secondary targets creating FFLs between BRC1 and the BTFs. The observation that 369 BRC1 is not induced by ABA treatments6 and is not co-expressed with ABA 370 networks in other tissues, supports the view that BRC1 is an ABA-independent, 371 master regulator of the ABA network, exclusively in buds. Perhaps BRC1 increases 372 chromatin accessibility of the BTFs to their targets in axillary buds, a tissue in which 373 otherwise BTFs might not bind their targets. Some BTFs (e.g. ABF3, ABI5, ATAF1) 374 require ABA-induced phosphorylation for full activation30. BRC1-mediated ABA 375 biosynthesis (by upregulation of NCED3 and BCH2) may be instrumental for this 376 post-translational activation, and for additional support of the ABA-related 377 responses. 378 The BRC1 transcriptional network is enriched in multi-output FFLs26 formed by 379 BRC1, BTFs and bona fide BRC1 targets. FFLs regulate reversible decisions and act 380 as sign-sensitive delay elements and persistence detectors. Some target genes may 381 preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted December 15, 2020. ; https://doi.org/10.1101/2020.12.14.394403doi: bioRxiv preprint require both BRC1 and a BTF for induction. These would show delayed responses 382 to BRC1, but no delay when BRC1 is downregulated. This would help filtering out 383 signals that transiently induce BRC1 expression, such as a brief shading. Only 384 persistent signals (e.g. a dense plant canopy) would lead to network induction. In 385 contrast, other genes could be induced by either BRC1 or a BTF. These targets 386 would have a rapid and persistent response even during a transient loss of BRC1 387 induction. BRC1NET is also enriched in multi-output FBLs formed by the BTFs - 388 which regulate each other- and jointly control downstream genes. Unlike FFLs, 389 FBLs provide memory of an input signal even after the signal is gone. They are 390 frequent in transcriptional networks that transduce signals into stable 391 developmental decisions. In BRC1NET, these motifs could ensure that 392 environmental or endogenous signals upregulating BRC1, can then be locked into 393 a steady state-response that stabilises bud dormancy. The final output will depend 394 on the prevailing network motifs, their sign (positive or negative), and additional, 395 BRC1-independent regulation of network components. It is plausible that the 396 degree of reversibility of bud dormancy (i.e. para- and ecodormancy vs 397 endodormancy) is in part determined by the prevalence of FFLs that allow delayed 398 responses and reversibility, or FBLs that promote more stable and lasting responses 399 of growth arrest. The BRC1NET is very robust against mutation, an essential feature 400 to maintain a high fitness state. In line with this, the network connectivity is only 401 mildly affected even after removal of up to 10 TFs (Supplemental Fig. 10). 402 In summary, during the evolution of angiosperms, BRC1/TB1 genes seem to have 403 recruited an ABA-related network comprising a collection of master regulators of 404 ABA signalling and their targets. This highly connected network, robust against 405 preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted December 15, 2020. ; https://doi.org/10.1101/2020.12.14.394403doi: bioRxiv preprint noise and mutation, constitutes a transcriptional cascade that in the first steps is 406 reversible and filters out transient signals (using FFLs), but becomes stable once 407 established (due to the FBLs). The flexibility of this network is essential for an 408 adaptive response and a correct carbon allocation in the plant under growth- 409 limiting conditions. This knowledge will prove fundamental to optimise plant 410 architecture and crop production. 411 412 Materials and methods 413 Plant materials and growth condition 414 Wild-type Arabidopsis thaliana plants were of the Columbia-0 (Col-0) ecotype. The 415 brc1-2 mutant was described10 and the GFP:BRC1ind lines have been generated 416 previously6. All lines were grown on 1/2MS medium 1% (w/v) agar plates under 417 long day conditions (16/8 light/dark cycle) at 21°C. 418 419 Promoter:GUS constructs 420 Genomic fragments of 1728, 1709, 1830, 2265, 1846, 540, and 1990 base pairs 421 comprising the BRC1 binding peaks upstream the ATG of ABI5, ABF3, GBF2, GBF3, 422 MYB3, ATAF1 and NAP, respectively, were amplified from genomic DNA using 423 primers in Supplemental Table 5 and Taq Phusion polymerase (New England 424 BioLabs), cloned in pGEM-t Easy (Promega), re-amplified by PCR with primers with 425 attB tails and BP-cloned into pDONR207 (Invitrogen). Promoter:GUS binary 426 vectors were generated by LR recombination using Gateway LR clonase II 427 (Invitrogen) according to manufacturer’s instructions into the destination vector 428 pGWB3. 429 preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted December 15, 2020. ; https://doi.org/10.1101/2020.12.14.394403doi: bioRxiv preprint 430 431 Generation of transgenic plants 432 Binary vectors were transformed into the Agrobacterium tumefaciens strain AGL-0. 433 Transgenic Arabidopsis plants were generated by agroinfiltration using the floral 434 dip method33. At least 10 T3 homozygous lines were generated and analysed for 435 each construct. 436 437 GUS staining 438 GUS staining was conducted as described in34. 439 440 RNA-seq and ChIP-seq material collection 441 For the RNA-seq experiment, plants of brc1-2 (control) and GFP:BRC1ind lines were 442 grown on plate. Ten-day old seedlings were induced with 5mL of 10µM estradiol 443 per plate. After an induction period of five hours, 1 g of tissue was harvested for 444 each of the three biological replicates. Total RNA was isolated using the InviTrap® 445 Spin Plant RNA Mini Kit according to the manufacturer’s protocol. TURBO™ 446 DNase was used to clean the RNA samples from DNA. Library preparation for 447 whole genome RNA sequencing was done using the Illumina Truseq Library 448 Preparation Kit. Library quality was evaluated using a Bioanalyzer and an RNA 449 Nano 6000 kit (Agilent). RNA concentrations were determined using the Xpose 450 ‘DSCVRY’ (Trinean). The libraries were then sequenced on the Illumina Hi-Seq 451 2500. 452 preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted December 15, 2020. ; https://doi.org/10.1101/2020.12.14.394403doi: bioRxiv preprint For the ChIP-seq experiment, plants of GFP:BRC1ind lines were grown on plate. Ten- 453 day old seedlings were induced with 5mL of 10µM estradiol per plate. After an 454 induction period of five hours, 1.5g of tissue was collected per sample. ChIP was 455 performed as described35, using µMACS Anti-GFP beads (Miltenyi). Input DNA was 456 used as control, which is isolated from the sonicated chromatin prior to 457 immunoprecipitation. 458 459 Identification of Genes Coregulated with BRC1-dependent genes 460 BRC1-dependent genes were obtained from4 Coregulation of the 181 UP and 57 461 DOWN was analysed by hierarchical clustering based on 15,200 publicly available 462 microarray experiments (Hcluster, ATTED-II36). Additional coregulated genes were 463 obtained with CoEx-Search (ATTED-II36). The expression of these genes was 464 validated in the original arrays 4. Only genes upregulated (positive fold change FC) 465 in WL+FR treated wild type dormant buds but not in WL treated or brc1 mutant 466 buds were included in the extended list of BRC1-dependent genes. 467 468 Functional Annotation of BRC1-dependent clusters 469 GO automated analyses of function prediction were carried out using a statistical 470 overrepresentation test (pval < 0.05) followed by Bonferroni correction for multiple 471 testing in the PANTHER classification system8. 472 473 Differential gene expression and gene ontology (GO) enrichment analyses 474 Libraries of three biological replicates were sequenced individually and analysed 475 using the Bowtie–Tophat–Cuffdiff (BTC) pipeline37. Differential gene expression was 476 preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted December 15, 2020. ; https://doi.org/10.1101/2020.12.14.394403doi: bioRxiv preprint determined for all estradiol induced GFP:BRC1ind samples, using induced brc1-2 477 samples as control. The cut-off was set at a false-discovery rate (FDR) <0.05 in all 478 analyses performed. The RNA-seq data is made available via NCBI and can be 479 accessed through GEO accession number GSE155028. The BINGO 3.03 plug-in38, 480 implemented in Cytoscape 2.8139, was used to determine and visualise the Gene 481 Ontology (GO)-enrichment categorization. A hypergeometric distribution statistical 482 testing method was applied to determine the enriched genes in combination with 483 the Benjamini and Hochberg multiple testing correction (FDR<0.05). 484 485 ChIP-seq data analysis 486 Libraries of three biological replicates of both the ChIP and Input samples were 487 sequenced after which Bowtie2 was used to map the reads to the Arabidopsis 488 genome 40. MACS2 was used for peak calling (the statistical detection of protein 489 binding sites in the DNA)41 using the Input samples as control. The ChIP-seq data 490 is made available via NCBI and can be accessed through GEO accession number 491 GSE155028. The function Distance2Genes, part of the R-package ‘CSAR’, was used 492 to determine the nearest gene42. These nearest genes are then presented in a list of 493 candidate target genes35 for further analysis. The analysis of read distribution of our 494 ChIP-seq Dataset was done by the R-package ChIPpeakAnno: 495 ‘assignChromosomeRegion’. As a cutoff 1kb of promoter region with an immediate 496 cutoff of 0.5 kb was used. All three biological replicates were initially analysed 497 separately. To confirm the high rate of reproducibility the overlap between the two 498 highest enriched samples and the overlap of all three samples was analysed. To 499 assess the reproducibility of the datasets, a pair-wise correlation analysis was done, 500 preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted December 15, 2020. ; https://doi.org/10.1101/2020.12.14.394403doi: bioRxiv preprint revealing Pearson correlation co-efficiencies ranging from 0,77 to 0,83 501 (Supplemental Fig. 3b) which supports the reproducibility and quality of our 502 experiments44. Because of the lower number of peaks in sample 3 (Supplemental 503 Fig. 3a), we decided to perform further analyses on the first two ChIP-seq datasets. 504 We approached these datasets by creating a file representing their overlap, which 505 was inferred by the distance between peaks in the individual samples. The average 506 peak size was ±300bp, so if the centre of a peak in one sample is no more than 507 150bp apart from the centre of a peak in the other sample, the peaks are assumed 508 to be in fact the same and were saved in a separate file for subsequent analyses. 509 This criterion resulted in a file containing roughly 60-70% of the peaks of sample 1 510 and 2. Since this merging of peaks may result in a slight shift of the actual binding 511 site, analyses looking specifically at the exact location of the binding site were done 512 on the individual datasets. 513 514 Analysis of peak-overlap between BRC1 and BTFs 515 We use as datasets the 4227 BRC1 peak locations and the peak locations and target 516 genes for 10 BTFs. We defined promoter regions as 3kb upstream region according 517 to TAIR. For each BTFs we performed the following analysis: we obtained the 518 promoter region of the BTF target genes; retrieve the peak locations of BTF 519 occurring in those promoter regions (at least one bp of the peak should be in the 520 promoter region); obtain the peak locations of BRC1 occurring in those promoter 521 regions; for each BRC1 peak occurring in those promoter regions, check if there is 522 a BTF overlapping peak (defined as having at least one bp in common). 523 524 preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted December 15, 2020. ; https://doi.org/10.1101/2020.12.14.394403doi: bioRxiv preprint Motif enrichment analysis 525 Motif discovery was carried out using MEME-ChIP45. This program consists of 526 several sub-programs that each perform a specific analysis. In the present study 527 MEME, DREME and CentriMO were used. MEME46 looks for overrepresented motifs 528 in a set of sequences compared to a background model of nucleotide frequencies. 529 The background model was generated using fasta-get-markov for the MEME suite, 530 using a second order Markov model, and a set of 207 randomly selected sequences 531 from Arabidopsis promoters as input. DREME15 looks for motifs occurring 532 significantly more often in the input set of sequences compared to a control set; as 533 control set, shuffled versions of the input set was used. CentriMO47 is used to test if 534 the motifs found are centrally enriched. In the present study, motifs were 535 investigated if they were in the top six motifs as defined by MEME-ChIP. Defaults 536 settings for DREME were used, including an E-value threshold of 0.05. Default 537 settings for MEME were used, with the following exceptions: the setting meme-mod- 538 anr was used, which assumes zero, one or multiple motif occurrences per 539 sequence. Minimum and maximum width for the motifs were set to 4- and 15bp, 540 respectively. 541 542 Network construction and visualization 543 The transcriptional network comprising the BRC1-dependent genes and the 544 potential targets of the BRC1-dependent TFs was constructed using the R package 545 igraph for network construction and analysis and ChIP-seq and DAP-seq data 546 available17,24. Nodes represent genes, and edges are drawn from genei to genej when 547 genei encodes a TF that binds to the promoter of genej. A network adjacency matrix 548 preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted December 15, 2020. ; https://doi.org/10.1101/2020.12.14.394403doi: bioRxiv preprint of order 947x947; number of rows (genes) x number of columns (genes); was 549 constructed by setting component (i,j) to 1 when genei encodes a TF that binds to 550 genej promoter. Otherwise it was set to 0. The function graph.adjacency with input 551 the previously described adjacency matrix and mode equal to directed was used to 552 generate our network comprising 947 nodes/genes and 10240 edges/potential 553 direct transcriptional regulations. A web app was developed using the R package 554 shiny to analyse this network (https://greennetwork.us.es/BRC1NET/) 555 556 Network motifs identification and analysis 557 A network motif is defined as a subgraph whose occurrence is significantly greater 558 in our network than its occurrence in random networks with the same topological 559 properties. The analysis was performed using the R package igraph. 1000 random 560 networks were generated with the same number of nodes as our transcriptional 561 network, namely 947 nodes. To capture the same topology as our network, 37 562 nodes were chosen randomly that would play the role of the TFs in our network. 563 Each of these nodes was connected to the same number of nodes as one of the TFs 564 in our network, but the targets were also chosen randomly. 565 566 Supplemental Information 567 Supplemental Information includes Supplemental Results, Supplemental figures 568 and Supplemental tables. 569 570 Acknowledgements 571 preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted December 15, 2020. ; https://doi.org/10.1101/2020.12.14.394403doi: bioRxiv preprint The research was supported by grants from the Dutch Scientific Organization 572 (NWO); (NWO-JSTP grant 833.13.008 and VENI 15060), Spanish Ministry of 573 Economy (MINECO) and fondos FEDER [grants BIO2014-57011-R, BIO2017-84363-R and 574 BIO2017-84066-R]. CT was a La Caixa fellow and Aitor Muñoz a FPI (MINECO) fellow. 575 We thank Juan Carlos Oliveros, Mónica Franch and F.J.M. van Workum for help with the 576 RNA- and ChIP-seq data analysis, Estíbaliz Bustos for amplification of BTF promoters and 577 Javier Paz-Ares and Desmond Bradley for constructive criticisms of the manuscript. 578 579 preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted December 15, 2020. ; https://doi.org/10.1101/2020.12.14.394403doi: bioRxiv preprint Genet. 8, 450–461 (2007). 716 27. Joanito, I., Chu, J. W., Wu, S. H. & Hsu, C. P. An incoherent feed-forward 717 loop switches the Arabidopsis clock rapidly between two hysteretic states. 718 Sci. Rep. 8, 1–16 (2018). 719 28. Mine, A. et al. An incoherent feed-forward loop mediates robustness and 720 tunability in a plant immune network. EMBO Rep. 18, 464–476 (2017). 721 29. Cutler, S. R., Rodriguez, P. L., Finkelstein, R. R. & Abrams, S. R. Abscisic acid: 722 Emergence of a core signaling network. Annual Review of Plant Biology 61, 723 (2010). 724 30. Fujita, Y., Fujita, M., Shinozaki, K. & Yamaguchi-Shinozaki, K. ABA- 725 mediated transcriptional regulation in response to osmotic stress in plants. J. 726 Plant Res. 124, 509–525 (2011). 727 31. Trigg, S. A. et al. CrY2H-seq: A massively multiplexed assay for deep- 728 coverage interactome mapping. Nat. Methods 14, 819–825 (2017). 729 32. Lumba, S. et al. A mesoscale abscisic acid hormone interactome reveals a 730 dynamic signaling landscape in Arabidopsis. Dev. Cell 29, 360–372 (2014). 731 33. Clough, S. J. & Bent, A. F. Floral dip: A simplified method for Agrobacterium- 732 mediated transformation of Arabidopsis thaliana. Plant J. 16, 735–743 (1998). 733 34. Chevalier, F., Iglesias, S. M. & Sánchez, Ó. J. Plastic Embedding of Stem 734 Sections. Bio-protocol 4, 1–7 (2014). 735 35. Mourik, H. Van, Muiño, J. M., Pajoro, A., Angenent, G. C. & Kaufmann, K. 736 in Methods in Molecular Biology 1284, 93–121 (2015). 737 36. Obayashi, T. et al. ATTED-II: A database of co-expressed genes and cis 738 elements for identifying co-regulated gene groups in Arabidopsis. Nucleic 739 preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted December 15, 2020. ; https://doi.org/10.1101/2020.12.14.394403doi: bioRxiv preprint Acids Res. 35, 4–6 (2007). 740 37. Trapnell, C. et al. Differential gene and transcript expression analysis of RNA- 741 seq experiments with TopHat and Cufflinks. Nat. Protoc. 7, 562–78 (2012). 742 38. Maere, S., Heymans, K. & Kuiper, M. BiNGO: A Cytoscape plugin to assess 743 overrepresentation of Gene Ontology categories in Biological Networks. 744 Bioinformatics 21, 3448–3449 (2005). 745 39. Shannon, P. et al. Cytoscape: a software environment for integrated models 746 of biomolecular interaction networks. Genome Res. 2498–2504 (2003). 747 doi:10.1101/gr.1239303.metabolite 748 40. Langmead, B. & Salzberg, S. L. Fast gapped-read alignment with Bowtie 2. 749 Nat. Methods 9, 357–359 (2012). 750 41. Zhang, Y. et al. Model-based analysis of ChIP-Seq (MACS). Genome Biol. 9, 751 (2008). 752 42. Muiño, J. M., Kaufmann, K., van Ham, R. C. H. J., Angenent, G. C. & 753 Krajewski, P. ChIP-seq Analysis in R (CSAR): An R package for the statistical 754 detection of protein-bound genomic regions. Plant Methods 7, 11 (2011). 755 43. Kaufmann, K. et al. Target genes of the MADS transcription factor 756 SEPALLATA3: integration of developmental and hormonal pathways in the 757 Arabidopsis flower. PLoS Biol. 7, (2009). 758 44. Bardet, A. F., He, Q., Zeitlinger, J. & Stark, A. A computational pipeline for 759 comparative ChIP-seq analyses. Nat. Protoc. 7, 45–61 (2012). 760 45. Machanick, P. & Bailey, T. L. MEME-ChIP: Motif analysis of large DNA 761 datasets. Bioinformatics 27, 1696–1697 (2011). 762 46. Bailey, T. L. & Elkan, C. Fitting a Mixture Model by Expectation Maximization 763 preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted December 15, 2020. ; https://doi.org/10.1101/2020.12.14.394403doi: bioRxiv preprint to Discover Motifs in Bipolymers. Proc. Second Int. Conf. Intell. Syst. Mol. 764 Biol. 28–36 (1994). doi:citeulike-article-id:878292 765 47. Bailey, T. L. & MacHanick, P. Inferring direct DNA binding from ChIP-seq. 766 Nucleic Acids Res. 40, 1–10 (2012). 767 768 preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted December 15, 2020. ; https://doi.org/10.1101/2020.12.14.394403doi: bioRxiv preprint preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted December 15, 2020. ; https://doi.org/10.1101/2020.12.14.394403doi: bioRxiv preprint preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted December 15, 2020. ; https://doi.org/10.1101/2020.12.14.394403doi: bioRxiv preprint AGI Name Description TF family At4g36740 ATHB-40 Homeobox-leucine zipper protein ATHB-40 HD-zip At5g66700 ATHB-53 Homeobox-leucine zipper protein ATHB-53 HD-zip At2g18550 ATHB-21 Homeobox-leucine zipper protein ATHB-21 HD-zip At4g01120 GBF2 G-box-binding factor 2 bZIP At2g46270 GBF3 G-box-binding factor 3 bZIP At4g34000 ABF3 ABA Responsive Element Binding Factor 3 bZIP At2g36270 ABI5 ABA Insensitive 5 bZIP At1g77450 ANAC32 NAC domain containing protein 32 NAC At1g01720 ATAF1/ANAC002 ATAF1/ANAC domain-containing protein 2 NAC At1g70000 MYBD MYB-like transcription factor MYB At2g40970 MYBC1 Homeodomain-like superfamily protein Homeodomain At5g13330 ERF113/Rap2.6L Ethylene-responsive transcription factor ERF113 ERF/AP2 At3g61630 CRF6 Ethylene-responsive transcription factor CRF6 ERF/AP2 At5g49700 AHL17 AT-hook Motif Nuclear Localized Protein 17 AT-hook At5g44260 TZF5 Tandem CCCH Zinc Finger protein 5 CCCH Zinc Finger AT5G62020 HSFB2A/HSF6 Heat Stress transcription factor B-2A HSF At1g19050 ARR7 Response Regulator 7 ARR At4g31800 WRKY18 WRKY transcription factor 18 WRKY At1g07090 LSH6 Light-dependent Short Hypocotyl 6 ALOG Table 1. BRC1-regulated Transcription factors (BTFs). List of bona fide BRC1 direct targets encoding TFs. In bold, those related to ABA signaling and/or response. preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted December 15, 2020. ; https://doi.org/10.1101/2020.12.14.394403doi: bioRxiv preprint preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted December 15, 2020. ; https://doi.org/10.1101/2020.12.14.394403doi: bioRxiv preprint preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted December 15, 2020. ; https://doi.org/10.1101/2020.12.14.394403doi: bioRxiv preprint preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted December 15, 2020. ; https://doi.org/10.1101/2020.12.14.394403doi: bioRxiv preprint preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted December 15, 2020. ; https://doi.org/10.1101/2020.12.14.394403doi: bioRxiv preprint