Gene family structure, expression and functional analysis of HD-Zip III genes in angiosperm and gymnosperm forest trees

Sequence analysis of HD-Zip III genes from conifer and hardwood trees

Four putative full-length HD-Zip III coding sequences were isolated from P. glauca by EST data mining, RT-PCR, and RACE cloning (Rapid Amplification cDNA End) with degenerate primers. Two class-IV sequences from P. abies have been previously reported and were denoted PaHB1 and PaHB2 [24]. Therefore, we designated the sequences that we isolated as PgHB3 [25] to PgHB6 (Additional file 1 Figure 1 HQ391914 to HQ391917). Predicted amino-acid sequences display the structural features of HD-Zip III, except that PgHB6 has a partially degenerated leucine zipper motif.

The Populus trichocarpa genome sequence [26] was reported to contain eight different HD-Zip III sequences, which are designated HB1 to HB8 [6]. HD-Zip III genes are distributed on seven of the nineteen poplar chromosomes (Additional file 1). We isolated full-length coding cDNA sequences for eight on the putative poplar HD-Zip III genes by RT-PCR, amplification, starting from the P. trichocarpa (Torr. & Gray) × P. deltoides (W. Bartram) hybrid clone H11-11 and from the P. tremula Minch × P. alba L. clone 717-1B4. For each of the eight cDNA clones, nearly perfect sequence identities were used to match the cDNA sequences with previously identified ESTs and genes predicted from the poplar genome [6], thus providing evidence that all of the predicted genes are expressed in Populus spp.

There are five HD-Zip III genes in the Arabidopsis genome belonging to the two major phylogenetic clades RVB and C8, each of which is divided into two subclades [27]. Floyd et al. (2006) [21] and Prigge and Clark (2006) [22] conducted phylogenetic investigations that included HD-zip III sequences from diverse plants, along with full-length and partial Pinus taeda L. cDNA sequences. They concluded that conifer HD-Zip III genes could be assigned to the two major angiosperm clades of C8 and RVB, but two of the conifer sequences were likely part of gymnosperm-specific clades. In this report, a neighbour-joining (NJ) tree [28] was constructed with complete amino acid sequences from several seed plants, including gymnosperms such as P. glauca and P. taeda, and angiosperms such as A. thaliana and P. trichocarpa, as well as lower plants such as the moss Physcomitrella patens (Hedw.) Bruch & Schimp. The resulting tree topology was consistent with previous reports; however, our data suggest that conifer sequences may in fact be uniquely represented in the C8 clade and absent in the RVB clade (Figure 1). The conifers that we analysed may thus have three C8 members, including sequences previously assigned to the RVB clade. The full-length P. glauca PgHB6 and the partial P. taeda PtaHD-34 and PtaHD-35 fell outside angiosperm clades and formed a monophyletic group, consistent with previous reports [21, 22]. Sequence similarity and tree topology clearly grouped the Populus sequences as four pairs of closely related paralogues, which is consistent with the ancestral salicoid genome-wide duplication and reorganisation described in modern Salicaceae [29].

HD-zip III transcripts accumulate during secondary vascular growth in Picea and Populus

Transcript accumulation was profiled in young P. glauca and P. trichocarpa × deltoïdes trees (refered as PtdHB) grown under controlled conditions by using RT-qPCR to compare steady mRNA levels in several organs and tissues (Figure 2). Transcripts of the four spruce sequences accumulated preferentially in the differentiating secondary xylem of stems (2X) and roots (R2X) and gave similar profiles overall (Figure 2A). PgHB3, PgHB4 and PgHB5 RNAs were also abundant in the differentiating secondary phloem (2P), and PgHB5 had the highest relative abundance in the young foliage (YL) (Figure 2A). The data suggested that the different transcripts differ substantially in abundance since the normalised number of RNA molecules varied by two orders of magnitude between the highest and lowest RNAs, i.e., PgHB3 and PgHB6, respectively. The aforementioned data are consistent with putative roles in vascular differentiation, with little indication of diversification between the gene sequences.

Compared with spruce, poplar HD-Zip III genes gave more diversified transcript accumulation profiles across the panel of organs and tissues, even within pairs of closely related paralogues (Figure 2B). The pair PtdHB1 and PtdHB2, which are close homologues of REVOLUTA, gave relatively similar profiles across the panel, except that PtdHB1 was less abundant in mature and old leaves than in developing tissues. Furthermore, PtdHB1 transcript abundance was two orders of magnitude higher than PtdHB2. The pair PtdHB5 and PtdHB6, closest homologues of Corona/AtHB-15, shared similar transcript profiles which varied strongly between the organs surveyed. Both were clearly most abundant in the developing secondary xylem (2X), but also accumulated in the apex and primary stem. On average, PtdHB5 was five to ten times more abundant than PtdHB6. The pair PtdHB7 and PtdHB8, which are the closest homologues of AtHB-8, gave dissimilar and even opposite transcript profiles. PtdHB7 transcripts were abundant in nearly all organs and lowest in the apex (A) and developing secondary xylem (2X), whereas PtdHB8 transcripts were most abundant in these same tissues (A, 2X). Transcripts of PtdHB3, which was a close homologue of PHV and PHB, largely accumulated in the apex and to a much lower degree than in other parts of the trees, especially the roots. Data are not reported for PtdHB4 because its amplification by RT-qPCR was not strong enough for reliable determinations.

