A Genomic Approach to Suberin Biosynthesis and Cork Differentiation.

Cork (phellem) is a multilayered dead tissue protecting plant mature stems and roots and plant healing tissues from water loss and injuries. Cork cells are made impervious by the deposition of suberin onto cell walls. Although suberin deposition and cork formation are essential for survival of land plants, molecular studies have rarely been conducted on this tissue. Here, we address this question by combining suppression subtractive hybridization (SSH) together with cDNA microarrays, using as a model the external bark of the cork-tree ( Quercus suber ), from which bottle cork is obtained. A SSH library from cork-tree bark was prepared containing 236 independent sequences; 69% showed significant homology to database sequences and they corresponded to 135 unique genes. Out of these genes, 43.5% were classified as the main pathways needed for cork biosynthesis. Furthermore, 19% could be related to regulatory functions. In order to identify genes more specifically required for suberin biosynthesis, cork ESTs were printed on a microarray and subsequently used to compare cork (phellem) to a non suberin producing tissue such as wood (xylem). Based on the results, a list of candidate genes relevant for cork was obtained. This list includes genes for the synthesis, transport and polymerization of suberin monomers such as components of the fatty acid elongase complexes, ABC transporters and acyltransferases, among others. Moreover, a number of regulatory genes induced in cork have been identified, including MYB, NAM and WRKY transcription factors with putative functions in meristem identity and cork differentiation.


INTRODUCTION
Land plants have evolved lipophilic barriers that protect the internal living tissues from dehydration, injuries and pathogens, and have evolved regulatory networks to adjust the barriers to the changing physiological and environmental conditions of the plant. Plant primary organs such as young stems and leaves, are protected by the cuticle, a lipophilic extra cellular polymer membrane composed of cutin and waxes. Secondary (mature) stems and roots, tubers and healing tissues are protected by cork, a tissue with multiple layers of cells that are dead at maturity. Key compounds for cork impermeability are suberin -a complex polymer comprising both aliphatic and aromatic domains-and associated waxes.
Cork is part of the plant constitutive defense system and contains secondary compounds such as triterpenoids and soluble phenylpropanoids that act on herbivores, microbes and fungi.
Cork or phellem, which is the technical term for cork, is formed by the phellogen (cork cambium). Cork formation involves proliferation and commitment of the phellogen derivatives, cell expansion and extensive deposition of suberin and waxes, and an irreversible program of senescence ending in cell death. The two best known and most studied examples of cork are the suberized skin of potato tuber (Sabba and Lulai, 2002) and the bark of the cork-tree (Quercus suber), from which bottle cork is obtained (Silva et al., 2005). In the cork-tree, the phellogen forms a continuous layer of cells that envelops the tree trunk and produces, each year, a 2-3 mm thick layer of almost pure cork that adheres to that of the previous year (Caritat et al., 2000). The chemical composition of the cork tree bark has been widely analyzed by chemical fractionation (Silva et al., 2005 review). Although the amounts of the different components can show significant variations (Pereira, 1988;Lopes et al., 2001), on average it contains 15% extractives (7.5% waxes and 7.5% tannins), 41% aliphatic suberin (referred to as "suberin" in cork-tree literature), 22% aromatic suberin (referred to as "lignin" in cork-tree literature), 20% polysaccharides, and 2% ashes (Pereira, 1988). The monomeric composition of the aliphatic (suberin) fraction has been analyzed by Holloway (1983) and by Bento et al. (1998) and that of the aromatic (lignin) fraction by Marques et al. (1999) and Lopes et al. (2000), among others. Waxes have been analyzed by Castola et al. (2005) and mainly consist of terpenes and sterols. Tannins, mostly ellagitannins, and other soluble polyphenols have been analyzed by Cadahia et al. (1998) and Varea et al. (2001). Many of the enzymatic activities necessary for cork biosynthesis may be inferred from the chemical composition. Cork is a site for four main secondary metabolic pathways: acyl-lipids, phenylpropanoids, isoprenoids and flavonoids. The acyl-lipids pathway is required for the biosynthesis of the linear long-chain compounds forming the aliphatic suberin domain, which share upstream reactions with cutin biosynthesis (Kolattukudy, 1981;Heredia, 2003;Kunst and Samuels, 2003;Nawrath, 2002). The phenylpropanoids pathway is needed for the biosynthesis of the cork aromatic Schreiber, 2005). Recently, a genome-wide study of the shoot epidermis of Arabidopsis (Suh et al., 2005) highlighted a series of new candidate genes relevant in cutin and wax synthesis. In lignin research, molecular genetic approaches have been widely used in the past. Global xylem transcript profiling have been reported for Arabidopsis (Ko et al., 2004;Ehlting et al., 2005) and for several tree-species (Whetten et al., 2001;Kirst et al., 2004;Andersson-Gunneras et al., 2006). For a number of genes, involvement in lignin biosynthesis was demonstrated by forward and reverse genetic approaches (Anterola and Lewis, 2002;Sibout et al., 2005;Abdulrazzak et al., 2006). An excellent review (Carlsbecker and Helariutta, 2005) summarizes the present knowledge of the molecular genetics of regulatory networks in xylem.
In order to reveal the genetic repertoire of cork cells and to identify genes likely to be related to suberin synthesis, we used a two-steps strategy. First, by means of suppression subtractive hybridization (SSH), a library of ESTs preferentially induced in cork was obtained. Then, these ESTs were printed on a microarray and subsequently used for a global comparison between a suberin producing (cork/phellem) and a non suberin producing (wood/xylem) tissue. Isolation of suberin genes in the corktree is particularly attractive because of its exceptional capacity to produce suberin. Peeling of the external bark from the cork-tree trunk allowed the harvesting of differentiating cork layers (Fig. 1A) and provided a highly enriched material for molecular investigations. In the following pages we present an initial analysis of the genomics of cork cells in cork-tree bark; as far as we know, this is the first global approach to cork and suberin molecular biology.

