extfieldsofstudy: Medicine Biology

text:
The RRM domain in GW182 proteins contributes to miRNA-mediated gene silencing

Proteins of the GW182 family interact with Argonaute proteins and are required for miRNA-mediated gene silencing. These proteins contain two structural domains, an ubiquitin-associated (UBA) domain and an RNA recognition motif (RRM), embedded in regions predicted to be unstructured. The structure of the RRM of Drosophila melanogaster GW182 reveals that this domain adopts an RRM fold, with an additional C-terminal α-helix. The helix lies on the β-sheet surface, generally used by these domains to bind RNA. This, together with the absence of aromatic residues in the conserved RNP1 and RNP2 motifs, and the lack of general affinity for RNA, suggests that the GW182 RRM does not bind RNA. The domain may rather engage in protein interactions through an unusual hydrophobic cleft exposed on the opposite face of the β-sheet. We further show that the GW182 RRM is dispensable for P-body localization and for interaction of GW182 with Argonaute-1 and miRNAs. Nevertheless, its deletion impairs the silencing activity of GW182 in a miRNA target-specific manner, indicating that this domain contributes to silencing. The conservation of structural and surface residues suggests that the RRM domain adopts a similar fold with a related function in insect and vertebrate GW182 family members.

INTRODUCTION
Proteins of the GW182 family are essential components of mRNA processing bodies or P-bodies (1)(2)(3). They interact with the Argonaute proteins and are required for miRNAmediated gene silencing in animal cells (3)(4)(5). GW182 proteins are characterized by the presence of two to three distinctive blocks of glycine-tryptophan repeats (referred to as N-terminal, middle and C-terminal GWrepeats) and two predicted structured domains: a central ubiquitin-associated (UBA) domain, and a C-terminal RNA recognition motif (RRM). Furthermore, a predicted unstructured, glutamine-rich (Q-rich) region lies between the UBA and the RRM domains [ Figure 1A; (1,6)].
The role of GW182 proteins in the miRNA pathway is well established in Drosophila melanogaster cells. Here, depleting the sole GW182 family member found in insects suppresses silencing of miRNA targets, irrespective of whether they are translationally repressed or directed to degradation (5)(6)(7)(8)(9). In these cells depleted of GW182 there is no decrease in the expression levels of AGO1 [the Argonaute protein dedicated to the miRNA pathway in D. melanogaster, (9)], suggesting that Argonaute proteins cannot silence miRNA targets in the absence of GW182. These and additional observations lead us to propose that target silencing by miRNAs is effected by a protein complex consisting minimally of an Argonaute and a GW182 protein (5,9).
The N-terminus of all members of the GW182 family, characterized by several GW-repeats, mediates the *To whom correspondence should be addressed. Tel: +49 7071 601 1350; Fax: +49 7071 601 1353; Email: elisa.izaurralde@tuebingen.mpg.de interaction with Argonaute proteins (6,18). How the other domains of GW182 proteins affect silencing activity is less well understood. For D. melanogaster GW182, we have shown that the N-terminal region containing GW repeats, together with the UBA domain and the Q-rich region is necessary and sufficient to localize GW182 to P-bodies (6). Furthermore, this protein fragment promotes the accumulation of AGO1 in P-bodies, a process that depends on the interaction with GW182 in D. melanogaster cells (6).
In C. elegans, AIN-1 and AIN-2 interact with the Argonaute proteins and function in the miRNA pathway, yet these proteins lack a UBA domain, the Q-rich region and the RRM domain. This raises the question: how do these domains affect the activity of vertebrate TNRC6A-C and insect GW182?
In this study, we characterized the RRM domain of D. melanogaster GW182. We show that this domain adopts an RRM fold in solution, with an additional C-terminal a-helix shielding the b-sheet surface, which in canonical RRMs is involved in RNA binding. This, together with the lack of general affinity for RNA and the absence of aromatic residues in the conserved RNP1 and RNP2 motifs, suggests that the GW182 RRM is not involved in RNA recognition. Rather, this domain may engage in protein-protein interactions through an unusual hydrophobic cleft exposed on its helical side. We also show that the RRM domain is not required for GW182 to interact with AGO1 or with miRNAs. Furthermore this domain is dispensable for the accumulation of both GW182 and AGO1 in P-bodies. Unexpectedly, however, the GW182 RRM contributes to silencing of a subset of miRNA targets.

MATERIALS AND METHODS
Cloning, expression and purification of the RRM domain of D. melanogaster GW182 The sequence encoding the RRM domain of D. melanogaster GW182 (gi:18447359; Trp1114 to His1198) was cloned into the pETM60 vector (derived from pET24-d; Novagen). The protein was expressed in the Escherichia coli strain BL21(DE3) Rosetta II at 20 C overnight. To uniformly label GW182 RRM with 15 N/ 13 C or 15 N, cells were grown in M9 minimal medium supplemented with 15 NH 4 Cl with or without 13 C 6 -glucose. The protein was purified by affinity chromatography using Ni-NTA (Ni-nitrilotriacetic acid) HiTrap chelating HP columns (GE Healthcare). The tag was then cleaved by overnight exposure to TEV protease. The protein was purified to homogeneity by two subsequent gel filtrations using a HiLoad 26/60 Superdex 75 preparative-grade column (GE Healthcare). The purity of the resulting protein, consisting of the cloned sequence plus an additional four residues at the N-terminus, was confirmed by sodium dodecyl sulfate-polyacrylamide gel electrophoresis. Samples for nuclear magnetic resonance (NMR) spectroscopy were prepared at 1 mM in PBS (pH 7) containing 1 mM dithiothreitol.

NMR measurement and structure calculation
All spectra were recorded at 298 K on Bruker AVIII-600 and AVIII-800 spectrometers. Backbone sequential assignments were completed using standard triple resonance experiments, implemented using selective proton flipback techniques for fast pulsing (19). Aliphatic sidechain assignments were completed by a combination of CCH-COSY and CCH-TOCSY experiments, while aromatic assignments were made by linking aromatic spin systems to the respective C b H 2 protons, in a 2D-NOESY spectrum. Stereospecific assignments and the resulting 1 rotamer assignments were determined for 34 of 57 prochiral C b H 2 protons and for the C g H 3 groups of all three valine residues. Assignments of 1 rotamers were also available for all four isoleucine residues and six of seven threonine residues. Assignment of 2 rotamers were made for all 4 isoleucine and 8 of 11 leucine residues.
Distance data were derived from a set of five 3D-NOESY spectra, including the heteronuclear edited NNH-, CCH-and CNH-NOESY spectra (20) in addition to conventional 15 N-and 13 C-HSQC-NOESY spectra, and a 2D-NOESY spectrum recorded on an unlabeled sample. Distance and dihedral angle restraints applied for the 82 high-confidence predictions found by the program TALOS (21) and 38 hydrogen bond restraints were derived as detailed elsewhere (22). Refinement was carried out by comparing all five of the experiments and back-calculated 3D-NOESY-HSQC and 3D-HSQC-NOESY-HSQC spectra (in house software). This process adjusted side-chain rotamers for several residues.
Structures were calculated with XPLOR (NIH version 2.9.4) using standard protocols with modifications as described (23). For the final set, 50 structures were calculated and 23 chosen on the basis of lowest restraint violations. An average structure was calculated and regularized to give a structure representative of the ensemble.