Over-expression of wild-type PtaHB1 and PtaHB7 genes in transgenic poplars

Transgenic poplar trees that overexpressed the complete coding sequence of PtaHB1 and PtaHB7 were obtained to investigate the potential roles of these HD-Zip III genes in tree development. The hybrid poplar clone INRA-717-1B4 (P. tremula × P. alba) was transformed using Agrobacterium with either one of the PtaHB constructs or an empty vector control (WT). Several hygromycin-resistant and GUS-positive lines were recovered and used to produce viable plants grown to an average height of 1.20 m in the greenhouse. All of the lines had transgene transcript accumulation levels which were significantly above levels detected for the INRA-717 endogene (Table 1). Interestingly, all of the lines overexpressing PtaHB1 (UBI::PtaHB1) had a visible external phenotype that was not seen in the controls (Figure 3), but no phenotype was observed upon over-expression of PtaHB7 (data not show).

Transgene construct	Gene	Mean log2 ratio	SD	p-value
UBI:PtaHB-1	PtaHB-1***	2.7320	0.4340	<0.001
	PtaHB-2	-0.7020	0.4580	0.0540
	PtaHB-3	-0.5330	0.4730	0.1500
	PtaHB-5*	-0.7920	0.4960	0.0430
	PtaHB-6	-0.5060	0.4840	0.1770
	PtaHB-7	-0.4920	0.4490	0.1470
	PtaHB-8	-0.5100	0.4600	0.1280
UBI:PtaHB-7	PtaHB-1*	0.5521	0.1869	0.0318
	PtaHB-2	0.4272	0.3613	0.2761
	PtaHB-3*	0.7675	0.1940	0.0148
	PtaHB-5*	0.7263	0.1578	0.0046
	PtaHB-6	0.6951	0.2727	0.0573
	PtaHB-7***	2.8634	0.2689	0.0008
	PtaHB-8	0.3194	0.1490	0.0838

