Temporal analysis of Arabidopsis genes activated by Eucalyptus grandis NAC transcription factors associated with xylem fibre and vessel development

Secondary cell wall (SCW) deposition in Arabidopsis is regulated among others by NAC transcription factors, where SND1 chiefly initiates xylem fibre differentiation while VND6 controls metaxylem vessel SCW development, especially programmed cell death and wall patterning. The translational relevance of Arabidopsis SCW regulation theory and the utility of characterized transcription factors as modular synthetic biology tools for improving commercial fibre crops is unclear. We investigated inter-lineage gene activation dynamics for potential fibre and vessel differentiation regulators from the widely grown hardwood Eucalyptus grandis (Myrtales). EgrNAC26, a VND6 homolog, and EgrNAC61, an SND1 homolog, were transiently expressed in Arabidopsis mesophyll protoplasts in parallel to determine early and late (i.e. 7 and 14 hours post-transfection) gene targets. Surprisingly, across the time series EgrNAC26 activated only a subset of SCW-related transcription factors and biosynthetic genes activated by EgrNAC61, specializing instead in targeting vessel-specific wall pit and programmed cell death markers. Promoters of EgrNAC26 and EgrNAC61 both induced reporter gene expression in vessels of young Arabidopsis plants, with EgrNAC61 also conferring xylem- and cork cambium-preferential expression in Populus. Our results demonstrate partial conservation, with notable exceptions, of SND1 and VND6 homologs in Eucalyptus and a first report of cork cambium expression for EgrNAC61.