Cork subtractive library
Suberin is a product of the secondary metabolism that is regulated in a tissue-specific manner. It was our intention to find candidate genes for cork and suberin biosynthesis, therefore, we chose as driver tissue for the SSH a fully undifferentiated tissue consisting of the proliferative mass obtained from Q.suber somatic embryo cultures (Fig. 1B). The proliferative mass is a translucent fully undifferentiated non vascularized tissue that develops in the hypocotyls of the recurrent somatic embryos. Cork tree somatic embryogenesis has been carefully characterized at anatomical and ultrastructural levels in previous works (Puigderrajols et al., 1996 and2001). Neither suberin deposition nor multilamellated cell walls could be detected in somatic embryos using optical and fluorescent microscopy techniques or by electron microscopy. The extraction of RNA from cork is difficult due to its high proportion of dead cells and phenols; however, using the protocol of Chang et al. (1993) that prevents oxidation of phenols, good quality RNA was obtained. A further mRNA purification step was added in order to reduce rRNA contamination and to increase SSH efficiency. Finally, the SSH products were cloned into PCR4-TOPO vector resulting in a library of 975 cDNA clones. Subtraction efficiency was checked by co-hybridizing tester (cork) and driver (embryo) RNAs on the cork-tree microarray described in material and methods section. This control experiment confirmed a high yield of subtraction efficiency: 96% of reliable ESTs (CV<0.3) were up-regulated in cork (FC>1.5); 75% of reliable ESTs showed signal intensity not significantly higher than background in embryo (data are shown in Table S2).
Single run sequencing of the library yielded 694 readable sequences longer than 100 bp. Of these, 579 grouped into 121 contiguous sequences (contigs) and 115 were single sequences (singletons).
Thus, in total, 236 independent sequences were obtained. Sequence redundancy [100 x (1-(contigs+singletons / readable sequences))] was of 66%. BLASTX analysis (Altschul et al., 1990) showed that 69 % of the independent sequences (71 singletons and 92 contigs) showed high sequence homology at the amino acid level to database sequences (e-value <e -21 ), and that 31 % (44 singletons and 29 contigs) could not be assigned to any gene ontology and were classified as no hits (e-value >e -10 ).
Sequences with the same GenBank entry were assumed to represent the same gene and, thus, the library was found to contain 135 unique genes (121 with assigned function and 14 with unknown function). The genes with assigned function were manually grouped into functional categories using NCBI, TAIR and published data. Categories were established taking into account the main metabolic and cellular processes leading to cork biosynthesis. Table I shows a selected list of genes with putative functions thought to be important for suberin biosynthesis and cork regulation grouped in functional categories.
The complete list is given in Table S1.
The relative contribution of the genes to the different categories is shown in Fig. 2. Acyl-lipids, Isoprenoids, Phenylpropanoids and Flavonoids, the four categories that represent the major pathways for the synthesis of cork chemical components, amounted to 43.5% of the genes. The Regulatory proteins category, which includes transcriptional regulation, signal transduction and regulated proteolysis related genes, amounted to 19% of the genes. The category Stress, that combines genes related to detoxifying enzymes and cell wall strengthening, amounted to 9.5%; and the category Unknown, which groups those genes with no assigned biological function, amounted to 10%. The genes not fitting to any of the above classes were grouped according to their annotations into two different categories named Miscellaneous and Others. The Miscellaneous category, which includes genes compatible with the main pathways leading to cork biosynthesis but whose substrates have not been characterized, amounted to 9%. The genes in Others are not further discussed in this paper. On the other hand, as can be observed in Tables I   and S1, the number of ESTs (N, redundancy) showed remarkable differences among the genes. In SSH libraries, although SSH should, in principle, decrease the frequency of abundant transcripts while increasing the probability of rare transcripts, genes both differentially and strongly expressed become over-represented (Ranjan et al., 2004). Therefore, with all precautions, we can hypothesize that genes showing a high redundancy are preferentially and strongly expressed in cork. This is the case of genes encoding long-chain acyl CoA synthase (LACS), hydroxycinnamoyl-coenzyme A/benzoyl transferase (HCBT) and cytochromes of the CYP72A subfamily, among others.

