Transcriptome profiling of Pinus radiata juvenile wood with contrasting stiffness identifies putative candidate genes involved in microfibril orientation and cell wall mechanics

Background The mechanical properties of wood are largely determined by the orientation of cellulose microfibrils in secondary cell walls. Several genes and their allelic variants have previously been found to affect microfibril angle (MFA) and wood stiffness; however, the molecular mechanisms controlling microfibril orientation and mechanical strength are largely uncharacterised. In the present study, cDNA microarrays were used to compare gene expression in developing xylem with contrasting stiffness and MFA in juvenile Pinus radiata trees in order to gain further insights into the molecular mechanisms underlying microfibril orientation and cell wall mechanics. Results Juvenile radiata pine trees with higher stiffness (HS) had lower MFA in the earlywood and latewood of each ring compared to low stiffness (LS) trees. Approximately 3.4 to 14.5% out of 3, 320 xylem unigenes on cDNA microarrays were differentially regulated in juvenile wood with contrasting stiffness and MFA. Greater variation in MFA and stiffness was observed in earlywood compared to latewood, suggesting earlywood contributes most to differences in stiffness; however, 3-4 times more genes were differentially regulated in latewood than in earlywood. A total of 108 xylem unigenes were differentially regulated in juvenile wood with HS and LS in at least two seasons, including 43 unigenes with unknown functions. Many genes involved in cytoskeleton development and secondary wall formation (cellulose and lignin biosynthesis) were preferentially transcribed in wood with HS and low MFA. In contrast, several genes involved in cell division and primary wall synthesis were more abundantly transcribed in LS wood with high MFA. Conclusions Microarray expression profiles in Pinus radiata juvenile wood with contrasting stiffness has shed more light on the transcriptional control of microfibril orientation and the mechanical properties of wood. The identified candidate genes provide an invaluable resource for further gene function and association genetics studies aimed at deepening our understanding of cell wall biomechanics with a view to improving the mechanical properties of wood.


Background
Wood cell (such as tracheids and fibres) initials are produced by the vascular cambium and subsequently undergo cell expansion, primary cell wall biosynthesis, secondary wall deposition, lignification, and finally programmed cell death [1,2]. In perennial woody plants, secondary xylem (wood) is derived from the annual activity of the vascular cambium, with wood laid down in different seasons and tree ages frequently having distinct mechanical properties [3][4][5], particularly in gymnosperms [5]. The mechanical properties of secondary xylem not only provide support for woody plants to maintain their shape, resist maturation stress (gravity), and respond to various environmental forces (wind, snow, etc); they also affect the suitability of wood for different commercial applications. A number of factors influence the mechanical properties of wood, including individual cell walls, anatomical structure, cell-cell adhesion and cell turgidity [6]. The mechanical properties of plant cell walls also play a crucial role in cell expansion [7,8], tissue or organ morphogenesis [9] and responses to various signals [10,11]. Understanding wood biomechanics at the cell and organ levels will provide unique insights into many biological processes in plants, such as cell wall biosynthesis and wood formation.
The mechanical properties of plant cell walls and organs are largely controlled by the architecture of the cytoskeleton [6,10]. Of the three main types of cytoskeleton polymers, microtubules are the stiffest while actin filaments are the least rigid [10]. Tethering of galactose residues in xyloglucans to cellulose microfibrils is essential for mechanical strength such as tensile properties of primary cell walls [6,12] and expansion relies on cooperation between specific expansins and XETs [13]. Previous studies have showed that cellulose microfibril orientation dominates mechanical properties in both primary and secondary cell walls of xylem [14][15][16][17], particularly microfibril angle (MFA) in the middle layer (S2) of secondary cell walls. In Eucalyptus, MFA alone has been estimated to account for 86% of the variation in wood stiffness with an additional 10% to be influenced by wood density [16].
The complex nature of the genetic control of wood properties was first showed by quantitative trait loci (QTL) mapping. Several QTLs have been identified that contribute to MFA and wood stiffness in pines [18][19][20], eucalypts [21][22][23] and poplar [24]. Recently several individual genes have been found to play roles in the regulation of microfibril orientation and mechanical properties of xylem cell walls. A b-tubulin gene has been found to influence cellulose microfibril orientation in wood cell walls of Eucalyptus [25]. Two fasciclin-like arabinogalactan proteins (FLA11 and FLA12) have been found to affect MFA and tensile stiffness in eucalypts and Arabidopsis by altering cellulose deposition and the integrity of the cell wall matrix [26]. Association studies have recently showed single nucleotide polymorphisms (SNPs) in a number of cell wall-related genes that influence MFA in loblolly pine (a-tubulin, COMT2, CCR1) [27], radiata pine (RAC13, SuSy) [28], white spruce (glycosyl hydrolase 10) [29] and Eucalyptus nitens (CCR) [30]. Furthermore, associations between allelic variation and other mechanical properties of wood have also been observed, including stiffness (COMT) and density (PAL1 and PCBER2) in radiata pine [28] and specific gravity (CAD and SAMS2) in loblolly pine [27]. While the research described above provides important insights into the molecular mechanisms of microfibril orientation and mechanical properties of wood cell walls, genes involved in the regulation of cell wall mechanics remain poorly characterised at the transcriptome level.
Radiata pine (Pinus radiata D.Don) is the most important conifer species in commercial forest plantations in Australia and several other countries and is grown primarily for structural timber. To gain further insight into the molecular control of microfibril orientation and cell wall mechanics, cDNA microarrays were used to investigate differential gene expression in radiata pine developing xylem with contrasting stiffness and MFA. Radiata pine trees with contrasting wood stiffness were selected in two progeny trials using acoustic velocity with an IML electronic hammer [31]. Stiffness, MFA and other wood property traits were further assessed in wood cores of the sampled trees using Sil-viScan 2 technology [32]. This high resolution technology can measure the properties of wood produced in different seasons and different years, and allowed us to compare gene expression and wood properties (such as MFA and stiffness) in a single season. The aim of this study was to identify putative candidate genes involved in the regulation of microfibril orientation and cell wall mechanics in wood cells.