through VND7 are vessel-specific, with VND6 specifically regulating metaxylem SCW deposition and VND7 regulating both meta-and protoxylem vessel formation [15][16][17][18][19][20][21][22][23][24] . Thus, VND6 is a key regulator of the reticulated and pitted wall patterning observed in secondary xylem vessels, the deposition of which is determined by the bundled microtubule structure of the cytoskeleton 25,26 . In woody angiosperms such as Populus the expression of NST and VND clade homologs appear to overlap somewhat, with both the SND1 homolog PtrWND1B (a Populus trichocarpa wood-associated NAC domain protein) and the VND7 homolog PtrWND6A being expressed in xylem and phloem fibres, while vessel-specific differentiation appears to be regulated by the exclusive expression of PtrWND6A in vessels 27,28 . High-resolution spatial transcript profiling in aspen from phloem through the cambium to the lignified xylem zone revealed biphasic expression peaks for SND1 homologs in phloem and early xylem SCW deposition, while VND6 homolog transcripts peaked either during xylem SCW deposition or its cessation 29 , suggesting a specialization of SND1 homologs in phloem and xylem fibre formation and VND6 homologs in vessel differentiation. However, dominant repression in Populus of either the SND1 homolog PtrWND2B or the VND7 homolog PtrWND6B resulted in significantly reduced xylem SCW deposition in both fibres and vessels 30 . Although this suggests less distinct roles for SWN-mediated regulation of fibre and vessel SCW formation in woody angiosperms, the question of partially overlapping versus distinct roles in secondary xylem development remains poorly resolved. For example, homologs of SND1 and VND6/VND7 in monocots (a lineage lacking secondary xylem derived from a vascular cambium) appear to be expressed indistinguishably in sclerenchyma fibres and vessels 31,32 , while in Norway Spruce (a woody gymnosperm lacking fibres and vessels) sufficiently distinct roles could be inferred for VND and NST homologs during xylogenesis 33 .
Xylem fibre and vessel differentiation is distinguished by differences in SCW thickness and patterning as well as the timing and rate of programmed cell death (PCD) and autolysis, which in the case of water-conducting vessel elements yields hollow lumens shortly after SCW deposition 6,34,35 . Congruent with their proposed functions in xylem vessel development, VND6 and VND7 strongly activate PCD genes in Arabidopsis 9,36 . PCD and autolysis, which are distinct biological processes, are initiated by Ca 2+ influx signals resulting from extracellular proteolysis by serine proteases 37,38 . Proteins currently known to be involved in vessel autolysis include XYLEM CYSTEINE PEPTIDASE 1 (XCP1), XCP2 and METACASPASE 9 (MC9) that together cause autolysis of the cytoplasm, the endonuclease ZEN1, XYLEM SERINE PEPTIDASE 1 (XSP1), and possibly ARABIDOPSIS THALIANA SUBTILASE 1.1 (SBT1.1) [39][40][41][42][43][44] . The later onset of PCD in fibres may be attributed to the weaker activation of PCD genes by SND1 and slower accumulation of the necessary enzymes to prompt cellular degradation 45,46 . Compared to fibres, metaxylem vessels exhibit extensive, pronounced pits in their SCWs that allow for the lateral movement of water to living fibres and xylem parenchyma 47,48 . Pits are established by the recruitment of microtubule disassembly proteins MIDD1 and Kinesin 13 A by pit localization of GTP-bound ROP11, a Rho GTPase that is localized to oval planar domains at the plasma membrane by IQD13, activated by ROPGEF4 and deactivated by ROPGAP3 [49][50][51] .
Despite what is believed to be substantial evolutionary conservation in SCW transcriptional regulation 52 , the prevalence and significance of lineage-specific adaptations to the SCW transcriptional network architecture among angiosperms are poorly understood. In Populus, for example, numerous co-orthologs of SND1 exist, of which alternative splicing of PtrWND1B in particular produces full-length and intron-retaining, truncated isoforms of the protein that can form non-functional heterodimers 27,[53][54][55] . Since Arabidopsis does not harbour known SND1 splice variants, this shows that SND1-mediated SCW regulation may differ between plant lineages. In the Eucalyptus woody model, a small number of studies have aided in our understanding of SCW regulation [56][57][58][59][60][61] . To date, our knowledge of the functions of SWN master regulators in Eucalyptus is limited to a few gene targets of EgWND1, an E. gunnii homolog of SND1 30 , with some evidence that an intron 2-retaining splice variant may also exist in E. grandis 27 . In previous work 62 , we identified putative E. grandis orthologs of SND1 (EgrNAC61, also orthologous to EgWND1) and VND6 (EgrNAC26) which we hypothesised regulate fibre and vessel SCW development, respectively.
While functionally characterized SCW TFs and their promoters are potential tools for the re-engineering of woody biomass, their crude manipulation generally results in undesirable phenotypes and hence refined approaches adopting synthetic and systems biology are favourable 63 . The development of high-precision multipartite synthetic gene constructs, however, assumes individual parts possess modularity (functional fidelity of discrete parts in different host organisms) and orthogonality (independence from and minimal interference by native endogenous processes), both of which are poorly demonstrated for SCW TFs. For instance, what gene targets might we expect to be activated when SCW TF homologs from diverged lineages are heterologously and ectopically induced, and how is this influenced when their endogenous promoters are exchanged with synthetic promoters? In this study, we functionally characterized E. grandis SCW master regulators through observation of their SCW network targets over time in the context of Arabidopsis mesophyll protoplasts, and their promoter activities in woody and non-woody models. Our results show that EgrNAC61 is a more potent activator of SCW genes while the VND6 homolog EgrNAC26 mainly regulates PCD and autolysis in vessels.

