A cross-species transcriptomics approach to identify genes involved in leaf development

Leaves are of fundamental importance to life on earth, representing the powerhouses of most food chains and are, ultimately, the source of energy that sustains humanity. Although much has been learnt about the biochemical and physiological functioning of leaf-level processes such as photosynthesis [1], still relatively little is known about why leaves are the size and shape they are and how they come to be so. There are genes known to affect meristomatic pattern formation (e.g. AS1 and WUS, KNOX and CLV, see [2] for a review of leaf development) and the rate of leaf primordia initiation [3] as well as genes that contribute to the determination of leaf length (ROT3 [4], LNG [5]) and width (AN, [6]), with less being known about the determination of leaf size (although [7] provide exciting insight). Although the two now classic examples of ROT and AN clearly affect leaf width and length, each will likely explain only a small proportion of the variation that exists for these traits, as numerous examples suggest that leaf size and shape are under complex polygenic control [8–17]. The vast majority of work identifying genes that function in the control of leaf size and shape has, to date, been conducted using the model species Arabidopsis thaliana, with most genes being identified through mutant screens, and, to the best of our knowledge, with no examples of evidence that such genes show variation between ecotypes of Arabidopsis with contrasting leaf shape or size. Beyond this, the work in [18] showed how conclusions drawn about gene function in Arabidopsis, which has simple leaves, may not hold true in species with more complex, dissected leaf forms.

We are interested in identifying genes that are functionally important in determining natural variation in leaf size, shape, and development. Forward genetic screens represent an important approach for identifying genes that have the potential to contribute to natural variation of phenotypic traits. However, such mutant screens are not directed by natural variation in phenotype and there is no reason to expect that genes identified through induced mutation (or other means of disrupting gene function) are those that have been acted upon by natural selection, which could render them interesting but unimportant in an ecological/evolutionary context. The opposite paradigm of reverse genetics offers a way of identifying functionally important loci but, in the case of QTL or association mapping, is time-consuming and often significantly limited by the biology of the system being studied – for example in forest tree species the generation of mapping populations is a very long-term commitment and the subsequent QTL (Quantitative Trait Loci) mapping resolution is of questionable value if the aim is to identify genes, or even causal polymorphisms, underlying a QTL (although association mapping largely overcomes this resolution limitation, for example see [19], but requires greater starting information to enable targeting of relevant loci). A third approach, which represents a refinement to a reverse genetics pipeline, is to identify genes whose expression pattern(s) suggest they might be functioning during a developmental process or response to stimuli/stress (i.e. functional genomics, or perhaps initially 'guilt-by-association') and to combine this information with QTL mapping, and approach that was termed 'genetical genomics' [20]. This proposed view of how to consider gene expression (or similarly peptide/metabolite abundance) in a genetical context is currently receiving considerable attention, extension and refinement of methodology [21–25]. Association mapping and genetical genomics both represent major motivations for the efficient identification of candidate genes, particularly where whole-genome assays cannot be utilised, such as in non-model systems or species such as aspens, where whole-genome SNP assays with adequate genomic resolution for association mapping due to rapid LD decay are still a long way from being available. In such cases well-conceived 'omics' studies may act as a magnet to help find the needles in the haystack that is the genome.

The use of microarrays is now a firmly established method, with transcriptomics being the most mature of the now multitude 'omics' fields. Expression microarrays allow the parallel profiling of most or all genes in a genome, and sensibly designed biological experiments allow the application of microarrays to identify patterns of gene expression associated with a trait of interest. Such traits can be temporal, spatial, or adaptive. As well as conducting individual experiments aimed at answering specific hypotheses, the concentration on MIAME (Minimum Information About a Microarray Experiment, [26]) compliance has enabled the use of meta-analysis approaches that aim to identify patterns across many experiments, or across species. Much effort has been put into the development and subsequent refinement of methods to isolate biologically meaningful information from the background noise [27–30].