Selection of sampled trees with contrasting wood stiffness
Measurements of acoustic velocity of standing trees in the two progeny trials showed that variation in juvenile wood stiffness, or the longitudinal modulus of elasticity (MOE) [33], had roughly normal distributions at the family level (Figure 1a, b). Continuous variation of wood stiffness in the two trials is a typical feature of quantitative traits. Overall average MOE in the Flynn and Kromelite trials were 4.66 and 3.72 GPa with standard deviations (SD) of 0.42 and 0.22, respectively. The lower MOE in the Kromelite trial may be due to the higher tree growth rate in that trial, as there is a negative correlation between MOE and tree growth in radiata pine [34,35].
In each trial five families with the highest MOE and five families with the lowest MOE were selected for further study (Figure 1c, d), thus between 4.0 and 9.1% of the total families were selected from each trial. The two groups of selected families had large differences in wood stiffness (Figure 1c, d). The families MPB215 and MPB003 were among the five low stiffness (LS) families selected in both trials, which is likely to be due to the moderate to high heritability of juvenile wood stiffness [35][36][37][38][39]. Two individuals with higher (or lower) MOE and above average growth rates were sampled from each selected family for further study (Figure 1e, f). Developing xylem tissues of these individuals were collected at three different seasons (spring, summer and autumn), and total RNAs were extracted for microarray experiments.