Results
EgrNAC26 and EgrNAC61 transcript splicing, transcriptional activity and spatial expression patterns. We PCR-amplified the EgrNAC26 and EgrNAC61 coding sequences from E. grandis developing secondary xylem cDNA for cloning. We found no evidence of alternative splice variants among the amplicons (Fig. S1). Similarly, analysis of RNA-seq reads from an independent developing secondary xylem sample 64,65 showed no significant coverage across any of the introns (Fig. S2). These results therefore do not support those of the previously observed EgrNAC61 intron 2 retention reported by Zhao et al. 27 . Since EgrNAC26 and EgrNAC61 are likely master regulators of SCW deposition in E. grandis based on homology, we next tested the ability of these proteins to activate transcription of the His reporter gene in two strains of S. cerevisiae. Both strains expressing either EgrNAC61 or EgrNAC26 grew on histidine-deficient media supplemented with up to 15 mM of the His reporter inhibitor 3AT (Fig. 1), while for the empty-vector control line significantly less growth was observed on media lacking histidine and no growth was observed on media supplemented with 3AT. For both strains, stronger reporter gene activation at 15 mM 3AT was observed for EgrNAC61 than EgrNAC26. These results demonstrate that EgrNAC26 and EgrNAC61 have strong intrinsic transcriptional activation activities.
We analysed promoter::GFP:GUS reporter fusions of the EgrNAC26 and EgrNAC61 promoter sequences to determine their heterologous expression patterns. In the vasculature of Arabidopsis seedling roots, which contain only xylem vessels, the metaxylem vessel marker VND6pro::GFP:GUS produced predominant GUS staining in the metaxylem vessels in the centre of the vasculature, with less blue staining in protoxylem vessels (Fig. 2a). In contrast, a 35Spro::GFP:GUS constitutive control showed staining in all root cell types. In most cases, GUS staining in EgrNAC26pro::GFP:GUS lines was not visible; in a few cases, faint expression was observed in the protoxylem vessels. However, GFP expression for this construct was clearly detected in the seedling hypocotyl vasculature (Fig. 2a). Previously, SND1 promoter-directed expression was observed predominantly in interfascicular and xylary fibres of inflorescence stems and hypocotyls of flowering Arabidopsis plants, and also in developing vessels undergoing SCW thickening 16,20 . In EgrNAC61pro::GFP:GUS plants, GUS staining was evident in protoand metaxylem vessels in the root vasculature (Fig. 2a). Analysis of GFP expression in the vasculature of seedling hypocotyls for the VND6, EgrNAC26 and EgrNAC61 promoters corroborates our findings of GUS staining in the seedling vessels (Fig. 2a).
We also explored the expression patterns of EgrNAC61 during wood formation through promoter::GUS analysis in transgenic poplar plants. GUS signal was strongest in the developing xylem region of the first internode, just before lignification of the primary xylem vessels distinguishable with phloroglucinol staining (Fig. 2b). In secondary xylem in the fifth internode, expression was observed throughout the cortex but was strongest in the xylem expansion zone, sandwiched by the vascular cambium and the lignification zone, as well as in the cork cambial region (phellogen). The latter is a novel report. GUS expression appeared to be visible in developing secondary vessel and fibre cells, forming a near-continuous ring of signal along the expansion zone (Fig. 2b, internode 5). Together, the Arabidopsis and poplar results support a role for EgrNAC26 in the regulation of vessel formation, and EgrNAC61 in the regulation of secondary xylem fibre and vessel deposition.