An interest in the application of clustering methods to microarray data, whereby groups of genes exhibiting similar expression profiles are identified, was borne from the regulon concept whereby it is expected that a single 'master regulator' (typically a transcription factor) will control the expression of many genes as an efficient way of initiating gene expression 'programs'. For example, leaf development can be considered as a progression through multiple, successive, modules of a developmental program running from primordia initiation to leaf maturity. We have previously shown that some such successive program modules can be identified by profiling gene expression of a field-grown aspen throughout the growing season in multiple years. This identified modules that could broadly be defined as 'cell division', 'cell elongation', and 'differentiation/maturation' [31]. Identifying such clusters (regulons/transcriptional modules, hereafter referred to as transcriptional modules) also allows inference about the function of unknown genes [32]. Clustering was first carried out at the level of an individual experiment using mainly hierarchical clustering [33], k-means clustering [34] and Self Organising Maps (SOMs, [35]) with other methods following. Interest has now extended to identifying transcriptional modules across more extensive collections of expression data [36–41] as well as to the application of network theory and analysis to reconstruct gene regulatory networks [42–47].

We have previously established UPSC-BASE, a database repository with integrated transcriptomics analysis tools [48] for our own Populus cDNA microarray platform, which now contains by-far the largest collection of transcriptomics data available for Populus. Here we show how this data can be mined using a meta-analysis and cross-validation approach (using external datasets and cross-experiment comparisons) to identify candidate genes involved in leaf development. Our aim was to answer the questions a) by comparing the expression profiles of contrasting developing tissue types, can we identify genes that are of particular importance in leaves? b) Using the available microarray data, can we identify transcription factors and 'leaf transcriptional modules'? c) Can we use the results to identify meaningful patterns associated with the process of leaf development? This work represents the first meta-analysis of Populus transcriptomics data and is one of only a small number of cross-experiment transcriptomics studies in plants.

After an early wave of enthusiasm for high-throughput transcriptomics, attention is shifting somewhat towards the next steps up the ladder of the central dogma, proteomics and metabolomics, and to the integration of information across all levels i.e. the systems biology approach. This shift in attention is partly due to the fact that control at the level of the transcriptome tells only part of the whole story in the plot of genotype to phenotype. However, transcriptional data is still important, telling which genes are being transcribed even if they are eliminated as characters from the plot at later stages by factors such as miRNA degradation. The criticism now commonly placed on transcriptomics (i.e. the potential for transcription not to represent protein activity) will, of course, hold equally true when proteomics data is examined in increasing detail; knowing which proteins are present tells you nothing about their half life, rate of replacement or side-branch decoration.

