Induction of PrMADS10 on the lower side of bent pine tree stems: potential role in modifying plant cell wall properties and wood anatomy

The molecular mechanisms underlying inclination responses in trees are unclear. In this study, we identified a MADS-box transcription factor differentially expressed early after inclination in the stems of Pinus radiata D. Don. PrMADS10 has a CDS of 582 bp and encodes a group II MADS-box transcription factor. We measured highest accumulation of this transcript on the lower side of inclined pine stems. In an effort to identify putative targets, we stably transformed Arabidopsis thaliana with a 35S::PrMADS10 construct. Transcriptome analysis revealed 1,219 genes differentially-expressed, with 690 and 529 genes up- and down-regulated respectively, when comparing the transgenic and wild-type. Differentially-expressed genes belong to different biological processes, but were enriched in cell wall remodeling and phenylpropanoid metabolic functions. Interestingly, lignin content was 30% higher in transgenic as compared to wild-type plants consistent with observed changes in gene expression. Differentially expressed transcription factors and phenylpropanoid genes were analyzed using STRING. Several MYB and NAC transcription factors showed interactions with genes of the phenylpropanoid pathway. Together, these results implicate PrMADS10 as a regulatory factor, triggering the expression of other transcription factors and genes involved in the synthesis of lignin.


Results
Sequence and phylogenetic analysis of PrMADS10. PrMADS10 full-length cDNA sequence was obtained using a partial EST sequence as template from the SSH library 29 , followed by 5′-and 3′-RACE-PCR. The sequence of 943 bp long contains 111 and 250 bp of 5′-and 3′-UTRs, respectively. PrMADS10 has a CDS of 582 bp, encoding a deduced protein of 193 amino acids and 22 kDa (pI 9.42; GenBank accession number, KM887510; Fig. 1A). The predicted PrMADS10 protein has the typical conserved structural features of MADS-box TFs and possesses a MIKC type protein structure (Fig. 1A). A phylogenetic analysis was performed with 48 MADS-box amino acid sequences, including proteins from poplar, pine and spruce. PrMADS10 is classified in group II (MIKC c ) according to MADS sub classification, and very close to orphans genes like AtSVP and StMADS16, genes that are expressed in vegetative tissue (Fig. 1B).
PrMADS10 transcripts are preferentially accumulated in pine stems and its protein has nuclear localization. RT-qPCR analysis shows a higher transcript accumulation of PrMADS10 in the lower side of inclined stems, reaching a maximal difference between upper and lower sides after 24 h (Fig. 2). A four-fold transcript increase was detected on the lower side of inclined seedling stems 24 h after the onset of the inclination response in comparison to non-inclined control. In addition, transcript levels were strongly downregulated in needles and roots after inclination.
To determine the subcellular localization of PrMADS10 tobacco leaves were agroinfiltrated 55 (Fig. 3). PrMADS10 fused to GFP co-localizes with SYTO83, a nuclear marker. This data reveals PrMADS10 is a nuclear localized protein, consistent with its predicted function as a transcription factor based on protein sequence.