Identification of early and late transcriptional targets of EgrNAC26 and EgrNAC61 in
Arabidopsis protoplasts. In order to identify early (likely direct) and late (likely indirect) gene targets of EgrNAC26 and EgrNAC61 in a heterologous background, we transiently overexpressed these genes in linkage with a dual GFP reporter cassette in Arabidopsis leaf mesophyll protoplasts. Total RNA was isolated 7 and 14 hours post-transfection (hpt) from three independent protoplast isolations per time point. Transfection efficiencies were approximately 30-50% (Fig. S3). Following RNA quality analysis (Fig. S4) and strand-specific RNA sequencing, between 17 and 36 million paired-end reads were obtained per library, most with >90% concordant mapping of read pairs to the Arabidopsis TAIR10 genome (Table S2)  . The overlap between EgrNAC26-OX and EgrNAC61-OX datasets was highly significant (P < 2.2 × 10 −16 ), only one transcript was alternately up-or down-regulated across time points, and DEGs common to EgrNAC26-OX and EgrNAC61-OX datasets were all up-or down-regulated concordantly, demonstrating robust and consistent activation of a subset of shared genes. Overrepresented biological processes were identified for EgrNAC26-OX and EgrNAC61-OX DEGs through Gene Ontology (GO) enrichment analysis. Significantly overrepresented terms among up-regulated DEGs in the EgrNAC26-OX and EgrNAC61-OX datasets were largely associated with SCW biosynthesis, particularly at 14 hpt. Common terms included secondary cell wall biogenesis, cellulose biosynthesis and (glucurono)xylan biosynthesis (Fig. 3b), supporting the hypothesized roles of these genes as SCW master regulators since the overexpression of the transgenes was sufficient to cause ectopic activation of these genes. Differences between the datasets provided insight into known and novel biological functions that were likely distinctly regulated by EgrNAC26 and EgrNAC61. Specifically, supporting evidence for the role of EgrNAC26 as a vessel-preferential master regulator was the enrichment of genes in the category "developmental PCD" at 7 hpt (Fig. 3b). In contrast, EgrNAC61 induced a good representation of cytoskeleton-related processes, a number of developmental processes such as trichome and epidermal differentiation, and interestingly several primary metabolic processes such as amino acid and pyrimidine metabolism in addition to secondary metabolism. These results are consistent with a role of EgrNAC26 and EgrNAC61 in the regulation of SCW synthesis.   (Tables 1, S3). Likewise, several lignin-specific genes were upregulated in EgrNAC61-OX at 14 hpt only, including 4CL3, CCoAOMT1, and HCT. Interestingly, far fewer SCW-related DEGs were unique to EgrNAC26-OX; noteworthy examples included LBD30/ASL19, VND2, TED7, a key cytoskeleton-associated protein MAP70-5 and the ROP11 activator ROPGEF4 (Table 1). No activation of the endogenous SND1 gene was observed, suggesting that the self-activation of SND1 observed in plasmid-based   promoter transactivation and electrophoretic mobility shift assays 68 did occur in the context of native mesophyll chromatin. Several SCW-related genes and TFs were upregulated in at least one time point for both EgrNAC26-OX and EgrNAC61-OX. However, we noted that apart from the abovementioned SCW-related genes and TFs being exclusively activated by EgrNAC61, this TF activated several key TFs and many SCW biosynthetic genes five-to ten-fold more strongly than EgrNAC26 by 14 hpt. These included TFs BLH6, C3H14, KNAT7, MYB46, SND2 and SND3, cellulose-associated CesA4, CesA7, CesA8 and COBL4, xylan-associated genes such as IRX9, IRX10, IRX15, GUX1, GUS2 and PARVUS, and lignin-associated genes such as CCoAOMT7, LAC4 and LAC17 (Table 1). This, despite a lower transfection efficiency for the EgrNAC61-OX construct (Fig. S3b). In contrast, although PCD genes were expressed in all datasets, BFN1, MC9, XCP1 and XCP2 were upregulated by as much as an order of magnitude more in EgrNAC26-OX compared to EgrNAC61-OX at 14 h (Table 1). Supported by the apparently stronger transactivation ability of EgrNAC61 in yeast (demonstrated by its persistent growth at 15 mM 3AT, Fig. 1), our results demonstrate that EgrNAC61 is a more potent transcriptional activator of the Arabidopsis SCW pathway compared to EgrNAC26, while the latter shows a similar trend for PCD-associated genes in particular.
To understand the temporal activation of SCW transcriptional network nodes following EgrNAC26 and EgrNAC61 induction in Arabidopsis protoplasts, we mapped the DEGs to a curated Arabidopsis SCW transcriptional network compiled from the literature 52 . Time-lapse activation of the known SCW transcriptional network nodes by EgrNAC61 and EgrNAC26 is represented in Figs 4 and 5, respectively. EgrNAC61 induced strong, early upregulation of second and third tier TFs as well as structural genes associated with cellulose, hemicellulose and lignin biosynthesis by 7 hpt and a few additional targets in all tiers by 14 hpt (Fig. 4). The effect of EgrNAC26 induction clearly differed from that of EgrNAC61, where strong activation of PCD-associated genes but weak activation of only a few SCW structural genes and second-or third-tier TFs were apparent at 7 hpt, with a better representation of mid-tier SCW-related TFs and biosynthetic genes by 14 hpt, although weakly activated compared to EgrNAC61-OX (Fig. 5). Given that EgrNAC26 and EgrNAC61 transfections were performed in parallel on a common preparation of protoplasts, these results demonstrate a weaker and possibly delayed onset of SCW structural gene activation for EgrNAC26 in comparison. Finally, we noted that two clusters of mid-and lower-tier TFs implicated in the SCW network 52 were unaffected by induction of either EgrNAC26 or EgrNAC61 (Figs 4 and 5, left). These clusters are not known to be activated by SWNs, consistent with our results and suggesting an independent route of regulation.