Total RNA from the same samples as those used for the microarray profiling: portion of the main undergoing primary growth (IT) for UBI::PtaHB1; scrapped secondary xylem (S2X) for UBI::PtaHB7 transgenic poplar. *p < 0.05, ***p < 0.001(Student's t test, N = 6); SD, one standard deviation.

Further characterisation of the PtaHB1 transformed trees showed that PtaHB1 transgene transcripts were five to eight times more abundant than the PtaHB1 endogene in the controls. The most obvious phenotype in these trees was their drooping leaves. The trees appeared to have a water-stress phenotype (Figure 3A) which was clearly not the case given that they were grown alongside perfectly healthy control trees. Upon closer inspection, it was evident that PtaHB1 over-expression resulted in altered petiole development, causing the leaves to hang downward. Other than the petiole, the leaves seemed to develop normally and to be perfectly healthy, with no indications of altered water relations. On average, the transgenic poplars had petioles that were 15% shorter, and the angle between the adaxial side of the leaf and the stem was 30% wider than those of control trees (Figure 3B). The increased angle and decreased length were statistically significant starting at the 10^th internode from the apex (where the first internode is the first leaf longer than 1 cm) (p < 0.05) (Figure 3C, D).

The vascular organisation of petioles from mature leaves was examined to further investigate the altered development. Cell wall autofluorescence associated with lignin accumulation was observed in transverse sections under UV-illumination, and clearly indicated that the distribution of fibres and vessels was altered in the transgenic trees (Figure 4A-B). The ratio of fibres (small lignified cells; Figure 4) to vascular elements (large lignified cells) was 0.80 in the transgenic trees, compared to 1.55 in the controls (p < 0.001).

Furthermore, quantitative determinations of fibre lengths in stems and petioles through FQA (Fiber Quantitative Analyser, Methods) showed that petiole fibres from immature (LPI 6) and mature leaves (LPI 21 and LPI 43) were slightly, but significantly shorter in the transgenic trees (p < 0.05, Figure 4D). Shorter fibre classes < 0.5 mm were over-represented in the transgenic trees, whereas longer fibre classes (from 0.75 to 1, from 2.0 mm and up) were significantly under-represented compared to the control trees (p < 0.05, Figure 4C). The fibre length classes from primary stems (internode between LPI 5 and LPI 6) followed a similar distribution pattern, but the mean fibre length was not significantly different between controls and transgenics. Secondary xylem fibres from the main stem, which were sampled from the same internodes as the petioles (LPI 6, LPI 21, LPI 43), did not differ between the transgenic and wild-type trees.

Effect of PtaHB1 over-expression on the transcriptome

Primary stem tissues from two control lines (two trees per line, n = 4 individual samples) and two PtaHB1 transgenic lines (two trees per line, n = 4 individual samples) were compared using a 3.4 K low redundancy cDNA microarray (GSE24703 for raw data on GEO database). A total of 48 transcripts that accumulated differentially were expressed with a false discovery rate (FDR, [30]) threshold set to 5% (q < 5.00) (SAM package release 3.0; [31]). Out of the 48 significantly misregulated genes, 8 transcripts were up-regulated and 40 were down-regulated in the transgenic trees (Table 1). The portion of the stem that we targeted in this analysis is also the part of the tree where petioles are actively developing and growing. Approximately one-third of the misregulated genes (14 out 48) had strong statistical support (q < 0.001). However, the fold-change of all genes identified was less than two (-1 < M < 1).

The predicted functions of the misregulated transcripts in the PtaHB1 transgenics were examined and separated into four categories: growth factor-related, cell wall-related, membrane trafficking, and general functions. The growth factor group included sequences related to brassinosteroid action, which are putative leucine-rich BAK1-like proteins (CN520805, CN519565). These genes were down-regulated and suspected to be involved in steroid signal transduction [34]. Genes for ethylene perception and response were also down-regulated (HO702822, CN522424, and HO702885). Cell wall-related sequences were an abundant category of down-regulated transcripts. Other sequences related to cell expansion and cell proliferation were down-regulated, including sequences encoding two fasciclin-like proteins (CN518490) [35] two glyoxalases (CN519263, CN521180) [36], a farnesylated protein (HO702768) [37] and an elongation factor 2 (CN524724). The down-regulated sequences also included a 4CL gene (CN522696) that is involved in the synthesis of G-lignin precursors, and which is consistent with the decrease in auto-fluorescent fibres [38]. The up-regulated sequences encoded transketolase-like proteins putatively involved in isoprenoid biosynthesis (CN523609) [39] and in decreasing cell proliferation in preparation to dormancy [40].

The impact of PtaHB1 and PtaHB7 transgene expression on the other HD-Zip III gene transcripts was investigated in the transgenic poplars (Table 1). In general, the UBI::PtaHB1 constructs led to decreased transcript accumulation of all other HD-Zip III genes. However, the number of RNA molecules was quite variable and the effect was significant only for PtaHB5 (Student's t test, mean log₂ ratio -0.7920, p = 0.0430). In UBI::PtaHB7 transgenic trees, the HD-Zip III genes had slightly increased transcripts levels, but only PtaHB1, PtaHB3 and PtaHB5 were significantly upregulated (mean log₂ ratio 0.5521, p = 0.0318; 0.7675, p = 0.0148 and 0.7263, p = 0.0046).

Accumulation of some, but not all HD-Zip III transcripts is linked to auxin in Poplar

Distinct HD-Zip III gene family evolution in gymnosperm and angiosperm trees

Gene sequences isolated from the moss P. patens with features typical of HD-Zip genes of class I, II, and III clearly indicate that they were acquired early in plant evolution [41]. The sequence analyses presented here (Figure 1) are consistent with the idea that HD-Zips have evolved through gene or genome duplications and potential gene losses [21, 22].

The phylogenetic tree we described (Figure 1) included four full length HD-Zip III cDNA sequences of P. glauca and was similar but not entirely congruent with the tree topology previously predicted with a Bayesian procedure that used full length and partial cDNA sequences from P. taeda [21, 22]. On one hand, previous authors have reported that gymnosperm HD-Zip III sequences could be assigned to both the C8 and RVB clades defined in angiosperm plants. On the other hand, they showed that gymnosperms also formed two independent clades not represented in angiosperms. Our results are consistent with the existence of gymnosperm clades with representatives from Pinus and Picea. These findings support the hypothesis that modern HD-Zip III family structure derives from four ancestral sequences, and that two of the ancestral sequences have been lost in angiosperms leading to clades C8 and RVB, whereas all four clades have potentially been retained in gymnosperms. However, our finding that the RVB clade lacked conifer sequences and the lack of a reference gymnosperm genome sequence led us to conclude that further analyses are needed to confirm whether or not gymnosperms are in fact represented in the RVB subclade.

Poplars have three more HD-Zip III sequences than Arabidopsis, which is consistent with the inferred genome evolution of the former [26]. The poplar sequences clearly formed four pairs of closely related paralogues. The salicoid plant lineage that gave rise to the family Salicaceae (including Populus spp.) appears to have undergone a relatively recent genome-wide duplication and reorganisation [26], whereas Arabidopsis is thought to have undergone genome size reduction [42]. These different evolutionary paths could have led to the loss of certain functions as well as neofunctionalisation or subfunctionalisation within the angiosperms.

Transcription profiles identify HD-Zip III putatively involved in vascular development

Delineating the potential role of HD-Zip III genes in regard with vascular development is aided by comparing RNA transcript accumulation in different organs, tissues and cell types, despite the overlapping profiles that may be observed within the family. Members of the C8 clade have been most strongly linked to vascular development and have not been implicated in leaf formation as such. In Arabidopsis, AtHB-15/CNA is expressed in procambial cells where it is involved in early initiation of vascular cells, and has been implicated in embryo polarity. The AtHB-8 gene product has been shown to promote the proliferation and differentiation of xylem cells. Its expression also localizes to pro-cambium cells, in addition to being modulated by auxin [16]. The transcripts corresponding to the three P. glauca sequences we assigned to the C8 clade were detected in all tissues but preferentially in differentiating secondary vascular tissues both in the stem and in the roots. This observation may represent evidence in support of the phylogenetic position of Picea sequences PgHB-3 to PgHB-5, along with several other gymnosperm sequences, in clade C8 rather than RVB. In Populus, there are four C8 sequences PtdHB5 to PtdHB8 with varied transcript accumulation profiles in vascular tissues. The accumulation of PtdHB5 and PtdHB-6 transcripts were also clearly preferential to secondary xylem tissues. In contrast, the paralogous sequences PtdHB7 and PtdHB8 have very dissimilar profiles and were distinctly not preferential to secondary xylem. These transcript accumulation profiles of PtdHB7 and PtdHB8 indicated that Populus C8 sequences may have undergone relatively recent neofunctionnalisation or subfonctionnalisation, compared to the pair of PtdHB5 and PtdHB6 which share the most similar expression patterns. Overall, it appears that gene duplications found in gymnosperm C8 clade, and even the more ancient duplications at the family level (PgHB6), have not led to strong diversification of expression profiles compared to that observed in angiosperms.

Lateral organ formation has been assigned to RVB clade that includes REV, AtHB-9 (PHB) and AtHB-14 (PHV). The closely related genes PHB and PHV are involved in leaf polarity, while REV has been implicated in several developmental processes, including vascular cambium identity and activity, as well as fibre differentiation. Two putative homologues of Arabidopsis REV genes have been detected in the genomes of Populus, Z. elegans, O. sativa and Z. mays L. [43, 44]. The functions of the Zinnia REV homologues appear to have diverged, with one being implicated in vascular development and the other in lateral organ formation [9]. In contrast, the Populus HB1 and HB2 have similar transcript patterns, except that HB2 transcripts accumulate more strongly in maturing leaves (Figure 2B). Arabidopsis may represent a unique case with a REV paralog potentially having been lost during ancestral genomic rearrangements [42], and resulting in a gain of function for the remaining REV sequence in developing xylem and leaves. Ko et al. (2006) [6] found that PtdHB1 was associated with secondary growth in poplar stems and hypothesised that HD-Zip III genes played a role in secondary xylem differentiation in trees. Our expression survey indicated that PtdHB1 transcripts are present at a similar level in the apex, primary stems, secondary xylem, and young roots.

Poplar HD-zip III genes play a role in fibre development

Constitutive over-expression of the poplar PtaHB1 gene in poplar led to greater transcript abundance corresponding to this gene, and resulted in shorter petioles and a wider angle between the stem and adaxial side of the petiole. The fibres with reduced lignification and shorter length suggested that development of primary xylem fibres was either impaired or delayed in the transgenic trees. Our hypothesis is that the increased angle between the petioles and the stem is caused by a delayed or incomplete fibre development relative to leaf expansion. Asynchronous development may cause the petioles to lack the necessary strength to support a fully expanded leaf. This phenotype bears a resemblance to that of the ifl-1 mutant (REV gene) in which interfascicular fibre development is impaired and inflorescence stems lack sturdiness [43]. Moreover, it has been observed that REV gain-of-function promotes xylem differentiation and accumulation, leading to dysfunctional vascular patterning [13]. It thus appears that the PtaHB1 over-expression phenotype is closer to the ifl-1 mutant phenotype than the REV gain-of-function phenotype. This result contrasts with our initial hypothesis that PtaHB1 over-expression might promote fibre differentiation.

Furthermore, the vascular phenotype observed in PtaHB1 transgenic poplars was not consistent with previous observations that over-expression of non-mutated HD-Zip III genes had no effect on plant development because of gene silencing by microRNAs. In Arabidopsis, all of the HD-Zip III family transcripts are targeted and negatively regulated by microRNAs (MiR165/166) found in multiple copies in the genome [6]. Moreover, the over-expression of the HD-Zip III sequences was shown to trigger the production of these microRNAs [45]. An hypothesis to explain our observations may be that HD-Zip III PtaHB1 over-expression triggered the accumulation of microRNAs that down-regulated the other members of the HD-Zip III family. Transcript accumulation of the other poplar HD-Zip III sequences provided evidence in support of this hypothesis. Interestingly, PtaHB7 over-expression did not appear to have a similar effect on other family members.

It has long been known that fibre differentiation or stem elongation is controlled by growth factors such as auxin, GA and brassinosteroids, and that overproduction of auxin decreases primary stem elongation in Arabidopsis [46]. The REV/IFL-1 gene of Arabidopsis has been implicated in polar transport of auxin [43]. Considering the known relationship between plant growth regulators and HD-Zip III genes, the delay in petiole elongation in PtaHB1 transgenics may be linked to a perturbation of growth regulation activity. For example, auxin accumulation may be shifted due to altered polar transport. Our observations in young wild-type poplar trees suggested a putative link between auxin and the expression of three different HD-Zip III genes (PtdHB1, PtdHB5 and PtdHB8), the transcript levels of which were affected by the application of an auxin transport inhibitor (Table 3). This observation is consistent with the fact that those gene transcripts are well represented in secondary xylem (Figure 2B). Plant decapitation (from the apex to LPI 3 inclusively) significantly down-regulated transcripts of PtdHB5 and had a smaller effect on PtdHB7 transcripts. These results are partially consistent with previous reports from Arabidopsis, where AtHB-8 was clearly modulated by auxin. Gene regulation may have evolved differently in poplar compared to Arabidopsis, reflecting the developmental and physiological differences between a woody perennial plant and an herbaceous annual.

The differentially expressed sequences identified in the PtaHB1 over-expressing poplars were classified into a few broad categories based upon putative functional assignments. Most of the transcripts were down-regulated and included sequences related to development pathways: brassinosteroid and ethylene growth regulators, ethylene perception and response, and putative steroid signal transduction proteins, in addition to cell wall-related and cell expansion or cell proliferation proteins. Many of the sequences have putative functions that appear consistent with developmental activities of IFL1/REV in Arabidopsis. They also appear to be consistent with the functions or pathways related to the phenotype observed in transgenic poplars.

Plant material and growth conditions

Picea glauca (Moench) Voss (white spruce) tissues were obtained from two sources. The gene isolation work used plantlets of clone Pg-653, which are produced by somatic propagation by the Canadian Forestry Service (Klimaszewska et al., 2001). For the transcript-level survey, samples were obtained from two wild field grown 33-year-old trees in a progeny trial that had been established near Québec City. Two hybrid poplar clones were used: Populus trichocarpa (Torr. & Gray) × P. deltoides (clone H11-11, for gene expression experiments) and Populus tremula L. × P. alba L. clone INRA clone 717 1-B4 (referred to as clone 717 for transgenic analysis). Rooted softwood poplar cuttings were produced in 25 cm³ pots protected in clear plastic bags, transferred to 3 L pots after five weeks, and maintained in a greenhouse. Both the spruce and poplar plants were grown in a greenhouse with 16 hours light per day, with a temperature regime of 22/17°C (day/night), and relative humidity of at least 70%. Natural daylight was supplemented with light from HQI-TS 400W/DH metal halogen lamps (Osram, Munich, Germany). Plants were fertilised weekly with 1 g L^-1 20/20/20 (N-P-K) and supplemented with calcium every two weeks.

Tissues for transcript profile survey

All tissues were frozen in liquid N₂ immediately upon removal from the tree and stored at -80°C until further use. Several organs/tissues were isolated from two 33-year-old field-grown white spruce P. glauca for transcript accumulation profiling. These included the terminal leader (1T), young needles from the upper crown (YN), differentiated secondary phloem (2P) and xylem (2X), as well as bark (B) tissues collected from three 30-40 cm bolts taken from the lower third of the main stem. The secondary phloem and xylem issues were scraped with scalpel immediately after removing the bark: 2P was scrapped from the exposed inner side of the bark and the 2X from surface of the exposed wood. This method relies on cleavage of tissues at the cambial zone and yields relatively pure tissue types, although the purity of 2X and 2P samples were not verified as such. Similarly, tissues from large roots, including differentiating xylem (R2X) and phloem (RPP; with phelloderm) were collected taken in a one-meter radius from the base of the stem. Each sample was kept separate for total RNA extraction. Tissues were also collected from two trees of clone H11-11 selected to be similar in size, i.e., average of 80 cm tall and possessing at least 25 leaves with a plastochron index (LPI) greater than zero [47]. LPI 0 denotes a leaf blade that is 1 cm long and is undergoing laminar expansion. The tissues consisted of: Apex LPI 0 without leaves (A); shoot tips up to and including LPI 1 (1T); young leaves up to LPI 3 (YF); mature leaves LPI 15-16 (ML); old leaves LPI 30-31 (OL), and differentiating xylem and phloem scrapped with a scalpel from the portion of the main stems exhibiting secondary growth (2X and 2P, respectively) between LPI 15 and LPI 20; and actively growing white coloured roots (R).

Auxin-related treatments in poplar

Two treatments were applied to young poplar trees (P. trichocarpa × P. deltoides H11-11, described above) averaging 80 cm in height in order to assess the potential effect of auxin on HD-Zip III transcript accumulation in poplar. Six randomly selected trees which possessed at least 25 leaves with a plastochron index (LPI) greater than zero were assigned to each treatment and to the control untreated group. The first treatment removed the shoot apex, which involved cutting off the top part of plants down to LPI 3 inclusively. The second treatment applied an inhibitor of polar auxin transport, N-(1-naphthyl) phthalamic acid (NPA, 1 mM for treatment and 0 mM for control) mixed with lanolin as a carrier. The mixture of NPA and lanolin was placed on the bark (after removing the layer of cuticle), entirely covering a segment of the stem of 2.5 cm centered at the internode LP15, and covered with paraffin film. Secondary xylem tissues were collected 72 h after the treatments were applied, from internodes LPI14 to 16 and consisted of whole stem segments without the bark. For the NPA treated trees, the sampling consisted of stems segments without bark 3 cm above and 3 cm below the region of NPA application (LPI15 was excluded). The transcript accumulation data (see below) obtained from these tissues were analysed using Student's t test applied to each gene separately comparing the treated and untreated plants. For the NPA treated trees, tissues obtained below and above the region of application were analysed separately and the data were averaged for each tree.

RNA extraction, sequence isolation and phylogenetic analyses

Total RNA extractions were carried out using the CTAB method of Chang et al. (1993) [48] from spruce and poplar tissues, which were stored at -80°C until used. The quantity and integrity of total RNA were evaluated by spectrophotometry (OD 260/280 ratio of 1.8:1 to 2.1:1), and by using a Bioanalyzer 2001 using an RNA 6000 Nano Kit (Agilent Technologies, Palo Alto, CA, USA) to achieve a 28S:18 S ratio of 1.2:1 to 2.2:1 and an RNA integrity number (RIN) above 7.

Partial HD-Zips III gene sequences were identified among EST (Expressed Sequenced Tags) obtained from the conifers P. taeda [49] and P. glauca [50] through BLAST analysis (Basic Local Alignment Search Tool). Completed coding sequences and untranslated 3' and 5' region (UTRs) were isolated by using PCR-cloning with degenerate primers or 3' RACE, 5' RACE (SMART RACE cDNA Amplification, BD Biosciences Clontech, Mountain View, CA, USA) or both cloning methods with mRNAs from needles or xylem. The cDNAs were cloned in pCR2.1 with the TA cloning kit (Invitrogen, Carlsbad, CA, USA), then sequenced with a 16-capillary genetic analyzer (ABI Prism 3130XL and an ABI Prism 3100XL, Applied Biosystems, Foster City, CA, USA). Public databases were searched by use of the BLAST algorithms to identify HD-Zip III sequences in the Poplar genome (JGI, P. trichocarpa v1.0, [26] and updated v2.0 http://www.phytozome.net/poplar), and the non-redundant NCBI database (nr; http://blast.ncbi.nlm.nih.gov/). Corresponding cDNAs were amplified from first-strand cDNA which had been derived from pooled total RNA extracted from developing secondary xylem, tips and young leaves (P. trichocarpa × P. deltoides and P. tremula × P. alba). The identity of each cDNA clone was confirmed by complete and partial sequence comparisons with the poplar genome and previously identified EST sequences [6].

The predicted HD-Zip III amino acid sequences from multiple species were aligned using ClustalW algorithm [51], which is included in the Bioedit software [52], and then refined manually. Phylogenetic analyses used the alignment of complete amino acid sequences via MEGA software [53]. A Neighbour-Joining phylogenetic tree [28] was constructed based on a Poisson correction model and pair-wise deletion algorithm (1,000 bootstrap replicates). The phylogenetic tree was rooted with three members of other HD-Zip class 1 and 4 sequences found in the Genbank database [54].

RT-qPCR procedures

Total RNA (2 μg) was DNase-treated and reverse-transcribed using the SuperScriptII First-Strand Synthesis System for RT-PCR following manufacturer's recommendations (Invitrogen, Carlsbad, CA, USA). The resulting first-strand cDNA was diluted 1:10 in deionised water before quantitative PCR determinations. An aliquot of cDNA equivalent of 20 ng of total RNA was used per 20 μL of PCR reaction. Amplifications were performed in a QuantiTect™PCR SYBR^® Green Kit (Qiagen, Mississauga, ON, Canada) with 0.3 μM of 5' and 3' primers, with a DNA Engine Opticon TM 2 System (MJ Research Inc., Ramsey, MN, USA) following the manufacturer's instructions. Primers were designed with Primer3 software [55], with a melting temperature (Tm) between 55°C and 62°C, and produced amplicons between 100 and 250 bp (Additional file 1). The poplar RT-qPCR primers were tested on the clones H11-11 and INRA-717. Amplification reaction efficiencies were between 90 and 105% for each primer pair. The thermal-cycling parameters were as follows: 95°C for 15 min; 45 cycles of 94°C for 15 sec, 55°C for 1 min and 72°C for 40 sec; followed by a melting curve analysis from 54°C to 95°C with increments of 0.1°C per step to verify the specificity of the amplification and presence of primer dimers. The number of target copies in each sample was determined from Ct values (cycle threshold values) using a linearised plasmid or purified PCR product to produce a standard curve, which was obtained by averaging values from several runs. Two to six independent sample assays were performed and each sample was loaded in duplicate. Results were normalised relative to the absolute RNA used in a single reaction and with the transcript level of a reference gene. The spruce data were normalised to the reference of transcript level gene encoding Elongation factor 1 alpha, EF1a, (AJ132534) [56]. The EF1a gene were strongly expressed and showed low variation between spruce tissues [57]. In poplar, the reference was a cell division cycle 2 gene homologue, CDC2 (AF194820) [58]. CDC2 cycle threshold values (Ct) averaged 18.5 + 0.5 (Standard Deviation), which is the experimental error range of the RT-qPCR cycler device. For statistical analyses of the RT-qPCR data from auxin-related treatments and transgenic experiments (Tables 1, 2, 3) a log₂ transformation was applied to the number of target copies.

Agrobacterium transformation and growth of transgenic poplars

Over-expression constructs were obtained by inserting the complete coding sequences of PtaHB1 and PtaHB7 from P. tremula × P. alba clone INRA-717 between the maize ubiquitin promoter [59] and a 35 S terminator into the pCambia1305.2 vector http://www.cambia.org. The resulting plasmids were then transferred into A. tumefaciens strain C58 pGV2260 [60].

For transformation, in vitro plantlets of the hybrid poplar clone INRA-717 were micropropagated on hormone-free 1/2 MS medium [61] supplemented with vitamin D3. Internodes from in vitro plantlets were co-cultivated with the engineered A. tumefaciens according to the method of Leple et al. (1992) [62], with the following modifications. After co-cultivation, the explants were decontaminated of A. tumefaciens with cefotaxim and transferred onto M2 medium containing cefotaxim alone [62]. Transgenic calluses were selected on M3 medium supplemented with hygromycin (10 mg L^-1) and emerging shoots were transferred for rooting to 1/2 MS medium without hygromycin, screened for positive X-glucuronidase activity [62], and assayed for poplar PtaHB1 and PtaHB7 mRNA accumulation by RT-qPCR. Shoots with elevated PtaHB1 and PtaHB7 mRNA levels relative to WT control were transferred to the greenhouse where they were grown under a photoperiod of 16 h of light, at 24/20°C (day/night). For each construct, four to six trees for each of six lines (independent transformation events), and 10 to 15 WT INRA-717 trees were propagated and carried forward together for phenotypic determinations.

Phenotype characterisations of transgenic poplars

Morphometric measurements were performed on two trees from three independent transformation events per line and six control trees (WT). Internode length, petiole length, and petiole angle were measured from LPI 1 to the bottom of each tree. Vascular tissue organisation was visualised on samples taken from fresh stem and petiole sections from LPI 20. Samples were fixed in 4% paraformaldehyde - 2% glutaraldehyde cacodylate buffer (pH 7.2) for 12 hours under vacuum, and stored at 4°C before embedding in paraffin blocks. Microtome cross-sections (10µm thick) were mounted on glass slides; after the paraffin was removed, the sections were rehydrated for observation under UV light on an Olympus BX51 microscope (Olympus, Montreal, QC, Canada). Fibre and vessel distributions in petioles were visualised by lignin autofluorescence and counted (75 to 150 cells per field imaged). The average fibre to vessel ratio was computed from four transgenic plants and four WT plants. Fibre Quality Analysis (FQA; OpTest Equipment, Hawkesbury, ON, Canada) used freshly debarked internodes, and intact petiole sections from LPI 6, LPI 21 and LPI 43. The samples were macerated in Franklin's solution until they were completely bleached and the fibres could be separated [63]. The FQA weighted, fibre length (lw) data (8000 individual fibres per petiole or internode) were obtained from three trees per transgenic line on six lines and on six control trees. For analyses, the data were placed into bins of 0.10 to 0.25 mm.

Statistical analyses of phenotypic and RT-qPCR data

Statistical treatment of phenotypic data compared transgenic poplars transformed empty vector (control lines) and those transformed with the PtaHB1 gene construct. Except where otherwise noted, the data were not transformed and each tree was considered an individual experiment unit. Normality was confirmed and Student's t tests were used for petiole length and angle data (Figure 3), with each internode tested separately in SYSTAT 13 (Cranes Software International Ltd., Chicago, IL). Student's t tests were also used to compare the proportion fibres to vessels within the petioles, determined by using UV microscopy (Figure 4A and 4B). The fibre length data (FQA) were analysed by comparing the proportion (%) of fibres in each length class (bin) (Figure 4C) with a one-way ANOVA, on each bin independently, and also by comparing the overall fibre length data (Figure 4D). Before analysis the normality of the data was confirmed, and ANOVAs were carried out using procGLM in SAS (Version 9.01, SAS Institute Inc., Cary, NC).

For RT-qPCR data comparing control and PtaHB1 transgenic lines (Table 1, 2) and comparing the impact auxin-related treatments (Table 3) normality and Student's t tests were applied to log₂ transformed numbers of transcript targets in SYSTAT 13 (Cranes Software International Ltd., Chicago, IL).

Microarray RNA profiling

A poplar 3.4 K cDNA microarray was prepared by the ARBOREA project and is described in Pitre et al. (2010) [64]. RNA transcript profiling was carried for each of the transgenic constructs, with four samples from two independent transformation events compared with four control samples from two independent transformation events (empty vector controls). Primary stem tissues (LPI 0 to 5) were used for profiling PtaHB1 over-expressing lines and secondary stem xylem tissues (LPI 21 to 43) were used for the PtaHB7. The microarray hybridization methods were as described in Pitre et al. (2010) [64]. Each microarray hybridisation used 1 μg of total RNA that was amplified using the SuperScript™ Indirect RNA Amplification System (Invitrogen, Carlsbad, CA, USA) and 5 μg of the aRNA were labelled with Alexa Fluor^®555 and 647 dyes (Invitrogen Carlsbad, CA, USA), for use in dye-swap experiments. The poplar 3.4 K microarray slides were pre-hybridised for 2 hours at 42°C in a solution containing 5× SSC, 0.1% SDS, 0.02% BSA (w/v), 0.01% herring sperm DNA (w/v), and 50% formamide. The slides were then washed twice in 0.1× SSC, once in water, rinsed in 2-propanol, and finally dried by centrifugation. The labelled targets (3.5 μL) were mixed with 52.5μL of hybridisation solution containing 5× SSC, 0.1% SDS, 0.01% herring sperm DNA (w/v), and 50% deionised formamide. The mixture was heat-denatured for 4 min at 95°C and cooled for 5 min on ice prior to hybridisation to the microarray. The microarray was then covered with a LifterSlip (Erie Scientific Company, Portsmouth, NH, USA) and placed in a hybridisation chamber II with increased depth (Corning, Lowell, MA, USA) and incubated for 12 h at 45°C in a model 1012 hybridisation oven (Shel Lab, Cornelius, OR, USA). After hybridisation, the slides were iteratively washed for 15 min in 2× SSC + 0.5% SDS, 0.5× SSC + 0.1% SDS and 0.1× SSC solutions at 45°C.

The slides were scanned using a ScanArray™Express scanner (Packard BioScience, Meriden, CT, USA) and the image files were analysed using QuantArray^® software (Packard BioScience, Meriden, CT, USA). Scan intensities were comparable between sets of slides for a given hybridisation. Data analysis was carried out using Bioconductor packages http://www.bioconductor.org distributed in R (R Development Core Team 2008). Median foreground intensity minus median background intensity was the response variable used for the statistical analysis. Data quality was assessed using graphical analysis tools in the marray and olin packages available in Bioconductor, and by assessment of within- and between-slide Pearson product-moment correlation coefficients (r), which were calculated both from the raw intensities and after normalisation. The composite normalisation method [65] was applied by using the two functions maNorm2 D and maNormLoess in the marray package [66]. We identified differentially expressed sequences with the SAM package release 3.0 [31] using a false discovery rate (FDR) [30] threshold set to 5% (q < 5.00). Data reported in the article are log₂-ratios of Alexa Fluor^®555/647, which are denoted M. Raw data are available in the Gene Expression Omnibus (GEO) database (accession number: GSE24703, http://www.ncbi.nlm.nih.gov/geo/).