Overexpression of PrMADS10 in Arabidopsis alters lignin and flavonoids accumulation. The
CDS of radiata pine MADS10 gene was isolated and introduced into pBI121 binary vector under the control of CMV35S promoter. Homozygous Arabidopsis T3 transgenic lines were obtained and used for performing morphological and transcriptional analyses. Four different PrMADS10 over-expressing lines were obtained (Suppl. Fig. 1). The evaluation of plant morphology and the weight of one hundred seeds indicated no changes between transgenic and control lines, however the length of rosette leaves were smaller in the transgenic lines (Suppl. Fig. 2). At the same time, lignin content in three of the transgenic lines was 30% higher than wild-type plants www.nature.com/scientificreports www.nature.com/scientificreports/ Networking interaction analysis in Arabidopsis overexpressing PrMADS10. MapMan 53,54 was used to visualize genes modulated by 35S::PrMADS10 in Arabidopsis. The terms for MapMan were assigned for 166 differentially-expressed genes; but due to a lack of matching MapMan terms, no functional assignment was done for the other 1,055 genes. The results obtained show a general view of different functional categories, which are being affected by 35S::PrMADS10, such as major and minor carbohydrate, amino acid, nucleotide, fermentation, lipids, secondary metabolisms and cell wall (Fig. 6). Metabolism-related genes are further divided as 56 genes for cell wall, 30 for lipid, 29 for secondary, 12 for amino acid, 5 for light, 9 for major carbohydrate, 7 for minor carbohydrate, and a few for glycolysis, fermentation, tricarboxylic acid (TCA) cycle, S-assimilation, nucleotide metabolism, tetrapyrrole synthesis, or mitochondrial electron transport. When the phenylpropanoid pathway was analyzed, several genes were positively-or negatively-regulated (Fig. 6).
Differentially expressed TFs and phenylpropanoid pathway genes were selected from the microarray assay to perform a MapMan classification. The list of TFs and phenylpropanoid pathway´s genes were introduced in STRING in order to determine network interactions (Fig. 7). Genes without interactions were erased from the Transcripts levels of PrMADS10 in young radiata pine seedlings after inclination stimuli. Stem samples were taken at different times of inclination either on the lower stem side (L) or upper stem side (U) of the stem. Roots and needles were obtained from the same inclined (I) seedlings. Ctrl (control) means noninclined seedlings. Data correspond to mean ± SE of three biological replicates and different letters indicate statistical differences (p ≤ 0.05). www.nature.com/scientificreports www.nature.com/scientificreports/ image. The interactions observed can be of different types, depending on the color of the line joining the spheres, and recognized as co-expression, co-localization and text-mining, and therefore it cannot be taken strictly as a direct protein-protein interaction. It can be observed that genes were clustered into different groups. TFs like MYB and NAC are included in one group (MYB50, MYB52, MYB42, NAC010, NAC073, MYB43, MYB83, HB16, MYB85, MYB32, NAC012 and NST2), and PAL, 4-coumarate-CoA ligase (4CL), and caffeoyl-CoA O-methyltransferase (CCoAOMT) are in a different one. Interestingly, seven genes from the phenylpropanoid pathway are grouped and interconnected. MapMan showed that these genes (PAL4, 4CL8, 4CL2, CCOAMT, 4CL5, CCR2; AT4G26220; UDP glucosyl transferase 72E1) are differentially expressed. The gene AT4G26220 has been reported as a probable CCoAOMT7 (amplified view of Fig. 7).
In a small network interaction, MYC4 and MYB76 are grouped to other TFs. At the same time, LHY is related to AGL24, TOC1 and CCA, and indirectly to DOF, PRR5 or PRR9. Finally, a zinc transporter protein correlates with WRKY70, which has direct interactions with WRKY38, WRKY51 and WRKY54.

