Identification of a new gene regulatory circuit involving B cell receptor activated signaling using a combined analysis of experimental, clinical and global gene expression data

INTRODUCTION

Diffuse large B-cell lymphoma (DLBCL), the most common lymphoid malignancy in adults, is derived from germinal centre or post-germinal centre B cells [1]. DLBCLs are clinically and molecular heterogenous [1-4]. Long-term disease-free survival is now a reality for at least 50 percent of patients. However, approximately 30 percent of patients with DLBCL, mainly those with advanced disease and those who relapse, will not respond appropriately to current treatments [4, 5]. A better understanding of the specific pathways underlying the malignant phenotype in DLBCL is urgently required. This would allow not only the development of novel targeted therapeutic approaches for defined subgroups of DLBCL, but also the identification of new biomarkers [6].

Activation of B lymphocytes is associated with rapid proliferation, somatic hypermutation and targeted DNA double strand breaks in the immunoglobulin (Ig) locus. This process is under tight spatial and temporal control achieved by extracellular signals sensed via B cell receptor (BCR), CD40, Toll-like receptors (TLR), B-cell activating factor (BAFF) receptor, interleukin-21 (IL21) receptor and cell intrinsic mechanisms [1, 7-13]. Molecular signatures suggest that specific signaling networks and survival mechanisms exist in DLBCL, caused either by extracellular signals from the lymphoma microenvironment or related pathway mutations within the lymphoma cells, or both [1, 13].

Based on the cell of origin two main molecularly defined subgroups of DLBCLs can be described: activated B cell (ABC) - like and germinal centre B cell (GCB) – like DLBCL [14]. This cell-of-origin-approach is driven by class labels. In parallel, signatures such as stroma, host response, OxPhos and BCR/proliferation or Jak/STAT and IKK have been described by comprehensive consensus clustering, and identify the lymphoma microenvironment as a defining feature [15-19]. The level of c-Myc activity also allows the subdivision of DLBCL [20-23]; a high c-Myc index or a high number of Myc positive lymphoma cells is associated with shorter survival [20-23], but the subentities defined by Myc expression only partly overlap with the ABC/GCB-like signatures. Although these studies predict activity of different oncogenic pathways for DLBCL subtypes the functional consequences of differential gene expression are still poorly understood [15, 18, 19, 24-27],

Recently we introduced guided clustering as a strategy that integrates experimental, clinical and global gene expression data to investigate aggressive non-Hodgkin lymphoma (NHL) [28-30]. This approach does not utilize the class labels from the clinical data (e.g. disease types or clinical outcomes) to derive a molecular signature, but instead is based on the underlying biology for example the activity of an entire pathway or pathway networks established by cell perturbation experiments. For example, in our recent proof of principle analysis we were able to captured quantitatively a link between BCL-6 gene regulation and TLR signaling. The combined analysis of experimental and global gene expression data can allow the development of functional hypotheses [29, 31]. Therefore, transcriptional modules (groups of coherent expressed genes) can be described as potential major oncogenic hubs that are conserved between patients.

The present study aimed to extract functional information on oncogenic pathway activities in distinct DLBCL cases based on gene expression from primary tumours and multiple in vitro interventions. To do this the guided clustering-driven identification of transcriptional modules was extended and multiple single in vitro interventions were combined in silico to identify transcriptional modules conserved between patients but dominantly affected by for example BCR signaling but not TLR, IL21, CD40L or BAFF. This is to enable the assessment of pathway specific activities in primary lymphoma cases.

A new regulatory circuit in B cells was identified. A group of dominantly suppressed genes upon BCR activation is involved in an overall diminished capacity of the cells to enter mitosis, leading to defects in metaphase as well as increased chromosomal aberrations. In a subgroup of GCB-like DLBCL with low Myc activity this new regulatory circuit is dominant. This regulatory circuit is nearly absent in BL but active in c-Myc^high DLBCLs. Our data supports the view that BCR signaling is context dependent and capable not only of promoting cell survival and proliferation but also delaying cell cycle progression thereby potentially increasing chromosomal aberrations. It further underpins the notion that defined pathways stimulated by microenvironmental factors activating the BCR are involved in DLBCL development and that these pathways might be of therapeutic relevance. Our analysis shows how guided clustering lead to the discovery of biomarkers for cancer stratification.

RESULTS

A combined analysis of experimental and tumour derived global gene expression data identifies a set of genes specifically suppressed by BCR activation

Ligands activating pattern recognition receptors, BCR, CD40, BAFF-receptors and IL21 receptor are well known mediators of signalling in B cells and important components of the GC B cell reaction. Furthermore, it is well known that elements of the corresponding signalling pathways are mutated in DLBCL [1, 7-13]. Thus, the signalling pathways activated by these factors represent promising candidates for the identification of oncogenic pathway signatures in DLBCL via guided clustering. To answer these questions, as a model cell line, BL2 was chosen. The criteria for their selection were: absence or low pathway activity, a strong signal induction by stimuli, and measurable global gene expression changes suitable for bioinformatic analysis as we have previously described [32].

Microarray data sets obtained from human transformed germinal centre B cells (BL2) stimulated with CD40L, BAFF, IL21, αIgM F(ab)₂ fragments and lipopolysaccharide (LPS) were processed as described previously, combined, and analysed by guided clustering using large-scale gene expression data from 175 DLBCL patients [28, 32]. The patients were selected from the MMML-cohort and are representative of non-mBLs without chromosomal MYC translocations [30]. Guided clustering was performed in the following way: the guiding datasets were obtained from in vitro stimulated BL2 cells and only genes driven dominantly by one stimuli, but not the others, included. These data sets were integrated with gene expression profiles of primary lymphoma material. Ten different gene clusters were identified characterized by increased or suppressed gene expression in experiments and concordantly expressed in lymphoma patients: CD40.1, CD40.2, IL21.1, IL21.2, BAFF.1, BAFF.2, BCR.1, BCR.2, LPS.1 and LPS.2 (Figure 1A, Table I). The suffix “.1” denotes genes mainly suppressed and “.2” those genes mainly activated (Table I, Supplementary Table S1). These clusters most likely represent surrogates of pathway activity dominated by one of the stimuli. To delineate so far undescribed biological outcomes the following experiments were focused on αIgM driven suppression of gene expression.