SilviScan profiling of microfibril angle and other wood traits
The acoustic velocity-based method using an IML hammer provides an overall measurement of wood stiffness (or MOE) of a standing tree. Wood cores of the 20 high stiffness (HS) and low stiffness (LS) trees selected in each trial were further analysed using SilviScan 2. This high resolution technology measured MFA, MOE and wood density along wood cores in earlywood and latewood of each ring formed in different seasons of each year. The three wood traits were compared in each ring between the HS and LS trees ( Figure 2). Annual growth rates as seen in the width of each ring were similar between the two groups (P-values ≤ 0.05), suggesting growth rate of the sampled trees should have a limited impact on comparisons of wood stiffness. In the first 3-4 rings (from the pith) the differences in MFA, MOE and wood density between the two groups were not statistically significant; while from ring 4 or 5 onwards MFA in each ring of HS trees was significantly lower and MOE significantly higher than that in LS trees (P-values ≤ 0.05) ( Figure 2). In contrast, wood density in these rings was not significantly different between the two groups. Thus, the 20 standing trees with contrasting stiffness based on acoustics were further confirmed by the SilviScan measurements. The SilviScan results also showed that the two groups of trees had contrasting stiffness and MFA in individual rings produced after 3 or 4 years. Thus, when developing xylem tissues were sampled for microarray experiments from the two plantations at year 6 and 7, respectively, the single ring produced in the two groups of sampled trees had contrasting stiffness and MFA. These provided a basis for the relatively precise analysis of wood traits (wood stiffness and MFA) and gene expression in the same ring. Thus, differential gene expression between the two groups of sampled trees is expected to be largely related to their contrasting stiffness and MFA rather than wood density. For comparisons between gene expression and wood formed in a particular growing season, SilviScan data was analysed in different seasons of each year. Three wood properties (MFA, MOE and density) were compared between HS and LS trees in both trials in earlywood (EW) and latewood (LW) (Figure 3 and Additional file 1). MFA showed a marked decline while MOE increased with tree age (up to 6 yrs) in the EW and LW of both HS and LS trees. Compared to the 10 trees with LS (IML-based) in both trials, the 10 HS trees (IML-based) had a significantly higher MOE (SilviScan) and lower MFA in the EW and LW of each ring (Pvalues ≤ 0.05). In the two groups of sampled trees, LW stiffness of each ring was generally greater than that of EW stiffness; while MFA of LW in each ring was consistently lower than that of EW in both trials ( Figure 3 and Additional file 1). Unlike MOE and MFA, wood density was not significantly different between the two groups No. of ring (from the pith) Density (kg/m3) Figure 2 Variation in mechanical and other wood properties between the two groups of sampled trees: individual rings. Ring width, modulus of elasticity (MOE), microfibril angle (MFA) and wood density were measured at a high-resolution in wood cores from the 20 sampled trees using SilviScan 2. Wood traits in each ring were compared between wood with contrasting stiffness and MFA, including ring width, MOE, MFA and density in Flynn (a) and Kromelite (b). Error bars represent the standard deviation of the mean value of each trait. Variation in both MOE and MFA between the two groups of sampled trees is statistically significant while no significant variation appears in ring width and density (P-values ≤ 0.05).
of trees in both EW and LW, except for a few rings produced in the earlier years in the Kromelite trial (Additional file 1). In summary, the two groups of sampled trees produced annual rings and both EW and LW with contrasting MOE and MFA, but similar wood density (particularly in rings produced in the later years). These results further suggested that differences in xylem gene expression in different seasons between the two groups of sampled trees would most likely be related to wood stiffness and MFA rather than density.
In order to obtain more reliable patterns of wood property variation between HS and LS trees, an additional 20 and 40 trees from the 10 HS and LS families in the Flynn and Kromelite trials, respectively, were measured by SilviScan and the data combined with SilviScan data from the 20 trees used for the microarray experiments. The magnitude of MFA variation between the HS and LS trees in both trials was similar in EW and LW of each ring, particularly in the outer rings closest to the sampling time ( Figure 4). In contrast, in both trials variation in MOE between these trees was typically higher in EW than in LW in each ring, suggesting that EW tissues formed in spring could be a major contributor to overall variation in wood stiffness. However, wood density had an opposite pattern compared to wood stiffness. In the Kromelite trial, the density difference between the HS and LS trees tended to be higher in LW than in EW of each ring; while in the Flynn trial greater density differences were also observed in LW in the two rings closest to the bark. No. of ring (from the pith)

Figure 3
Variation in mechanical and other wood properties between the two groups of sampled trees in the Flynn trial: earlywood and latewood in each ring. Microfibril angle (MFA), modulus of Elasticity (MOE) and wood density in wood cores from the 20 sampled trees were measured by SilviScan 2. Wood variation between the two groups of sampled trees was compared in earlywood (EW) (a) and latewood (LW) (b) of each ring for MOE, MFA and density. Error bars represent the standard deviation of the mean value of each trait. Variation in MFA and MOE is statistically significant (P-values ≤ 0.05) but it is not statistically significant for wood density variation in all rings.