Differential screening between phellem and xylem by microarray hybridization
Cork (phellem) and wood (xylem) tissues have common features: both originate from secondary (cambial) meristems and both synthesize aromatic polymers. However, only cork tissue produces suberin and associated waxes. Therefore, in order to identify genes mostly related to suberin synthesis, the cork ESTs from our library were printed on a microarray (see for details the material and methods section) and hybridized to cork and wood tissue. For hybridization, cork and wood RNA was obtained from field grown cork-trees during the vegetative season when both cambial layers are in full activity.
Three independent cork-trees (biological replicates) were sampled and a dye swap for each biological sample (technical replicates) was performed. Microarray data were lowess-normalized to account for intensity-dependent differences between channels. After normalization, dye swap replicates showed no strong deviations from linearity ( Fig. 3A) proving low dye bias. The comparison between the biological and technical replicates showed a high degree of inter-array reproducibility with Pearson's correlation coefficients ranging from 0.95 to 0.98 (Fig. S1). In order to select those genes with good evidence of being differentially expressed we used a Volcano plot (B/M, odds versus ratio) (Fig. 3B) and established a cut-off of B>3. Fold change (FC) values for genes with B>3 are given in Tables I and S1; all data can be found in Table S2.
The great majority of the library genes were up-regulated in cork (B>3, FC>2). Regarding the main cork biosynthetic pathways, genes within the Acyl-lipids and Isoprenoids categories, more relevant for suberin and wax biosynthesis, showed much higher FC values than genes within the Phenylpropanoids and Flavonoids categories, more relevant for aromatic compounds biosynthesis. The fact that most genes in the Phenylpropanoids category were up-regulated in cork could indicate that specific paralogs are induced in this tissue. This hypothesis is supported, for instance, by two 4CL (4-coumarate: CoA ligase 1)-coding genes, with only one paralog differentially expressed in cork. Quite the opposite, the two paralogs coding HCBT (N-hydroxycinnamoyl/benzoyltransferase) are both strongly cork up-regulated (FC = 37 and FC = 32, respectively). Since such a strong cork induction suggests a specific role in cork synthesis, HCBT could be a key enzyme for synthesis of phenylpropane derivatives characteristic of suberin such as feruloyltyramine. Most genes of the Acyl-lipids category were strongly up-regulated in cork. This applies to genes possibly involved in the synthesis of suberin monomers such as the ω hydroxylase CYP86A1 or the β-ketoacyl-CoA synthase (KCS), and enzymes that catalyze ester bonds such as glycerol-3-P acyltransferase (GPAT). The putative lipases/esterases, including GDSL-motif putative lipases, were highly cork up-regulated, and, although the lipase function of these proteins has not been proved in plants, a possible lipase role in cork cannot be discarded. Moreover, interestingly, the highest FC values within the functional categories were shown by genes of the Miscellaneous category. This is a remarkable result taking into account that the Miscellaneous category contains genes encoding enzymes, such as cytochrome P450s, transporters and one putative acyltransferase, which may catalyze reactions important in the biosynthesis of suberin or other cork chemical components.
With regard to the Stress, Regulatory proteins and Others categories, changes of the FC within each category showed diverse behaviors. Two genes involved in regulated proteolysis (ubiquitin/26S proteasome regulatory subunit, FC=91; cysteine proteinase, FC=81), were the two most phellem upregulated genes in the library. It should also be noted that some transcription factors (MYB, FC=60; WRKY, FC=28 and NAM, FC=16) and some signal transduction genes (protein kinase, FC=31; calcium binding annexin, FC=27) exhibited high FC values. On the other hand, all genes of the Unknown category showed strong cork up-regulation, a fact pointing to possible phellem specific functions for these genes.

Validation of differential expression of genes by RT-PCR
We used the RT-PCR with incremental cycle numbers to validate the cork to wood gene expression ratios measured by the microarray. The transcript abundance was analyzed for six relevant candidate genes having moderate to high FC values. The genes selected for validation were: hydroxycinnamoylcoenzyme A/benzoyl transferase (HCBT), ferulate 5-hydroxylase (F5H), long-chain acyl CoA synthase (LACS), palmitoyl-acyl carrier protein thioesterase (PATE/FAT), WRKY transcription factor (WRKY) and phytosulfokine receptor (PSKR). As control, the transcript levels of three constitutive genes (actin, ACT; elongation factor, EF; and Cap Binding Protein, CBP) were measured in order to verify that equal amounts of cDNA were used for both tissues (Fig. 4A). Gene-specific oligonucleotides (Table S3) were used in PCR reactions containing equal amount of both cork and wood cDNAs as templates. Products of incremental cycle numbers were subsequently analyzed. The difference in cycle numbers required for equal amplification of the corresponding PCR product in cork and wood respectively was used to estimate levels of differences in expression within the two tissues. The cDNA was obtained from the three biological replicates used in the microarray hybridization. The possible contamination by genomic DNA was excluded using actin primers specifically designed to differentiate genomic DNA from cDNA.
Results of RT-PCR with incremental cycle numbers (Fig. 4B) confirmed the differential expression of all six selected genes. Amplification products of their transcripts exhibited differences of three to nine cycles, which correspond to the cork-to-wood gene expression ratios measured by the microarray.