Next the 137 DLBCL patient samples used for guided clustering were sorted according to an individual index value derived from BCR.1. This index reflects the extent to which BCR.1 genes are expressed within an individual sample. A high BCR.1 index indicates strong down-regulation of the BCR.1 genes. Finally, lymphoma samples were sorted according to their individual index (Figure 1B). This results in a spectrum of lymphoma samples ranging from low to high BCR.1 index scores. A substantial number of analysed patient samples is characterized by a corresponding strong BCR.1 signal.

Molecular functions, biological processes, cellular components and pathways of BCR.1 genes were identified by gene ontology (GO) based gene set enrichment analyses (Table II, Supplementary Table S2). Among the top 100 BCR.1 genes, the most significantly enriched gene sets corresponded to biological processes which included cell cycle checkpoints, microtubule-based processes, microtubule cytoskeleton organization, spindle organization, mitotic chromosome segregation or sister chromatid segregation as well as response to DNA damage stimuli, organelle organization or metabolic processes (Table II). In Figure 1C BCR.1 index genes associated with cell cycle regulatory functions are presented. Genes are grouped according to the corresponding cell cycle phases in which they are involved, demonstrating that more than 40% of cell cycle regulatory genes within this gene cluster belong to G2 and M cell cycle regulation.

Verification and validation of gene expression after αIgM treatment reveals the downregulation of BCR.1 index genes

Table III summarises the fold changes in expression of genes from the BCR.1 cluster involved in cell cycle regulation that were suppressed by αIgM treatment of BL2 cells.

Of the 288 genes comprising the BCR.1 index PLK1, AURKA, NEK2, BUB1B, CDC20, MPS1, BIRC5, HMMR, TACC3 and KIF14A were selected for validation and verification experiments by qRT-PCR in two independent cell lines [33-40]. These genes were selected because they are mitotic entry and progression regulating genes. BL2 and Ramos cells were stimulated with αIgM for 3hrs and gene expression changes analysed by qRT-PCR (Figure 2A). In line with our microarray data, qRT-PCR analyses show suppression of all selected BCR.1 genes by BCR activation. Furthermore, the expression of these genes was analysed over time following αIgM treatment (Figure 2B). Data were taken from a global gene expression analysis in BL2 cells. Cells were stimulated by αIgM and harvested at different time points starting from 30 min after adding αIgM until 480 min. RNA of three independent biological experiments was hybridized to Human ST1.0 microarrays. Most of the genes show a comparable time dependent suppression after αIgM treatment within the first 3-5hrs. Interestingly, some of the genes started to return to their basal expression level within 8hrs, suggesting a temporary effect on cell cycle regulation or corresponding negative feedback mechanism when cells are treated with αIgM (Figure 2B).

Next, we investigated whether BCR activation affects the same genes in lymphoma precursor cells in the same way as in cultured lymphoma cells. Therefore, the expression of cell cycle regulatory genes affected by αIgM treatment in BL2 cells was analysed in human tonsillar CD10⁺ B cells. Data were taken from a global gene expression analysis of CD10⁺ B cells from nine different donors (Table III, Supplementary Table S3). Although there are single BCR.1 genes that were not suppressed by BCR activation in human tonsillar CD10⁺ B cells, most of the analysed genes are suppressed by BCR stimulation as observed for BL2 cells. Thus, based on the common suppression of BCR.1 genes in cultured B lymphoma cells and in CD10⁺ human B cells, we conclude that this transcriptional regulation represents a common functional pathway of activated BCR signaling.

In addition, cell lines derived from DLBCL, U2932, HT, SUDHL5 and SUDHL6, were analysed for the presence of the BCR.1 index and whether it is affected by αIgM treatment. The BCR.1 index was detectable in HT and U2932 cells (Supplementary Figure S1A). Stimulation of DLBCL cell lines with αIgM led to a strong downregulation of BCR.1 genes (index increase) with the exception of SUDHL6 (Supplementary Figure S1B), indicating that down-regulation of BCR.1 index genes is also observed in DLBCL.

A prolongation of G2/M transition in αIgM stimulated human B cells

Since many genes of the BCR.1 gene cluster are involved in the regulation of mitotic processes, we next studied if BCR signaling had an inhibitory effect on mitosis. Cell cycle progression was monitored using flow cytometry to define the percentage of the diploid and tetraploid sets of chromosomes [41, 42]. In asynchronously growing BL2 and Ramos cells, αIgM treatment led to a delay in G2/M as revealed by an increase in the percentage of cells with a 4N DNA content (Figure 3A). To obtain a better functional insight into the effects of BCR activation on cell cycle regulation, Ramos cells were synchronized at the G1/S transition. After release from thymidine block, cells where left untreated or stimulated by αIgM. Cell cycle progression was monitored using flow cytometry at defined time points as shown in Figure 3B. Within 4hrs both untreated or αIgM treated cells entered the G2 phase of the cell cycle. While untreated cells progressed through mitosis within 5-6hrs, αIgM treated cells were delayed by approximately 1-2hrs. Figure 3C shows the percentage of the 2N and 4N DNA content. Analyses of the phosphorylation of Histone 3 as a mitotic marker (pHistoneH3) showed that pHistone-H3 peaks after 5hrs in control cells in contrast to αIgM treated cells, which peaks at 6.5hrs (Figure 3D). After 8hrs no phosphorylation remains. Interestingly, the phosphorylation of this mitotic marker is generally decreased in cells activated via the BCR. To quantify the effect, the mitotic index of thymidine synchronized BL cells was analysed by measuring the percentage of pMPM3 positive cells. In αIgM treated cells the number of pMPM3 positive cells was lower compared to untreated cells (Figure 3E) further supporting the view that the activation of the BCR signaling pathway leads to a delay in G2 and/or lag in progression through the mitotic phase of the cell cycle. This is in agreement with the observed reduced expression of cell cycle regulators from the BCR.1 gene cluster. Immunoblot analysis for Aurora-A, Aurora-B, Tpx2 or Mad2 and the phosphorylated forms of Aurora-A and -B revealed corresponding changes in protein levels and phosphorylation patterns (Figure 3F). In line with a delay in mitotic progression, activation (phosphorylation) of the key mitotic kinases, Aurora-A and Aurora-B, was prolonged in BCR-activated BL cells (Figure 3F). Untreated cells enter M phase 4-6 hrs after release of the cell cycle block. After 8 hrs most of the cells have left mitosis as shown by a peak of TPX2 levels at 6hrs and the sharp decline of TPX2 after 8hrs (Figure 3F). However, cells that have been activated via BCR enter mitosis at the latest after 8hrs and need 10hrs to complete mitosis. The prolongation of the G2 phase as well as the observed delay of entry into mitosis could be an indication that BCR activation suppresses mitotic processes.

