Characterization of microRNAs Expressed during Secondary Wall Biosynthesis in Acacia mangium

MicroRNAs (miRNAs) play critical regulatory roles by acting as sequence specific guide during secondary wall formation in woody and non-woody species. Although thousands of plant miRNAs have been sequenced, there is no comprehensive view of miRNA mediated gene regulatory network to provide profound biological insights into the regulation of xylem development. Herein, we report the involvement of six highly conserved amg-miRNA families (amg-miR166, amg-miR172, amg-miR168, amg-miR159, amg-miR394, and amg-miR156) as the potential regulatory sequences of secondary cell wall biosynthesis. Within this highly conserved amg-miRNA family, only amg-miR166 exhibited strong differences in expression between phloem and xylem tissue. The functional characterization of amg-miR166 targets in various tissues revealed three groups of HD-ZIP III: ATHB8, ATHB15, and REVOLUTA which play pivotal roles in xylem development. Although these three groups vary in their functions, -psRNA target analysis indicated that miRNA target sequences of the nine different members of HD-ZIP III are always conserved. We found that precursor structures of amg-miR166 undergo exhaustive sequence variation even within members of the same family. Gene expression analysis showed three key lignin pathway genes: C4H, CAD, and CCoAOMT were upregulated in compression wood where a cascade of miRNAs was downregulated. This study offers a comprehensive analysis on the involvement of highly conserved miRNAs implicated in the secondary wall formation of woody plants.