Differential gene transcription in wood with distinct stiffness and MFA
Transcripts differentially accumulated in radiata pine developing xylem with contrasting wood stiffness and MFA were identified in the Flynn trial in spring (Flynn-EW) and autumn (Flynn-LW) and the Kromelite trial in summer (Kromelite-LW). Of 3, 320 xylem unigenes present on the cDNA microarrays, 3.4% were differentially regulated in Flynn-EW, including 46 and 66 unigenes preferentially transcribed in HS (low MFA) and LS (high MFA), respectively ( Figure 5a). In contrast, more unigenes (8.9%, 2.6 times) were differentially transcribed in the same trees in autumn (LW), including 151 unigenes preferentially transcribed in HS and 144 in LS trees ( Figure 5a). However, the largest proportion of differentially accumulated transcripts was found in trees (14.5%) sampled in summer (LW) in the Kromelite trial. Thus, LW tissues had 3-4 times more differentially accumulated transcripts than EW tissues, suggesting that expression of genes regulating wood stiffness and microfibril orientation is impacted by the season. While the number of differentially transcribed unigenes was 3-4 times higher in LW than in EW (Figure 5a), the greatest differences in wood stiffness (MOE) in rings formed in the sampling years were observed in EW laid down in spring ( Figure 4). This suggests that the number of differentially transcribed unigenes alone does not explain the seasonal pattern of MOE variation. Comparisons of differentially accumulated unigenes in the Flynn trial showed that only 10 (or one) are common to both EW (spring) and LW (autumn) in juvenile wood with HS and low MFA (or LS and high MFA) (Figure 5b, c). Similarly, only eight (or 13) unigenes were common in both EW (spring) in Flynn and LW (summer) in Kromelite. In contrast, considerably more differentially regulated unigenes were conserved in LW tissues sampled in autumn (Flynn) and summer (Kromelite), including 43 unigenes for HS (low MFA) and 45 for LS (high MFA) wood ( Figure  5b, c). Taken together, these results suggested that distinct sets of genes control microfibril orientation and cell wall mechanics in different seasons.
Microarray expression of 10 differentially transcribed unigenes (eight from HS and low MFA, and two from LS and high MFA) selected in the Flynn-LW experiment were validated by RT-MLPA using the same RNA samples as for the microarray experiments. Transcript accumulation measured by the RT-MLPA method was relatively consistent with the microarray results for all 10 validated genes ( Figure 6), particularly in the ranking of expression magnitude, suggesting the microarray experiments in this study were sufficiently reliable for the identification of candidate genes that may influence juvenile wood stiffness and MFA.

Functional analysis of differentially accumulated transcripts
A total of 307 and 446 xylem unigenes (after redundant unigenes were removed) preferentially transcribed in juvenile wood with HS (low MFA) and LS (high MFA), respectively, were identified in the three microarray experiments. More than 64% of the unigenes had homologs in the UniProt known protein [40] and RefSeq [41] databases (tblastx, E-value < 1e-5) ( Table 1). Functional analysis showed that 52-59% of the unigenes had been assigned a gene ontology (GO) term, while the remaining 41-48% were functional unknowns (Table 1). Among the differentially transcribed unigenes with assigned GO terms, the majority were classified as molecular functions (86-93%) or biological processes (81-82%), with fewer as cellular components (61-64%) ( Table 1). Slightly more GO terms were assigned to unigenes preferentially transcribed in wood with HS and low MFA (59.3%) compared to LS trees with high MFA (52%), and more genes with molecular functions and/or involved in cellular components are preferentially transcribed in LS wood with high MFA (Table 1).
Further comparisons of lower level GO terms assigned to the differentially transcribed unigenes showed distinct differences between juvenile wood with contrasting stiffness and MFA in 17 sub-categories (Table 2). In terms of cellular components, genes involved in vacuolar membrane and proton-transporting two sector ATPase complexes had higher representation in HS wood with low MFA; whereas genes with functions in external encapsulating structure were preferentially transcribed in LS wood with high MFA. Genes more abundantly transcribed in trees with HS and low MFA included more genes with different molecular functions, such as binding activities (actin, copper ion, heme, tetrapyrrole and transition metal ion), translation elongation factors and transmembrane transporters. In contrast, only two sub-categories of molecular functions (coenzyme binding and carbon-oxygen lyase activity) had higher representation in juvenile wood with LS and high MFA. In the HS wood (low MFA), genes involved in translational elongation, transport (ion, electron and cation), stress responses and ribonucleotide biosynthetic process were more highly transcribed; whereas LS wood with high MFA had increased transcription of genes involved in different metabolic and developmental processes and RNA processing.
Putative candidate genes for microfibril orientation and wood stiffness To identify putative candidate genes involved in the regulation of wood stiffness and microfibril orientation, differentially accumulated transcripts identified in a single  Juvenile wood included two types: wood with high stiffness (HS) and low microfibril angle (LM), and wood with low stiffness (LS) and high microfibril angle (HM). b Refer to the matches in the UniProt known protein and RefSeq databases (tblastx, E-value < 1e-5). c CC-cellular component; MF-molecular function; BPbiological process. Percentage (%) of CC, MF and BP were calculated using the total unigenes assigned with GO terms.
microarray experiment were examined in the other two microarray experiments. In the two experiments with LW tissues, 46.5% of transcripts identified in Flynn-LW (autumn) were confirmed in Kromelite-LW (summer), and 30.7% transcripts identified in Kromelite-LW (summer) were detected in Flynn-LW (autumn). In contrast, fewer genes identified in Flynn-EW (spring) were confirmed in Flynn-LW (19.4%) and Kromelite-LW (28.7%). This is likely to be caused by differential gene expression in wood tissues formed in different seasons. A total of 51 and 57 candidate genes were identified with differential transcription in juvenile wood with HS (low MFA) and LS (high MFA), respectively, in at least two of the three microarray experiments (Additional file 2). Among these candidate genes 43 were functional unknowns based on the UniProt known protein and RefSeq databases. The remaining candidate genes for HS (and low MFA) with relatively clear functions are listed in Table 3 (29 genes for HS -low MFA) and Table 4 (34 for LS -high MFA). Of the 29 candidate genes for HS (low MFA) eight genes are involved in secondary cell wall formation, including four lignin genes (C4H, C3H, CAD and PCBER) and four secondary wall cellulose synthases (PrCesA1, 3, 7, 11) [42]. Other genes involved in cell wall synthesis included sucrose synthase (SuSy) and genes involved in actin filament (actin depolymerising factor and actin) and microtubule development (tubulin beta-3). In contrast, only three of the 34 candidate genes preferentially transcribed in LS (high MFA) wood had a clear function in cell wall development (cellulase, pectate lyase, peroxidase) and none were specifically involved in secondary cell wall synthesis.