Therefore, defects in mitosis were investigated and quantified. Synchronized BL cells were released from cell cycle block as described above in the presence or absence of αIgM. 2hrs before harvesting, cells were treated with the proteasome inhibitor MG132 to arrest cells in metaphase. In Ramos cells treated with MG132 alone, metaphase chromosomes were found to be properly aligned at the equatorial plane. Abnormal metaphases were increased in αIgM treated cells (Figure 3H), suggesting that BCR signaling causes chromosome missegregation and subsequent accumulation of karyotype abberations.

To test this, a chronic activation of BCR was mimicked by culturing Ramos cells over 3 weeks with daily addition of αIgM F(ab)₂ fragments. In Figure 4 A/B cell viability and cell doubling data are presented from these experiments. A decrease of cell viability within the first 10 days was accompanied by a very low cell doubling. This fits well with earlier observations of increased cell death in BL cells after in vitro BCR activation monitored as cell cycle arrest within G1 [43]. However, after two weeks both cell proliferation and viability were increased. Analysis of the karyotypes of chronically treated Ramos cells revealed additional chromosomal aberrations including translocations (t(15;21)(q14?;q22)) as well as additions (add(16)(q24)) (Figure 4). These structural abnormalities were not observed in untreated BL cells. One explanation for the induction of additional chromosomal abberations could be that chronic exposure of B cells to antigens and thus BCR activation can lead to subsequent cell cycle changes predisposing cells for the acquisition of subsequent chromosomal aberrations.

Gene suppression is negatively associated with c-Myc activity after B cell receptor activation

It was shown recently that in BL cell lines αIgM treatment is associated with a PI3K dependent decrease in MYC expression [32]. Therefore, the involvement of c-Myc in BCR.1 gene regulation, particularly on cell cycle regulators and on cell cycle deregulation induced by αIgM was evaluated.

First, a patient-based comparison between the BCR.1 index with our recently established Myc-index was performed [20]. The correlation coefficients of these two indices were calculated in gene expression profiles of 219 B-NHL including molecularly defined BLs (mBLs), non-mBLs and intermediate lymphoma (molecularly unclassified B-NHLs) [30]. The correlation coefficient is -0.76 (Figure 5A). Thus, the BCR.1 index inversely correlates with the c-Myc index. Figure 5A shows both index values per sample. Therefore, mBL cases are characterized by a high c-Myc but a low BCR.1 index, whereas non-mBLs and intermediate lymphoma are more heterogeneous with a lower c-Myc and higher BCR.1 index. Lymphomas with a high BCR.1 signal can be characterized by a very low c-Myc activity and therefore can be grouped as BCR.1^high/Myc^low NHL. The majority of non-mBLs and intermediate lymphoma are characterized by corresponding higher BCR.1 but lower c-Myc activity relative to mBLs. Comparison of the BCR.1 index with BCL6 regulated genes was performed in patient samples (Supplementary Figure S3) [28, 44]. Although a good separation of mBLs from intermediate and non-mBLs was observed, the index correlation (-0.1) was not as strong as for c-Myc. Therefore, the c-Myc protein levels were scored on tissue microarrays (TMA) of 99 lymphoma cases taken from the MMML-cohort described above [30]. In parallel, the BCR.1 index for these cases based on their gene expression profiles was calculated. This analysis revealed that cases with higher numbers of c-Myc positive lymphoma cells were significantly more likely to express BCR.1 genes (low BCR.1 index) (p=0.0012) (Figure 5B). This finding underscores the negative correlation of c-Myc and BCR.1 indices and supports the definition of a BCR.1^high Myc^low NHL subgroup (mainly DLBCL).

To assess whether the cell cycle delay after BCR activation might be the result of the BCR induced inhibition of c-Myc, we investigated whether the inhibition of c-Myc alone can induce G2-phase prolongation of thymidine synchronized Ramos cells. Double thymidine block synchronized Ramos cells were released from cell cycle block and treated with the chemical Myc inhibitor 10058-F4, αIgM, or both. Within 4hrs all cells entered the G2 phase of the cell cycle. While untreated cells entered the G1 phase of the cell cycle within 5-6hrs αIgM, 10058-F4 and αIgM/10058-F4 treated cells were delayed (Figure 5C). The prolonged G2 phase in 10058-F4 treated cells was comparable to the effect of αIgM treatment alone and thus supports the inverse correlation of BCR.1 and the Myc index shown in Figure 5A/B. The strongest delay was observed by the combination of BCR activation and c-Myc inhibition. Therefore, αIgM treatment together with 10058-F4 shows an additive effect on G2 phase prolongation. This suggests that the BCR effect on the cell cycle regulation is not solely the result of inhibiting MYC expression. In Figure 5D the percentage of the diploid and tetraploid sets of chromosomes is shown.

The role of aberrant c-Myc activity on cell cycle and BCR.1 gene regulation was analysed by qRT-PCR as shown for BUB1B and other BCR.1 genes in BL2 and Ramos cells treated with 10058-F4, αIgM or both (Figure 5E, Supplementary Figure S4). The expression of some BCR.1 genes is dependent on aberrant c-Myc activity, but BCR activation and c-Myc inhibition act in an additive manner.

In addition we tested by chromatin-immunoprecipitation (ChIP) whether c-Myc binding to a known E-box containing region in the first intron of the BUB1B gene is altered upon BCR activation [45, 46]. The gene locus of the BUB1B gene was enriched by ChIP of c-Myc in unstimulated cells (Figure 5F). Importantly, following BCR activation this enrichment was abolished. As positive control a ChIP of acetylated histone H3 was used as a marker for active transcription. Therefore, a critical cell cycle regulator is suppressed in response to suppression of MYC gene expression by αIgM. However, the observed additive effect of BCR activation together with 10058-F4 treatment also indicates that additional pathways are involved through BCR activation which could mediate the suppression of the BCR.1 genes.

The BCR.1 gene module characterizes individual NHL and is more active within GCB-like DLBCL

