An Arabidopsis Gene Regulatory Network for Secondary Cell Wall Synthesis

Summary The plant cell wall is an important factor for determining cell shape, function and response to the environment. Secondary cell walls, such as those found in xylem, are composed of cellulose, hemicelluloses and lignin and account for the bulk of plant biomass. The coordination between transcriptional regulation of synthesis for each polymer is complex and vital to cell function. A regulatory hierarchy of developmental switches has been proposed, although the full complement of regulators remains unknown. Here, we present a protein-DNA network between Arabidopsis transcription factors and secondary cell wall metabolic genes with gene expression regulated by a series of feed-forward loops. This model allowed us to develop and validate new hypotheses about secondary wall gene regulation under abiotic stress. Distinct stresses are able to perturb targeted genes to potentially promote functional adaptation. These interactions will serve as a foundation for understanding the regulation of a complex, integral plant component.

Plant cell shape and function are in large part determined by the cell wall. Almost all cells have a primary wall surrounding the plasma membrane. Specialized cell types differentiate by depositing a secondary cell wall upon cessation of cell elongation. In addition to providing mechanical support for water transport and a barrier against invading pathogens, the polymers contained within the wall are an important renewable resource for humans as dietary fibre, as raw material for paper and pulp manufacturing, and as a potential feedstock for biofuel production. Secondary cell walls account for the bulk of renewable plant biomass available globally.
The secondary cell wall consists of three types of polymer-cellulose, hemicelluloses and lignin-and is found in xylem, fibres and anther cells. Cellulose microfibrils form a main load-bearing network. Hemicelluloses include xylans, glucans, and mannans. Lignin is a complex phenylpropanoid polymer that imparts 'waterproofing' capacity as well as mechanical strength, rigidity and environmental protection. Despite the importance of the plant secondary cell wall, our knowledge of the precise regulatory mechanisms that give rise to these metabolites is limited. The expression of cell wall-associated genes is tightly spatiotemporally co-regulated 1,2 . However, the pervasive functional redundancy within transcription factor families, the combinatorial complexity of regulation, and activity in a small number of cell types render functional characterization from single gene experiments difficult. A model of master regulators has been proposed with NAC domain and homeobox HD-ZIP Class III (HD-ZIPIII) transcription factors initiating cell specification and secondary cell wall synthesis in Arabidopsis thaliana. In this model, VASCULAR-RELATED NAC DOMAIN6 (VND6) and VND7 are sufficient but not necessary to regulate xylem vessel formation; additionally, the HD-ZIPIII transcription factor PHABULOSA (PHB) also regulates vessel formation, and acts in a highly redundant manner with four other HD-ZIPIII factors 3 . In anthers, two NAC domain transcription factors, NAC SECONDARY WALL THICKENING1 (NST1) and NST2, are sufficient to drive the secondary cell wall biosynthetic program, but act redundantly 4 . Thus, regulation of this process is highly redundant and combinatorial. However, no comprehensive map of interactions has been developed at cell-type-resolution over time, nor have upstream regulators been identified. We therefore chose to pursue a network-based approach to comprehensively characterize the transcriptional regulation of secondary cell wall biosynthesis.