Discussion
In this study, we explored the temporal transactivation dynamics of Arabidopsis genes (Brassicales) by homologs of SCW master regulators of fibre and metaxylem vessel formation from E. grandis, a commercially important hardwood. Our results provide (1) insight into the possible roles of EgrNAC26 and EgrNAC61 SWNs in fibre and vessel SCW regulation in eucalypts, (2) evidence of considerable conservation in SWN gene targets between Eucalyptus as a representative of the basal eudicot lineage the Myrtales which divergence from the Brassicales over 100 million years ago 69 , and (3) a demonstration of the regulatory dynamics and specific targets of orthologous master regulators associated with fibre and vessel differentiation in a heterologous context typical of a synthetic biology strategy.
Few studies have explored the genome-wide gene targets of Arabidopsis SND1 and VND6/VND7 9,24,36,70 . These previous investigations revealed preferential activation of PCD-associated genes such as XCP1 and XCP2 during vessel differentiation by VND6 and VND7. In woody angiosperms, only one study has reported the genomic targets of a fibre-associated SCW master regulator, in this case a time-course experiment of poplar stem differentiating xylem protoplasts transfected with Ptr-SND1-B1, a homolog of SND1 from P. trichocarpa 71 . Ours is the first study to probe the global transcriptional targets of EgrNAC26 and EgrNAC61, albeit in a heterologous system (Arabidopsis). One advantage of our approach is that the transfections for each construct were performed in parallel from a single protoplast preparation, allowing for a direct comparison of the gene targets and temporal activation dynamics of each SWN candidate. Two limitations should be noted: first, mesophyll protoplasts may not reflect the accessible chromatin landscape or contain the same co-factors as xylem mother cells. However, SWNs and their homologs are known to be particularly powerful SCW master regulators, sufficient  Table 1. Secondary cell wall-associated genes regulated by EgrNAC26 and EgrNAC61. *Inf refers to instances where the FPKM was 0 in the empty-vector control, but was larger than 1 in the test sample. hpt, hours post transfection.
for ectopic tracheary element transdifferentiation in many cell types 20,[72][73][74] . The inherent transcriptional transactivation properties of EgrNAC26 and especially EgrNAC61 in yeast (Fig. 1) and the extensive and significant enrichment of SCW-associated biological processes in all EgrNAC26-OX and EgrNAC61-OX DEG datasets   Fig. 3b) support this statement. Thus, we can infer that EgrNAC26 and EgrNAC61 activate homologs of genes similarly in protoplasts as they do in developing xylem. Second, since the use of a heterologous system may yield Arabidopsis-specific targets and cannot detect novel gene targets unique to E. grandis, validation of the Eucalyptus orthologs of the target genes in this study is recommended. Given the predominant fibre-and vessel-preferential SCW regulatory functions of SND1 and VND6 respectively in Arabidopsis, we expected to observe similar cell type associations for their putative E. grandis orthologs EgrNAC61 and EgrNAC26. In contrast to this hypothesis, we observed a substantial overlap in promoter activity between EgrNAC61pro and EgrNAC26pro in during Arabidopsis vascular development, namely in seedling vessels (Fig. 2). The EgrNAC26 promoter-driven reporter expression in Arabidopsis primary xylem vessels was substantially weaker, however. It is possible that the EgrNAC26 promoter construct, which included only 1,223 bp upstream of the predicted transcription start site due to a large 2,422 bp 5′ UTR region containing an embedded intron, lacked a critical cis-element for stronger expression, or that the promoter is regulated differently in E. grandis.
In comparison with our promoter analysis, expression of the aspen SND1 homolog peaked in phloem and early-onset SCW deposition in the xylem from cryosectioned xylogenesis developmental series 29 , consistent with our observation for developing xylem although it is unclear why no EgrNAC61pro activity was observed in phloem fibres (Fig. 2b). The aspen series 29 did not include cork cambium, but we here report strong GUS induction by EgrNAC61pro in this tissue layer, while EgrNAC26pro was not assessed. Cork (phellem) consists of lignified and suberized cells that undergo PCD and have cell walls in which a guaiacyl-rich lignin backbone predominates 75,76 . Phloroglucinol staining revealed small amounts of lignin in the phellem of the poplar transgenics (Fig. 2b, bottom panel). Since EgrNAC61 activates the phenylpropanoid pathway (Figs 3, 5), we postulate that it has been co-opted to activate phellem lignification in the cork cambium too.
In addition to its homology to VND6, several lines of evidence point to a vessel-specific SCW regulatory role for EgrNAC26 in the transient expression experiment: (i) EgrNAC26 transfection uniquely induced vessel-specific regulators of SCW deposition, among them LATERAL ORGAN BOUNDARIES 30 (LBD30/ ASL19), VND2 and TRACHEARY ELEMENT DIFFERENTION-RELATED 7 (TED7). LBD30/ASL19 and LBD18/ASL20 are targets of VND6 that participate in a positive feedback loop during vessel differentiation 77 , and VND2 is a vessel-specific SCW master regulator that interacts with VND7 15,23,74 . TED7 is a vessel-specific protein with an unknown role in tracheary element SCW formation, acting with TED6 and interacting with the SCW-related CesA complex 78 . While TED6 was significantly upregulated by both EgrNAC26 and EgrNAC61 at 14 hpt (Table 1), the exclusive activation of TED7 at the earliest time point by EgrNAC26 suggests a vessel-related specialization. (ii) EgrNAC26 activated proteins associated with SCW pit formation in vessels, uniquely among them RHO GUANYL-NUCLEOTIDE EXCHANGE FACTOR 4 (ROPGEF4) and MICROTUBULE-ASSOCIATED   (Table 1). ROPGEF4 is essential for pit formation through its highly localized activation of the small GTPase ROP11, which was a common target of EgrNAC26 and EgrNAC61 in addition to the microtubule-associated protein MIDD1 and the ROP11-deactivating protein ROPGAP3 79 . MAP70-5 is a tracheary element-specific microtubule-implicated protein, the abundance of which has a profound effect on the frequency of pitted, reticulated or smooth wall patterns observed in transdifferentiated cells 80 . Thus, its strong upregulation at 14 hpt in EgrNAC26-OX may re-structure the cytoskeleton for vessel-specific microtubule-directed SCW deposition. (iii) EgrNAC26 transfection induced stronger activation of the chief genes associated with autolysis of tracheary elements during PCD, namely XCP1, XCP2 and MC9 46 , as well as BIFUNCTIONAL NUCLEASE 1 (BFN1), the closest homolog of the Zinnia endonuclease ZEN1 known to degrade nuclei in tracheary elements 41 ( Table 1). The exclusive enrichment of "programmed cell death" as a biological function term in the early induction (7 hpt) of EgrNAC26 (Fig. 3b) strengthens this observation. A ribonuclease with a currently unknown role in xylem PCD, RNS3, was rapidly induced following transfection by EgrNAC26 and EgrNAC61 (Table 1), while expression of suspected but non-validated xylem PCD genes XSP1 and SBT1.1 were not affected in this experiment.
We also observed what appears to be a much weaker and possibly delayed response in SCW target activation by EgrNAC26 compared to EgrNAC61. Interestingly, the expression of EgrNAC26 in protoplasts was substantially lower than that of EgrNAC61, despite the use of the same constitutive promoter and the fact that the transfection efficiency for EgrNAC26 was somewhat higher (Figs S3b, S5a). We cannot currently explain this phenomenon. However, several promoter transactivation experiments in Populus and Arabidopsis showed that NST clade SWNs activate mid-level regulators PtrMYB021 (i.e. MYB46), PtrMYB020 (MYB83), PtrMYB128 (MYB103), PtrNAC156 (SND2), PtrKNAT7 (KNAT7) and PtrMYB28 (MYB58/MYB63) more strongly than VND clade SWNs 13,28,30,55 , a trend that also emerged in our time-course experiment where SCW structural genes were preferentially activated by EgrNAC61 much earlier and to a stronger level than for EgrNAC26 (Figs 4, 5). Since fibre cells exhibit thicker SCWs than vessels, stronger activation of SCW structural genes in fibres is developmentally plausible, and transcript levels of SCW cellulose, xylan and lignin genes appear higher in laser micro-dissected fibres compared to vessels in Populus in support of this 81 . An enrichment for primary metabolic processes in EgrNAC61-OX (Fig. 3b) may reflect a metabolic shift toward carbon investment in this carbon sink which is in agreement with previous observations of several amino acid pathways being upregulated during in vitro tracheary element transdifferentiation 82 . Fibre-specific mutant phenotypes for EgrNAC61-preferential targets such as SND2 and SND3 13,83,84 suggest that EgrNAC61 is an important fibre SCW regulator, but the clear reporter gene expression conferred by its promoter in Arabidopsis vessels and all cell types in the developing xylem in poplar plants (Fig. 2) suggests a dual role in SCW formation in fibres and vessels, and possible other xylem cells.
Given the combined evidence of vessel-specific specialization of EgrNAC26 and preferential PCD activation, versus vessel-and fibre-preferential expression of EgrNAC61 and strong activation of SCW structural genes and TFs following induction in Arabidopsis protoplasts, we favour the hypothesis that EgrNAC61 is the chief fibre SCW master regulator that plays a supporting role in vessel development, while EgrNAC26 complements ScIEntIFIc REPORTS | (2018) 8:10983 | DOI:10.1038/s41598-018-29278-w EgrNAC61 in vessels with PCD and SCW pit-related regulatory functions. The potential role of EgrNAC61 in regulating cell wall deposition in cork cells is an interesting topic for further investigation.
Plant transformation and GUS analysis. The pBGWFS7 constructs (VND6pro::GFP:GUS, EgrNAC26pro::GFP:GUS and EgrNAC61pro::GFP:GUS) were transformed into Arabidopsis thaliana using the Agrobacterium-mediated floral dip method 90 and the transformed lines were selected on half-strength MS agar medium using Glufosinate-ammonium (20 mg/l). Transgenic Arabidopsis seeds were surface-sterilized, plated on nonselective nitrate-enriched media and allowed to germinate for 2-3 weeks. For poplar transformation, the EgrNAC61pro::GUS construct in the pMDC162 backbone was introduced into poplar hybrid Populus alba × P. tremula 717-1B4 by co-cultivation of Agrobacterium solution with leaf discs followed by shoot regeneration on selective media containing 30 mg/l Hygromycin. Transgenic poplar plants were propagated using greenwood stem cuttings and young developing stems (internodes 1 and 5) were harvested. GUS staining was performed as described by Spokevicius et al. 91 , excluding acetone treatment, but adding 0.5 mM potassium ferrocyanide and 0.5 mM potassium ferricyanide. After incubation at 37 °C overnight, seedlings were destained in 70% ethanol and visualized with differential interference contrast and brightfield microscopy.