Discussion
The mechanical properties of wood are complex traits Tensile stiffness of xylem cells is largely determined by the angle of cellulose microfibrils in the middle layer (S2) of the secondary cell wall [14,15,43]. Wood stiffness has been shown to be under moderate to strong genetic control [35][36][37][38][39] and a number of QTLs for MOE or its major component traits (wood density and MFA) have been identified in a range of tree species [18][19][20][21][22]24]. The present study showed that wood stiffness had continuous variation with a nearly normal distribution at the family level, thus providing additional evidence that wood stiffness is a typical quantitative trait which is likely to be controlled by variation in numerous genes of small individual effect. Over 100 candidate genes were differentially transcribed in juvenile wood with contrasting stiffness and MFA, providing further molecular evidence that wood stiffness and MFA are complex traits controlled by numerous genes. The 3, 320 xylem unigenes on the microarrays were derived from relatively abundant xylem transcripts in radiata pine [42] and they may only represent 30-40% of genes moderately to highly transcribed during wood formation [44]. Therefore, the total number of genes influencing wood stiffness and MFA in radiata pine may be up to 2-3 times more than the number identified in this study. The complexity of wood stiffness and MFA was also demonstrated by the relatively low magnitude of changes in transcript accumulation observed between juvenile wood with contrasting stiffness and MFA. In the three microarray experiments, the unigenes showing differential accumulation had average expression ratio  values of 1.50 (Flynn-EW, spring), 1.46 (Flynn-LW, autumn) and 1.32 (Kromelite-LW, summer), respectively. These values are much lower than those observed in comparisons between EW and LW [45], and between juvenile and mature wood [46] using the same cDNA microarrays and experimental protocols. The smaller changes in transcript accumulation suggested that genes involved in the regulation of wood stiffness and microfibril orientation may individually have small effects.