Mapping the secondary cell wall synthesis regulatory network
To systematically map this regulatory network at cell-type-resolution, we used a combination of high-spatial-resolution gene expression data 5 and the literature 1,6 to identify fifty genes implicated in xylem cell specification. These included transcription factors and enzymes involved in cellulose, hemicellulose and lignin biosynthesis that are expressed in root xylem cells (Supplementary Table 1; Methods). Selection of both developmental regulators and downstream functional genes allowed us to interrogate upstream regulatory events that determine xylem specification and differentiation associated with secondary cell wall synthesis. Promoter sequences were screened using an enhanced yeast one hybrid (Y1H) assay against 467 (89%) of root-xylem-expressed transcription factors 7 . Protein interactions were identified for 45 of the promoters (Supplementary Table 2). The final network comprises 242 genes and 617 protein-DNA interactions ( Fig. 1a; http://gturco.github.io/ trenzalore/stress_network). Thirteen of the transcription factors have been previously identified as having a role in xylem development or secondary cell wall biosynthesis. Six of the transcription factors were previously shown to bind to these promoters and a further nine of the protein-DNA interactions were implied in gene expression studies, that is, without demonstrating direct binding 6,[8][9][10][11] . These interactions represent independent validation of our approach (Supplementary Table 2; Extended Data Fig. 1). All together, the network contains 601 novel interactions, although false negatives and false positives are a component of all network approaches 12 .
Our Y1H approach revealed a highly interconnected regulatory network. On average, each cell wall gene promoter was bound by 5 transcription factors from 35 protein families with over-representation of AP2-EREBP, bHLH, C2H2, C2C2-GATA and GRAS gene families (Supplementary Table 3). Our network now adds an additional layer of gene regulation with novel factors upstream of VND6 and VND7 and supports feed-forward loops 9,11,13 as an overarching theme for regulation of this developmental process with a total of 96 such loops (Fig. 1a, b).
To organize the network, we employed a power graph compression approach to condense the network into overlapping node sets with similar connectivity. Protein-DNA interactions (edges) between proteins and promoters (nodes) in the original network were replaced by 'power edges' between overlapping 'power nodes' 14 . A power edge exists between suites of transcription factors that bind to the same set of promoters. Using this approach, 24 power edges were observed (Supplementary Table 4; Fig. 1c). Some sets could be distinguished on the basis of target gene function. For instance, one power edge connects 16 transcription factors with promoters of two lignin genes, 4CL1 and HCT, while another power edge connects three transcription factors with genes related to cellulose and hemicellulose biosynthesis such as CESA4, CESA7, IRX9, COBL4 and GUX2.

Testing interactions predicted by the network
Using our network, we hypothesized that E2Fc is a key upstream regulator of VND6, VND7 and secondary cell wall biosynthesis genes. This hypothesis is based on our finding that E2Fc bound to 23 promoters including those of VND6, VND7 and MYB46, and cellulose-, hemicelluloseand lignin-associated genes (Fig. 2a). VND7 and MYB46 are also known to bind to the promoters of many of these genes as well 9,13,15 , creating a suite of feed-forward loops. E2Fc is a known negative regulator of endoreduplication 16,17 . Before terminally differentiating, xylem cells elongate and likely undergo endoreduplication before secondary cell wall deposition. E2Fc can act as a transcriptional repressor [16][17][18] as well as a transcriptional activator [19][20][21][22] and here we report both. E2Fc activated VND7 expression in a dose-dependent manner ( Fig. 2b and Extended Data Fig. 2a, b)

RESEARCH ARTICLE
low E2Fc levels resulted in VND7 repression, whereas moderate E2Fc levels resulted in VND7 activation (Extended Data Fig. 2b). This dynamic regulation was also observed in an E2Fc-knockdown line 23 , where transcript abundance of VND6 and VND7 was significantly increased (Fig. 2c). Based on our results, we propose that E2Fc acts in a complex, concentration-dependent manner to regulate gene expression either as an activator or a repressor. Coincident with the repression observed in E2Fc-knockdown lines, ectopic patches of lignin were observed near the root-shoot junction using phloroglucinol staining (Fig. 2d). A significant increase in crystalline cellulose in the knockdown line was observed using an Updegraff assay (Fig. 2e). All five HD-ZIPIII transcription factors, including REVOLUTA (REV), PHB, and PHAVOLUTA are jointly necessary for xylem cell specification and secondary wall synthesis 3 . We found that VND7 bound REV and PHB promoters in yeast. VND7 has been to shown to act as a transcriptional activator 9 or as a repressor when complexed with VNI2 24 . With a dexamethasone-inducible version of VND7 25 , transcript levels of REV and PHB were significantly decreased by 2.5-fold following induction (Fig. 3a). The REV transcription factor bound to the promoter of the lignin biosynthesis gene PHENYLALANINE AMMONIA LYASE4 (PAL4). In a rev-5 loss-of-function mutant, PAL4 significantly increased in transcript abundance (Fig. 3b) and transient induction of REV by a glucocorticoid receptor fusion 26 resulted in a decrease of PAL4 expression (Fig. 3c). Taken together, these data suggest that E2Fc can activate VND7 expression in a dose-dependent manner, while VND7, possibly in concert with VNI2, can repress REV expression, and REV can repress expression of PAL4. This series of interactions predicted by the network model and tested by perturbation analyses ensures that activation of VND7 and coordination of lignin biosynthesis is tightly regulated.
We next sought to identify all transcription factors that potentially regulate secondary cell wall biosynthesis genes, not just in root xylem cells but also in above-ground cell types including xylary fibres, interfascicular fibres and anthers. Many of the biosynthetic genes downstream of the key NAC domain transcription factors act in both the root and the shoot 9 . To expand the network, we used Y1H to screen multiple smaller promoter fragments of a subset of promoters included in the root xylem network, including genes associated with cellulose, hemicellulose and lignin biosynthesis against a library of 1,664 full-length Arabidopsis transcription factors (Supplementary Tables 5, 6). We observed a total of 413 interactions that included proteins from 36 of the 75 protein families tested (Supplementary Table 7; Fig. 1d; http://gturco. github.io/trenzalore/secondary_cell_wall). We found an over-representation of AP2-EREBP, bZip, ZF-HD, MYB and GeBP families (Supplementary Table 8). Each promoter interacted with an average of 38 different proteins, generating even more possibilities for combinatorial, redundant or condition-specific gene regulation. Like the root xylem network, previously reported protein-DNA interactions were observed in this screen including MYB46 and MYB83 binding the promoters of CESA genes (Supplementary Table 7) 8,27 . Since most of these interactions were novel, a subset was additionally validated. Transient expression of AIL1, MYB83, MYB54, NAC92, NST2 and SND1 caused a significant increase in CESA4::LUC activity in tobacco (measured by luciferase activity), indicating binding and activation of the CESA4 promoter (Fig. 4a). We further tested three regions of the CESA4 promoter with two NAC family proteins, SND1 and NST2 (Fig. 4b, c), using an in vitro electrophoretic mobility shift assay (EMSA). Extracts of Escherichia coli expressing either glutathione-S-transferase-conjugated NST2 (GST:NST2) or GST:SND1 in the presence of a CESA4-2pr promoter probe produced DNA species with retarded mobility (Fig. 4b, c). We also observed binding of the CESA7, CESA8 and KOR promoter fragments with the NST2 protein and CESA8 with the SND1 protein (Extended Data Fig. 3). These interactions between NST2 and CESA4, CESA8, and KOR promoters were further confirmed in planta by chromatin immunoprecipitation (ChIP). An antibody to green fluorescent protein (GFP) was used to immunoprecipitate NST2 protein from extracts of 35S::NST2::GFP plants. The complex was significantly enriched for fragments from the CESA4, CESA8 and KOR promoters (Fig. 4d). The tracheary element-regulating ciselement (TERE, CTTNAAAGCNA) is a direct target of VND6 28,29 . A perfect TERE is present in the CESA4 promoter (CTTGAAAGCTA) and TERE-like sequences are present in CESA8 (CTTCAATGTTA) and KOR (CTTGAAAATGA). Taken together, these data clearly demonstrate that the expression of CESA4 and other secondary cell wall genes is mediated by the direct binding of the NAC-domain binding transcription factors NST2 and SND1 to the target gene promoters via the TERE.

Abiotic stress can co-opt the xylem regulatory network
Having generated a gene regulatory network supported by in vivo and in vitro approaches, we sought to test if the model could allow us to predict responses under abiotic stress perturbation. Co-opting a developmental regulatory network is likely a key mechanism to facilitate adaptation in response to stress. Thus, we hypothesized that stress responses are likely integrated into the gene regulatory network that determines xylem cell specification and differentiation and that we can predict the exact genes that these stresses manipulate within our network.
We first identified genes within the network whose expression was altered specifically in the root vasculature in response to salt, sulphur, iron and pH stress 30,31 and nitrogen influx 32 . Genes within the root xylem secondary cell wall network were significantly differentially regulated in response to sulphur stress, salt stress and iron deprivation (Supplementary Table 9). Substantial overlap was observed between iron AT5G15710 control in rev-5 relative to wild-type. c, PAL4 expression relative to UBC21 control following one hour dexamethasone (Dex) treatment of 35S:REV:GR relative to untreated. *P , 0.05 for panels b and c, n 5 2 biological replicates with 3 technical replicates. All panels show data as means 6 s.d., with P calculated from Student's t-test.

ARTICLE RESEARCH
deprivation and salt stress gene responses and was further characterized (Fig. 5a). We filtered the xylem network to include only genes differentially expressed in salt or iron, creating stress-specific sub-networks (Extended Data Fig. 4). Previously, we determined that key developmental transcription factors have significantly more upstream regulators compared to other genes 33 . In response to iron deprivation, REV has the most upstream regulators, while in response to salt stress, VND7 and MYB46 have the most upstream regulators.
On the basis of these data from the iron-deprivation sub-network, we hypothesized that REV plays a key role in regulating secondary cell wall development in response to iron deprivation. To additionally determine directionality and sign (activation or repression) in the network, we constructed a network of 16 key nodes using the consensus network from four unsupervised and one supervised network inference method. REV was also predicted to be an important regulator of lignin biosynthesis gene expression in response to iron deprivation using these methods (Extended Data Fig. 5). First, to test the modelgenerated prediction that lignin biosynthetic gene expression is altered in response to iron deprivation, we measured phenylpropanoid-related gene expression. An increase in 4CL1, PAL4 and HCT gene expression was observed (Fig. 5b). Additionally, iron deprivation stress altered the timing and spatial distribution of the 4CL1 transcript (Fig. 5c). These expression changes are accompanied by an increase in fuchsin staining, indicative of increased phenylpropanoid deposition (Extended Data Fig. 6b). Expression in a rev-5 loss-of-function mutant in irondeficient conditions revealed a REV-and stress-dependent influence on CCoAOMT1, PAL4 and HCT expression (Fig. 5d), thus validating our model predictions.
In the high-salinity sub-network VND7 and MYB46 contain the most upstream regulators (Extended Data Fig. 4). VND7 and MYB46 expression is greatly increased in roots in response to salt stress, but lignin biosynthetic gene expression is unaltered (Fig. 5e; Extended Data Fig. 6a). In corroboration with this hypothesis, the network model constructed using the described in silico methods also predicts VND7 and MYB46 as main regulators in response to salt stress but not iron deprivation (Extended Data Fig. 7), and indeed this was observed with an expansion of the domain of VND7 expression after salt treatment but not iron deprivation (Fig. 5e, f; Extended Data Fig. 6c). In conjunction with this ectopic increase, we observed an additional strand of metaxylem in roots exposed to high salinity (Fig. 5g).

Discussion
Owing to functional redundancy among regulators of secondary cell wall biosynthesis, transcription factors have largely eluded identification by loss-of-function genetic screens. Our network approach has identified hundreds of novel regulators and provided considerable insight into the developmental regulation of xylem cell differentiation. The network, which includes a cell cycle regulator, is comprised of many feed-forward loops that are likely to ensure robust regulation of this process (Fig. 5h). Accordingly, we revealed that perturbation at distinct

RESEARCH ARTICLE
nodes changes the network subtly, including phenylpropanoid biosynthesis in response to iron deprivation and ectopic xylem cell differentiation in response to salt stress (Fig. 5h). We anticipate that these findings will be instrumental in biotechnology and in our understanding of cell fate acquisition.
Online Content Methods, along with any additional Extended Data display items and Source Data, are available in the online version of the paper; references unique to these sections appear only in the online paper.

METHODS
Yeast one-hybrid (Y1H) protein-DNA interaction assays. The root vascularexpressed transcription factor collection is described in ref.  Table 5) were arrayed in 96-well plates and plasmids were prepared using the Promega Wizard SV 96 plasmid purification DNA system according to manufacturer's recommendations. Root secondary cell wall gene promoters (2-3 kb of upstream regulatory region from the gene's translational start site, or the next gene, whichever comes first) were cloned and recombined with reporter genes according to ref. 33. Promoter sequences and primers used are described in Supplementary Table 1 Table 6.
Root bait promoters were screened against the stele-expressed transcription factor collection using the Y1H protocol as previously described 7 . The 1,663 transcription factor library was transformed into each yeast strain and the b-galactosidase activity was determined as described in ref. 38, but in 384-well plates. Positive interactions were visually identified as incidence of yellow caused by the presence of ortho-nitrophenyl cleavage from colourless ortho-nitrophenyl-b-D-galactoside by b-galactosidase. The DNA bait strains were similarly tested for self-activation before screening by not transforming with prey vectors in the presence of thiamine. All interacting transcription factors were assembled into a cell wall interaction library and the screen was repeated to confirm the results and each clone was sequenced to reconfirm identity. Statistical analysis for protein family enrichment. Enrichment was determined using the hypergeometric distribution online tool (http://stattrek.com/). The population size is the number of transcription factors in the xylem transcription factor collection while the successes within the population is the number of transcription factors within that transcription factor family in the xylem. The number of successes in the sample was the number of proteins belonging to that family, and the number in the sample is the total number of transcription factors within the network. The A. thaliana transcription factor list is as described in ref. 7. Power graph compression approach. The power graph compression was performed using the algorithm as previously described 14 . Plant material. The E2Fc RNAi line is described in ref. 23 and was verified by quantifying E2Fc transcript abundance relative to the Col-0 control using an E2Fc primer compared to an ACTIN control primer (Supplementary Table 1). VND7::YFP lines are described in ref. 39. The VND7 glucocorticoid induction line is described in ref. 9. The rev-5 loss-of-function mutant was described in ref. 40. Cloning and insertion of the 4CL1 promoter into a pENTR p4-p1R donor vector was performed according to ref. 33 (for sequence, see Supplementary Table 1). The promoter was then recombined into binary vector pK7m24GW,3 along with pENTR 221 ER-GFP::NOS. The resulting 4CL1::GFP vector was transformed into Agrobacterium strain GB3101. Col-0 plants were then transformed using the floral dip method. Plant growth conditions. All plants were grown vertically on plates containing 13 Murashige and Skoog salt mixture, 1% sucrose, and 2.3 mM 2-(N-morpholino) ethanesulphonic acid (pH 5.8) in 1% agar. NaCl plates were made by adding 140 mM NaCl to this standard media. Iron control and deprivation media were made according to ref. 30. Plants grown on stress media (iron or salt) were first germinated on nylon mesh placed over control media for four days before transferring mesh with seedlings to iron deprivation or NaCl plates. Plants used for RNA isolation were also grown on nylon mesh placed over the agar to facilitate the collection of root material 5 . Determination of crystalline cellulose. Roots of 7-day-old plants were harvested and lyophilized. Six to ten plates of seedlings grown at the same time on the same media were pooled to make a single biological replicate. Crystalline cellulose was measured according to ref. 41. After hydrolysis of non-cellulosic polysaccharides from an alcohol insoluble residue wall preparation with the Updegraff reagent (acetic acid:nitric acids:water, 8:1:2 v/v), the remaining pellet was hydrolysed in 72% sulfuric acid. The resulting glucose quantity was determined by the anthrone method 42 . Phloroglucinol staining. Five day after imbibition seedlings to be stained with phloroglucinol were fixed in a 3:1 95% ethanol:glacial acetic acid solution for 5 min. Samples were then transferred to a solution of 1% phloroglucinol in 50% HCl for 1-2 min. Whole seedlings were then mounted in 50% glycerol on slides and viewed using an Olympus Vanox microscope. Images were captured with a PIXERA Pro-600ES camera. Confocal laser scanning microscopy. Confocal laser scanning microscopy was carried out on a Zeiss LSM700. Cell walls were stained using propidium iodide as previously described 30 . Transient protein-DNA interaction detection in tobacco. b-glucuronidase. For transient transactivation expression assays, the VND7, GAL4, and/or CyclinB1 promoters were cloned into pGWB3 to generate GUS (b-glucuronidase gene) fusion reporters for E2Fc transcriptional activity. The E2Fc effector vector 43 (in PYL436) was provided by S. D. Kumar (UC Davis, CA). The effector and reporter constructs were transformed into Agrobacterium tumefaciens strain GV3101 and co-infiltrated with the p19 silencing inhibitor into 3-weeks-old Nicotiana benthamiana leaves at A 600 nm 0.6:0:6:1, respectively. Leaves were harvested 3 days after agro-infiltration and homogenized in GUS extraction buffer (50 mM Na 2 PO 4 pH 7, 10 mM Na2-EDTA, 0.1% SDS, 0.1% Triton TX-100 and 10 mM b-mercaptoethanol). Quantitative MUG fluorescent assay for GUS determination was performed using 100 mg of protein/sample in 500 ml of GUS assay buffer (1 mM 4-methyl umbelliferyl b-Dglucuronide, SIGMA, in Extraction Buffer). Samples were covered in aluminium foil and incubated at 37 uC. Reaction was stopped at different time points by transferring 50 ml to a tube with 450 ml of Stop Buffer (0.2 M Na 2 CO3). 4-methylumbelliferone fluorescence was determined using a Infinite 200 Pro-series reader (excitation at 365 nm, emission at 455 nm). Luciferase (Fig. 2). Overnight cultures of Agrobacterium (GV3101, D 600 nm 5 0.6) carrying VND7 promoter fused to luciferase (LUC) and 35S::E2Fc were prepared in infiltration medium (2 mM Na 3 PO 4 , 50 mM MES, 0.5% glucose, 100 mM acetosyringone) at D 600 nm 5 0.1. Subsequently, cultures containing VND7::LUC and 35S::E2Fc at respective ratios of 1:0, 1:0.5, 1:1, 1:2, 1:5, and 1:10 were spot-infiltrated into 6-7-week-old Nicotiana benthamiana leaves. To prevent gene silencing, Agrobacterium strain carrying the pBIN19 suppressor from tomato bushy stunt virus was included in each of the combinations 44 . The LUC activity was inspected at 72 to 96 h post infiltration using CCD camera (Andor Technology).
Luciferase imaging of VND7::LUC was performed as previously described with modifications 45 . Briefly, tobacco leaves were cut off after 3 days of transient transformation and sprayed with 1 mM luciferin (Promega) in 0.01% Tween-80, then were imaged using an Andor DU434-BV CCD camera (Andor Technology). Images were acquired every 10 min for 12 pictures. Luciferase activity was quantified for a defined area as mean counts per pixel per exposure time using Andor Solis image analysis software (Andor Technology). Statistical analyses were performed using two-tailed Student's t-tests. The difference was considered significant if P , 0.05. Luciferase (Fig. 4). A vector system was created to generate a single vector with the CaMV 35S constitutive promoter (35S) fused to a transcription factor, a promoter fragment fused to the firefly luciferase reporter gene, and 35S fused to the Renilla luciferase reporter gene. The constitutively expressed Renilla gene served as a control to normalize for transformation efficiency. This system includes one destination vector pLAH-LARm and three entry vectors pLAH-TF, pLAH-PROM and pLAH-VP6435T using MultiSite Gateway Pro Technology (Invitrogen) to simultaneously clone three DNA fragments (Extended Data Fig. 8). To develop the expression vector, promoter fragments and transcription factors were cloned, using the BP system (Invitrogen), into pDONR-P3-P2 and pDONR-P1-P4 to create pLAH-TF and pLAH-Prom, respectively. PacI-digested pMDC32 was ligated with the 2.427 kb pFLASH fragment following HindIII and SacI digestion to yield pLAR-L with the firefly luciferase (LUC) reporter gene. The 3 kb pRTL2-Renilla HindIIIdigested fragment was inserted into SacI-digested pLAH-L to create pLAR-LR with both firefly LUC and Renilla luciferase (REN) genes. To generate pLAH-LAR, a SpeI-digested PCR fragment containing the AmpR gene amplified from pDEST22 was ligated with SpeI-digested pLAR-LR. To add the minimal CaMV 35S fragment (Mini35S) before the LUC reporter gene, the gateway cassette ccdB/CmR of pLAR-LAR was replaced by a HindIII-digested PCR fragment Mini35S-ccdB-CmR amplified from pMDC32 using specific primer pHindIII-Rv and primer Mini35S-attR2. The final destination vector is referred to as pLAH-LARm.
The protein coding regions of select transcription factor genes were amplified. Each amplified fragment was recombined with pDONR-P1-P4 vector by performing BP reactions to produce pLAH-TF. Target promoter fragments were amplified from A. thaliana genomic DNA using appropriate primers with attB3 and attB2 sites (Supplementary Table 10). Each amplified fragment was cloned into pDONR-P3-P2 vector by performing BP reactions to produce pLAH-PROM. A third pDONR RESEARCH ARTICLE vector (pLAH-VP64Ter) was designed to create a carboxy-terminal fusion of the strong transcription activation domain VP64 to the transcription factor followed by the 35S transcription terminator (35St). A PCR fragments containing VP64 region and 35S terminator was amplified from pB7-VP64 using specific primers with attB4r and attB3r sites (Supplementary Table 10) into pDONR P4r-P3r to produce pLAH-VP6435T. Finally, the fully functional expression vector was generated by Gateway LR cloning of destination vector and the three entry clones: pLAH-LARm, pLAH-TF, and pLAH-VP64Ter (Extended Data Fig. 7).
Agrobacterium tumefaciens strain GV3103 (MP90) carrying expression constructs were grown in Luria-Bertani media with rifampicin and ampicillin and suspended in infiltration buffer 10 mM MES, pH 5.7, containing 10 mM MgCl 2 and 150 mM acetosyringone. The cultures were adjusted to a D 600 nm of 0.8 and incubated at room temperature for at least 3 h before infiltration. The cultures were hand infiltrated using a 1 ml syringe into 3-to 4-week-old N. benthamiana leaves. Leaf samples were harvested 36 h after infiltration and assayed for luciferase activity according to manufacturer instructions using the Dual-Luciferase Reporter Assay Systems (Promega). Approximately 100 mg of tissue was frozen in liquid nitrogen and homogenized using a Retsch Mixer Mill MM400 for 1 min at 30 Hz. Ground tissue was then thawed in lysis buffer (0.1 M HEPES, pH 7.8, 1% Triton X-100, 1 mM CaCl2 and 1 mM MgCl 2 ) at 25 uC for 15 min. Then 50 ml of Luciferase Assay Reagent II was added to 10 ml aliquots of the lysates to measure firefly luciferase activity, 1,000 ms integration time, using a Spectra Max M5/M5e plate reader to measure total light emission. Firefly luciferase activity was quenched with 50 ml of Stop & Glo Reagent, which contains Renilla luciferin substrate, also measured, 100 ms integration time, as total light emission. An expression vector containing part of the coding sequence (1X/1Y) of the b-glucuronidase reporter gene rather than a transcription factor gene was used for baseline measurement of firefly luciferase activity. To estimate relative transcription factor affinity with each promoter fragment, three biological replicates of transcription factor expressing vectors were compared to the average results for the GUS expression vector. First, dividing firefly luciferase activity by Renilla luciferase activity normalized the transformation efficiency of each infiltrated leaf sample. Relative binding of the transcription factors to the promoter bait sequences was determined relative to the GUS control using a Student's t-test in R v2.11.0. Electrophoretic mobility shift assays. To express recombinant NST2 or SND1 protein, coding sequence was cloned and fused to glutathione S-transferase tag in the pDONR211 vector and then transferred into pDEST15 (Invitrogen). E. coli strain BL21-AI (Invitrogen) transformed with pDEST15-GST:NST2 were grown in liquid media to a D 600 nm of 0.4, treated with 0.2% L-arabinose to induce expression overnight and harvested by centrifugation the following day. Cells were treated with 1 mg ml 21 lysozyme on ice for 30 min in minimal volume of 13 PBS buffer and lysed by sonication. Cell lysates were clarified by centrifugation and incubated with 100 ml of glutathione Sepharose beads (GE Healthcare) for 30 min at 4 uC with rotation. The beads were transferred to a column, washed with 10 volumes of 13 PBS. Protein was eluted in 100 mM Tris-HCl pH 8.0, 100 mM NaCl and 3 mg ml 21 glutathione buffer and purified protein was resuspended in 50% glycerol and stored at 280 uC.
Three overlapping probes were generated for CESA7, CESA8 and KOR promoters using the same oligonucleotides described in Supplementary Table 1, whereas three probes were generated for CESA4 using the following primers: CESA4pr-1fwd, CACCGGGCCTTTGTGAAATTGATTTTGGGC; CESA4pr-1rev, TGTA TTTCTACTTTAGTCTTAC; CESA4pr-2fwd, CCAGATTTGGTAAAGTTTAT AAG; CESA4pr-2rev, GTGTCATAAGAAAGCTTCAAG; CESA4pr-3fwd, TCTT ATGACACAAACCTTAGAC; CESA4pr-3rev, ACACTGAGCTCTCGGAAGC AGAGCAG. Reactions were carried out in binding buffer (10 mM Tris, pH 7.5, 50 mM KCl, 1 mM DTT, 2.5% glycerol, 5 mM MgCl 2 , 0.1% IGEPAL CA-630, and 0.05 mg ml 21 calf thymus DNA). Following the addition of 150 ng of protein from the GST purification eluate, reactions were incubated at room temperature for 30 min. Protein-DNA complexes were separated from the free DNA on 1% agarose/13 TAE gels at 4 uC. The agarose gels were stained with ethidium bromide and bands visualized under ultraviolet light. For the titration of promoter DNA with NST2 protein, CESA4 promoter fragment-2 DNA and KOR promoter fragment-1 DNA in 30 ng were titrated with increasing amounts of NST2 protein: 25, 50, 150, 300, and 600 ng. Binding reaction and the separation of protein-DNA complexes were carried out as described above. Chromatin immunoprecipitation of NST2. Chromatin immunoprecipitation was conducted as described in ref. 46 with the following modifications. Roughly 5 g (fresh weight) whole stems from six-week-old Arabidopsis were harvested and crosslinked for 15 min under vacuum in crosslinking buffer (10 mM Tris, pH 8.0, 1 mM EDTA, 250 mM sucrose, 1 mM PMSF and 1% formaldehyde). Technical replicates containing approximately 1.5 mg DNA were resuspended in 800 ml SII buffer, incubated with 2 mg anti-GFP antibody (ab290, Abcam) bound to Protein G Dynabeads (Invitrogen) for 1.5 h at 4 uC and then washed five times with SII buffer. Chromatin was eluted from the beads twice at 65 uC with Stop buffer (20 mM Tris-HCl, pH 8.0, 100 mM NaCl, 20 mM EDTA and 1% SDS). RNase-and DNasefree glycogen (2 mg) (Boehringer Mannheim) was added to the input and eluted chromatin before they were incubated with DNase-and RNase-free proteinase K (Invitrogen) at 65 uC overnight and then treated with 2 mg RNase A (Qiagen) for 1 h at 37 uC. DNA was purified by using Qiagen PCR Purification kit and resuspended in 100 ml H 2 0. Quantitative PCR reactions of the technical replicates were performed using Quantifast SYBR Green PCR Kit (Qiagen), with the following PCR conditions: 2 min at 95 uC, followed by 40 cycles of 15 s at 95 uC, 15 s at 55 uC and 20 s at 68 uC. Primers used in this study are listed in Supplementary Table 4. Results were normalized to the input DNA, using the following equation: 100 3 2 (Ct input 2 3.32 2 Ct ChIP) . Quantitative RT-PCR. Primers for qRT-PCR were designed to amplify a 100 bp region (or a 400 bp region for REV, PHB, and PHV transcripts due to sequence similarity) on the 39 end of each transcript 33 . Primer sets used for qRT-PCR are listed in Supplementary Table 1. Each plate was considered a biological replicate and Columbian and reference genotypes were plated on the same plate. Five days after imbibition, total RNA was extracted from seedling roots using an RNeasy Kit (QIAGEN). cDNA was synthesized by treatment with reverse transcriptase and oligo(dT) primer (SuperScript III First-Strand Synthesis System; Invitrogen). qRT-PCR was performed in an iCycler iQ Real-Time PCR Detection System (Bio-Rad) using the Bio-rad iQ SYBR green Supermix. Gene expression was measured between wild-type and mutant pairs across at least two biological replicates with three technical replicates using the D 2 DC T method 30 . VND7 induction experiments. VND7-VP16-GR 9 plants were grown vertically on sterile mesh placed on top of MS media with sucrose. Five days after imbibition, seedlings were transferred, with the mesh, to MS media containing 10 mM dexamethasone and roots were collected for qRT-PCR (RNeasy Kit; Qiagen) after 0, 1, 2, 3, or 4 h on dexamethasone (n 5 3). As a positive control, upregulation of MYB46 expression was confirmed using qRT-PCR. Nitrogen influx, salt stress, iron deprivation, sulphur stress, pH stress analysis. The data sets used contained mean expression values for each gene in both control and treatment, and a q value for each gene indicating the significance of the hypothesis that the expression values of control and treatment are drawn from distributions with the same means. These data sets were filtered to extract only those genes whose q value was # 0.01 and whose fold change between mean expression values was $ 1.5 in either direction. Fisher's exact test was used to test whether the number of such genes is overrepresented in the xylem cell specification and differentiation gene regulatory network. Gene regulatory network inference. Expression data 30 were used, after normalization with the mmgMOS method from the PUMA R package 47 . The supervised regulatory interactions network was constructed using SIRENE 48 . The directionality of the interactions is defined by the protein-DNA interactions from Y1H data. The interaction sign (activation or repression) is derived by Pearson's correlation coefficient for each protein-DNA interaction. The analysis performed was categorized as (1) supervised tier Ia, network inferred with SIRENE with the provided Y1H gene regulatory connections and the corresponding gene expression profiles (16 genes, 4 transcription factors); (2) supervised tier Ib, an additional three verified connections from the supervised tier Ia and unsupervised tier I were considered in the inference. The unsupervised regulatory interaction network was constructed using the consensus from four different gene regulatory network inference methods, GENIE3 49 , Inferelator 50 , TIGRESS 51 and ANOVerence 52 . The data used were the same as the supervised TIERIa network. The default parameters were used in all methods and a rank-based method was used to build the consensus network as in ref. 53.