To further study the functional relevance of the BCR.1 gene module in specific subgroups of B-NHL and specifically in DLBCL, we compared by the gene expression profiles of 389 primary lymphoma samples [30, 47]. Figure 6A shows that mBLs (colored in red) are characterized by a low or absent BCR.1 index. This was also observed for samples characterized by MYC-aberrations. About 20% of analysed non-mBLs show a strong coherently low expression of BCR.1 genes and therefore are characterized by a strong corresponding BCR signal. This is in good agreement with the observation presented in Figure 1B within the „training cohort“ and the observed inverse correlation of the BCR.1 and c-Myc indices (Figure 5A). A comparison of ABC- and GCB-like DLBCLs in this analysis is difficult as mBLs are labelled also as GCB-like lymphomas. Therefore, an additional analysis was performed (Figure 6B). Comparing 281 non-mBL DLBCL samples, we observed a significantly higher BCR.1 index in GCB-DLBCL (n=159) compared with ABC-like DLBCL (n=122) (p=9.35e−06). This supports the view that more GCB-like DLCBL are characterized by a specific coexpression of BCR.1 genes than are ABC-like DLBCL.

DISCUSSION

Gene expression based molecular tumour signatures have been identified in aggressive NHL. In DLBCL specific oncogenic signalling networks have been deduced from these signatures. These networks are driven by mutations or the microenvironment acting onto the lymphoma cells [1, 13, 48, 49]. However, the functional effects revealed by these gene expression signatures, are poorly understood. A better understanding could improve of lymphoma treatments. Therefore, the present study aimed to gain further functional insights into oncogenic pathways active in aggressive NHL. The guided clustering approach was extended to identify transcriptional modules that are conserved between patients and in vitro stimulated cells but dominantly activated or inhibited by for example BCR signaling but not by other well known mediators of signalling in B cells. We have provided a comprehensive genome-wide resource for the functional exploration of molecular and cellular regulation in B cells. Ten new transcriptional modules were identified (Figure 1).

Of the ten transcriptional modules we identified, the BCR.1 module comprises genes, which are suppressed by BCR activation. It was shown that these genes are downregulated by BCR activation in a number of different cell lines including BL2, Ramos, U2932, HT and SUDHL5. We conclude that in patients BCR or related pathways are involved in the described way of gene coregulation. Furthermore, the BCR.1 module is strong in a substantial number of GCB-like DLBCL but not in mBL and supports the hypothesis, that in some DLBCL BCR signaling is involved in the downregulation of cell cycle regulators. Because in ABC-like DLBCLs the BCR.1 signal is attenuated compared to GCB-like DLBCLs (Figure 6B), it would be interesting to determine if the coactivity of classical Jak/STAT and NF-κB signaling characteristic of ABC-like DLBCL can enhance the BCR.1 gene expression independent of c-Myc [19, 50].

Remarkably, BCR.1 genes are enriched for those involved in cell cycle regulation (Figure 2 and Table III). The functional outcome of this gene module is a delayed entry to and progression of mitosis most probably as a result of reduced expression of the identified cell cycle regulators (Figure 3). After continuous BCR activation some B cells gain additional genetic aberrations (Figure 4) as a direct or indirect consequence of the observed changes in G2/M cell cycle phase delays. By large-scale follow-up analyses of pairs of primary and relapsed aggressive NHL it can now be tested which role BCR.1 gene expression is playing to affect increases of chromosomal aberrations or disease progression. However, this is also an important observation as the genes suppressed by BCR activation are also found to be downregulated in isolated human tonsillar CD10⁺ B cells. Therefore, it is proposed that these genes are functionally relevant in both transformed and in reactive B cells.

There is growing evidence of a role for BCR signaling in the pathogenesis of different subtypes of NHL. In BL a tonic BCR signal is observed whereas in DLBCL both tonic and chronic active BCR signaling seem to be important [26, 27, 51-53]. As well as supporting B cell proliferation, survival or differentiation, BCR activation can also induce G1 cell cycle arrest, apoptosis and increased cell death in different B cells [54-59]. Our observations suggest that these suicide pathways can be, at least partially, attributed to a delay in G2/M with associated mitotic defects potentially leading to subsequent chromosomal aberrations. Our data also provide a plausible mechanism to explain the link between inflammation and the development of lymphomagenic mutations; especially in DLBCL with a high BCR.1, chronic inflammation and/or chronic antigen driven activation of oncogenic pathways with no or low c-Myc activity. The inverse correlation of the Myc-index and the BCR.1 index suggests that in lymphomas with a high MYC expression corresponding BCR.1 genes are highly expressed (low BCR.1 index, high gene expression) and thus most likely indicates an absence of a strong antigen stimulation.

The newly identified oncogenic network present in a substantial number of GCB-like DLBCL will likely define more rational treatment targets in this heterogeneous disease. Whether this has implications for stratifying NHL patients for molecular-based prognostication and for targeted therapy has to be investigated in new clinical trials using R-CHOP and/or currently introduced targeted drugs for Btk, PI3K and NF-kB pathways or BET-domain containing protein inhibitors [25, 48, 51, 60].

MATERIALS and METHODS

Cell culture

BL2, Ramos, SUDHL5, SUDHL6, U2932 and HT (all obtained from the DSMZ Braunschweig, Germany) cell cultivation as well as purification and sorting of primary human tonsillar B cells was performed as previously described [20, 61] [32]. A detailed protocol can be found in the “Extended View Methods Section”. For stimulation studies, cells were cultured as described previously [32]. Cells chronically activated by αIgM were grown for up to 28 days and αIgM was added every 24hrs. Cells were harvested using corresponding inhibitors of phosphatases and proteases and RNA isolated using the RNeasy Plus Mini Kit (Qiagen). Cell doubling and viability was determined by counting cells in the presence of trypan blue.

Synchronisation of BL cells was performed using thymidine treatment as recently described [62]. Cell cycle analysis was performed on the basis of analysing the DNA content in the nuclei of the cells by propidium iodide staining followed by flow cytometry. For a detailed description see the “Expanded View Methods Section”.

Gene expression analysis

qRT-PCR was performed using SYBR green. ΔCt values were normalised to ß2m and abl expression and ΔΔCt values calculated. Oligonucleotides used are summarized in Supplementary Table S4. For whole genome microarrays Human Genome U133A 2.0 plus Arrays (Affymetrix) or Human ST1.0 Arrays (Affymetrix) as indicated was performed according to manufacturer’s recommendations by the TAL (UMG, Germany) For further details of gene expression analyses see the “Expanded View Methods Section”.

Cytogenetic analysis