Cytoskeletal genes are implicated in stiffness and MFA
There is growing evidence that cortical microtubules, the dynamic arrays of alpha-and beta-tubulins in the cytoskeleton, play a key role during the crystallization of cellulose microfibrils [43] by directing their orientation during deposition in the wall [47]. In this study a betatubulin gene was preferentially transcribed in HS (low MFA) wood of radiata pine, providing molecular evidence that tubulins have functional roles in regulating MFA in the secondary xylem of gymnosperms. Several alpha-and beta-tubulins have previously been observed to be highly expressed in poplar xylem and tension wood with significantly reduced MFA [48]. Allelic variation in an alpha-tubulin was significantly associated with MFA in loblolly pine [27]. In Eucalyptus nitens a beta-tubulin gene was preferentially expressed in xylem from the upper sides of branches which had low MFA [49] and its over-expression in transgenic eucalypt xylem tissue directly influences MFA [25]. However, different beta-tubulin proteins have been observed to be expressed in both compression and opposite wood of maritime pine [50]. This study showed genes encoding actin and actin depolymerising factor (ADF) that were more strongly transcribed in juvenile wood with HS (low MFA) in radiata pine. ADF plays an important role in regulating the optimum balance between unpolymerised actin molecules and assembled actin filaments [51]. Actin filaments are one of the three major polymers in the cytoskeleton, but they are much less rigid compared to microtubules [10]. Actin microfilaments have been observed to align with cortical microtubules which in turn aligned with the orientation of cellulose microfibrils in cultured cotton fiber cells [52]. The actin and ADF genes identified in this study may affect wood stiffness and MFA by influencing interactions between actin filaments and microtubules.

Secondary cell wall genes are implicated in stiffness and MFA
Candidate genes involved in secondary and primary cell wall formation may contribute to variation in MOE and MFA observed between trees. Many candidate genes preferentially transcribed in HS (low MFA) wood had clear roles in secondary cell wall formation. This included four secondary cellulose synthases (PrCesA1, 3, 7, 11) [42] and several genes involved in the biosynthesis of lignin (C4H, C3H, CAD, PCBER). In addition, transcripts (contigs) preferentially transcribed in HS (low MFA) juvenile wood in a single microarray experiment included many other secondary cell wall genes, such as 4CL, CCoAOMT, chitinase-like, SAMS, PPBG, AGP4, GRP2, PRP, etc (data not show). In contrast, none of the candidate genes transcribed preferentially in LS wood (high MFA) have clear roles in secondary wall development. In fact, cell wall genes identified in juvenile wood with LS (higher MFA) were involved in cell division (cyclin-like F-box) and primary wall formation (pectate lyase, peroxidase, cellulase and ovule/fiber cell elongation protein). Taken together, increased secondary cell wall synthesis appears to be a fundamental process in the formation of HS (low MFA) wood.
Crystalline cellulose microfibrils are synthesised by cellulose synthase (CesA) complexes in the plasma membrane then extruded into the external matrix [53]. This study showed four secondary wall PrCesA genes and a SuSy gene that were transcribed more highly in HS wood with reduced MFA. Three other secondary cellulose synthases were present on the microarrays but were not differentially regulated (PrCesA5, 6,8), suggesting that certain members of the CesA family may have specific roles associated with modification of cell wall stiffness and microfibril orientation. Previous studies showed that an over-expressed SuSy gene in poplar increased cellulose content, secondary wall thickness and wood density, but did not affect microfibril orientation [54]. However, a recent study identified an allelic variant in a radiata pine SuSy gene that was significantly associated with MFA [28].
Lignification is a distinct feature of secondary xylem development. Many lignin pathway or related genes were preferentially transcribed in HS (low MFA) wood, including C4H, C3H, CAD and PCBER, as well as several other genes that were identified in individual microarray experiments (4CL, CCoAOMT, and SAMS). More transcription of these genes in HS wood (with lower MFA) could increase lignification of secondary walls, resulting in higher density as seen in HS wood ( Figure  2, 3 and Additional file 1). The lignified matrix in secondary walls can also generate compression resistance which additionally contributes to tensile stiffness of cells [43]. Antisense 4CL transgenic poplar produced low stiffness wood with reduced lignin content [55]. Associations of lignin genes with MFA have been observed in Eucalyptus (CCR) [30] and radiata pine (PAL, COMT and PCBER) [28], suggesting lignin genes may also be involved in the regulation of microfibril orientation. Previously, several genes involved in lignin biosynthesis have been found to be differentially regulated in the formation of tension wood [49,56] and compression wood [57] in which stiffness and MFA have been drastically altered. Taken together, candidate genes involved in lignin biosynthesis may influence both MFA and wood density through their key roles in secondary wall deposition.
Hormone signalling genes associated with wood stiffness and MFA It is well known that wood formation is largely regulated by different hormones. In the poplar xylem transcriptome approximately 2% of ESTs are involved in hormone biosynthesis [58]. Higher ethylene levels have been observed in compression and tension wood formation during which MFA is significantly altered [59], and an ACC oxidase protein involved in ethylene biosynthesis was more strongly expressed in compression wood [11]. Increased transcription of a gene encoding an ethylene-forming enzyme in high MFA wood and an ethylene responsive element binding factor (ERF) in low MFA wood was observed in this study. These data suggested that ethylene levels may affect MFA and wood stiffness via the transcription factor, ERF. The present study also identified a cytokinin-binding protein gene (CBP) that was more abundantly transcribed in HS wood (low MFA) in both summer and autumn. CBP has an essential role in cytokinin signal transduction pathways [60] and extra supply of cytokinin increases the formation of secondary xylem with higher lignification and thicker cell walls [61]. Thus, hormone signalling genes could be involved in the regulation of cellulose and lignin biosynthesis, and their differential expression may produce wood with contrasting stiffness and MFA. Identification of ethylene-and cytokinin-responsive genes in this study suggests their possible co-regulatory roles in juvenile wood formation leading to distinct stiffness and MFA. Interestingly, phytochrome, a lightresponsive gene, was more abundantly transcribed in HS (low MFA) wood in radiata pine. Phytochrome regulates a large number of genes involved in hormone signalling or enzymes involved in cell wall modification [62]. However, the role of the phytochrome gene in the regulation of microfibril orientation and deposition remains unclear.