Here we have used an extensive set of transcriptomics data from a diverse range of sources, covering different microarray platforms, species of Populus as well as Arabidopsis to identify genes that are actively transcribed during the process of leaf development. A similar starting point to ours was previously reported for Arabidopsis, where it was found that different tissue types could be differentiated on the basis of gene expression data [56]. However, only four leaf-specific markers were identified. We believe that Populus represents a more suitable model for answering these questions as the degree of tissue differentiation is greater than in Arabidopsis (contrast the trunk of a Populus tree to the stem of Arabidopsis), a fact that should increase the resolution with which tissue-specific expression patterns can be discriminated. Here we were also able to clearly separate tissue types on the basis of gene expression but with greater within tissue-type resolution than was possible using Arabidopsis. Figure 1 shows that leaf samples form a unique group that is distinct from other tissue types. A second group was found for 'wood' samples, including old roots (which were most similar to bark), and a third, less distinct, group containing young roots and flowers. When considering only those gene models predicted to have transcription factor activity, we could see further differentiation within the three main groups (Figure 1B). EST library representation (Figure 2) and functional enrichment analysis (Figure 3) of gene models that were highly representative of leaves in comparison to other tissues identified apical shoot and young leaf EST libraries and GO categories for photosynthesis and chloroplasts (as one would expect) but also a number of categories likely to be involved in leaf development. As functional enrichment was tested for in a list of genes that are more representative of leaves than other tissues, these genes are not simply general markers of high cell division/expansion but are more specifically involved in those processes primarily in leaves. We were interested primarily not in identifying genes that are highly representative of leaves per se but more specifically, genes that are involved in leaf development or that could serve as developmental markers. Our own interest in this approach is in the identification of candidate genes for association mapping in the SwAsp collection of Swedish aspens [62]. We also wanted to develop a set of developmental marker genes as we have previously shown that the magnitude of transcriptional regulation during development is greater than that induced by changes in the weather for a field-grown aspen [31] and that the induced transcriptional response to a treatment can be masked by developmental alterations induced by the treatment [63]. The availability of developmental markers would allow rapid comparisons to be made between leaves thought to be at the same developmental state to confirm whether like is being compared with like. As many more genes than can be sensibly used for downstream applications were identified in the OPLS method, we took the approach of concentrating on transcription factors and of identifying transcriptional modules within our list of candidate genes as a means of selecting a subset of genes that would maximise the information content. We were able to identify transcriptional modules (Figure 4A) that appear to be robust (Figure 4B) and for which broad biological function could be assigned to. For each module, we plotted the expression of the most highly connected gene across a range of microarray data sets in both Populus and Arabidopsis (Figure 5E–F). As can be seen, the green and blue modules appear to be important in early stages of leaf development, the red module in later stages (differentiation rather than cell production), and the turquoise module important in later development and in mature leaf functioning. These patterns of expression make sense when considered alongside the functional classification for each module: The green (protein biosynthesis) and blue (DNA replication and structure) are highly expressed in young leaves that are involved primarily in cell division, the red (organelle/intracellular) module then having a later peak as cells begin to expand and differentiate, with the latest peak in the turquoise (photosynthesis/chloroplast) module, which peaks and remains high as differentiation nears completion and expansion reaches an end, marking the point where leaves switch from sinks to sources. This progression through development is also revealed by the connectivity seen between modules in Figure 4B, with there being a progression in developmental time from the bottom (green and blue) to the top (turquoise) of the diagram. The orange (secondary metabolism) module has a more variable expression profile (data not shown) and is also the least inter-connected of the network modules (Figure 4). Although a regulatory role cannot easily be ascribed to each of these genes (Table 2 – one would expect the most highly connected gene to be the best candidate for being the module regulator gene), for some a regulatory role does seem plausible. For example, the most highly connected gene in the blue module is estExt fgenesh4 pm.C LG V0721, the ortholog of AT1G76540, a cyclin-dependant kinase involved in regulation of the G2/M transition of the mitotic cell cycle [64]. Most of the genes identified in Table 2 do not have reported under- or over-expression phenotypes and are currently of undetermined function. There is currently little information on how transferable and conserved such complex developmental transcriptional modules identified from transcriptional data are likely to be, especially in plants.

Due to their key role in controlling developmental programs and their ability to influence the expression of a number of downstream targets (trans effects), we concentrated our analysis on transcription factors. Transcription factors of interest were selected using OPLS and EST library over-representation tests. This combined approach provides evidence that the C2C2-YABBY, CCAAT-HAP3 and 5, ZF-HD, and MYB families of transcription factor families are of particular importance in leaves, with C2C2-YABBY being the most highly represented family. YABBY transcription factors are known to be involved in the process of abaxialisation [59–61, 65–67] and ectopic expression of YABBY genes is sufficient to cause the abaxialisation of adaxial epidermal tissues [59–61]. Five gene models in the C2C2-YABBY family were identified by the OPLS analysis of the UPSC-TC experiment with one of those (grail3.001817701) also identified by the EST library over-representation test. Three of the five gene models are orthologs of AT2G45190 (AFO/FIL/YAB1) and the other two of AT2G26580 (YAB5). Although no functional evidence is available for these genes in Populus yet, homologous Arabidopsis genes have been chacterized. [61] report that there was no observable leaf phenotype in the fil single mutant but suggest this is the result of redundant activities of other members of the YABBY family. [60] report aberrant leaf phenotypes for 35S over-expression lines of FIL, with strong over-expression apparently being lethal and weaker over-expressing lines showing the formation of wrinkled leaves. The expression of YAB5 was recently suggested to be negatively modulated in the adaxial domain by AS1 and AS2, two genes critical for the development of properly expanded leaves [68]. Two of the identified transcription factors have mutant phenotypes reported in TAIR. Eugene3.00021070 is the ortholog of AT4G37740, Growth-Regulating Factor 2, the mutants of which have smaller leaves [69], and which is one of the predicted targets of the ptc-MIR396 miRNA family that we found to be over-represented for target sites in the OPLS leaf genes. EstExt fgenesh4 pm.C LG VI0283 is the ortholog of ASYMETRIC LEAVES 1 (AS1, AT2G37630), which is homologous to the Antirrhinum PHANTASTICA (PHAN) and maize ROUGH SHEATH2 (RS2) genes [70–72] and which acts as a negative repressor of class I KNOTTED1-like homeobox (KNOX, in particular KNAT1 and KNAT2) genes and that shows genetic interaction with SHOOT MERISTEMLESS (STM; [70]) and BREVIPEDICELLUS (BP; [18]). as1 mutants have highly lobed leaves, occasionally forming ectopic shoots on leaves [70, 73] show that AS1 and AS2 together with ERECTA (ER) are involved in the establishment of adaxial-abaxial polarity (with ER appearing to protect the AS1/AS2 pathway from heat stress; [74]). [75] show interaction of AS1 with CUP-SHAPED COTYLEDON (CUC) 1, a gene involved in shoot apical meristem formation. It would therefore appear that AS1 is a key and central player in leaf formation and pattern development. Through a combination of the various analyses undertaken in this study, it would appear that miRNA family ptc-MIR396 and the GRF gene family may represent interesting targets for follow-up studies.