Introduction
Recently, tremendous interest has been devoted to better understand the molecular basis of wood formation. Wood formation involves several events, from the differentiation and division of the xylem cells until the formation of heartwood tissues [1]. Elucidation of the mechanisms involved in producing secondary walls will help plant biologists to optimize the production of environmentally cost-effective and renewable alternative resources as wood has been widely used for pulping, timber industry and biofuel industry [2]. In the last decades, genomics studies in model species have revealed various wood associated genes playing an independent and additive role during secondary wall formation [3,4].
Genomics studies in herbaceous and woody species have identified pathway genes playing an interconnected role in lignin biosynthesis. These complexities reveal groups of transcriptional factors playing an independent and additive regulatory role in lignin biosynthesis [5]. Although important progress has been made into the identification of transcriptional factors that bind the AC elements of monolignol biosynthetic pathway genes, little is known about the transcriptional regulators that are responsible for turning on the secondary wall biosynthetic program [2]. This was hampered by the overlapping roles in the transcriptional regulation of lignin biosynthetic pathway with the biosynthesis of other secondary wall components [6]. Such a multilayered regulatory network affects multiple target genes and also points to the existence of switches that can be tuned to regulate different cell wall pathways in new spatial and temporal patterns [7].
In model plant species, studies have validated the roles of miRNAs in regulating lignin biosynthetic pathway genes. Transcriptional regulation alone is insufficient to provide an efficient fine-tuning of target gene expression in certain cells or tissues. Several studies have suggested that miRNAs play important roles in plant growth and development. More than 21,643 miRNAs have been identified from 168 species (miRBase release 18, November 2011; http://microrna.sanger.ac.uk;http://cgrb.orst. edu/smallRNA/db/). This large numbers of miRNAs identified using cloning, direct sequencing and computational prediction suggest that non-coding RNA play a prominent role in the regulatory networks of plant and animal development. Identification of a cascade of highly conserved miRNA across plant species suggests a common mechanism regulating vascular development in all plants [8]. Small RNAs like miRNAs act as posttranscriptional regulators by regulating the meristem cell differentiation and tissue patterning and have been shown to spatially regulate gene expression at various developmental stages in plants [9,10].
Deep sequencing of small RNAs has suggested that certain miRNA families are indeed conserved across diverse plant species. Although important progress has been made in characterizing all miRNAs in plants, limited information to justify their roles in various tissues and developmental stages has hampered our understanding of their roles. miRNAs have been reported to play important roles in plant growth and development [11]. In this study, we detected the expression of highly conserved amg-miRNAs in various wood forming tissues across different genotypes in A. mangium. Expression patterns of the highly conserved amg-miRNAs in different tissues suggested amg-miR166 as the potential regulatory sequence in vascular tissues. Subsequently, mapping of the amg-miR166 cleavage site via 59RACE and 39RACE were employed to isolate the gene targets. Bioinformatics pipeline analysis and expression studies in high and low lignin tissues were undertaken to unravel the regulatory roles of this amg-miR166. Our study also suggests miRNAs may play an important regulatory role in monolignol biosynthetic pathways despite the complexity of miRNA interaction in secondary wall formation.

Results
Expression profile of highly conserved miRNA families in various tissues in A. mangium Despite the importance of small RNAs in plant growth and development, information on the roles of small RNA in Acacia mangium is still in its infancy. Unlike Arabidopsis and poplar, potential roles of highly conserved miRNA families in various wood forming tissues in A. mangium are unclear. From the expression profiles of small RNAs generated between low lignin Am54 and high lignin Am48 [12], we have identified six highly conserved miRNA families that are potentially important regulatory sequences in secondary wall formation (Table S1). We report here the expression level of these six miRNAs in various wood forming tissues across different genotypes of four years A. mangium using quantitative real time -PCR.
Altogether, we have identified six amg-miRNA families with potentially important regulatory roles in secondary wall biosynthesis from their expression profile in high lignin compression wood and low lignin tension wood across three different genotypes ( Figure 1). Secondary cell walls are composed of cellulose, hemicellulose and lignin. This novel discovery is consistent with the findings reported earlier across three 2-year-old A. mangium genotypes: NSW22, AVW22 and WMH16 [12]. Hence, these highly conserved amg-miRNA families may be important regula-tors in secondary wall formation since genotype and age differences have minimal impacts on their expression patterns.
Our findings suggest the involvement of a well conserved regulatory network in secondary wall forming cells with some miRNA families exhibiting dual roles. This is evident as amg-miR172, amg-miR156, amg-miR168, amg-miR159 and amg-miR394 exhibited more or less similar expression between phloem and xylem wood tissues (Figure 1). Although the five amg-miRNAs showed differences between compression and tension wood (Figure 1), the presence of these amg-miRNA in phloem tissues indicates their potential involvement in phloem and xylem development. Only amg-miR166 exhibited strong downregulation in xylem tissues which suggests its roles in xylem formation and lignin regulation for further validation.
Identification of key regulatory sequences in vascular tissue development through quantitative real time-PCR To improve our understanding on the complexity in the regulatory network in vascular tissues, we used quantitative realtime PCR technology to track the changes in miRNA expression in various wood forming tissues. We selected one miRNA with highest expression among the members in each family from six different families (Table S1) and examined their expression in leaf, phloem, compression wood and tension wood. We compared the ct value of the selected miRNA member from each family.
Of the six different miRNA members, we found that amg-miR166 has the highest ct value in compression wood and tension wood (Table 1). Ct values are inversely proportional to the amount of target nucleic acid in a sample. Every increase in 3.33 ct value represented a tenfold decrease in concentration of the investigated target nucleic acid. Phloem tissues showed a -1000 increase in amg-miR166 concentration as compared to xylem tissues ( Figure 1) implying the strong downregulation of this miRNA from phloem to xylem tissue. Moreover, near absence of this miRNA in compression wood and tension wood as compared to other amg-miRNAs using quantitative real time PCR requires in-depth studies to verify the role of this miRNA in xylem formation and lignin regulation in A. mangium. This expression pattern is similar to previous studies reporting miR166 playing an important role in vascular development and lignin regulation [8,12,13,14,15,16,17].

Experimental and binding energy characterization of genes regulated by miR166
To strengthen our findings and to understand how miRNAs regulate their target genes [11,18,19], we performed 39-and 59-RACE using various tissues. Pooled xylem tissues (from compression wood and tension wood) were used to obtain different members of HD-ZIP III regulated by amg-miR166. We were unable to obtain any HD-ZIP III members via 59RACE mapping of amg-miR166 cleavage site in leaf and phloem tissues. RACE PCR products revealed nine mRNA are indeed the gene target of amg-miR166 ( Figure 2). All the amg-miR166 binding site are located at the 59UTR of the characterized HD-ZIP III mRNA. Using psRNA Target Analysis, we found that all our HD-ZIP III mRNA will be cleaved by amg-miR166 with scores of 3.0 and below regardless of the polymorphism in the sequences ( Figure 3). 59 RACE mapping of miRNA cleavage site revealed nine different members of HD-ZIP III transcription factors were successfully characterized in this study. We obtained five novels HD-ZIP III mRNA out of the nine members ( Figure S1). HD-ZIP III (1) shared 91% in nucleotide level and 89% in protein level to Glycine max HD-ZIP III (REVOLUTA) transcription factors with e-value of 0. HD-ZIP III (2) shared 92% in nucleotide and protein level to Medicago truncatula HD-ZIP III (ATHB15 group) transcription factors with e-value of 0. HD-ZIP III (3) shared 91% in nucleotide and 92% in protein level to Medicago truncatula HD-ZIP III (ATHB15 group). HD-ZIP III (4) shared 89% in nucleotide and 88% in protein level to Populus trichocarpa class HD-ZIP III (ATHB8 group). Blast analysis at nucleotide level between HD-ZIP 2 and HD-ZIP 3 revealed that both of them share 87% similarity in nucleotide and 92% in protein level (data not shown). Although Arabidopsis genome contains five class III HD-ZIP genes (ATHB8, PHAVOLUTA/ATHB9, PHABULOSA/ATHB14, CORONA/ATHB15 and REVOLUTA/IFL1) which are targeted by miR166, our study revealed only three different class III HD-ZIP are present in A. mangium. In summary, blast analysis at

Real-time PCR validation of HD-ZIP III transcription factors
To pinpoint the role of each HD-ZIP III transcription factor in wood forming tissues, we analyzed their expression level in various wood forming tissues. Preferential expressions in each of the different transcription factors in various tissues suggest a potential role in xylem formation. Interestingly, the four different HD-ZIP III transcription factors were abundantly expressed in compression wood and their expression was downregulated in tension wood, tissue characterized by lower lignin content. They were weakly expressed in phloem tissues and absent in leaf tissues. From the expression data of these four different HD-ZIP III as presented in Figure 4, we suggest that amg-miR166 play a pivotal role in maintaining the abundancy of the various members of HD-ZIP III transcripts through their miRNA cleavage site and also via epigenetic silencing. Hence, the study presented here strongly suggests vascular tissue development is very complex and require the involvement of various members from a particular family ( Figure 3) which are in agreement with other studies [8,9,15].

Sequence comparison of conserved and novel miRNA target
As both miRNAs and their corresponding target site are highly conserved across a wide range of species [20], information on the homologous miRNA and mRNA dataset will provide useful information on conserved regulatory roles even in different species [21]. To prove whether the homologous miRNA and mRNA in other species have a similar sequence conservation and binding energy with our species, we isolated the mRNA target of amg-miR166 using 59RACE and 39RACE approach. We chose amg-miR166 in this study because this miRNA showed the greatest difference in its expression between phloem and xylem tissue compared to the other miRNAs. Moreover, it had the lowest abundance in compression wood and tension wood among the various amg-miRNA families chosen for validation in this study. This amg-miR166 is the first plant miRNA identified and is found in all plant species. Bioinformatics analysis using psRNA Target Analysis Server between amg-miR166 with HD-ZIP III from ten different species indicated that variations in the miRNA: mRNA target sequences only occur at the 39 end of miRNA ( Figure 5). We verified here the four isolated HD-ZIP III mRNA (HD-ZIP III (1)(2)(3)(4) share high similarity with other species at nucleotide and protein levels ( Figure S1). In addition to these four different HD-ZIP III mRNA, we also obtained five novel amg-miR166 targets cleavage site (HD-ZIP III (5-9) ( Figure S1).
Expression of C4H, CCoAOMT and CAD genes in high lignin compression wood and low lignin tension wood Previous investigation has identified six genes involved in lignin biosynthetic pathway in Acacia mangium x Acacia auriculiformis hybrid: caffeate O-methyltransferase (COMT), cinnamoyl-coenzyme A reductase (CCR), cinnamyl-alcohol dehydrogenase (CAD), phenylalanine ammonia-lyase (PAL), cinnamate 4-hydroxylase (C4H) and caffeoyl-coenzyme A O-methytransferase (CCoAOMT) [22]. C4H, CCoAOMT and CAD are three enzymes that play an important role in the beginning, middle and the end of the lignin biosynthetic pathway [3]. C4H, CCoAOMT and CAD have been isolated and characterized in A. mangium x A. auriculiformis hybrid [23]. In this study, downregulation of these six highly conserved miRNAs in high lignin compression wood as presented in Fig. 1 was accompanied by the upregulation of the C4H, CCoAOMT and CAD genes (Figure 4). Generally lignin content in compression wood in the different A. mangium trees was higher than in tension wood ( Table 2). As the expression of major pathway genes involved in secondary wall biosynthesis are highly coordinated and could be switched on and off [16,24,25,26,27,28,29,30,31,32,33,34], upregulation and downregulation of miRNA affects not only lignin content but also cellulose and/or hemicellulose levels.

Characterization of amg-miR166 Precursors
Although miR166 has been well studied in various species, no studies have explored the biogenesis of this miRNA in A. mangium. Biogenesis of miRNA involves several enzymes and processes [35,36,37,38,39]. In order to address this question, we first obtained the precursor structure of all the isoforms belonging to amg-miR166 family ( Figure 6). Interestingly, the precursor structures varied greatly among all the members of this family. Jones-Rhoades and Bartel [40] reported that miRNA: miRNA* portion of the pre-miRNAs are more conserved than other parts of the precursor's structure even within the same miRNA family. In agreement with Jones-Rhoades and Bartel [40], four out of the five amg-miR166 precursors' structures in this study shared more conserved portion in their miRNA: mRNA* structure.
Several studies in plant genome have proven that similar mature miRNA may arise from multiple precursors. As shown in miRBase, certain isoforms of mature miR166 in Arabidopsis and Poplar were encoded by more than five different precursor structures. Since Acacia genome is not available yet, we are unable to determine the complete precursor structure of this amg-miR166 family. As structure of miRNA precursors may influence the rate of miRNA production [40], identification of the complete set of all amg-miR166 precursor structures is important for the selection of the right amiRNA precursor that can best express mature miRNA. The precursor structure of amg-miR166 obtained in this study fulfilled several parameters; have negative folding free energies [41], MFEI threshold above 0.7 [42], predicted mature miRNA sequence differed from their homologous miRNA by only (0-3) mismatches [42], no large break or loop in the miRNA sequence [42] and range from 90 nucleotides to more than 210 nucleotides ( Figure 6). The precursor structure of the different amg-miR166 members obtained can be used for the development of artificial miRNA (amiRNA) gene construct for silencing of HD-ZIP III transcription factors.

Discussion
Regulation of class III HD-ZIP transcription factor via posttranscriptional and transcriptional gene silencing The first plant genes identified to be the target under miRNA regulation are HD-ZIP III transcription factors [43]. Bioinformatics analysis using psRNA server indicated all the five different members belonging to the amg-miR166 family in our study will bind to the nine different HD-ZIP III transcripts and direct its cleavage (Figure 3). We demonstrate that most of HD-ZIP III transcripts in compression wood and tension wood are safe from low abundance of amg-miR166 cleavage. No direct evidence of HD-ZIP III transcripts cleavage by amg-miR166 in leaf and phloem tissues at the expected site was found. However, the  presence of high level of CpG island in the HD-ZIP III transcripts ( Figure S1) could contribute to high level of methylation which have been associated with low rates of transcription of miR166 targets [44]. Hence, amg-miR166 was suggested to play a critical regulatory role via RNAi and epigenetic silencing in controlling gene expression in secondary wall formation in A. mangium. We show that expression of amg-miR166 was corroborated with the expression of HD-ZIP III transcripts in this study (Figure 1 and 4). For instance, upregulation of miR166 in leaf was accompanied with strong downregulation of their transcripts and similar expression also reported in phloem, compression wood and tension wood in this study.
HD-ZIP III transcription factors play an important role in xylem development. Overexpression of ATHB8 (subgroup of HD-ZIP III) gene resulted in precocious proliferation of xylem tissue while absence of vascular phenotype in their knockout cannot be explained by redundancy among the HD-ZIP III members [8,15]. Overexpression of sense and antisense ATHB15 (subgroup of HD-ZIP III) transcripts resulted in moderately dwarfed and severely dwarfed transgenic plants which may due to the suppressed vascular system [9]. ATHB15 negatively regulate cell differentiation and is critical for vascular development [9]. Xylem formation are no longer detected in root when all the five HD-ZIP III transcription factors are knocked out as this group of transcription factor drives the de novo xylem formation [45,46]. In agreement with Carlsbecker et al. [45] and Cano-Delgado et al. [46], our study suggest that amg-miR166 may play a potential role in xylem formation as the expression of this miRNA family was strongly downregulated in xylem tissue, accompanied by the upregulation of HD-ZIP transcription factors (Figure 1 and 4). In the March issue of Plant & Cell Physiology, Zhong and Ye [15] reported that overexpression of miR165 leads to the downregulation of all the five HD-ZIP III transcript genes. Although HD-ZIP III transcripts are regulated by miR165 and miR166 [15], our deep sequencing study reported the absence of miR165 member in low lignin Am54 and high lignin Am48 [12].

Spatial and temporal expression of highly conserved miRNA families
Our results revealed the involvement of a cascade of miRNAs in the regulatory network in secondary wall formation in A. mangium. Experimental analysis suggests that temporal expression of amg-miR168, amg-miR159, amg-miR172, amg-miR156 and amg-miR394 in compression wood and tension wood indirectly assisting the role of amg-miR166 (Figure 1). These six amg-miRNA members demonstrated obvious differences in the intensity between high lignin compression wood and low lignin tension wood. However, only amg-miR166 demonstrated obvious differences in the intensity between phloem and xylem tissues, suggesting its roles in xylem development ( Figure 1). As flowering plays a pivotal role in xylem expansion [47], differential expression of the miRNA members in this study annotated with flowering properties (amg-miR159, amg-miR156 and ang-miR172) suggest its complexity in plant growth and development. Secondary wall formation is very complex as it involves various groups of regulatory network with some families playing dual roles by regulating the synthesis of cellulose and lignin simultaneously. These findings suggest that the five different amg-miRNA play dual roles in secondary wall formation.
Lu et al. [11] reported that a large number of highly conserved miRNA families exhibit contrasting species-tissue-specific expression patterns in Arabidopsis and Populus despite of sequence conservation between these two species. In this study, amg-miR156, amg-miR159, amg-miR172, amg-miR394, amg-miR168 and amg-miR166 have the most conspicuous expression in tension wood compared to compression wood ( Figure 1). These miRNA families were also abundantly expressed in tension wood compared to compression wood [11]. Although miR156 family showed tissue specific expression level in Populus, Arabidopsis and in this study, however, this miRNA family was expressed to a similar level in all tissues in apple [13]. Expression of ptr-miR159 was significantly higher in compression wood than in tension wood with no preferential differences in expression between phloem and xylem tissue [11]. In apple stem, this miRNA family was weakly expressed in xylem tissue compared to phloem tissue [13].
We hypothesized that miRNAs play an independent and additive role during secondary wall formation with the expression of some miRNAs indirectly affecting the expression level of other downstream miRNAs. It is evident from the expression pattern of various miRNA members chosen for validation in this study ( Figure 1). Spatial and temporal expression of these miRNA members in leaf, phloem, compression wood and tension wood tissue is likely to fine-tune the expression of target transcripts beyond a certain level (Figure 1, 3 and 4). For instance, amg-miR166 was highly expressed in leaf and phloem tissue and only mild expression was detected in xylem tissue. These differential expression patterns of amg-miR166 in different tissues suggest downregulation of this amg-miR166 family being important for HD-ZIP III transcription factors playing a role in xylem development and indirectly regulate the expression of various pathway genes involved in lignin biosynthesis. Studies by Ohashi-Ito et al. [17] revealed that overexpression with a mutation in the START domain of Ze-HB12 (HD-ZIP III, zinnia homologue of Arabidopsis REVOLUTA) triggered the upregulation of four genes encoding enzymes involved in lignin monomer synthesis. Du et al. [16] reported overexpression of a miRNA-resistant poplar HD-ZIP III (POPCORONA, poplar homologue of ATHB15) resulting in the up-regulation of genes involved in cellulose biosynthesis while lignification and metabolite related genes (putative laccase, cinnamoyl CoA reductase, chalcone synthase, putative pectin methyltransferases and pectinesterase) are downregulated. Overexpression of ATHB8 in Arabidopsis promotes vascular cell differentiation and accompanied with an increased production of lignified tissues [8]. Our study suggest the involvement of other highly conserved amg-miRNAs which indirectly regulate the expression of various pathways genes involved in secondary wall biosynthesis from their expression patterns in compression wood and tension wood. Expression of lignin pathway genes are reported to be regulated by various transcription factors [16,17,48,49,50,51,52].

MicroRNA plays an important regulatory role in secondary wall formation
Lignin formation involves complex interaction between various groups of transcription factors with the monolignol pathway genes. Various groups of transcription factors as MYB, NAC, LIM and HD-ZIP III were reported to independently regulate the expression of various pathway genes in lignin biosynthesis [6,8,9,52,53]. Studies in models species have validated the roles of highly conserved miRNAs with transcription factors target in regulating the monolignol biosynthetic pathway genes which suggest that the regulatory mechanism in vascular development and lignin biosynthesis might be conserved in diverse plants species [54]. This effect was seen very clearly in the expression profile of various amg-miRNA families in various tissues (Figure 1), low lignin Am54 and high lignin Am48 (Table S1).
Better understanding of the biological functions of miRNAs requires knockdown of their mRNA target. Various knockdown studies have revealed that both miRNA and their target genes exhibited negative correlation [55,56,57] and this implies that expression of miRNA is indeed an indicator of their target abundance. In line with this, our analysis proves that a single miRNA can regulate the stability of several different transcripts from a particular family (Figure 3). In vitro assay of Arabidopsis miRNA gene target via bioinformatics approach indicated that distantly related transcripts can have the same miRNA target site demonstrating expression of miRNA being best to reflect the regulatory events in secondary wall formation ( Figure S2). In Arabidopsis, all the class III HD-ZIP members display overlapping and distinct roles in vascular development [58,59,60,61] although all these different members were regulated by miR166. Interestingly, 59RACE mapping of our amg-miR166 target cleavage site revealed this miRNA family has nine different transcripts.
Although MYB are the most investigated transcription factor in the regulation of lignin biosynthetic pathway genes, antagonistic and additive roles of the various members belonging to this family have added additional complexity to the transcriptional regulation of monolignol biosynthetic pathway. For instance, AtMYB58, AtMYB63 and AtMYB85 are identified as the transcriptional activators of lignin biosynthesis in vessels and fibers [2,62] while AtMYB103, in cellulose biosynthesis [2]. These antagonistic roles were complicated with the discovery of other MYB members like AtMYB52, AtMYB54, and AtMYB69 that regulate the entire cellulose, lignin and xylan biosynthesis [2]. In case of miR159 family, we showed that six different MYB members are regulated by amg-miR159 ( Figure S2).
Our analysis on the expression of these various amg-miRNA families and indirect in vitro assay via bioinformatics approach suggests miRNA being involved in fine tuning the expression of its transcripts via RNAi and epigenetic silencing. To gain more insight into this mechanism, we have shown in our RT-qPCR the miRNA members and their corresponding HD-ZIP III gene targets were inversely expressed in various tissues (Figure 1 and 4). Downregulation of these amg-miRNA members in compression wood is critical for transcriptional regulation of various pathway genes involved in secondary wall biosynthesis. This was evident from expression of C4H, CCoAOMT and CAD genes in compression wood, a region enriched with lignin content in this study. Complexity in the regulatory network in monolignol biosynthetic pathway Wood (secondary xylem) formation involves various stages of division zone, expansion zone, maturation zone and programme cell death zone [63]. Various transcription factors are implicated in the regulation of cell division, differentiation, expansion and patterning of vascular tissue during secondary wall formation [64,65,66,67,68,69]. This effect was seen very clearly in the expression profile of various amg-miRNA families in various tissues (Figure 1), low lignin Am54 and high lignin Am48 [12]. Moreover, this observation raise important questions on the complexity involved in lignin regulation as various groups of short small RNAs with very different expression level were reported between low lignin Am54 and high lignin Am48 (data not shown).
For comparison, Oh et al. [70] showed that transcription factors were highly expressed in xylem tissue which explains the events where secondary xylem formation occurs. Xylem-abundant MYB proteins are involved in the transcriptional regulation of secondary xylem formation from their upregulation in the xylem tissue compared to the bark tissue [70,71]. In contrast to class I KNOX genes in negatively regulating the differentiation of cambial cells, transcription factors like NAC domain promote differentiation of vessel elements [6,72,73]. Hence, it is obvious that regulatory network in monolignol biosynthetic pathway involves a complex interaction of various groups of transcription factors.
Since the lignin biosynthetic pathway genes are more or less highly conserved across all plant species, it is likely a common mechanism regulating vascular tissue development in all plants [10]. This is evident from studies reported on the roles of miR166 in vascular differentiation and development via HD-ZIP III genes; a possible general rule in all vascular plants [9]. In addition, miR166 and its target genes HD-ZIP III are ancient and the most well studied small RNA in plants [8,58,59,60,74,75]. In agreement with many other studies, our findings suggest amg-miR166 play a potential roles in xylem development. Moreover, members of this amg-miR166 family regulate different HD-ZIP III transcripts during xylem development. In addition to amg-miR166, our study suggests secondary wall biosynthesis may involve other highly conserved miRNA families like amg-miR394, amg-miR159, amg-miR156, amg-miR172 and amg-miR168.

Inconsistency in findings on alteration of the pathway genes
Although much attention has been devoted into isolation and characterization of all the genes involved in monolignol biosynthetic pathway, many fundamental questions about their regulation remain unanswered. Different laboratories reported conflicting results on alteration of the pathway genes even within a similar species. These complexity raise questions on our understanding in the regulatory network involved in monolignol biosynthetic pathway as various groups of transcription factors have been shown to regulate the pathway genes.
System biology has documented that altering the pathway genes often modifies its interaction with the environment [76,77]. These were true from the contrasting reports on the effects of downregulation of 4CL via antisense RNA on the growth performance. Hu et al. [25] reported greenhouse evaluation of suppressed 4CL activity exhibiting up to a 45% reduction in lignin and these were compensated with substantially enhanced leaf, root and stem growth in transgenic poplar. Subsequent studies found no growth enhancement using transgenes suppressing the expression of 4CL in Chinese white poplar (P. tomentosa) [78] and greenhouse grown aspen [79,80,81]. However, 4CL downregula-tion in poplar grown under field condition resulted in considerable variation in productivity [82]. In agreement with Pilate et al. [83] and Leple et al. [84], Voelker et al. [82] reported that high level of lignin reduction observed in approximately one-third of the transgenic events led to reduced growth and serious physiological abnormalities.
Since different research groups reported inconsistent findings on the effects of downregulation of pathway genes even within same model species, our study suggest better understanding of the roles of miRNA is critical as miRNA were reported to regulate various transcripts with overlapping roles, which in return regulate the monolignol biosynthetic pathway genes. Thus, we believe that the findings presented in this study will add more blocks to the network and expand the knowledge on the underlying regulatory details as the major pathways in lignin biosynthesis are more or less conserved across various species.

Plant materials
Plant materials for gene expression analysis and lignin content determination were used from three 4-year-old trees of A. mangium from Plot A, Plant Biotechnology Centre, Universiti Kebangsaan Malaysia, Bangi, Malaysia. Three seed sources (NSW19, NSW20 and ERC22) were used in this investigation. Four different types of tissue-leaf, phloem, compression wood and tension wood-were used in this investigation. Compression wood was taken from the lower portion of the bend and tension wood was collected from the upper portion of the bend. Phloem tissues were taken from the outer layer of the inner bark tissues of the main stem.

Quantitative real time PCR of selected miRNAs
Quantitative real time PCR was performed using iCycler iQ TM Real Time PCR Detection System (Biorad, USA). From the expression of 12 highly conserved amg-miRNAs between low lignin Am54 and high lignin Am48 [12], we have identified six highly conserved amg-miRNAs (Table S1) play an important regulatory roles in secondary wall formation. Expression analysis using semi-qPCR and rt-qPCR indicated the remaining six amg-miRNAs are lowly expressed and no obvious differences were observed between high lignin compression wood and low lignin tension wood ( Figure S3 for semi-qPCR, method unpublished). From here, one miRNA from each family with strong differences in expression between low lignin Am54 and high lignin Am48 (Table S1) were selected for further validation in various tissues using iQ TM SYBR Green Supermix (Biorad, USA) among three different individuals. Total RNA containing small RNA was isolated from four different tissues using miRVana miRNA Isolation kit (Ambion, Austin, TX, USA) and treated with RNase-free DNase (Qiagen, Germany). Isolated Total RNA was analyzed for their integrity using RNA 6000 Nano kit (Agilent Bioanalyzer, USA) and Nanodrop ND-100 Spectrophotometer. Only Total RNA with RIN value above 7.5 was selected for the first strand cDNA synthesis. Total RNA containing small RNA was reverse transcribed using miScript Reverse Transcription kit (Qiagen, Germany). 5.8 S rRNA was used to normalize the cDNA concentration in various tissues. In each reaction, 15 ml consisting of 6 ml of 26SYBR Green Supermix (Biorad, USA), 1 ml of cDNA template equivalent to 100 pg total RNA, and 2 ml of 106 miScript Universal Primer (Qiagen, Germany) and 5 pmol of forward primers. PCR amplification protocol was as follows: 95uC for PCR initiation activation step, 45 cycles consisting of 15 s at 94uC, 30 s at 57uC and 30 s at 70uC and final extension of 4 min at 70uC.

Lignin content determination
The lignin content of compression and tension wood tissues from three different A. mangium trees were determined using pyrolysis gas chromatography mass spectrometry (GC-MS) method. Each analysis was done in triplicates. Dried wood samples were ground to fine powder using a electric blender. Extractive content in wood meal was removed by refluxing them with 2 volumes of toluene and 1 volume of ethanol for 8 hours at 6 volts. Later on, the wood meal was refluxed with ethanol for 4 h to remove any residue of toluene. Pyrolysis GCMS was performed using GCMS G1374A (Agilent, USA). Compounds were identified by their mass spectra and retention time through comparison with those reported in literature [85,86]. Total lignin was calculated by adding all the S, G and H units.
miRNA: mRNA target prediction miRNA: mRNA duplexes of all the isoforms belonging to amg-miR166 were characterized using psRNA Target Analysis Server [87]. In these analyses, scoring systems for mismatches in the miRNA: mRNA were as follow; 0 score assigned to a complementary pair, G:U wobble assigned 0.5 score, non-G:U wobble mismatches assigned 1 score and finally 2 scores assigned to each bulged nucleotide [40,87,88,89]. miRNA: mRNA duplex with scoring of 3.5 or less were considered high confidence [40,87,88,89]. Nine different HD-ZIP III mRNAs isolated using 59RACE were characterized for their binding energy with the five different members of amg-miR166. In addition, HD-ZIP III sequences of ten different species were downloaded from National Center for Biotechnology Information (NCBI) and used to predict their binding energies with amg-miR166. psRNA Target Analysis and binding energy were conducted with default settings. Subsequent miRNA:mRNA analysis was done using various members with overlapping roles from Arabidopsis with each different highly conserved amg-miRNA (amg-miR159, amg-miR156, amg-miR168, amg-miR172 and amg-miR394).

and 59 Rapid Amplification of cDNA Ends (RACE) for mapping of amg-miR166 cleavage site
Aliquots of total RNA from compression wood and tension wood from three 4-year-old individuals of A. mangium, pooled prior to first strand cDNA synthesis. First strand cDNA was constructed using SMART TM RACE cDNA Amplification Kit (CLONTECH Laboratories, Inc. USA) following protocol instruction. 39 RACE was first employed to obtain the gene target of amg-miR166. In the 39RACE PCR amplification, primer sequence flanking the amg-miR166 target site was designed (Table S2). About 20 clones were randomly selected and were sequenced. To amplify the 59 ends of the target miRNA, two amplifications were performed using the 39RACE sequence (Table S2). The outer primer sequence of 39RACE was used in the first 59RACE PCR amplification. The 59RACE amplified products were purified and re-amplified using the inner 39RACE sequence in the following 59RACE PCR amplification (Table S2). In addition, amg-miR166 sequence ID 256(21) was used to obtain the novel HD-ZIP III (HD-ZIP (5-9) in 59 mapping of amg-miR166 cleavage site. The amplified products were gel purified, cloned into pGEM-T vector, sequenced and analyzed. To obtain the full lengths of the highly conserved HD-ZIP III mRNA (HD-ZIP (1-4), forward and reverse primers from each of the different HD-ZIP III mRNA were designed from the partial RACE HD-ZIP III sequence. First strand cDNA was synthesized using SMART TM RACE cDNA Amplification Kit (CLONTECH Laboratories, Inc. USA) and was used to obtain the full length of the different HD-ZIP III sequences in the three different A. mangium trees. 59 RACE mapping of the cleavage site were then employed in each of trees using leaf, phloem, pooled compression wood and tension wood tissues.

RT-qPCR of HD-ZIP III transcription factors and key lignin genes
Quantitative real time -PCR was performed using iQ TM SYBR Green Supermix (Biorad, USA) in iCycler iQ TM Real Time PCR Detection System (Biorad, USA) in a final volume of 15 ml. PCR primers (Table S3) were designed from the RACE sequence. The PCR protocol employed was the same as in quantitative real time PCR analysis of the selected miRNAs. Actin gene was included as an internal control for normalization in cDNA amounts used. The reactions were performed in triplicate for each run with three biological replicates. Melt curve analysis was conducted to determine the amplification specificity and the amplified product was further analyzed in 2% agarose. Relative normalized expression of the four different HD-ZIP III and C4H, CCoOAMT and CAD genes in leaf, phloem, compression wood and tension wood were compared.

Precursor amg-miR166 isoforms characterization
We designed primers flanking the conserved region of mature miRNA and the opposite arm from precursor structure sequence available in miR Base (http://microrna.sanger.ac.uk/). Arabidopsis and Poplar miRNA precursor databases from miR Base were used in the primer design. In addition, we used all five different amg-miR166 sequences as forward and reverse sequences to obtain a variety of precursors structure as miRNA and miRNA* often have quite perfect matches. PCR was conducted using DNA from A. mangium NSW19 with amplification protocol as follows: 95uC for PCR initiation activation step, 35 cycles consisting of 45 s at 94uC, 30 s at 50uC and 30 s at 70uC and final extension of 4 min at 70uC. Primers used are listed in Table S4. For each PCR, 20 ml reaction volumes consist of 1 U/ml i-Taq DNA polymerase (intron), 2 ml 106 PCR buffer (intron), 2 ml dNTP mixture (intron), 1 ml of DNA equivalent to 100 ng and 2 ml mixture of 10 mM forward and reverse primer respectively. PCR products were gel purified, cloned into pGEM-T vector and sequenced. The sequences obtained were analyzed for their ability to fold into hairpin structures using Mfold 3.1 [90,91] with default settings.

Conclusion
Our study critically demonstrated the roles of highly conserved miRNAs in secondary wall formation. From the expression pattern of six different highly conserved amg-miRNA families in various tissues, only amg-miR166 was strongly downregulated in xylem tissue compared to phloem tissue across six different genotypes. This phenomenon suggests that amg-miR166 might play a potential role in xylem development and indirectly regulate the genes involved in lignin biosynthesis in Acacia mangium. Moreover, this miRNA family was found in all vascular plants. Since different research groups have reported inconsistent findings on the effects of downregulation of pathway genes even within same model species, this study provides a better understanding of the overlapping roles of miRNA in regulating the pathway genes involved in secondary wall formation. Table S1 The 6 highly conserved plant miRNA families with strong differences in the expression level in each of the isoforms between low lignin Am54 and high lignin Am48. (DOC)