DISCUSSION
We report a collection of cork genes potentially important for cork biosynthesis and differentiation based on sequence homology and microarray comparison. This list includes a set of genes possibly involved in the biosynthesis, transport and polymerization of suberin that, in general, agrees with the biosynthetic pathway suggested by Kolattukudy (2001), Bernards (2002) and Franke et al. (2005). The list also contains a number of putative cork regulatory genes which might be of particular interest considering the lack of knowledge in this field. Finally, a number of genes with unknown function strongly induced in cork appeared in this study. Although the present work does not prove the involvement of the candidate genes in cork differentiation, on the basis of this study direct experimental approaches can be designed.
In the two following sections we discuss the putative roles of a set genes potentially relevant for suberin biosynthesis and the regulation of cork differentiation, considering available data on homologous (best BLASTX hit) and related genes.

Candidate genes for suberin
Synthesis of aliphatic monomers. Dihydrolipoamide S-acetyltransferase and biotin carboxyl carrier protein (BCCP) are enzymes involved in de novo fatty acid (FA) biosynthesis, a step necessary for the synthesis of long chain fatty acids (LCFA) and very long chain fatty acids (VLCFA). VLCFAs are precursors of waxes and some suberin monomers. Their up-regulation in cork may be due to a higher demand for acyl chains in this tissue. Palmitoyl-acyl carrier protein thioesterase (FAT/PATE) and longchain acyl-CoA synthetase (LACS) are enzymes involved in the export of FA and LCFA from the chloroplast, a required step for the synthesis of VLCFA. Genes encoding FAT and LACS enzymes are required for normal wax and cuticle development (Bonaventure et al., 2003;Schnurr et al., 2004).
Cytochromes P450 catalyze several key reactions in the synthesis of the aliphatic monomers. CYP86A1 is a hydroxylase that has capacity for ω -hydroxylation of C12-C18 FA (Benveniste et al., 1998).
Members of the CYP86A family, LACERATA and ATT, are involved in cuticle and cutin synthesis (Wellesen et al., 2001 andXiao et al., 2004). Lipoxygenases (LOX) catalyze the addition of molecular oxygen to polyunsaturated FA and have been related to epoxydation of cutin monomers by the lipoxygenase/peroxygenase pathway (Blee and Schuber, 1993;Lequeu et al., 2003). FA elongase complexes (FAE) in the endoplasmic reticulum catalyze chain elongation required for the biosynthesis of VLCFA. Two genes encoding condensing enzymes of these complexes are present in our library, a β -ketoacyl-CoA synthase (KCS) and a β -ketoacyl CoA reductase (KCR). These two genes show, respectively, high similarity to At5g43760 (91% similarity) and At1g67730 (79% similarity), which have been identified as possible candidates for cuticular wax biosynthesis (Costaglioli et al., 2005). Todd et al., 1999;FIDDLEHEAD, Yephremov et al., 1999, Pruitt et al., 2000CUT1, Millar et al., 1999;CER6, Fiebig et al., 2000) and KCR (Dietrich et al., 2005)  Assembly of the aliphatic polyester. The polymerization of the suberin glycerol polyester is poorly understood, but is thought to take place in the apoplast by esterification of the carboxyl groups of α,ωdiacids and ω -hydroxyacids with the alcohol groups of either glycerol or ω -and mid-chain hydroxy fatty acids (Kolattukudy, 2001). Our library contains two putative glycerol-3-P acyltransferase (GPAT) that catalyze the formation of ester bonds. Although the known plant GPATs were previously thought to be involved in membrane and storage lipid synthesis, the analysis of the transcriptome of the stem epidermis suggested a role of these enzymes in the synthesis of aliphatic polyesters (Suh et al. 2005) which has now been confirmed by analysis of the gpat5 mutant (Beisson et al 2007).
Synthesis of aromatic monomers. Genes encoding phenylpropanoid-related enzymes which act upstream in the synthesis of aromatic compounds are represented in the library. Most of them catalyze biosynthetic reactions leading to the phenylpropanoid precursors of the aromatic part of suberin (see Bernards, 2002). Some of these enzymes, such as phenylalanine ammonia-lyase (PAL) and cinnamate 4-hydroxylase (C4H), are key enzymes for the regulation of the phenylpropanoids pathway. Other enzymes, such as F5H (ferulate 5 hydroxylase), are thought to be important for the synthesis of the hydroxycinnamic acids typical of suberin (Bernards and Lewis, 1998). Two genes, both strongly upregulated in phellem, code for hydroxycinnamoyl/benzoyltransferases (HCBTs), a family of enzymes capable of catalyzing the synthesis of N-hydroxycinnamoyl amides such as feruloyltyramine (Yang et al., 1997). Feruloyltyramine is found in potato wound suberin (Negrel et al. 1996) and some evidence supports the hypothesis that it is also a component of the aromatic suberin in Q. suber (Marques et al., 1999). HCBTs belong to the BAHD acyltransferase superfamily (D'Auria, 2006), discussed below. On the other hand, the biosynthesis of hydroxycinnamic acid amides and their subsequent polymerization in the plant cell wall is generally accepted as a plant defense response (Facchini et al, 2002).
Assembly of the aromatic polymer. Despite the fact that peroxidase-mediated coupling was proposed for the assembly of the aromatic monomers in suberizing cells and that several suberin-associated class III peroxidases have been described (Kolattukudy, 2001;Bernards, 2002), no class III peroxidase was found in our cork library. However, cell wall peroxidase activity could be developed by other proteins present in the library like the apoplastic annexin (classed in the Regulatory proteins category), which has peroxidase activity (Gorecka et al., 2005). The extracellular H 2 O 2 needed for the peroxidase catalytic activity could be provided by the copper-containing amine oxidase (Rea et al., 2004).
Conversely, the present study underscores three genes coding for laccases, which are extracellular oxidases with capacity for coupling phenylpropanoids. A role for laccases in lignin synthesis has been proposed (Kiefer-Meyer et al., 1996;LaFayette et al., 1999;Ranocha et al., 2002;Ehlting et al., 2005) and a similar function in the synthesis of the aromatic part of suberin can be hypothesized (Liang et al., 2006).
Linkages between aliphatic and aromatic units. Esterified ferulates are very important suberin monomers (Adamovics et al., 1977;Bernards and Lewis, 1992;Lofty et al., 1994). Some members of the BAHD acyltransferase superfamily catalyze the esterification of hydroxycinnamates with fatty acids (Lotfy et al., 1995). BAHD acyltransferases is a large superfamily of enzymes showing in vitro high catalytic versatility and wide substrate specificity (D'Auria, 2006 classed in Miscellaneous shows high homology at the amino acid level with several members of the BAHD superfamily (e-value <e -30 ). Interestingly, CER2/Glossy2, a gene whose knockout leads to a wax mutant, encodes a BAHD member (Kunst and Samuels, 2003;D'Auria, 2006). Conversely, as previously discussed, the HCBT genes that also belong to this superfamily are more probably related to synthesis of aromatic amides.
In summary, our cork library contains a set of structural enzymes that are probably good candidates for the synthesis of the aliphatic and aromatic monomers of the suberin, and also provides putative candidates for the assembly of the polymer. One interesting observation is the relative importance of the cytochrome P450 superfamily in the cork library. The abundance of P450s monooxygenases and of oxidases reflects the high complexity of synthesizing the cork polymeric matrix and indicates that cork cell metabolism must generate reactive oxygen species (ROS) in high amounts.
This corresponds to previous observations showing that cork cells suffer from fairly high oxidative stress (Pla et al., 1998;Pla et al., 2000). Here, we found that ascorbate peroxidase (APX), a crucial enzyme for H 2 O 2 detoxifying in plant cells (Davletova et al., 2005), was strongly up-regulated in cork.

Candidate genes for regulation of cork formation
Only very limited knowledge is available about the hormonal control of cork formation. The ethylene forming enzyme aminocyclopropane-carboxylate (ACC) oxidase was highly expressed in cork and wood, without showing significant differences between both tissues. Although the possible role of ethylene in these tissues is unclear (Andersson-Gunneras et al., 2003;Lulai and Suttle, 2004), it could be a common regulator in cork and wood. Phytosulfokines (PSKs) are peptide hormones which induce cell dedifferentiation and re-entrance into the cell cycle (Matsubayashi et al., 2002). The presence of a cork up-regulated PSK receptor-kinase suggests a possible role of PSKs in phellem regulation. On the other hand, enzymes acting on lipid catabolism such as lipoxygenases (LOX1) and putative lipases (GDSLmotif proteins) could be involved in the synthesis of wound hormones such as jasmonic acid (Wasternack et al., 2006). Interestingly, the library contains a highly phellem up-regulated annexin. This gene shows 89% similarity to At1g35720, an Arabidopsis annexin that senses the Ca 2+ signal elicited by ABA and transmits it downstream in the signaling pathway (Lee et al., 2004).
Phellogen derivatives undergo very rapidly phases of cell division, cell expansion, bulk suberin deposition and cell death marked by the complete autolysis of the cells. Regulated proteolysis is required during programmed cell death (PCD) and for the switch from one developmental phase to another, a process that requires removing pre-existing regulatory networks (Sullivan et al., 2003;Vierstra, 2003). Genes encoding regulated proteolysis, such as a cysteine protease and a Ub/26S proteasome system, were highly induced in cork. Cysteine proteases have been involved in programmed cell death (Rojo et al., 2004). Moreover, cysteine proteases and Ub/26S proteasome genes are among the most expressed genes during fiber cell death (Moreau et al., 2005).
We have found five transcription factors related to meristem identity that could play a key role in the maintenance of the phellogenic identity of cells or promote their differentiation into phellem cells.
One of them is a R2R3 MYB transcription factor involved in axillary meristem identity in Arabidopsis (Muller et al., 2006). Another one is an AS1 (Asymmetric Leaves 1)-interacting protein, which could be important for establishing cell fate in leaf development (Phelps-Durr et al., 2005). Two ESTs encode proteins of the NAM (No-Apical-Meristem) family (Souer et al., 1996). Finally, a SQUAMOSA/APETALA1 transcription factor, which is a floral meristem identity gene (Mandel and Yanofsky, 1995) was also highly up-regulated. On the other hand, the library also contains a WRKY transcription factor that is highly induced in cork. WRKY transcription factors are a large family of plant-specific regulators that mainly control senescence, stress and defense responses (Eulgem et al., 2000). WRKY factors modulate gene expression by binding to W-boxes of some stress induced genes, including P450s (Mahalingam et al., 2003;Narusaka et al., 2004).
In conclusion, a number of interesting regulatory candidate genes for cork regulation has been identified, although much more work is needed to elucidate their function.

Plant material and tissue harvesting
Cork (phellem) and wood (xylem) tissues were harvested from 15-20-year-old field-grown cork trees (Quercus suber L) at Peratallada (Girona, Spain) during the growing season. External bark (cork bark) was removed and, using sterile scalpels, the exposed phellem tissue could be harvested. Thus, fractions rich in differentiating phellem were obtained (Fig.1A). Wood was obtained after removing the internal bark (secondary phloem) and fractions enriched in differentiating xylem were harvested as described above. Harvested samples of cork and wood were immediately frozen in liquid nitrogen and stored at -80 ºC. In order to prevent genetic and environmental variability, both cork and wood samples cohybridized in the array were obtained from the same tree specimens.
As a source of somatic embryos, a cork tree recurrent embryogenic line maintained in a medium free of plant-growth regulators was used (Puigderrajols et al., 1996). Macronutrients were those from SH medium (Schenk and Hildebrant, 1972), and micronutrients and vitamins were from MS medium (Murashige and Skoog, 1962), including 3% (w/v) sucrose. The culture was solidified with 0.6% (w/v) agar and the pH was adjusted to 5.7. The cultures were incubated in a growth chamber at 25ºC and a 16/8 hr photoperiod at 50 µmol m -2 s -1 was provided by cool-white plus Grolux fluorescent lamps. Subcultured cotyledonary stage embryos showing signs of secondary embryogenesis were picked, their hypocotyledonary region dissected, immediately frozen in liquid nitrogen and stored at -80ºC.

RNA extraction and ds cDNA synthesis
Total RNA was extracted from cork and wood tissue as described by Chang et al. (1993) and from somatic embryo as described by Martin et al. (1993). Remaining traces of genomic DNA were removed by TURBO DNAse treatment (Ambion) and a further purification step was performed with the CleanUp protocol of RNeasy Plant Mini Kit (Qiagen). RNA quality and purity were checked by formamideformaldehyde denaturing agarose gel electrophoresis, spectrophotometry using Nanodrop TM and Bioanalyzer 2100 (Agilent). In order to obtain polyA + RNA, 150 µg of total RNA were purified using Oligotex mRNA kit (Qiagen). The quality and yield of the mRNA was checked by denaturing agarose gel electrophoresis and RiboGreen ® RNA quantitation Reagent (Molecular Probes). Double-stranded cDNA was synthetized from 300 ng of mRNA using primers from the PCR-Select cDNA Subtraction kit (Clontech, USA), Superscript II, E.coli DNA ligase, E.coli DNA polymerase, RNase H and T4 DNA polymerase from Invitrogen, all according to manufacturer's recommendations.

cDNA subtractive library construction (SSH)
The cork tree phellem subtractive library was made applying the suppression subtractive hybridization (SSH) technique (Diatchenko et al., 1996;Diatchenko et al., 1999) and using the PCR-Select cDNA Subtraction kit (Clontech, USA). To isolate ESTs preferentially induced in phellem, cDNA from cork tissue as tester and cDNA from somatic embryo as driver were used. Two rounds of subtractive hybridization and PCR amplification according to Clontech instructions were performed. In order to improve the yield of long ESTs, products from five secondary PCRs were pooled and size-selected on an agarose gel electrophoresis in three fractions: 150 to 800 bp, 800 to 1200 bp and 1200 to 2000 bp.
Each cDNA fraction was concentrated and purified using Minelute PCR Purification Kit (Qiagen), cloned into the pCR ® 4-TOPO ® T/A cloning vector and transformed in TOP10 competent cells (Invitrogen). Kanamycin resistant colonies were picked and grown overnight in liquid LB medium containing kanamycin at 37ºC and 320 rpm. Finally, glycerol was added to grown cultures to 15% final concentration and clones were stored at -80ºC in 96-deep-well plates.

Sequencing and EST analysis
Plasmids were isolated from overnight-grown bacterial cultures using a standard alkaline lysis protocol with SDS in 96-well format. Inserted fragments were amplified by PCR using M13 oligonucleotides. (CAGGAAACAGCTATGAC), 0.02 U µ l -1 DyNAzyme (Finnzymes), 2.5mM MgCl 2 and 50 to 100 µ g plasmid template and using PTC220 Multicycler (Dyad). PCR was done for 1 cycle at 95ºC for 2 min, 35 cycles at 95ºC for 45s, 50ºC for 2 min, 72ºC for 1 min and an additional cycle at 72ºC for 6 min.
Liquid handling steps in plasmid preparation as well as for PCR set up were carried out using a liquid handler robot RSP 200 (Tecan). All PCR products were separated on agarose gels. Gel images were analyzed in order to verify the amplicon length and quality using GelMaster software (Bajla et al., 2005). Single-run sequencing was carried out by MWG-Biotech AG using ABI capillary sequencers.
Raw sequences and base confidence scores were obtained from chromatogram files using the basecalling program Phred (Ewing and Green, 1998;Ewing et al., 1998). Poly-A-tails were removed with Trimest, a tool from the EMBOSS package. Vector masking and trimming as well as adaptor removal were performed with cross match (http://www.phrap.org). Sequences with less than 100 bases after trimming were discarded. Remaining sequences were assembled using StackPack (SANBI, http://www.sanbi.ac.za/Dbases.html) with the standard configuration (for initial clustering, a similarity cut-off value of 96 matching bases over a window size of 100 base pairs compared in words of length 6). Assembled sequences were annotated by using the Stackpack Annotation Module SAM (Universite Bordeaux 2, http://cbi.labri.fr). With SAM, the consensus sequences and singletons were blasted locally against the following databases: SWISSPROT (ftp://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/swissprot.gz), NR (ftp://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/nr.gz) and NT (ftp://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/nt.gz). EST sequences with an e-value >10 -10 were considered as no hit and excluded from further analysis. Functional assignment was conducted manually by analyzing the BlastX-and BlastN-annotated ESTs, where SWISSPROT hits were favoured over NR hits and NR hits were favoured over NT hits.

Construction of cDNA microarray
The cDNA inserts of 975 clones from the cork oak phellem SSH library genes were amplified as described in sequencing and EST analysis, and printed on a microarray. Along with cork ESTs, 2298 poplar ESTs (Dejardin et al., 2004) and 2096 oak ESTs and 1167 potato ESTs from the PICME repository (http://www.picme.at) were added and used to normalize intensity data (Table S2) was set to 45%. After printing, slides were scanned, cross-linked at 65 mJ and stored at room temperature in the dark.

Differential screening by microarray hybridization
For hybridization, 20 µ g of total RNA of each tissue were reverse transcribed using Superscript II (Invitrogen) and dNTP mix containing aminoallyl-dUTP (Clontech). Purification of cDNA, coupling of fluorescent dye and probes purification were performed using Atlas TM Glass Fluorescent Labelling Kit (Clontech) following manufacturer's recommendations. We used Cy3 and Cy5 Mono-Reactive Dye packs (Amersham). After labelling, the probes were quantified using a Nanodrop TM spectrophotometer.
Microarray hybridization was performed using a humidity hybridization chamber (In-Slide-Out, Boeckel Scientific). Briefly, slides were prehybridized with 5x SSC, 0.1% (w/v) SDS, 0.1 mg/ml BSA for 1h at 42ºC; followed by a wash in 0.1 x SSC for 5 min at 25ºC, two washes in water for 30 s, and dried by centrifugation at 1600 rpm for 2 min. Then, slides were hybridized overnight (18 h) at 37ºC in DIG Easy Hyb (Roche) containing labelled probes, yeast tRNA (0.1 mg/ml) and salmon sperm DNA (0.1 mg/ml). After hybridization, slides were washed 3 times in 0.1% (w/v) SDS, 1x SSC for 10 min at 50ºC, and washed 4 to 6 times in 1x SSC. Then, slides were dried by centrifuging at 500 rpm for 5 min and immediately scanned using the DNA Microarray Scanner (Agilent) at 10µm resolution and 100% laser intensity and PMT settings.
Microarray images were quantified using GenePix 6.0 (Axon) software. Only spots with signal intensities twice above the local background, not saturated and not flagged as absent by GenePix, were considered reliable and used for subsequent analysis. Extracted intensities were subtracted from the local background and the log 2 -ratios were normalized in an intensity dependent fashion by print-tip lowess. For phellem to xylem hybridization, normalized log 2 -ratios were scaled between arrays to make all data comparable. Statistically significant differences in gene expression were determined by computing a Bayesian statistic using all log 2 -ratios from replicate hybridizations. ESTs were considered as differentially expressed when their Bayesian statistic B was higher than 3. All quantitative and statistical analyses were performed using MMARGE tool, a web implementation of the Limma package in the R environment. Relative fold change data were recalculated using xylem as the reference tissue.
For phellem to embryo hybridization, ESTs showing log 2 -ratio coefficient of variation (CV) lower than 0.3 (calculated among duplicated spots) were considered as reliable.

RT-PCR with incremental cycle numbers
To verify the difference of expression levels between cork and wood, equal amounts of cDNA were used for PCR with gene specific oligonucleotides in 100 µl reactions. For each tissue sample, single stranded cDNA was synthesized from 1 µg total RNA using the Superscript II (Invitrogen) in a 20 µ l www.plantphysiol.org on August 28, 2017 -Published by Downloaded from Copyright © 2007 American Society of Plant Biologists. All rights reserved. reaction. Then, cDNA was 2.5 fold diluted and 1 µ l of diluted cDNA was used as template. Primers were designed with Primer3 software (Rozen and Skaletsky, 2000), sequences are shown on Table S3. PCR conditions were 2.5 mM MgCl 2 , 0.2 mM each dNTP, 0.5 µ M each forward and reverse oligonucleotide and 0.05 U µ l -1 Eurotaq polymerase (Euroclone). Thermal cycling parameters were as follows: 1 cycle at 95ºC for 3 min and 36 cycles at 95ºC for 30 s, 57ºC for 30 s and 72ºC for 1 min, using a T-Gradient termocycler (Biometra). Aliquots of 10 µ l were taken every three cycles from cycle 12 to 36 and analyzed by agarose gel electrophoresis stained with ethidium bromide as described (Lê et al., 2005). Control amplifications were performed with specific primers to actin (ACT), elongation factor (EF) and Cap Binding Protein (CBP). Target amplifications were performed with specific primers to hydroxycinnamoyl-coenzyme A/benzoyl transferase (HCBT), WRKY transcription factor (WRKY), long-chain acyl CoA synthase (LACS), ferulate 5-hydroxylase (F5H), phytosulfokine receptor (PSKR) and palmitoyl-acyl carrier protein thioesterase (FAT). In order to discard possible genomic DNA contaminations we designed the actin primers complementary to two exons separated by an intron of 94 bp using Populus genomic DNA sequence.

Figure S1: Comparison between biological and technical replicates.
Pearson's correlation coefficients are given. For each gene, EST GenBank accession numbers, best BlastX in GenBank and in TAIR and gene description are given. The putative functions were assigned based on the highest BlastX score match (evalue <10 -20 ). The cork/wood expression ratio is given as FC (Fold Change) for genes with B>3. N is the number of ESTs present in the library. Each EST was spotted twice on the microarray. Annotation and GenBank data concerning to Quercus suber SSH library are given. Spots are considered absent if signal intensity <150; marginal if 150<signal intensity<250; present if 250<signal intensity<65000; and saturated if signal intensity>65000. The cork to wood hybridization was performed using three biological replicates with dye swap. Statistical B, log 2 -ratio (wood as a reference), cork and wood signal intensities and median of intensities (A) are given. The cork to embryo hybridization was performed using same RNA samples     Table I. Selected list of the more relevant candidate genes for suberin biosynthesis and cork regulation grouped in functional categories. The complete list of genes is reported on Table S1. For each gene the EST GenBank accession numbers and the putative molecular function are given. The putative functions were assigned based on the highest BlastX score match (e-value <10 -20 ) and not always are supported by biochemical data in plants. The cork/wood expression ratio is given as FC for genes with good evidence of being differential expressed (B>3). N is the number of ESTs in the library.    R  O  S  s  c  a  v  e  n  g  i  n  g  E  E  7  4  3  6  5  9  A  s  c  o  r  b  a  t  e  p  e  r  o  x  i  d  a  s  e  (  A  P  )  ,  p  u  t  a  t  i  v  e  A  t  3  g  0  9  6  4  0  6  7  1   R  e  g  u  l  a  t  o  r  y  p  r  o  t  e  i  n  s   T  r  a  n  s  c  r  i  p  t  i  o  n  f  a  c  t  o  r  E  E  7  4  3  8  7  0  M  Y  B  -r  e  l  a  t  e  d  t  r  a  n  s  c  r  i  p  t  i  o  n  f  a  c  t  o  r  A  t  5  g  5  2  6  6  0  6  0  1  T  r  a  n  s  c  r  i  p  t  i  o  n  f  a  c  t  o  r  E  E  7  4  3  6  8  0  R  2  R  3  -M  Y  B  t  r  a  n  s  c  r  i  p  t  i  o  n  f  a  c  t  o  r  A  t  3  g  4  9  6  9  0  3  1  T  r  a  n  s  c  r  i  p  t  i  o  n  f  a  c  t  o  r  E  E  7  4  3  8  0  9  W  R  K  Y  t  r  a  n  s  c  r  i  p  t  i  o  n  f  a  c  t  o  r  A  t  2  g  4  6  1  3  0  2  8