As a means of indicating how generally-applicable and relevant our finding are, we were interested to see what degree of overlap would be found between the UPSC-TC and Pt-TC experiments if the same OPLS analysis was carried out on the Pt-TC experiment. We did not expect a substantial overlap as the two experiments have some key differences: The UPSC-TC experiment samples tissues to a greater resolution and, for the leaf samples, includes a set of leaves that are considerably younger than the young leaf samples in the Pt-TC experiment (see GEO for sample details of the Pt-TC experiment). However, 452 (40.5%) of genes were found in common between the UPSC-TC and Pt-TC experiments, which suggests that the analysis method is robust and applicable beyond the samples contained in the original experiment (data not shown). This also suggests that our identified genes are likely to be functionally important within the SwAsp aspen trees that we wish to use for future association mapping, which was a major motivation for this undertaking rather than simply selecting candidate genes on the basis of published forward genetics screens.

As an example of the potential uses of the information presented by this work, we examined the collocation of identified candidate genes and transcription factors to QTL for leaf development that we had mapped in previous work [9] and to those additional QTL available in the PopGenIE Browser. This allowed us to identify two collocating transcription factors, one of which (AFO) has functional support for a role in leaf development in Arabidopsis [59–61]. We have also shown that multiple members of a miRNA family identified as targeting a number of genes within the OPLS leaf gene set additionally collocate to QTL for leaf development, raising the potential for miRNAs being the underlying genetic loci of a QTL.

There are at least six publications containing QTL data for leaf related traits in the mapping population considered here alone (and a greater number if one considers all results published in any Populus mapping population), but in no case (our own publication included) would this collocation analysis have been possible using only the published data. We would therefore strongly suggest that a MIAME-equivalent standard should be established for the publication of QTL results to ensure that map, QTL, and phenotype data are always made publicly available as a requirement for the publication of QTL results, especially with the increasingly availability of eQTL results in addition to more classical phenotypic QTL data as this will ensure that such cross-species, comparative genomics studies are achievable in the same way that they are in the transcriptomics field thanks to standards compliance.

We have shown that a large, diverse collection of microarray data can be used in a combined-analysis approach to identify candidate genes involved in the processes of leaf development. Network-clustering analysis on a set of genes that are highly-representative of leaves identified modules with distinct patterns of expression associated with leaf development and functional enrichment analysis of genes within each module revealed that the pattern of expression between the modules makes biological sense in the context of their likely functional roles. Our cross-experiment and cross-species validation approach showed that identified genes, particularly the transcription factors, have patterns of expression that suggest they represent a robust set of candidate genes involved in leaf development. The analysis approach used identified genes with proven functional roles in leaf development in Arabidopsis and other species, suggesting that our method is biologically robust and the results meaningful. This suggests that the identified genes of as-yet unknown function are sensible targets for further targeted functional characterisation or as candidates for association mapping, as well as providing a set of genes that can serve as robust markers of leaf developmental age.