DNA constructs and transfection of S2 cells
Luciferase reporters and plasmids for expression of miRNAs and of N-HA protein fusions were described before (6,7,24). GW182-ÁRRM was generated by site-directed mutagenesis using the plasmid pAc5.1B-N-HA-GW182 as template, and the Quick Change mutagenesis kit from Stratagene.
Transfections of S2 cells were performed in 6-well plates, using Effectene transfection reagent (Qiagen). The transfection mixtures contained 0.1 mg of firefly luciferase (F-Luc) reporter plasmid, 0.4 mg of the Renilla transfection control and 0.5 mg of plasmids expressing miRNA primary transcripts or the corresponding vector without insert. When indicated, 0.025 mg of plasmids expressing recombinant proteins were cotransfected as indicated in the figure legends.
Firefly and Renilla luciferase activities were measured 3 days after transfection using the Dual-Luciferase Reporter Assay System (Promega). Total RNA was isolated using TriFast (Peqlab Biotechnologies).

Immunofluorescence
Three days after transfection, S2 cells were allowed 15 min to adhere to Poly-D-Lysine-coated coverslips, washed once in serum-free medium and fixed with 4% paraformaldehyde in PBS for 15 min, followed by 5 min incubation in methanol at À20 C. After fixing, cells were washed in PBS, permeabilized for 5 min with PBS containing 0.5% Triton X-100 and washed again with PBS. Cells were stained with affinity-purified anti-Tral antibodies (25) diluted 1:250 in PBS containing 1% BSA. Alexa-594coupled goat anti-rat secondary antibody (Molecular probes, catalog number A11007) was used at a dilution of 1:1000. HA-tagged proteins were detected with a monoclonal anti-HA antibody (Covance Research Products, catalog number MMS-101P) diluted 1:1000 in PBS containing 1% BSA. Alexa Fluor 594-coupled or Alexa Fluor 488-coupled goat anti-mouse secondary antibody (Molecular probes, catalog numbers A11001 and A11005, respectively) was used at a dilution of 1:1000. Cells were mounted using Fluoromount-G (Southern Biotechnology Associates, Inc.). Images were acquired using a Leica TCS SP2 confocal microscope.

Data Bank accession numbers
The chemical shifts and coordinates of the structure have been deposited in the BioMagResBank (accession code 16206) and PDB (accession code 2WBR).

Solution structure of the RRM domain of D. melanogaster GW182
The RRM domain of D. melanogaster GW182 (GW182 RRM, residues S1116 to H1198) adopts an RRM fold in solution comprising the typical b 1 a 1 b 2 b 3 a 2 b 4 topology, where a four-stranded antiparallel b-sheet is flanked on one side by two a-helices [a1 and a2; Figure 1B; (26,27)]. The domain also contains a third, C-terminal a-helix (a3), which covers the face of the b-sheet that in canonical RRMs is involved in RNA-binding ( Figure 1B). A similar C-terminal a-helix lies in the same location in a number of other RRM domains [reviewed in (27)]; of those, the N-terminal RRM domain of U1A [ Figure 1C; (26,28,29)], the RRM3 domain of U2AF 65 [ Figure 1D; (30)] and the quasi RRM1 domain of hnRNP F [ Figure 1E; (31)] are well characterized and superpose with root mean square deviation (r.m.s.d.) values of 2.28, 2.38 and 2.42 Å , respectively (for alignable residues). GW182 RRM also contains a short b-strand (b3 0 ) in loop L5, which forms a small b-hairpin (b-turn of type II) with the N-terminal portion of b-strand b4 ( Figure 1B).
The ensemble of 23 lowest energy NMR structures (Supplementary Data 1) has a good r.m.s.d. (calculated over 80 structured residues) of 0.39 Å for backbone atoms and of 0.83 Å for all non-hydrogen atoms ( Table 1). The restraint violations are also very low, with the final set having, on an average, three violations of distance restraints >0.09 Å per structure and one dihedral restraint violation >1 Å .