Metaphases were analysed using chromosome banding analysis and 24-color FISH. Chromosome preparation was performed using standard cytogenetic protocols and modified chromosome banding technology (GAG// Giemsa bands of acetic Saline Giemsa). M-FISH (24 color FISH) was done according to the manufacturer’s protocol (MetaSystems). Karyotypes were classified according to the International System for Human Cytogenetic Nomenclature (ISCN) [63].

Chromatin immunoprecipitation

BL2 cells were stimulated with 1.3µg/ml αIgM F(ab)₂ fragments for 3hrs or left untreated as control. Cells were sedimented and resuspended in PBS containing 1.42% formaldehyde and incubated for 15min at room temperature. Chromatin-Immunoprecipitation was conducted as described previously using 2µg anti c-Myc (clone N-262, Santacruz), 2µg anti acH3 (clone 06-599, Millipore), 2µg αIgG control (ab46540, abcam) antibodies. A detailed protocol is presented in the “supplemental Material and Methods section”.

Case selection and classification for immunohistochemistry

127 B-cell lymphoma samples were obtained from the files of the lymph node registry Kiel with ethical approval and classified according to the WHO classification using standard histological, immunohistochemical and molecular criteria [64]. In addition, 20 specimens obtained within the framework of the network project “Molecular Mechanisms in malignant Lymphoma” (MMML) were analysed. The latter cohort was diagnosed and classified based on the previously reported molecular signatures [30]. All lymphomas were analysed using tissue-micro-arrays (TMA) containing 2 cores of 1 mm or 0.6 mm (for the MMML specimen) for each case (ethical approval Ref. Nr. 1/1/05).

Bioinformatic analyses and analysed datasets

The identification of pathway activation clusters was performed using the Guided Clustering algorithm [28]. Guided Clustering identifies genes specific for one of the five stromal stimuli (BCR, CD40L, BAFF, IL-21, LPS). Within this framework the stimulated BL2 cells act as guiding data whereas the cohort of Myc-negative DLBCL represent the guided dataset [30]. The definition of the required binary vector depends on the respective stimulation with the numerical value 1 for the stimulus and 0 for the remaining stimuli as well as for the controls. For each stimulation two clusters were extracted, showing an opposite regulation. The respective genes in each cluster may represent a potential surrogate marker for pathway activity in lymphoma samples. To determine the extent of gene cluster activity in lymphoma patients, one index was calculated by gene cluster and per lymphoma sample. Expression values from those genes assembling the individual cluster were used to calculate a single representative value (the index) for each sample by fitting a standard additive model with independent gene and sample effects using Tukey’s median polish procedure as described [65]. The index values of genes inhibited by a stimulus (.1 modules (eg BCR.1)) were multiplied by minus one to enable the interpretation of the indices as absence or presence of stimulation. Hence, down-regulation of these genes leads to a high index value (stimulation is present) whereas a low index value indicates the absence of any regulatory effect due to stimulation. Gene set enrichment analysis (GSEA) of ranked gene list was performed using the Java implementation of GSEA obtained from http://www.broadinstitute.org/gsea/. GSEA was conducted in the mode for pre-ranked gene lists on the C2 set of curated gene signatures from the Molecular Signature Database (MSigDB). Genes are attributed to one ore more biological processes, molecular function and cellular location with respect to the Gene Ontology as well as pathways in the KEGG database [66, 67]. Hypergeometric testing was performed to test for non-random overrepresentation of module genes in each of these terms based on the hypergeometric distribution model. This GO/KEGG analyses was implemented in customized inhouse scripts.

The primary data are available from GEO (http://www.ncbi.nlm.nih.gov/geo/) - raw data for gene expression changes for LPS, CD40L, BAFF, IL-21 and αIgM stimulated BL2 cells are available under accession no. GSE42660 and GSE29700 [28, 32], data for CD10⁺ human tonsillar B cells (GSE71724), pathway inhibition (GSE68761) and gene expression kinetics (GSE71721) are summarized as SuperSeries record GSE71725.

ACKNOWLEDGEMENTS

We would like to thank S. Hengst for excellent technical support in some of the presented experiments. This work was supported by grants from the Deutsche Krebshilfe (Network “Molecular mechanisms of malignant lymphoma”), BMBF network on Medical Systems biology “HaematoSys” (BMBF-FKZ 0315452H) (D.K., R.S., W.K.), BMBF eBio “MMML-MycSys” (BMBF-FKZ 0316166E) (D.K., R.S., L.T., W.K.) and the Deutsche Forschungsgemeinschaft (GRK1034), Bloodwise (formerly Leukaemia and Lymphoma Research; M.V., P.G.M.). We thank all other members of the consortia “HaematoSys” and “MMML-MycSys” for their helpful discussions preparing this manuscript.

CONFLICTS OF INTEREST

The authors declare no competing financial interests.

Author contributions

Contribution:A.Sc., A.St., N.W., F.v.B., M.F., E.H. did most of the experiments with K.S. and M.V. contributing to specific experiments as cytogenetic analysis, FISH and preparation/characterization of CD10⁺ tonsillar B cells and W.K., S.Z. performing immunohistochemistry; K.M., C.K. and M.E. performed bioinformatics analysis; R.S. designed experiments, performed bioinformatics analysis and interpreted data, wrote the manuscript; H.B., P.G.M., L.T. designed experiments, final approval of the manuscript; A.Sc., D.K. designed the research, analysed, interpreted data, and wrote the manuscript.

REFERENCES

1. Shaffer AL, 3rd, Young RM and Staudt LM. Pathogenesis of human B cell lymphomas. Annu Rev Immunol. 2012; 30:565-610.

2. Klein U and Dalla-Favera R. Germinal centres: role in B-cell physiology and malignancy. Nat Rev Immunol. 2008; 8:22-33.

3. Coiffier B. State-of-the-art therapeutics: diffuse large B-cell lymphoma. J Clin Oncol. 2005; 23:6387-6393.

4. Martelli M, Ferreri AJ, Agostinelli C, Di Rocco A, Pfreundschuh M and Pileri SA. Diffuse large B-cell lymphoma. Critical reviews in oncology/hematology. 2013; 87:146-171.

5. Kubuschok B, Held G and Pfreundschuh M. Management of diffuse large B-cell lymphoma (DLBCL). Cancer treatment and research. 2015; 165:271-288.

6. Simon R. Genomic biomarkers in predictive medicine: an interim analysis. EMBO molecular medicine. 2011; 3:429-435.

7. Compagno M, Lim WK, Grunn A, Nandula SV, Brahmachary M, Shen Q, Bertoni F, Ponzoni M, Scandurra M, Califano A, Bhagat G, Chadburn A, Dalla-Favera R and Pasqualucci L. Mutations of multiple genes cause deregulation of NF-kappaB in diffuse large B-cell lymphoma. Nature. 2009; 459:717-721.

8. Kim SJ, Lee SJ, Choi IY, Park Y, Choi CW, Kim IS, Yu W, Hwang HS and Kim BS. Serum BAFF predicts prognosis better than APRIL in diffuse large B-cell lymphoma patients treated with rituximab plus CHOP chemotherapy. Eur J Haematol. 2008; 81:177-184.

9. Ngo VN, Young RM, Schmitz R, Jhavar S, Xiao W, Lim KH, Kohlhammer H, Xu W, Yang Y, Zhao H, Shaffer AL, Romesser P, Wright G, Powell J, Rosenwald A, Muller-Hermelink HK, et al. Oncogenically active MYD88 mutations in human lymphoma. Nature. 2011; 470:115-119.

10. Zhu X, Hart R, Chang MS, Kim JW, Lee SY, Cao YA, Mock D, Ke E, Saunders B, Alexander A, Grossoehme J, Lin KM, Yan Z, Hsueh R, Lee J, Scheuermann RH, et al. Analysis of the major patterns of B cell gene expression changes in response to short-term stimulation with 33 single ligands. J Immunol. 2004; 173:7141-7149.

11. Denepoux S, Razanajaona D, Blanchard D, Meffre G, Capra JD, Banchereau J and Lebecque S. Induction of somatic mutation in a human B cell line in vitro. Immunity. 1997; 6:35-46.

12. Henriquez NV, Floettmann E, Salmon M, Rowe M and Rickinson AB. Differential responses to CD40 ligation among Burkitt lymphoma lines that are uniformly responsive to Epstein-Barr virus latent membrane protein 1. J Immunol. 1999; 162:3298-3307.

13. Pasqualucci L. The genetic basis of diffuse large B-cell lymphoma. Curr Opin Hematol. 2013; 20:336-344.

14. Alizadeh AA, Eisen MB, Davis RE, Ma C, Lossos IS, Rosenwald A, Boldrick JC, Sabet H, Tran T, Yu X, Powell JI, Yang L, Marti GE, Moore T, Hudson J, Jr., Lu L, et al. Distinct types of diffuse large B-cell lymphoma identified by gene expression profiling. Nature. 2000; 403:503-511.

15. Monti S, Savage KJ, Kutok JL, Feuerhake F, Kurtin P, Mihm M, Wu B, Pasqualucci L, Neuberg D, Aguiar RC, Dal Cin P, Ladd C, Pinkus GS, Salles G, Harris NL, Dalla-Favera R, et al. Molecular profiling of diffuse large B-cell lymphoma identifies robust subtypes including one characterized by host inflammatory response. Blood. 2005; 105:1851-1861.

16. Polo JM, Juszczynski P, Monti S, Cerchietti L, Ye K, Greally JM, Shipp M and Melnick A. Transcriptional signature with differential expression of BCL6 target genes accurately identifies BCL6-dependent diffuse large B cell lymphomas. Proc Natl Acad Sci U S A. 2007; 104:3207-3212.

17. Lenz G, Wright G, Dave SS, Xiao W, Powell J, Zhao H, Xu W, Tan B, Goldschmidt N, Iqbal J, Vose J, Bast M, Fu K, Weisenburger DD, Greiner TC, Armitage JO, et al. Stromal gene signatures in large-B-cell lymphomas. N Engl J Med. 2008; 359:2313-2323.

18. Kloo B, Nagel D, Pfeifer M, Grau M, Duwel M, Vincendeau M, Dorken B, Lenz P, Lenz G and Krappmann D. Critical role of PI3K signaling for NF-kappaB-dependent survival in a subset of activated B-cell-like diffuse large B-cell lymphoma cells. Proc Natl Acad Sci U S A. 2011; 108:272-277.

19. Lam LT, Wright G, Davis RE, Lenz G, Farinha P, Dang L, Chan JW, Rosenwald A, Gascoyne RD and Staudt LM. Cooperative signaling through the signal transducer and activator of transcription 3 and nuclear factor-{kappa}B pathways in subtypes of diffuse large B-cell lymphoma. Blood. 2008; 111:3701-3713.

20. Schrader A, Bentink S, Spang R, Lenze D, Hummel M, Kuo M, Arrand JR, Murray PG, Trumper L, Kube D and Vockerodt M. High myc activity is an independent negative prognostic factor for diffuse large B cell lymphomas. Int J Cancer. 2012; 131:E348-361.

21. Walther N, Ulrich A, Vockerodt M, von Bonin F, Klapper W, Meyer K, Eberth S, Pukrop T, Spang R, Trumper L and Kube D. Aberrant lymphocyte enhancer-binding factor 1 expression is characteristic for sporadic Burkitt’s lymphoma. Am J Pathol. 2013; 182:1092-1098.

22. Horn H, Ziepert M, Becher C, Barth TF, Bernd HW, Feller AC, Klapper W, Hummel M, Stein H, Hansmann ML, Schmelter C, Moller P, Cogliatti S, Pfreundschuh M, Schmitz N, Trumper L, et al. MYC status in concert with BCL2 and BCL6 expression predicts outcome in diffuse large B-cell lymphoma. Blood. 2013; 121:2253-2263.

23. Dave SS, Fu K, Wright GW, Lam LT, Kluin P, Boerma EJ, Greiner TC, Weisenburger DD, Rosenwald A, Ott G, Muller-Hermelink HK, Gascoyne RD, Delabie J, Rimsza LM, Braziel RM, Grogan TM, et al. Molecular diagnosis of Burkitt’s lymphoma. N Engl J Med. 2006; 354:2431-2442.

24. Caro P, Kishan AU, Norberg E, Stanley IA, Chapuy B, Ficarro SB, Polak K, Tondera D, Gounarides J, Yin H, Zhou F, Green MR, Chen L, Monti S, Marto JA, Shipp MA, et al. Metabolic signatures uncover distinct targets in molecular subsets of diffuse large B cell lymphoma. Cancer Cell. 2012; 22:547-560.

25. Chapuy B, McKeown MR, Lin CY, Monti S, Roemer MG, Qi J, Rahl PB, Sun HH, Yeda KT, Doench JG, Reichert E, Kung AL, Rodig SJ, Young RA, Shipp MA and Bradner JE. Discovery and characterization of super-enhancer-associated dependencies in diffuse large B cell lymphoma. Cancer Cell. 2013; 24:777-790.

26. Davis RE, Ngo VN, Lenz G, Tolar P, Young RM, Romesser PB, Kohlhammer H, Lamy L, Zhao H, Yang Y, Xu W, Shaffer AL, Wright G, Xiao W, Powell J, Jiang J-K, et al. Chronic active B-cell-receptor signalling in diffuse large B-cell lymphoma. Nature. 2010; 463:88-92.

27. Schmitz R, Young RM, Ceribelli M, Jhavar S, Xiao W, Zhang M, Wright G, Shaffer AL, Hodson DJ, Buras E, Liu X, Powell J, Yang Y, Xu W, Zhao H, Kohlhammer H, et al. Burkitt lymphoma pathogenesis and therapeutic targets from structural and functional genomics. Nature. 2012; 490:116-120.

28. Maneck M, Schrader A, Kube D and Spang R. Genomic data integration using guided clustering. Bioinformatics. 2011; 27:2231-2238.

29. Bentink S, Wessendorf S, Schwaenen C, Rosolowski M, Klapper W, Rosenwald A, Ott G, Banham AH, Berger H, Feller AC, Hansmann ML, Hasenclever D, Hummel M, Lenze D, Moller P, Stuerzenhofecker B, et al. Pathway activation patterns in diffuse large B-cell lymphomas. Leukemia. 2008; 22:1746-1754.

30. Hummel M, Bentink S, Berger H, Klapper W, Wessendorf S, Barth TF, Bernd HW, Cogliatti SB, Dierlamm J, Feller AC, Hansmann ML, Haralambieva E, Harder L, Hasenclever D, Kuhn M, Lenze D, et al. A biologic definition of Burkitt’s lymphoma from transcriptional and genomic profiling. N Engl J Med. 2006; 354:2419-2430.

31. Bild AH, Yao G, Chang JT, Wang Q, Potti A, Chasse D, Joshi MB, Harpole D, Lancaster JM, Berchuck A, Olson JA, Jr., Marks JR, Dressman HK, West M and Nevins JR. Oncogenic pathway signatures in human cancers as a guide to targeted therapies. Nature. 2006; 439:353-357.

32. Schrader A, Meyer K, von Bonin F, Vockerodt M, Walther N, Hand E, Ulrich A, Matulewicz K, Lenze D, Hummel M, Kieser A, Engelke M, Trumper L and Kube D. Global gene expression changes of in vitro stimulated human transformed germinal centre B cells as surrogate for oncogenic pathway activation in individual aggressive B cell lymphomas. Cell Commun Signal. 2012; 10:43-64.

33. Gadde S and Heald R. Mechanisms and molecules of the mitotic spindle. Curr Biol. 2004; 14:R797-805.

34. Fry AM, O’Regan L, Sabir SR and Bayliss R. Cell cycle regulation by the NEK family of protein kinases. J Cell Sci. 2012; 125:4423-4433.

35. Lioutas A and Vernos I. Aurora A: working from dawn to dusk in mitosis. Cell Cycle. 2014; 13:499-500.

36. Nigg EA. Mitotic kinases as regulators of cell division and its checkpoints. Nat Rev Mol Cell Biol. 2001; 2:21-32.

37. Gruneberg U, Neef R, Li X, Chan EH, Chalamalasetty RB, Nigg EA and Barr FA. KIF14 and citron kinase act together to promote efficient cytokinesis. J Cell Biol. 2006; 172:363-372.

38. von Schubert C and Nigg EA. Polo-like kinases. Curr Biol. 2013; 23:R225-227.

39. Mita AC, Mita MM, Nawrocki ST and Giles FJ. Survivin: key regulator of mitosis and apoptosis and novel target for cancer therapeutics. Clin Cancer Res. 2008; 14:5000-5005.

40. Kinoshita K, Noetzel TL, Pelletier L, Mechtler K, Drechsel DN, Schwager A, Lee M, Raff JW and Hyman AA. Aurora A phosphorylation of TACC3/maskin is required for centrosome-dependent microtubule assembly in mitosis. J Cell Biol. 2005; 170:1047-1055.

41. Holtick U, Vockerodt M, Pinkert D, Schoof N, Sturzenhofecker B, Kussebi N, Lauber K, Wesselborg S, Loffler D, Horn F, Trumper L and Kube D. STAT3 is essential for Hodgkin lymphoma cell proliferation and is a target of tyrphostin AG17 which confers sensitization for apoptosis. Leukemia. 2005; 19:936-944.

42. Stolz A, Ertych N, Kienitz A, Vogel C, Schneider V, Fritz B, Jacob R, Dittmar G, Weichert W, Petersen I and Bastians H. The CHK2-BRCA1 tumour suppressor pathway ensures chromosomal stability in human somatic cells. Nat Cell Biol. 2010; 12:492-499.

43. Mackus WJ, Lens SM, Medema RH, Kwakkenbos MJ, Evers LM, Oers MH, Lier RA and Eldering E. Prevention of B cell antigen receptor-induced apoptosis by ligation of CD40 occurs downstream of cell cycle regulation. Int Immunol. 2002; 14:973-982.

44. Ci W, Polo JM, Cerchietti L, Shaknovich R, Wang L, Yang SN, Ye K, Farinha P, Horsman DE, Gascoyne RD, Elemento O and Melnick A. The BCL6 transcriptional program features repression of multiple oncogenes in primary B cells and is deregulated in DLBCL. Blood. 2009; 113:5536-5548.

45. Menssen A, Epanchintsev A, Lodygin D, Rezaei N, Jung P, Verdoodt B, Diebold J and Hermeking H. c-MYC delays prometaphase by direct transactivation of MAD2 and BubR1: identification of mechanisms underlying c-MYC-induced DNA damage and chromosomal instability. Cell Cycle. 2007; 6:339-352.

46. Seitz V, Butzhammer P, Hirsch B, Hecht J, Gutgemann I, Ehlers A, Lenze D, Oker E, Sommerfeld A, von der Wall E, Konig C, Zinser C, Spang R and Hummel M. Deep sequencing of MYC DNA-binding sites in Burkitt lymphoma. PLoS One. 2011; 6:e26837.

47. Klapper W, Stoecklein H, Zeynalova S, Ott G, Kosari F, Rosenwald a, Loeffler M, Trümper L, Pfreundschuh M and Siebert R. Structural aberrations affecting the MYC locus indicate a poor prognosis independent of clinical risk factors in diffuse large B-cell lymphomas treated within randomized trials of the German High-Grade Non-Hodgkin’s Lymphoma Study Group (DSHNHL). Leukemia. 2008; 22:2226-2229.

48. Emadali A, Rousseaux S, Bruder-Costa J, Rome C, Duley S, Hamaidia S, Betton P, Debernardi A, Leroux D, Bernay B, Kieffer-Jaquinod S, Combes F, Ferri E, McKenna CE, Petosa C, Bruley C, et al. Identification of a novel BET bromodomain inhibitor-sensitive, gene regulatory circuit that controls Rituximab response and tumour growth in aggressive lymphoid cancers. EMBO molecular medicine. 2013; 5:1180-1195.

49. Scott DW and Gascoyne RD. The tumour microenvironment in B cell lymphomas. Nat Rev Cancer. 2014; 14:517-534.

50. Ding BB, Yu JJ, Yu RY, Mendez LM, Shaknovich R, Zhang Y, Cattoretti G and Ye BH. Constitutively activated STAT3 promotes cell proliferation and survival in the activated B-cell subtype of diffuse large B-cell lymphomas. Blood. 2008; 111:1515-1523.

51. Roschewski M, Staudt LM and Wilson WH. Diffuse large B-cell lymphoma-treatment approaches in the molecular era. Nature reviews Clinical oncology. 2014; 11:12-23.

52. Chen L, Monti S, Juszczynski P, Daley J, Chen W, Witzig TE, Habermann TM, Kutok JL and Shipp MA. SYK-dependent tonic B-cell receptor signaling is a rational treatment target in diffuse large B-cell lymphoma. Blood. 2008; 111:2230-2237.

53. Wilson WH, Young RM, Schmitz R, Yang Y, Pittaluga S, Wright G, Lih CJ, Williams PM, Shaffer AL, Gerecitano J, de Vos S, Goy A, Kenkre VP, Barr PM, Blum KA, Shustov A, et al. Targeting B cell receptor signaling with ibrutinib in diffuse large B cell lymphoma. Nat Med. 2015; 21:922-926.

54. Lin CK, Zou HY, Kaptein JS, Yen CF, Kalunta CI, Nguyen TT, Park E and Lad PM. Anti-IgM-induced growth inhibition and apoptosis are independent of ornithine decarboxylase in Ramos cells. Exp Cell Res. 1997; 237:231-241.

55. Eeva J, Postila V, Matto M, Nuutinen U, Ropponen A, Eray M and Pelkonen J. Kinetics and signaling requirements of CD40-mediated protection from B cell receptor-induced apoptosis. Eur J Immunol. 2003; 33:2783-2791.

56. Mimori K, Kiyokawa N, Taguchi T, Suzuki T, Sekino T, Nakajima H, Saito M, Katagiri YU, Isoyama K, Yamada K, Matsuo Y and Fujimoto J. Costimulatory signals distinctively affect CD20- and B-cell-antigen-receptor-mediated apoptosis in Burkitt’s lymphoma/leukemia cells. Leukemia. 2003; 17:1164-1174.

57. Chaouchi N, Vazquez A, Galanaud P and Leprince C. B cell antigen receptor-mediated apoptosis. Importance of accessory molecules CD19 and CD22, and of surface IgM cross-linking. J Immunol. 1995; 154:3096-3104.

58. Carey GB, Donjerkovic D, Mueller CM, Liu S, Hinshaw JA, Tonnetti L, Davidson W and Scott DW. B-cell receptor and Fas-mediated signals for life and death. Immunol Rev. 2000; 176:105-115.

59. Yamashita Y, Miyake K, Miura Y, Kaneko Y, Yagita H, Suda T, Nagata S, Nomura J, Sakaguchi N and Kimoto M. Activation mediated by RP105 but not CD40 makes normal B cells susceptible to anti-IgM-induced apoptosis: a role for Fc receptor coligation. J Exp Med. 1996; 184:113-120.

60. Nagel D, Vincendeau M, Eitelhuber AC and Krappmann D. Mechanisms and consequences of constitutive NF-kappaB activation in B-cell lymphoid malignancies. Oncogene. 2014 Dec 11; 33:5655-65.

61. Vockerodt M, Morgan SL, Kuo M, Wei W, Chukwuma MB, Arrand JR, Kube D, Gordon J, Young LS, Woodman CB and Murray PG. The Epstein-Barr virus oncoprotein, latent membrane protein-1, reprograms germinal centre B cells towards a Hodgkin’s Reed-Sternberg-like phenotype. J Pathol. 2008; 216:83-92.

62. Stolz A and Bastians H. A novel role for Chk2 after DNA damage in mitosis? Cell Cycle. 2010; 9:25-26.

63. Simons A, Shaffer LG and Hastings RJ. Cytogenetic Nomenclature: Changes in the ISCN 2013 Compared to the 2009 Edition. Cytogenet Genome Res. 2013; 141:1-6.

64. Swerdlow SH, Campo E, Harris NL and al e. (2008). WHO Classification of Tumours of Haematopoietic and Lymphoid Tissues. In: Atlas of Hematologic Neoplasms 4th ed Lyon, France, IARC Press.

65. Beier CP, Kumar P, Meyer K, Leukel P, Bruttel V, Aschenbrenner I, Riemenschneider MJ, Fragoulis A, Rummele P, Lamszus K, Schulz JB, Weis J, Bogdahn U, Wischhusen J, Hau P, Spang R, et al. The cancer stem cell subtype determines immune infiltration of glioblastoma. Stem cells and development. 2012; 21:2753-2761.

66. Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, Davis AP, Dolinski K, Dwight SS, Eppig JT, Harris MA, Hill DP, Issel-Tarver L, Kasarskis A, Lewis S, Matese JC, et al. Gene ontology: tool for the unification of biology. The Gene Ontology Consortium. Nat Genet. 2000; 25:25-29.

67. Kanehisa M and Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000; 28:27-30.

68. Smyth GK, Michaud J and Scott HS. Use of within-array replicate spots for assessing differential expression in microarray experiments. Bioinformatics. 2005; 21:2067-2075.