Discussion
The importance in flower development was the initial interest to study MADS-box genes in plants 33 . However, a MADS-box gene (PdPI) was isolated from Populus deltoides, and it is expressed during flower development and in different vegetative organs. The evidence suggests that PdPI plays multiple functions in the development of this species 57 . In radiata pine, these TFs are expressed in male and female strobiles, for example PrMADS3 (AGL6-like) is transcribed in the primordium of the acule 47 . Even though evidence defines roles for MADS-box TFs in flowering and other reproductive processes, their expression in roots, leaves and stems 33 , like in the case of PrMADS10, suggest another type of function in these other organs, perhaps expressed in response to stress.
In a first effort to understand the role that this MADS-box gene plays in radiata pine, the full length of PrMADS10 was obtained by RACE. Deduced amino acid sequence was aligned with MADS-box proteins from different species available in Genbank database. PrMADS10 contains the typical MADS-box conserved domains: MADS-box, Intermittent, K-box and C-terminal 58 . According to Shore and Sharrocks 43 , a conserved sequence of 56 amino acids is characteristic in this domain, where 16 of which are identical in all family members. Structural features can also be observed in the MADS-box domain 44 , and the region of DNA responsible for protein interactions and dimerization 43 .
Phylogenetic analysis of PrMADS10 showed a distribution of sequences grouped into two main clades, with the sequences of interest being clustered in clade II. In this clade, all proteins have a MIKC-like structure 59,60 . www.nature.com/scientificreports www.nature.com/scientificreports/ In the phylogenetic tree, PrMADS10 is closely related to PaMADS1 from P. abies (Fig. 1B). PaMADS1 is expressed mainly in female and male pine cones, controlling the development of these tissues 61 . In addition, protein sequences of interest are also associated with a sub-clade classified as similar to SVP, AGL and MADS, such as PtSVP (P. trichocarpa) and AtSVP (A. thaliana), AtAGL24 (A. thaliana), StMADS16 (S. tuberosum) and PtMADS1 (P. tomentosa). PtSVP is a protein from the SVP (Short Vegetative Phase) family, and is very closely-related to AtSVP, which is a negative floral regulator 62 . This sequence is conserved in both angiosperms and gymnosperms. Similarly, AtAGL24 promotes the identity of the inflorescence in Arabidopsis whose expression occurs at the apex of the buds at the time of floral transition 63 . Most of the genes near to PrMADS10 were named orphan genes, because StMADS11, StMADS16, AtAGL24, and AtSVP are expressed in vegetative tissues, like vascular cambium region 64 , and having the same expression pattern than PrMADS10. Thus, the phylogenetic evidence suggests that PrMADS10 could be involved in floral development, but not observable in one-year-old seedlings. Nevertheless, this sequence was isolated from the lower side of pine stems, shortly after tilting 29 . This implies that its role in pine is related to the loss of verticality, modulating gene expression in vegetative tissues rather than in floral development.
PrMADS10 was differentially expressed 2.5 hours after tilting and accumulation of transcripts was preferentially observed on the lower side of bent young pine stems 29 . The result suggests that expression of PrMADS10 is temporally and spatially regulated in pine stem sections. The change in expression levels for this transcription factor in pine stems may be important for control of metabolic pathways that modify cell wall structure and properties which impacts on wood anatomy. Different molecular events related to cell wall modifications in response to verticality loss, including changes in a large number of TFs have been reported 1,29 . The work of Allona 1 and Ramos 29 were complementary, as both report new genes that respond to tilt, albeit there is no mention concerning the differential expression of MADS-box TFs in the former study.
Several different MADS-boxes have been described in radiata pine; PrMADS4 to 9 were detected in vegetative shoots, floral organs and roots, with greater expression in young floral tissue compared to adult tissue 46 . PrMADS1 to 3 could be involved in regulating female reproductive structures 47,48 and cone development 49 . The expression of MADS-box genes during flower development and vegetative organs have also been reported in P. deltoides 57 . In contrast, Cseke 44 reported PTM5 which is specific to vascular tissue and expressed in differentiated primary-secondary xylem and phloem. Therefore, although MADS-box genes are associated with flowering, other tissues like stem can be the tissue where this transcription factor can be expressed, which seems to be the case of PrMADS10, and more likely involved in stems modulating the response to inclination. www.nature.com/scientificreports www.nature.com/scientificreports/ The AraGene microarray shows that 35S::PrMADS10 differentially regulates the expression of 1,219 genes when heterologously-expressed in Arabidopsis. PrMADS10 induces the expression of 690 different genes and down-regulates that of 529. Aswath and Kim 33 state that the constitutive expression of MADS-box genes into tobacco or Arabidopsis plants has proved to be a useful tool to analyze gene functions, as shown by the characterization of SAG1 from the conifer black spruce (Picea mariana), an homologous to AGAMOUS 65 . The overexpression of SAG1 produces homeotic floral conversion in transgenic Arabidopsis (from sepal to carpel and petal to stamen). On the other hand, overexpression of a MADS-box from P. tomentosa (PtAP3) in tobacco plants causes a fast growth and early flowering phenotype 66. The 35S::PrMADS10 transgenic lines of Arabidopsis did not show any discernible phenotypic changes, including flowering time, even though the transgene was over-expressed by at least 20,000 fold in all 4 lines obtained. The first 50 most-overexpressed genes were those involved in phosphate starvation, as well as, TFs and sucrose synthase 4 (SUS4) ( Table 1). Unexpectedly, genes for phosphorous starvation were the most differentially-expressed. Pi was not a limiting factor in our assays, yet some reports indicate that Pi starvation is under the control of TFs required to maintain phosphorus homeostasis, which in turn is affected by environmental stress [67][68][69][70][71] .
Two genes, a multigenic AtSUS isoforms (SUS1-6) that synthesize monosaccharides for cellulose and starch biosynthesis 71 , and rhamnose biosynthesis 3 (RHM3), providing rhamnose for the synthesis of the cell wall 72 , are also up-regulated by PrMADS10. Basic chitinase is also up-regulated, which is a protein leading to www.nature.com/scientificreports www.nature.com/scientificreports/ inhibition of seedling growth 73 . In addition, the gene microtubule-associated protein 70-5 is positively-regulated by PrMADS10. This gene is essential for defining SCW polymers and is expressed in the cellular cortex of wood-forming cells, influencing the pattern of SCW thickenings in tracheary elements 74 . These results suggest  www.nature.com/scientificreports www.nature.com/scientificreports/ that PrMADS10 positively regulates several genes involved in wood formation, such as secondary metabolism, sugar synthesis or remodeling of cell wall components, as shown in the MapMan analysis (Fig. 7).
The accumulation of lignin is increased by 30% in the T3 transgenic lines, concomitant with the differential expression of genes involved in phenylpropanoid pathway, where a few genes from the pathway are positively-regulated by PrMADS10. Interestingly, Pinus taeda MYB4 has been reported as a positive regulator in lignifying cells 75 . Moreover, MYBs have also been related to cell wall metabolism, and have been shown to be the main master switches of SCW biosynthesis in different species 13,76,77 . More recently, the biosynthesis of cellulose/ xylan/lignin is triggered by MYB46/83 and MYB85 78 . Both MYBs (MYB83 and MYB85) genes are up-regulated in our experiment, which indicates that PrMADS10 can induce the expression of both TFs, thus enhancing the accumulation of lignin. A full network of protein-protein relationships reported by STRING indicates that these TFs are also associated with NAC. Moreover, the network correlates the interaction between MYBs and NACs with phenylpropanoid genes, such as CCoAOMT, 4CL2, 4CL5, PAL4 and 4CL8. The co-expression of these groups of genes was reported when transcriptomic analysis was performed in Arabidopsis stem during development 79 . Additionally, in woody plants like eucalyptus and poplar, the over-expression of MYB216 increases the accumulation of transcripts of 4CL5 and PAL4 80,81 . The knock-out for PtMYB156 in poplar showed an increment in the accumulation of 4CL5 transcripts indicating a selective regulation of phenylpropanoid genes 82 .
TFs like NAC and MYB are preferentially co-expressed in xylem in Picea glauca 83 and Pinus pinaster 84 , suggesting that both integrate a network which takes part in the development of complex wood traits, and as regulator for the synthesis of phenylalanine. In our assay, NAC10, NAC073 and MYB52, MYB42, MYB85 TFs were differentially expressed, and also showed relationships with other NACs and MYBs in the STRING network. Coincidently, several MYB TFs were characterized from radiata pine seedlings and are differentially-expressed on both sides of an inclined stem, yet preferentially on the upper side 85 . Whether NAC is modulating the expression of MYB and NAC genes, or if PrMADS10 could play a central role in this interaction in response to bending in pine, is an open question that requires further analysis.
The most regulated genes in the compression side of adult radiata pine tree trunks are those like: cell division, cellulose biosynthesis, lignin deposition and microtubules 86 . Genes like tubulin beta, sucrose synthase, proline-rich protein, and pectin lyase-like were differentially-expressed in 35S::PrMADS10 Arabidopsis, as www.nature.com/scientificreports www.nature.com/scientificreports/ well as, in adult pine trees 86 . Also, tubulin alpha chain, the cellulose synthase gene family, sucrose synthase and expansin were preferentially-expressed in Pinus taeda differentiated xylem 87 . These genes are active in some way during the process of cell wall remodeling, and are also differentially-expressed in 35S::PrMADS10 transgenic Arabidopsis.
Genes from the phenylpropanoid pathway are also differentially-regulated in 35S::PrMADS10 plants. The increase in COMT transcript levels was similar to that observed for PAL. If COMT levels increase, the amounts of synapyl alcohol should be modified, which is the last step for lignin (syringyl) synthesis. Similarly, if COMT and CAD transcript levels increase, they could modify the amounts of synapyl and p-coumaryl alcohols, the final step for the synthesis of lignin (syringyl and p-hydroxyphenyl, respectively). This is corroborated by previous studies, which suggest that the level of C3H transcripts could be determinant in the flow of metabolites towards lignin production, and the transcripts of CCoAOMT, CCR and CAD are modulated according to the metabolic demand 88 . Similar findings were reported in inclined radiata pine stems 29 . The net result of increasing the level of transcripts of these genes is an increment in flux through the monolignol portion for the lignin biosynthesis pathway 75 . Patzlaff 75 reported that PtMYB4 recognizes AC elements of the PAL promoter and genes encoding lignin biosynthesis enzymes are altered after overexpressing the PtMYB4 gene in transgenic tobacco plants. In addition, the deposition of lignin is increased, and spreads to other types of cells that normally do not lignify, showing that PtMYB4 was sufficient to induce lignification in a heterologous system. In contrast, overexpressing PtrMYB3 from populus in Arabidopsis led to an increase in the deposition of the three main polymers of cell wall 89 . A 30% increment in lignin was observed in 35S::PrMADS10 lines, suggesting stimulation towards the lignin biosynthesis pathway. Besides, when one-year old pine seedlings were tilted, lignin accumulation and wall thickening were observed after 15 days 90 . This suggests that remodeling of the cell wall could initially be regulated by PrMADS10 after tilting.