Isolation and transfection of Arabidopsis leaf mesophyll protoplasts. A. thaliana ecotype
Columbia (Col-0) seedlings were grown for four weeks at 22 °C under a 10 h photoperiod prior to mesophyll protoplast isolation from leaves as previously described 93 . Multiple plants were pooled for protoplast isolations and de-bulked to form three independent biological replicates. Large-scale transfections with pUC-35S-EgrNAC26-35S-GFP, pUC-35S-EgrNAC61-35S-GFP, and pUC-35S-NOS-35S-GFP (empty vector) constructs were conducted using 4 × 10 5 protoplasts for each of three technical replicates per biological replicate (n = 3) to account for transfection variation and incubated for 7 h and 14 h at room temperature in the dark. Technical replicates were bulked and harvested through centrifugation at 500 × g. Of the harvested cells, 10 μl was used to determine transfection efficiency (at 14 hpt only) and the remaining cells were flash-frozen and stored at −80 °C.
RNA-seq analysis. Total RNA was isolated using the Qiagen RNeasy Mini Kit (Sigma-Aldrich®, South Africa) and RNA quality was determined using an Experion instrument (Bio-Rad). Stranded 100 nt paired-end (PE100) mRNA sequencing was performed at Beijing Genomics Institute (BGI, Hong Kong) on the Illumina HiSeq4000 platform using fragment sizes of up to 200 bp. Raw RNA-seq reads passing initial quality filtering by BGI were further assessed for quality with FastQC 94 . The RNA-seq reads were formatted with FastQGroomer 95 prior to mapping to the A. thaliana genome (TAIR10) using TopHat 2 96 . The mapped RNA-seq reads were assigned to gene models and the fragments per kilobase of transcript per million reads mapped (FPKM) values for each mapped transcript was calculated in Cufflinks 97 using the A. thaliana TAIR10 annotation. The list of genes generated by Cufflinks was compiled into a single list with all nine libraries and Principle Component Analysis (PCA) was applied in R 98 using the prcomp function 99 . Differentially expressed genes (DEG) were identified using the mapped RNA-seq reads in Cuffdiff 97 , with FPKM ≥1, false discovery rate (FDR) >0.05 and fold change ≥2 or ≤0.5.
The BiNGO plugin in Cytoscape 100,101 was used to determine overrepresented Gene Ontology (GO) terms in the up-and down-regulated EgrNAC26 and EgrNAC61 DEG datasets. A hypergeometric test with Benjamini and Hochberg FDR correction with a P-value cut-off of 0.05 was applied with the whole A. thaliana annotation as a reference set.
Data availability. The RNA-seq data reported in this study is freely available at the NCBI Sequence Read Archive (https://www.ncbi.nlm.nih.gov/sra), accession SRP117192.