GW182 RRM lacks RNA binding properties
Canonical RRM domains are characterized by two conserved sequence signatures: Figure 2A; reviewed in (27)]. These strands provide aromatic side chains to the surface of the b-sheet (positions 3, 5 in RNP1 and position 2 in RNP2, bold characters). In GW182 RRM, the RNP1 (QGIALCKY) and RNP2 (LLLKNL) motifs deviate from the consensus signature in that the aromatic residues that usually play a key role in nucleic acid binding (27) are replaced by aliphatic residues (L1120 in b1, I1154 and L1156 in b3; Figure 2A and B). Furthermore, these residues together with W1118 (in b1) contribute to hydrophobic interactions with the C-terminal a-helix (residues V1192 and I1195 in a3; Figure 2B, red dots). The C-terminal a-helix is further stabilized by a hydrogen bond between the side chains of H1198 in a3 and Y1149 in b2 ( Figure 2B). All of these residues are either conserved or conservatively substituted in vertebrate TNRC6A-C proteins (Figure 2A), indicating that in GW182 RRMs the interaction of helix-3 with the b-sheet surface is a conserved feature (Figure 2A and B).
Given these tight aliphatic methyl-methyl contacts (349 Å 2 total buried surface), we consider it unlikely that helix a3 would move away from the observed position, like it does when the U1A RRM domain binds RNA [ Figure 1C; (28,29)]. Even if it did, the exposed b-sheet surface would be hydrophobic and would not expose any classical RNAbinding features ( Figure 2B).
GW182 RRM is therefore similar to the quasi RRM domains (qRRM1 and qRRM2) of hnRNP F, where the a3 helix is equally unlikely to liberate the b-sheet surface [ Figure 1E; (31)]. Nevertheless, the qRRM1 domain of hnRNP F binds RNA. This RNA binding occurs at an alternative surface of the domain that is positively charged and includes the equivalent of the b3 0 /b4 hairpin loop as well as loops L1 and L3 ( Figure 1E). Mutagenesis and NMR-shift perturbation experiments revealed both basic and aromatic side chains are involved in RNA binding [ Figure 2A, magenta residues: R16 (b1), W20 (L1), R52 (L3), R75-H80-R81 and Y82 (hairpin b3 0 -b4), and F86 (b4); (31)]. In the GW182 RRM none of these residues is conserved (Figure 2A). Furthermore, the GW182 RRM lacks any positively charged surface patch that would be indicative of RNA binding ( Figure 1B versus 1E, lower panels), arguing against an hnRNP F-like mode of RNA recognition.
In agreement with the assumption that GW182 RRM is unlikely to bind RNA either in a canonical or in the quasicanonical mode, no detectable RNA-binding activity was observed in vitro. Indeed, the GW182 RRM did not detectably interact with an oligomeric RNA consisting of 15 or 30 uracyls (poly(rU) 15,30 , as shown by analytical size exclusion chromatography measurements (data not shown). We also tested RNA binding by NMR, which detects low-affinity interactions (34); however, no binding (i.e. chemical shift changes) was detected, even when up to an equimolar amount of poly(rU) 15 or poly(rU) 30 was added to purified GW182 RRM (Supplementary Data 2). Although the possibility that GW182 RRM binds a specific RNA sequence motif cannot be ruled out, the lack of general affinity for RNA together with the absence of critical RNA-binding residues suggest that this domain does not bind RNA.

GW182 RRM may interact with other proteins via an unusual hydrophobic cleft
Several studies showed that some RRMs mediate protein-protein interactions, in addition to or instead of recognizing RNA [reviewed in (27)]. This is the case, for instance, for the RRM domains of the exon-junction complex protein Y14 (35), the nonsense-mediated mRNA decay effector UPF3 (36) and the C-terminal RRM (RRM3) of human U2AF 65 , which also bears a C-terminal a-helix (a3) across its b-sheet surface [ Figure 1D; (30)].
While Y14 and UPF3 bind their protein partners via the b-sheet surface generally used by RRMs to bind RNA (35,36), the RRM3 of U2AF 65 recruits the N-terminal peptide of Splicing Factor 1 (SF1) via its a-helical surface. The C-terminal end of the SF1 peptide inserts into a hydrophobic pocket between a-helices a1 and a2 and the adjacent loops L1 and L5 [ Figure 3C; (30)].
Interestingly, GW182 RRM has an even larger hydrophobic cleft between the a-helices and adjacent  (30). The C-terminal a-helix is colored in orange. Side chains of conserved aromatic residues in RNP1 and RNP2 motifs are shown as yellow sticks. Lower panels: corresponding surface representations colored according to the electrostatic surface potential (blue and red colors indicate positive and negative potential, respectively, ramped from À5 kT/e to +5 kT/e using the Adaptive Poisson Boltzmann Server tool (32). All structure representations were done in PyMol (http://www.pymol.org).
loops, which is rather unique among RRM domains ( Figure 3A-D). Compared with U2AF 65 RRM3 the space is liberated primarily by the orientation of the b3 0 /b4 hairpin (loop L5) such that side chains do not reach to contact residues on the opposing a-helix a1 ( Figure 3A). The entire cleft is lined by aliphatic side chains [ Figure 2A, residues in red: L1121 (b1), L1133 and L1136 (a1), L1172 (a2), L1177 (b3 0 ) and I1182 (b4)], creating a hydrophobic patch of considerable size, which is usually not exposed to the solvent ( Figure 3A). Each of these residues is conserved in vertebrate TNRC6s (Figure 2A). It is therefore likely that this large hydrophobic cleft represents a conserved site for protein-protein interactions, which could be relevant for miRNAmediated gene silencing.

The GW182 RRM is dispensable for the interaction with AGO1 and miRNAs
To investigate the role of the RRM domain in GW182 function, we generated a protein mutant lacking the complete domain (GW182-ÁRRM) and tested its interaction with AGO1 and miRNAs. To this end, wild-type and mutant HA-epitope-tagged GW182 were transiently expressed in S2 cells and immunoprecipitated from cell lysates using anti-HA antibodies. We then tested for co-immunoprecipitation of endogenous AGO1 and miRNAs by western and northern blot, respectively. We observed that GW182-ÁRRM interacted with endogenous AGO1 as efficiently as wild-type GW182 ( Figure 4A, lanes 6 versus 5). Similarly, GW182-ÁRRM co-immunoprecipitated endogenous bantam ( Figure 4A, lanes 6 versus 5). Together these results indicate that the RRM domain of GW182 is not required for the interaction with AGO1 and miRNAs.

GW182 RRM is dispensable for P-body localization
In previous studies, we showed that an N-terminal fragment of GW182 containing the N-terminal GW-repeats, the UBA domain and the Q-rich region, was necessary and sufficient to direct the protein to P-bodies (6). This protein fragment was also necessary and sufficient to promote the accumulation of AGO1 in P-bodies, a process that depends on GW182 (6). As expected from these previous studies (6), GW182-ÁRRM localized to cytoplasmic foci, which correspond to endogenous P-bodies as judged by the staining with anti-Trailer hitch (Tral) antibodies [ Figure 4B and C; (3,25)]. GW182-ÁRRM also promoted the accumulation of AGO1 in P-bodies ( Figure 4D-F), indicating that neither P-body localization nor the AGO1 interaction is impaired in this mutant, in agreement with the immunoprecipitation assays and with previous reports (6,8).

GW182 RRM contributes to target silencing by miRNAs
To investigate whether the RRM domain of GW182 plays a role in silencing, we established a complementation assay. In this assay, miRNA function is monitored using reporters in which the F-Luc Open reading frame (ORF)

. T Q G T A L I R Y S T K Q E A A K A Q T A L H M C V L G N T
. . is cloned upstream of 3 0 -UTRs of D. melanogaster genes regulated by miRNAs. In particular, we used the previously characterized reporters F-Luc-CG10011 and F-Luc-CG3548 (silenced by miR-12) and F-Luc-Nerfin (silenced by miR-279 and miR-9b). We depleted endogenous GW182 using dsRNAs targeting the 5 0 -and 3 0 -UTR sequences of GW182 mRNA. miR-12 directs the F-Luc-CG10011 and F-Luc-CG3548 reporters to degradation, whereas miR-279 and miR-9b silence the F-Luc-Nerfin reporter mainly at the translational level ( Figure 5A-D). Depleting endogenous GW182 suppressed silencing of these reporters and lead to a 3-to 8-fold increase of F-Luc activity [ Figure 5E-H; (6)(7)(8)(9)]. We then tested whether silencing could be restored by expressing either wild-type GW182 or the mutant lacking the RRM. These proteins were expressed from vectors including only the GW182 ORF, and thus were insensitive to the dsRNAs used to deplete the endogenous protein.
We found that in cells depleted of endogenous GW182, expressing wild-type GW182 restored silencing to levels comparable with those observed in control cells ( Figure 5E-H). In contrast, the GW182 mutant lacking the RRM failed to efficiently restore silencing of either the F-Luc-CG10011 by miR-12, or of the F-Luc-Nerfin reporter by miR-279 ( Figure 5E and G). The GW182-ÁRRM mutant fully restored silencing of the F-Luc-CG3548 reporter and of the F-Luc-Nerfin reporter when silencing was triggered by miR-9b ( Figure 5F and H). The expression levels of wild-type and mutant GW182 proteins were comparable ( Figure 4A). Together these results indicate that the RRM domain of GW182 is not essential, but contributes to silencing in a manner that is specific to miRNA-target pairs. This is so, independently of whether the target is silenced at the translational level or directed to degradation.

DISCUSSION
Proteins of the GW182 family are essential components of the miRNA pathway. They are required both to repress translation and to enhance degradation of miRNA targets (3,9). The insect and vertebrate members of this protein family are characterized by two domains predicted to be structured: a UBA and an RRM domain (1,2,4,6). Here, we show that the RRM domain of D. melanogaster GW182 adopts a classical RRM-fold, but lacks RNA-binding features. This domain is not essential for silencing, but rather contributes to silencing of a subset of miRNA-target pairs.

The GW182 RRM lacks RNA-binding features
The RRM is an ancient protein domain of about 90 amino acids; it is present in all kingdoms of life and particularly prominent in eukaryotes. The core of the domain consists of one four-stranded anti-parallel b-sheet packed against two a-helices (26,27). Originally identified as an RNAbinding module, this domain can mediate protein-protein interactions in addition to, or instead of binding RNA [reviewed in (27)]. RRM-RNA interactions often involve two RNA bases stacking with the ring of conserved aromatic residues in the RNP motifs I and II; additional contacts may vary amongst different RRMs. In GW182 RRM, the conserved aromatic residues in the RNP1 and RNP2 motifs are substituted with aliphatic residues. This, together with the lack of detectable RNA-binding affinity, the absence of regions with positive electrostatic surface potential and the packing of the C-terminal a-helix against the b-sheet surface, suggest that GW182 RRM does not bind RNA.
A C-terminal a-helix masking the canonical RNAbinding surface was observed in other RRMs, including the N-terminal RRM of U1A (29), the CstF-64 RRM (37) and the qRRM1 of hnRNP F (31). Despite the a-helix, however, these domains do bind RNA.
In the U1A RRM, the C-terminal a-helix moves into an open conformation upon RNA binding (29), whereas in CstF-64 the a3 helix unfolds in the presence of RNA (37). In both cases, conserved aromatic residues in the RNP1 and RNP2 motifs become exposed and available for basestacking interactions with RNA. The tight aliphatic interactions between helix a3 and the b-sheet surface of GW182 RRM makes it unlikely that this helix would swing away from the b-sheet as observed in U1A RRM. Finally, even if the C-terminal a-helix of GW182 would adopt an open conformation or unfold, the exposed b-sheet surface of GW182 lacks key aromatic and positively charged residues that could mediate RNA binding.
The C-terminal a-helix of the qRRM1 of hnRNP F does not change conformation in the presence of RNA (31). Nevertheless, as mentioned above, the qRRM1 of hnRNP F still binds RNA through aromatic and positively charged residues that are not on the b-sheet surface but are part of a short b-hairpin and loops L1 and L3 (31). These residues are not present in GW182 RRM, indicating this domain is also unlikely to bind RNA in the quasicanonical mode.

GW182 RRM: a protein interaction domain?
A survey of available structures of RRM domains in complex with other proteins revealed that RRM-protein interactions are very diverse, with no common features emerging [reviewed in (27)]. For instance, both Y14 and UPF3 RRMs interact with their binding partners (MAGOH and UPF2, respectively) through their b-sheet surfaces; however, the interaction interface between Y14 and MAGOH is mainly hydrophobic, whereas, the UPF3-UPF2 interaction is mediated by charged residues (35,36). Other RRMs use the helical surface for protein interactions. For instance, U2AF 65 RRM3 interacts with SF1 through residues located on helices a1 and a2 (30). In other cases, RRM-protein interactions are mediated by extensions or variations of secondary structural elements (27).
It is currently unclear how the GW182 RRM would interact with putative protein partners; but similar to U2AF 65 RRM3, the structure of the GW182 RRM reveals a hydrophobic cleft on the side of the b-sheet opposite to where helix a3 lies. This cleft is likely to engage in interactions with as yet unidentified partners.

Conservation of the GW182 RRM
The GW182 RRM is highly conserved (e.g. 58% identity between D. melanogaster GW182 and human TNRC6A), suggesting that in insect and vertebrate members of the GW182 protein family this domain adopts a similar fold with a similar function. In particular, the occlusion of the putative RNA-binding surface by the C-terminal a-helix and the presence of a hydrophobic cleft on the helical face are structural features conserved in RRM domains of GW182 proteins. Despite this conservation, the RRM is not essential for the silencing activity of GW182. Indeed, a GW182 protein lacking this domain can partially or fully rescue silencing in cells depleted of endogenous GW182. Nevertheless, it is still possible that this domain is important for silencing of only a specific subset of miRNA targets.
The conservation of the RRM domain in insects and vertebrates contrasts with the absence of this domain in the C. elegans proteins AIN-1 and AIN-2 (4,6). One possible explanation is that this domain plays a critical role in silencing a subset of miRNA targets only in insects and vertebrates. Alternatively, the C. elegans AIN-1 and AIN-2 may function in conjunction with additional partners that compensate for the lack of the RRM. Finally, the vertebrate and insect GW182 proteins may act in processes unrelated to silencing, which require the RRM domain. Indeed, inhibiting GW182 function by antibody injection in D. melanogaster syncytial stage embryos, leads to abnormal nuclear divisions, a phenotype that might or might not be a consequence of inhibiting the miRNA pathway (38).

Role of the GW182 RRM domain in silencing
In previous studies we showed that silencing by miRNAs is effected by AGO1-GW182 complexes (9). These complexes repress translation and/or enhance degradation of miRNA targets (3,9). GW182 RRM contributes to both of these mechanisms, as a GW182 protein that lacks the RRM is impaired in mediating both translational repression and mRNA decay. Although we cannot rule out that deleting the RRM domain affects GW182 folding, we consider this possibility unlikely because GW182-ÁRRM still interacts with AGO1 and miRNAs, localizes to P-bodies and silences a subset of miRNA targets. This raises the question: what step of silencing is affected by the deletion of the RRM domain? One possibility is that this domain contributes to target recognition by RISC, either by interacting with proteins bound to the mRNA target or by interfering with their binding, thereby increasing target accessibility. Alternatively, this domain may contribute to silencing downstream of target binding. For instance, the RRM may facilitate the interaction of RISC with either components of the translation machinery or general mRNA decay factors, and thus contribute to the effector step of silencing.


text:
ExprAlign - the identification of ESTs in non-model species by alignment of cDNA microarray expression profiles

Background Sequence identification of ESTs from non-model species offers distinct challenges particularly when these species have duplicated genomes and when they are phylogenetically distant from sequenced model organisms. For the common carp, an environmental model of aquacultural interest, large numbers of ESTs remained unidentified using BLAST sequence alignment. We have used the expression profiles from large-scale microarray experiments to suggest gene identities. Results Expression profiles from ~700 cDNA microarrays describing responses of 7 major tissues to multiple environmental stressors were used to define a co-expression landscape. This was based on the Pearsons correlation coefficient relating each gene with all other genes, from which a network description provided clusters of highly correlated genes as 'mountains'. We show that these contain genes with known identities and genes with unknown identities, and that the correlation constitutes evidence of identity in the latter. This procedure has suggested identities to 522 of 2701 unknown carp ESTs sequences. We also discriminate several common carp genes and gene isoforms that were not discriminated by BLAST sequence alignment alone. Precision in identification was substantially improved by use of data from multiple tissues and treatments. Conclusion The detailed analysis of co-expression landscapes is a sensitive technique for suggesting an identity for the large number of BLAST unidentified cDNAs generated in EST projects. It is capable of detecting even subtle changes in expression profiles, and thereby of distinguishing genes with a common BLAST identity into different identities. It benefits from the use of multiple treatments or contrasts, and from the large-scale microarray data.

Background
Transcript screening investigations have traditionally been led by sequence analysis of cDNA clone collections to define the identity of hybridisation probes included on microarrays for expression profiling [1]. Despite this, all eukaryotic EST collections contain large proportions of transcripts (~50%) that remain unidentified by unattended BLAST protocols. Some of these may represent new, undiscovered protein-coding or non-protein-coding transcripts [2][3][4]. Others may arise from untranslated regions of coding sequence RNA, which being non-conserved fail to align to reference databases. Finally, some may be concatenated constructs generated artefactually during the production of cDNA libraries.
We have experienced these kinds of difficulty in our analysis of ESTs from the common carp, Cyprinus carpio L., a well-used model species for research into environmental responses [5], and which is the subject of a substantial aquaculture interest for both food and ornamental uses. The common carp genome is widely thought to have become duplicated within the previous 12-15 Mya, and many duplicate paralogs are retained [6][7][8] to complicate the analysis. We originally generated a medium-scale collection of ~13.5K directional, cDNA clones from multiple tissues [9], though this has more recently been increased [10]. 9,202 directional EST were assembled into 6,033 transcriptional units. Of these only 3,252 were BLASTidentified leaving 2,701 as unclassified, many of which displayed interesting expression properties in response to a range of chronic stress treatments [9].
Additional information regarding the identity of ESTs may come from the comparison of expression profiles of one microarray probe with another since different probes arising from the same gene should have very highly correlated profiles whilst probes with an identical BLAST identity but arising from different members of a gene family might present divergent expression profiles. Either way, co-expression indices can be used as evidence in seeking an identity for a BLAST-unidentified cDNA clone, and can separate putative isoforms. To explore the limits of expression profiling, and the extent to which dissimilar but coregulated genes may confound the process, we have accumulated data from a very large number of microarray hybridisations, including RNA from all of the major organs of common carp exposed to a range of environmental stresses, including chronic cooling [9], chronic hypoxia [11] and starvation/refeeding protocols. This large dataset represents a substantial data resource that can be used to suggest gene identity through correlation analysis. Here we describe the Expression Alignment (ExprAlign) technique for assigning a putative gene identity, which, following the pioneering work of Kim et al. [12,13], is based on the clustering of gene expression profiles [14][15][16][17]. This resolves a number of issues relating to the identification of probes that were unidentified by conventional unattended BLASTx procedures, including those from untranslated regions of transcripts.

ESTs resource and common carp microarray data
We used the EST resources from carpBASE 2.1, which was constructed by the EST analysis package EST-ferret 2.1 http://legr.liv.ac.uk. This comprised the 13,349 directional cDNA clones already described [10], of which 9202 were 5' end sequenced, BLASTx identified and annotated with gene ontology, KEGG and CDD terms.
The cDNA microarray used in this work has been described in [9] and [11], and comprised 13440 PCRamplified cDNA probes, including blanks and standards. The raw expression data has been deposited in ArrayExpress E-MAXD-1 and E-MAXD-10, respectively. The gene expression data used in this analysis comprised 707 common carp RNA samples, hybridised to 1414 cDNA microarrays, all using a reference-based, dye-swap design against a common reference using dye-swap, and with 4fold or greater biological replication. These experiments were conducted with ethical approval and corresponding personal and project licences of the Home Office, U.K 189 RNA samples were generated from the study of chronic cold stress [9], including samples for time-course after transfer from a preconditioning temperature of 30°C to 23°C, 17°C, and 10°C. Tissues examined were brain, gill, heart, intestine, kidney, liver and skeletal muscle. 414 RNA samples were used in a time-course study of hypoxia stress [11] in 4 tissues (brain, heart, intestine, liver and skeletal muscle), conducted at two temperatures (17°C and 30°C). 104 samples were from a study of starvation over a 6-week time course during which the animals were starved and then re-fed (I. Hardevig and A.R. Cossins, unpublished). The 'starvation' data included samples for liver and skeletal muscle only. Normalisation and statistical processing of the resulting data has been previously described [9,11] using established techniques [18].

Computing and filtering Pearson's correlation coefficients
The Pearson product-moment correlation coefficient (r) was calculated for every pair of array probes. To speed this up a programme called "CORR" http://legr.liv.ac.uk was written in C, resulting in completion in approximately 40 minutes using a single Linux machine. The Receiver Operating Characteristic (ROC) curves were implemented as described in Additional File 1, Section A, to optimize the threshold of the correlations [19]. The resulting criterion was used to filter the r-values, and the selected relationships were stored in a spreadsheet consisting of a few thousand gene-pairs (the rows), with coefficient scores.
The optimized threshold for filtering Pearson correlation coefficients was determined by plotting the sensitivity (True Positive, P + ) of the comparison against the selectivity (False Positive, P -). The sensitivity P + indicates the probability of the observed true positives at a threshold and the selectivity Pshows the probability for the observed true negatives at a threshold. On the other hand, the probability of the missed true positives was given by (1 -P + ) and the probability of the missed true negatives by (1 -P -). The total probability of missing the true positives and true negatives was given by E = (1 -P + ) + (1 -P -), this providing the optimal threshold. The calculations for P and E are detailed in the Additional File 1, Section A.

Visualising expression alignments
The VxInsight package [20][21][22] was implemented for clustering gene expressions using the Pearson's correlation coefficients and to visualise the alignments. The VxInsight package consists of three parts: VxOrd, VxInsight and VxImport. VxOrd implements the force-directed ordination algorithm to assign the X, Y coordinates for each gene on a 2-dimensional surface based on the clone pair similarities of correlation coefficients. Then the coordinate maps were converted to the 3-dimensional mountain terrains in VxInsight. Finally, VxImport loaded the gene annotation from carpBASE 2.1 into VxInsight for biological interpretations. Figure 1 illustrates the stages used in ExprAlign to generate a landscape depiction of the expression profile, starting from the microarray intensity data and ending with the production of the colour-rendered 3D landscape.

Post ExprAlign analysis
Groups of cDNA clones included within landscape features were assessed using BLASTx, searches as described previously [10] followed by mapping of FASTA sequences onto the zebrafish genome v8 as described by Christoffels et al. [23], or onto common carp fosmid clones BX571686 and BX571725 [24]. Details of Pearsons correlation coefficient of closely related clones, and of mountain allocations, have been included as part of the clone entry in carpBASE2.1.

Correlation analysis
Using normalised data generated from 707 RNA preparations from the cold and hypoxia datasets, each containing 13,440 cDNA probes from multiple tissues we calculated the Pearson correlation coefficient [25] between every pair of probes represented on the microarray. This required 180 million calculations, covering different combination of stress and time, as well as all different tissues. We then used a Receiver Operating Characteristic (ROC) procedure [19] to optimize the threshold that minimises the representation of false positives at 0.858. The coefficients below this value were discarded leaving 30868 correlation coefficients representing 3039 gene probes for construction of the landscape.

Network analysis using 3D landscape for visualisation
We then created a network representation of the r-values that exceeded the ROC threshold, which was ordinated and visualized using VxInsight [20]. This clustered groups of gene probes with high r-values between them using a force-repulsion model and represented them as features in a 3-dimensional landscape metaphor, as an aid to easy interpretation (see Figures 2a and 2b). The height of each landscape feature was an indication of the number of gene probes contained within it, and large-scale features ('mountains', 'hills') can include a number of smaller features ('hillocks'). The 3039 gene probes previously identified comprised 1192 possessing a BLASTx identity whilst the remaining 1847 were unclassifiable. Big mountains were located around the centre of the landscape while smaller mountains were positioned away from the centre. The distance between each feature, both large and small, was representative of the r-values that connect the features. Thus a small distance between landscape features Flow diagram depicting the order of events (in boxes) and informatic tools and packages used (italics) in the ExprAlign pipeline Figure 1 Flow diagram depicting the order of events (in boxes) and informatic tools and packages used (italics) in the ExprAlign pipeline.
indicates that each cluster of gene probes was more similar to each other than if they were more distantly linked.
To assess the significance of the topographical patterns revealed by this process, we permuted the expression table by shuffling the values for all probes across all arrays. The resulting r-values were uniformly low with only 3 pairs of probes with r > 0.858, and 33870 with r > 0.25. We used the latter data to generate a landscape shown in Figure 2c, which showed no structuring. This indicates that the landscape features evident in Figure 2a do not result from random effects in the underlying data, and represent biologically significant outcomes.

Dependency of landscape features upon scale and diversity of datasets
To determine whether the structure of the gene co-expression landscape and the clustering of genes was influenced by the particular dataset used in their construction, we recalculated the landscape using only the data from the 386 arrays generated for the cooling experiment [9]. In this case 4236 microarray probes exceeded the ROC threshold and were ordinated onto the landscape map, comprising 1776 identified and 2460 unclassifiable gene probes. 46 mountains (labeled 'CE', 'cold expression') were generated of which 22 were identified.
We then compared the landscape features for both the GE and CE datasets (provided in reduced form in Figure 3, and shown in full in Additional File 1, Figure S1) using a matrix in which cells contained the number of gene probes that were included in corresponding landscape features. We found 21 highly similar mountain-pairs between GE and CE mountains, indicated by greyed cells. However, we noted some differences between the two landscapes. Thus, the clones in mountain CE209 were separately located in GE142 and GE103; clones in CE178 were separately located in GE111 and GE16b; and the clones in CE149 were separately located in GE142 and GE116. Feature CE26 (see Additional File 1, Figure S1) was linked to three GE features, namely GE35 (apolipoprotein A-I), GE17 (transferrin) and GE13 (transferrin). Moreover, CE119 (creatine kinases) was linked to GE57 and GE24 (both labelled as creatine kinases); and CE121 (parvalbumins) were broken down to GE16 and GE14b (both parvalbumins). The splitting of CE features into two GE features might indicate the separation of isoforms.
We conclude that the discrimination between closely linked gene clusters was significantly affected either by the scale or by the diversity of the data used in their construction. We then tested which of these was the important factor by randomly discarding 50% of the arrays included in the GE dataset. The resulting landscape contained 2444 clones in 27 mountains (labelled 'RE', randomised expression) of which 12 were identified. We then compared the contents of landscape features of the RE dataset with the original GE dataset (summarised in Figure 4 and displayed in full in Additional File 1, Figure S2). This shows that 17 of the RE mountains were related to single features on the GE landscape, and that the 'identity' of these features was conserved. However, 11 of the RE mountains were linked to multiple GE mountains; for example, RE237 and RE174 were each linked to 3 different features in the GE landscape. Thus, whilst the main features of the landscape were robust to reductions in the amount of array data some were sensitive to increasing data complexity, as produced by undertaking more diverse experimental treatments.
A landscape representation of the co-expression profile for the GE analysis of common carp genes Figure 2 A landscape representation of the co-expression profile for the GE analysis of common carp genes. (a) An overview landscape. The white spots, P1 to P8, represent the locations of different parvalbumin isoforms. (b) shows a magnified part of (a) showing labels for each feature displaying the number of clones within that feature. 'Identified' mountains were indicated in blue text and 'unknown' mountains in red. (c) Terrain map derived from permuted and thus random-shuffled array data.

Relating unidentified or unclassifiable clones to identified genes
To test the significance of the representation of identified genes within each of the landscape features we computed p-values for the GE dataset for each gene identity using Fisher's exact test. Table 1 lists the p-values, all of which lay between 10 13 and 10 101 , indicating highly significant enrichment of the identified gene within the mountain.
Secondly, we examined individual clusters either by attended curation of the ESTs followed by BLASTx or by mapping clones onto the zebrafish genome. Thus, for mountain GE13c, only 4 clones recovered a BLASTx hit for the common carp desaturase-2, 2 of them following a frame shift correction. A further 6 were confirmed as belonging to this gene using BLASTn against the 2 available common carp fosmid clones but in the untranslated region, two of them being very short sequence reads. The three remaining clones failed to provide a sequence, but were inferred as being the same gene. Mountain GE10d contained only unidentifiable sequences (thus, not included in the listing of identified mountains in Table 1), all of which belonged to the same carpBASE2.1 contig (1127-2) and which mapped to the same zebrafish location (Chromosome 23: 24767100). This lies within an unprocessed pseudogene (locus i.d. OTTDARG00000030773) so whilst the identity and status of this transcriptional product remains unclear, it is revealed as being expressed in particular common carp tissues and in response to particular environmental influ-Matrix to compare the distribution of selected gene identi-ties in GE and CE landscape co-expression features Figure 3 Matrix to compare the distribution of selected gene identities in GE and CE landscape co-expression features. Additional File 1, Figure S1 comprises all features. The labels for each column or row indicates the number of genes within that landscape feature. The GE features are listed in Table 1. The numbers contained within individual cells indicate the number of microarray probes common to the linked GE and CE features. The denser colour represents a greater level of the agreement.
Matrix for agreement and difference of global mountains and random mountains Figure 4 Matrix for agreement and difference of global mountains and random mountains. This figure includes selected landscape features whilst Additional File 1, Figure S2 comprises all features. Other details are as described in the legend to Figure 3. Column 1 represents the number of gene probes included within the landscape feature indicated in Figure 2, and this is also used to denote the identity. Column 2 indicates how many of these were BLAST identified as indicated in column 3. The number of the best-represented identified gene in a mountain must be >2, and its percentage of clones in mountain must be over 20%. P represents the significance (Fishers Exact test) of the gene being over-represented in the mountain comparing to the whole clone set included within carpBASE 2.1.
ences. Mountain GE10 contained 5 ESTs that were identified by BLASTx as myoglobin, and 5 ESTs that were frame shift corrected also to BLASTx as myoglobin. Finally, in mountain GE113, 50 of the 55 identified genes (i.e. ~91%) were BLASTed as glyceraldehyde 3-phosphate dehydrogenase. The p-value for this identity (6e-101) was very highly significant. However, 5 ESTs showed homology to different genes including microglobulin, adenylate kinse, 60S ribosomal protein, and a glycoprotein precursor. The remaining 58 clones lacking sequence data were thus inferred as having the same identity as the predominant gene.
For large mountains, in particular, we found multiple identified genes being represented, making it difficult to ascribe a single identity to unidentified clones. In some cases we have generated 3-5 sub-clusters using the Kmeans technique, each of which displayed a more singular gene identity with which to infer an identity for the unidentified or un-sequenced ESTs (see below).
Using these methods, we inferred a tentative identity for a total of 522 unidentified ESTs, which represented 17% of the 3039 unidentified ESTs. These are listed in Additional File 1, Table S1. From this it is evident that many clones failed to generate sequence data during the automated sequence analysis, yet they provided informative amplicons for array fabrication.

Discrimination of gene isoforms by ExprAlign
We have explored the extent to which ExprAlign supplements sequence alignment in distinguishing isoforms with gene families. Figure 5 (and Additional File 1, Section B) shows for a number of co-expression mountains the changes in transcript expression of constituent clones as a conventional heatmap across all three experimental conditions.
Aldolase A and B were separately resolved into two neighboring mountains: GE17b containing only S341 (isoform A) and GE33 containing only S488 (isoform B). The close proximity of these mountains in Figure 2 was indicative of the similarity of their expression patterns (see Figure 5), yet ExprAlign cleanly resolved even these subtly different expression patterns in precise agreement with the sequence alignments. Figure 5 shows that resolution between isoforms depended mainly on expression differences in the chronic cold-exposure experiment, and in just 3 tissues; brain, gill, heart and kidney. Thus, aldolase A was up-regulated not only in hypoxia muscle (17°C) but also in hypoxia liver, hypoxia brain, hypoxia intestine (30°C), cold brain, cold gill, cold heart, cold intestine and cold kidney. Aldolase B was up-regulated in cold-intestine, cold heart, hypoxia muscle 17°C, hypoxia intestine 30°C and hypoxia liver. Previous work [29] compared the absolute level of the gene expression of aldolase isoforms, but our study focused on stress-regulated changes.

(iii) Transferrin variant A
Crucian carp (C. auratus) has two transferrin variants [33], whilst silver crucian carp (C. auratus gibelio) has 5 and white crucian carp (C. auratus cuvieri) has 3 [34]. We found that ESTs blasting to transferrin variant A were located in two mountains, GE13 and GE17 that were located in two separate locations but in the same sector of the landscape. Whilst overall the expression profiles (Figure 5) were similar differences were noted in intestinal and cardiac tissues exposed to hypoxia. However, sequence alignment of the ESTs in these two mountains failed to demonstrate any sequence clustering which related to the separate expression mountains. The lengths of the GenBank accessions for carp transferrin variants (AF457152, EU71532-EU71535) were over 2200 bp whilst our ESTs were ~400 bp of which only 200 bp could be used for sequence alignment. So it is most likely that the mountains contain different isoforms and that the sequence domains linked to these expression differences lay outside of the sequenced stretches.

(iv) Apolipoproteins
The apolipoprotein gene family were represented in five different mountains, namely GE11 (containing apolipoprotein C-I precursor, apoC-I), GE14 (apolipoprotein Eb precursor, apo-Eb), GE35 (apolipoprotein A-I, apoA-1), GE125 (14-kDa apolipoprotein, apo-14 kDa) and GE79 (mixed apolipoprotein & its precursor). Pufferfish apoA-I was expressed mainly in liver but apo-14 kDa was mainly Heatmaps for identified GE mountains incorporating data for 707 RNA preparations comprising stress treatments in up to 7 different tissues Figure 5 Heatmaps for identified GE mountains incorporating data for 707 RNA preparations comprising stress treatments in up to 7 different tissues. The numbers on the right side indicates the mountain names which corresponds to the number of gene probes. K1, K2, ʜ, Kn shows the sub-clusters generated by the K-means clustering technique for each mountain. B -brain, G -gill, H -heart, I intestine, K kidney, L liver, and M -muscle. Red and green indicates up-and down-regulated gene expressions respectively. Heatmaps for additional identified GE mountains are shown in Additional File 1, Figure S3.
expressed in liver and less abundantly in brain [35]. We show in common carp that some members of this gene family displayed hypoxia responses notably in the intestine (GE35, Figure 4). By contrast, the features GE11 and GE14 showed stress responses of apoC-I and apo-Eb genes in brain and intestine, respectively. The transcripts of the latter gene were also observed in gill.

Discussion
Approximately 46% of the singletons and assembled contigs in our common carp EST project failed to yield an identity using the unattended BLASTx procedure [10], some of which represent non-overlapping, or 3' sequences of identified genes. Other clones failed to yield a sequence on automated analysis yet provided suitable hybridisation probes on the microarray. Given these special problems and complications in this species due to a recently duplicated genome, the sometimes incoherent nature of subsequent gene losses [24], and the divergent tissue-and response-specific expression patterns so generated [9], we have explored how expression alignment techniques might complement the more usual sequence alignment methods to assign an identity for an otherwise unidentifiable sequence.
Our approach was based on the idea that expression profiles for non-overlapping probes derived from a same gene should be highly correlated when tested across a range of experimental treatments, and this should enable unidentified clones to be identified by comparison with identified clones. Similarly, comparing expression profiles for cDNA microarray probes possessing the same BLASTx identity, offers a means of testing their common identity, given that they may represent unrecognised isoforms or variants of a given gene. Thus, combining alignment of sequence data with that from gene expression data offers a useful means of improving the quality of gene identification, and for discriminating isoforms or members of gene families whose separate identify may not easily be made evident using conventional methods.
For this work we chose to include all available cDNA clones on the microarray, resulting in up to 80 clones per contig, and to gather expression profiles from a wide range of major organs and tissues, exposed to a range of experimental treatments. Consequently, the carp array included the substantial repetition of some genes and this provided greater support for the identified gene clusters. Our approach was based on the comparison of the expression profiles of pairs of probes using Pearsons correlation coefficients which were used to create a network linking genes together on the basis of their shared expression properties. The VxInsight algorithm uses a force-repulsion mechanism to gather the distributed gene networks into discrete clusters, which are then presented in an easy-tounderstand landscape metaphor.
We show that the resulting landscape features, and the associated clusters were robust, first, because permuting and randomising the expression values generated neither high correlation coefficients nor landscape features, and second, because the form of the clusters are largely retained when using different scales of array data from small to large. We show that datasets that contain a wider range of experimental treatments and tissues can fragment the gene clusters into smaller forms, each with a distinctive character. Thus, the exact level of discrimination achieved depends upon the diversity of the data used in its construction, with extra experimental treatments offering additional changes in the expression relationships between genes, thereby refining the resulting correlations.

Gene identification using ExprAlign
Many of the resulting landscape features or mountains generated by VxInsight contained predominantly just one kind of BLAST-identified gene, and we show that there is substantial enrichment of these genes within the features compared to chance alone. Thus unidentified probes within that mountain were also tentatively labelled with that gene identity. Using this approach we were able to impute an identity to 522 unidentified clones in the GE landscape, which represented ~17% of all unidentified clones on the map. The validity of this assignment can be tested by the attended analysis of the clone sequences. This was achieved for mountain GE10 which possessed 5 different probes identified as myoglobin by BLASTx, and another 5 probes lacking an identify. Closer inspection of the corresponding sequences, and manual attempts at alignment, were subsequently able to demonstrate that all of the unidentified ESTs were also myoglobin, including additional examples of a unique brain-specific isoform [11]. This indicates the limitations of unattended BLAST analysis, and where appropriate, the need for manual verification of clones assigned with an ExprAlign identity. A similar result was obtained for parvalbumin [32].

Separation of isoforms using ExprAlign
ExprAlign has also proved useful in separating clones that have been assigned a common BLASTx identity but which have distinctive expression profiles. If clones possessing the same gene BLAST identity were indeed sourced from the same transcriptional start site, then they should display identical expression patterns. On the other hand, if the clones were representative of different isoforms with distinctive expression properties, then they would occupy different features on the expression landscape, or could be distinguished using additional clustering techniques. We show that these expectations are largely met for a series of test cases, as described above. As a further example, in the case of fatty acid-binding protein (mountain GE17c, not shown), we identified 13 cDNA clones all possessing the same BLASTx sequence alignment, but which displayed two contrasting expression profiles. In this particular case, the differences were very subtle and limited to just one tissue (liver) responding to just one treatment condition (hypoxia at 17°C). Separation required application of Kmeans clustering to the genes contained within the landscape feature. This again indicates that the level of func-tional dissection possible by ExprAlign depends to a significant extent on the diversity of treatments for which array data is generated.

Conclusion
By leveraging large-scale microarray data, the ExprAlign approach offers a practical aid to gene identification, and for the discovery of novel isoforms or variants of known genes. It bears comparison with the more conventional sequence alignment techniques in that both depend upon an association of co-expression properties, or sequence alignment, respectively, with a set of known standards. Both are limited by the available data, and both have limitations. We show that ExprAlign identifies subtle differences in expression properties which may not be evident from sequence similarity, particularly if the sequence data is limited to the 800 bp provided by a typical single pass read or if clones fail to generate sequence in automated procedures. Finally, it reveals ESTs possessing the same BLAST identity but which have distinctive expression profiles.