Materials and Methods.
One year-old half-sib Pinus radiata D. Don (radiata pine) seedlings were grown in a local nursery from seeds obtained from an open-pollinated population (half sibs). Seedlings of around 30 cm tall were maintained at 20 °C, and following the protocol established previously 29,91 . Nine inclined seedlings were collected at each sampling time; their stems were cut into different parts along the longitudinal axis into lower and upper sides, pooled, immediately frozen in liquid nitrogen and stored at −80 °C until RNA extraction. Additionally, 'non-inclined' seedlings were sampled as control. Previous inclination procedure was followed 29 .
Gene cloning and vector construction. Pinus radiata MADS10 (PrMADS10) was cloned from a cDNA library constructed from the stems of young seedlings exposed to inclination at different times 29 . The PrMADS10 ORF was amplified and cloned following the RACE strategy, and using RNA from inclined stems (BD Smart RACE cDNA Amplification Kit, Clontech, USA). The www.nature.com/scientificreports www.nature.com/scientificreports/ following primers: PrMADS10-F 5′-CAGATCTCCGTCGGCAGTTAAAGGAAC-3′ and PrMADS10-R 5′-GGCTGCGAAGATAACCCTAGATGCAAG-3′ were used to obtain 3′ and 5′ sequences, respectively. PCR conditions for amplification was: 30 cycles each of 1 min at 94 °C, 1 min at 66 °C, and 3 min at 72 °C. A final 20-extension step of min at 72 °C was performed.
For the construction of the PrMADS10-GFP fusion, the PrMADS10 CDS was amplified with PrMADS10FL-F 5′-CACCATGGCCGGCGAGAAAAGAAAGAT-3′ and PrMADS10FL-R 5′-AATCTTTGATTCGGACGACTGT-3′ primers, without the stop codon and inserted into pENTR/D-TOPO vector (Invitrogen). Subsequently, the insert was transferred to the plant binary vector pKF7WG2 by recombination using Gateway ® LR Clonase TM II Enzyme Mix kit (Invitrogen) following manufacturer's instruction. The 35 S::PrMADS10-GFP construct was used for subcellular localization studies. All clones used were confirmed by sequencing. Subcellular localization in tobacco leaves. Gateway ® LR Clonase TM II Enzyme Mix kit (Invitrogen) was used to perform a 35S:PrMADS10-GFP construction. Manufacturer's instruction was followed for recombination.
After two days of cultures Agrobacterium was collected and centrifuged at 6000 x g for 10 min at 4 °C. The LB medium was supplemented with three antibiotics: gentamycin (100 μg/ml), rifampicin (10 μg/ml) and spectinomycin (50 μg/ml). Pellet was re-suspended in distilled and sterile water. Young tobacco leaves were infiltrated on the abaxial side, and were analyzed after three days post-infiltration. Syto ® 84 Orange Fluorescent Nucleic Acid stain (Thermo Scientific) was used to label the nucleus. Subcellular localization of PrMADS10 in transient transformed leaves were analyzed through tissue visualization under a confocal fluorescence microscope (Carl Zeiss Confocal microscopy LMS 700) employing phase contrast image. Bar represent 10 μm.