Conclusions
Naturally occurring variation in microfibril orientation and mechanical properties of wood was correlated with distinct xylem transcriptomes, particularly in wood synthesized late in growing seasons (latewood). Genes involved in cytoskeleton development and secondary cell wall formation (cellulose synthesis and lignin pathway) were preferentially transcribed in secondary xylem with reduced microfibril angle and higher stiffness. In contrast, a few genes with a role in cell division and primary wall synthesis were more highly transcribed in wood with low stiffness and higher MFA. Many genes with unclear functions were also differentially regulated in wood with altered microfibril orientation and distinct mechanical properties. The identified candidate genes are a valuable resource for future transgenic studies and association analyses aimed at improving the mechanical properties of wood through manipulating cellulose microfibril orientation.

Field trials and selection of sampled trees
Two radiata pine (Pinus radiata D.Don) field trails planted at Flynn, Victoria (38°14' S, 146°45' E) and Kromelite, South Australia (37°50' S, 140°55' E) were used in this study. The two trials included 250 and 110 half-or full-sib families, respectively, and 16 full-sib families were common to both trials. At the time of sampling for wood property measurement and microarray experiments, the cambial ages (ring numbers at breast height) of the sampled trees in the Flynn and Kromelite trials were 7 and 6 yrs, respectively, both of which were producing juvenile wood. Tree diameters at breast height (DBH) were measured in both trials and overall wood stiffness of standing trees was assessed using an IML hammer (Instrument Mechanic Lab, Inc. Kennesaw, GA, USA) in the methods described previously [36]. The IML hammer provides a measurement of acoustic velocity (V) of an impulse to travel through the trunk of standing trees. Wood stiffness (MOE) was calculated based on the IML measurement: MOE = V 2 × 1000 [31].
Five families each with the highest and lowest MOE were selected at each site, on the basis of the IML-based measurement of standing trees. To maximize the difference between HS and LS trees in the comparisons, two individuals with higher (or lower) MOE were further selected in each of the five families with highest (or lowest) MOE. To account for the influence of growing conditions, all selected individuals had a straight bole with growth rates similar to or higher than the family average. The HS and LS trees being compared were also grown in a similar environment, such as no adjacent dead trees, similar soil and same slope direction. A total of 20 trees (5 families, 2 trees per family) in each trial were selected for microarray analysis.

Measurement of microfibril angle and other wood traits
Wood cores sampled from 40 trees at the Flynn site and 60 trees at the Kromelite site were analysed at a highresolution by SilviScan 2 [32,63]. These trees were selected from the HS and LS families identified above and included the trees used for microarray analysis. A wood core (12 mm in diameter) was drilled at breast height from each sampled tree. Cores were further trimmed to produce a 2 mm thick and 7 mm wide strip. MFA, wood density and MOE were measured across the wood strips. Wood trait variation based on SilviScan was compared between the two groups of sampled trees using average values in EW and LW of each ring, and P-values were calculated to indicate statistical significance.