Stable arabidopsis transformation.
Columbia ecotype (Col-0) of Arabidopsis thaliana (L.) Heynh plants were transformed using floral dip method 92 . Germinated seeds were placed in vessels containing rock wool and embedded in a hydroponic medium. The plants were maintained at 25 °C in the growth chamber with a long day photoperiod regime (16 h light/8 h dark). In the T2 generation, lines with single T-DNA insertion site were selected based on the segregation of resistant and sensitive seedlings to glufosinate ammonium containing medium and verified by PCR with specific primers. Then, T3 homozygous lines for PrMADS10 construct were selected and used in the analyses.

RNA extraction and quantitative RT-PCR (RT-qPCR).
Total RNA was extracted from radiata pine seedlings following the procedure described by Le Provost 93 . Integrity of RNAs was checked on agarose gels stained with GelRed (Biotium Inc.), and their concentration determined by a ND-1000 UV spectrophotometer (Nanodrop Technologies, Montchanin, DE, USA). cDNA synthesis was performed using First Strand cDNA Synthesis Kit (Fermentas Life Science, Glen Burnie, MD, USA).
Total RNA was isolated from 35S::PrMADS10 transgenic Arabidopsis, using the SV total RNA isolation system (Promega). Primers for quantitative real time-PCR (RT-qPCR) were designed using Beacon Designer v 2.0 software (Premier Biosoft, Palo Alto, CA, USA). All primers used in this work are listed in Table 1. YBR Green/ ROX quantitative PCR (qPCR) Master Mix (2 × ; Fermentas Life Science) was used for all qPCR quantifications in a final volume of 20 μL following the manufacturer's protocol. All experiments were run on a real-time Mx3000P PCR detection system (Stratagene, Cedar Creek, TX, USA). The instrument was set to measure SYBR green dye fluorescence at the end of each cycle. Initial primer concentrations were 250 nM for all reactions, and the cDNA template for each sample was synthesized using 1 μg of DNase-treated total RNA using a first-strand cDNA synthesis kit (Fermentas Life Science) according to the manufacturer's instructions. The first-strand RT reaction product was diluted ten-fold, and 2 μL was used for each qPCR reaction. The cycle threshold (Ct) line was determined manually as the point where the R2 value for the standard curve reached its highest point 94 . Standard curves were determined in duplicate reactions from the dilution series of each amplicon. qPCR determinations were run in duplicate and values of each sample corresponded to a mean ± SE of three biological replicates. A melting curve analysis was performed for each set of primers in order to avoid non-specific amplification. The expression levels were normalized with the stable expression level of three housekeeping genes (Suppl. Table 1). The overexpression of the PrMADS10 gene was obtained from 3 biological replicates. Relative expression was calculated using primers for AtFbox, AtUbi10 95 and AtPP2 as normalizing genes. Data were analyzed using the methods derived from the algorithm of Vandesompele 96 . Two-way ANOVA-LSD post hoc was used to determine the main effects of inclination and time of inclination exposure effect for each gene using Statistica for Windows (v. 7.0; StatSoft, Tulsa, OK, USA). Significant differences were inferred at P ≤ 0.05. Differences in PrMADS10 transgenic expression, lignin and anthocyanin quantification was analyzed by one-way ANOVA-LSD post hoc and the significant differences inferred at P ≤ 0.05.
Gene expression using Affymetrix ATH1 microarrays. For microarray hybridizations, total RNA was processed using the GeneChip one-cycle target-labeling kit (Affymetrix). Biotinylated cRNA was synthesized from 5 μg of total RNA from Arabidopsis stems (three months old) using the Affymetrix IVT kit according to the manufacturer's instructions. cRNA was used to hybridize ATH1 GeneChip expression microarrays. Three un-transformed plants were used as controls and four transgenic lines were considered as biological replicate. Affymetrix data were normalized in R (http://www.r-project. org/) using RMA 97 . For detecting differentially regulated genes, normalized log2-transformed data were analyzed using Rank product statistic, as described before [98][99][100] . The data was processed with bioinformatics tools available at the VirtualPlant web site (http://www.virtualplant.org). (2019) 9:18981 | https://doi.org/10.1038/s41598-019-55276-7 www.nature.com/scientificreports www.nature.com/scientificreports/ For the heatmap, Expander 7.11 (http://acgt.cs.tau.ac.il/expander/) was used. K-means of 3 clustering were obtained from microarray data. The clusters were obtained with 50 max iteration used distance metric with Pearson correlation. The GO pie charts were obtained using PANTHER14.1 DB and graphic genic onthology was obtained using GOrilla 56 .
Microarray data. Raw signal intensity values were first normalized with RMA method using affy package in R language 101,102 and probes were mapped to Locus ids of Arabidopsis genome. Genes with at least 2-fold change and p-value < 0.05 in Rank product analysis were considered as differentially expressed. Microarray data was analyzed using K-means clustering from Expander 7.11 (http://acgt.cs.tau.ac.il/expander/).

Functional classification based on MapMan.
Gene expression data in a context of metabolic overview was visualized using MapMan (version 3.6.0RC) software 53,54 . MapMan uses a plant-specific ontology that classifies genes into well-defined hierarchical categories, denominated BINs and suBINS 103 . The diagram shows positively and negatively-regulated genes in red and blue, respectively. The data set obtained from microarray analysis (AraGene-1_0-ST) where compared to Ath_AGI_TAIR9_Jan2010 integrated in MapMan from TAIR.

STRING interaction network.
String is a database where known and predicted direct (physical) interactions, as well as, indirect (functional) interactions can be established based on co-expression, co-localization, or text-mining and others [104][105][106] . Differentially-expressed genes related to transcription factors and those genes involved in the synthesis of lignin were picked and the web server was interrogated in order to uncover potential protein-protein association networks. The database was interrogated for the last time on 28 th March (2019).
Anthocyanin and lignin content. Leaf and stem samples from three months-old Arabidopsis were ground in liquid nitrogen and extracted overnight in 1.0 mL of 1% (v/v) HCl in methanol at 4 °C following a previous report 107 . Relative anthocyanin levels were determined by measuring the absorbance at 530 nm of the aqueous phase 108 . Lignin was extracted as described by Campbell and Ellis 109 . Samples from whole Arabidopsis stems were diluted in 1 M NaOH (1/3, v/v) and hydrolyzed. A colorimetric assay was performed using thioglycolic acid (Sigma-Aldrich), and the absorbance was measured at 280 nm. The results were expressed as μg lignin per gram of fresh weight (FW).