Collection of developing xylem tissues
Developing xylem tissues of selected trees were collected in the Flynn trial in spring (October) and autumn (April), when earlywood (EW) and latewood (LW) respectively, was being synthesized and from the Kromelite trial in summer (late November) when LW was being formed. The xylem tissues were scraped at breast height with a sharp chisel after removal of the bark. In the Flynn trial EW and LW tissues were collected on opposite sides of the trunk from the same trees. To avoid the presence of compression wood, developing xylem was collected from the tree trunk perpendicular to the prevailing wind direction. Fresh xylem tissues were immediately placed into liquid nitrogen in the field, and then stored at -80°C until RNA extraction.

RNA isolation and microarray experiments
Total RNA was extracted using a modified CTAB method [64]. Radiata pine cDNA microarrays contained 6, 169 xylem ESTs assembled into 3, 320 unigenes (986 contigs and 2, 334 singletons) [42,45]. These unigenes were derived from genes moderately to highly transcribed in radiata pine wood formation [42] and may account for about 30-40% of the genes expressed in wood formation based on white spruce and poplar [44]. Probe synthesis, microarray hybridization and data normalization were performed in methods described previously [45]. Raw expression values from the 12 microarrays were deposited in the NCBI GEO database (accession number GSE23020).
Transcript accumulation was compared in juvenile wood with HS (low MFA) and LS (high MFA) using three microarray experiments: the Flynn trial in spring (EW) and in autumn (LW), and the Kromelite trial in summer (LW). In each experiment, total RNA samples extracted from the 10 HS wood were pooled into two bulks (five trees each, one tree per family) for comparisons with RNA pooled into two bulks of five LS wood. This pooling approach could partly account for the genetic variation among the different genotypes and provided two biological replicates for each experiment. A dye swap was performed for each biological replicate.

Microarray data analysis and functional annotation
Wood stiffness is a quantitative trait [35][36][37][38][39] and is likely to be influenced by changes in expression of numerous genes. Some of these changes are expected to be small according to the multigenic model, in particular, the expression of transcription factors involved in the regulation of wood stiffness and microfibril orientation. Furthermore, gene expression usually changes in different biological replicates with various genetic backgrounds. Thus, a slightly lower threshold of transcript accumulation (20% change) and a moderately stringent average P-value (≤ 0.05) were used in this study to identify transcripts differentially accumulated in wood with distinct stiffness and MFA. Microarray expression of selected transcripts was validated using RT-MLPA (see below). Additional confidence in the identification of candidate genes was achieved by confirmation of the differentially accumulated transcripts in at least two of the three microarray experiments.
Identified differentially accumulated transcripts were functionally annotated using gene ontology (GO) terms [65]. Comparisons between transcripts differentially accumulated in wood with contrasting MFA and stiffness were performed according to the assigned GO term categories. GO terms differentially represented in HS (lower MFA) and LS (high MFA) wood (P-values ≤ 0.05) were identified using the "library comparison" function of the Bio301 system, an automated EST sequence management and functional annotation system [66].

Validation of differentially accumulated unigenes
Microarray expression of differentially accumulated transcripts identified in the Flynn trial in autumn (Flynn-LW) were validated using the reverse transcriptase-multiplex ligation dependent probe amplification (RT-MLPA) method [67]. Eight genes preferentially transcribed in wood with HS and lower MFA, including cellulose synthase 3 (PrCesA3), chitinase-like, protective protein for beta-galactosidase (PPBG), cinnamate 4hydroxylase (C4H), phytochrome, PrCesA11, CCoAOMT and transmembrane serine protease 9 (TSP9), and two genes preferentially transcribed in wood with LS and higher MFA (FadR and T12C24.11) were included in the validation (Additional file 3). Developing xylem tissues collected for the microarray experiment (Flynn-LW) were also used as the starting materials for the RT-MLPA validation. Eight biological replicates were used with four technical replicates. Differential transcription of each gene was summarized as a mean log-2 ratio from the 32 replicates.
Approximately 400 ng of DNase treated total RNA was reverse transcribed into first strand cDNA using the ImProm-II Reverse Transcription System (Promega, WI). The cDNA was hybridized at 60°C overnight with bulked RPO (right probe oligo) and LPO (left probe oligo) designed for the validated genes (Additional file 3). Ligation and PCR amplification were performed with SALSA D4 primer. Individual gene fragments were separated from the mixed PCR products using a CEQ™ 8000 Genetic Analysis System (Beckman Coulter, CA) and relative gene expression levels were determined using the built-in software.