Identification of Differentially Expressed Genes in Human Salivary Gland Tumors by DNA Microarrays¹

Pleomorphic adenoma is the most common benign tumor of major salivary glands. It shows a marked histological diversity with epithelial, myoepithelial, and mesenchymal components in a variety of patterns (1). Cytogenetic analysis have recognized, in addition to a subgroup with normal karyotypes, recurrent chromosomal translocations, with breakpoints at 8q12, 3p21, and 12q13–15 (2) and corresponding to PLAG1 (3), β-catenin (CTNNB1; Ref. 4), and HMGIC (5) genes, respectively. The second most common benign lesion of the salivary glands is Warthin’s tumor. It consists of abundant lymphoid stroma and bilayered oncocytic epithelium that forms cystic or papillary structures. The epithelial component represents the neoplastic proliferation of salivary ducts that have been entrapped within lymph nodes associated with the salivary gland.

Mucoepidermoid carcinoma is the most common malignant, locally invasive tumor of the salivary glands, especially of the parotid gland. Some markers of tumor progression, invasiveness, and prognosis are p21(Kip1), a cyclin-dependent kinase inhibitor, the oncoproteins Bcl-2 and Bax, the tumor suppressor gene product p53, terminal deoxynucleotidyltransferase-mediated nick end labeling staining, and the cell cycle antigen Ki-67 (6, 7). A rare form of malignant tumor, arising in pleomorphic adenoma of the salivary glands, is clear cell carcinoma. This tumor shows definite areas characteristic of pleomorphic adenoma combined with areas revealing evidence of malignancy.

The emerging technology of cDNA microarrays provides the ability to comparatively analyze mRNA expression of thousands of genes in parallel (8). Several studies have already demonstrated the usefulness of this technique for identifying novel cancer-related genes and for classifying human cancer at the molecular level (9–11). In this study, we have used a cDNA microarray containing 19,200 different human probes corresponding to different Unigene clusters to characterize benign and malignant salivary gland tumors. The main aim of the study was to identify gene expression alterations that are either common to most of the tumors studied or specific to clinical features of the tumors.

Samples were derived from surgically resected parotid tissue from five patients with pleomorphic adenoma, four with Warthin’s tumor, one with clear cell carcinoma, and two with mucoepidermoid carcinoma. Informed consent was obtained from patients to use their operated specimens for research purposes. Part of each sample was frozen in liquid nitrogen and stored at −80°C until RNA extraction, and the remaining part was used for histopathological analysis. Tumor paired control non-tumor tissues were obtained from a clinically unaffected site and were histologically normal.

Total RNA was extracted from the biopsies by using RNAzol. Ten μg of total RNA were used for each sample. cDNA was synthesized by using Superscript II (Life Technologies, Inc.) and amino-allyl dUTP (Sigma). Monoreactive Cy3 and Cy5 esters (Amersham Pharmacia) were used for indirect cDNA labeling. A pool of eight non-tumor salivary glands RNA was labeled with Cy3 and used as control against the Cy5-labeled tumor cDNA. Human 19.2 K DNA microarrays, containing ESTs ³ corresponding to at least 15,448 different Unigene clusters, were used (Ontario Cancer Institute). ⁴ One hundred μl of the sample and control cDNA in DIG Easy hybridization solution (Roche) were used in a sandwich hybridization of the two slides constituting the 19.2 K set at 37°C overnight. Washing was performed three times for 10 min with 1× SSC, 0.1% SDS at 42°C and three times for 5 min with 0.1× SSC at room temperature. Slides were dried by centrifugation for 2 min at 2000 rpm. Hybridized slides were scanned using the GenePix 4000A scanner (Axon Instruments, Foster City, CA). Arabidopsis RNA was used as a reference for RNA labeling.

Images were analyzed by using GENEPIX PRO 3.0 (Axon Instruments). Spots showing no signal or obvious defects were accordingly flagged by GENEPIX PRO and excluded from the analysis. Local background was subtracted from the remaining spots, and then spots with fluorescence intensity in both channels of <1000 were stringently flagged as absent. The background level for a typical array was of 67.5 (SD, 145.2) for the Cy5 channel and of 182.7 (SD, 305.4) for Cy3 channel. Maximum fluorescence intensity is at about 65,500 in both channels, and arrays were scanned ensuring that the <1% of the total number of spots measured on the slide had median intensities close to maximal. Ratios of net fluorescence from the Cy5-specific channel to the net fluorescence from the Cy3-specific channel were calculated for each spot, representing tumor mRNA expression relative to the corresponding non-tumor salivary gland tissue. A normalization factor was estimated from ratios of median by GENEPIX PRO. Cy5:Cy3 expression ratios were then log-transformed (base 2) and normalized, for each different array, by adding the log₂ of the respective normalization factor to the log₂ of the ratio of medians for each spot within the array, so that the average log-transformed ratio equaled zero. The threshold for significant RNA expression changes (3.0-fold; i.e., ∼1.5 on the log₂ scale) was established as three times the SD of an assay, where the same RNA sample was independently retrotranscribed and labeled with both cyanines. DNA spots present in at least 75% of the arrays (5,681 spotted cDNAs) and with expression ratios higher than the above-defined threshold, in at least one array, were selected for the clustering analysis (4,434 spotted cDNAs).

To select the probes that distinguish between two sample groups (e.g., pleomorphic adenoma versus Warthin’s tumor), we developed a program performing two statistical tests, t test and log likelihood. Log likelihood was defined as the difference between the mean of each probe in the two groups. To perform multiple t testing correction, the Benjamini and Hochberg procedure, as modified by Storey and Tibshirani (12), was applied. Briefly, it controls the pFDR, defined as the proportion of probes expected to be identified by chance (assuming probes are independent) relative to the total number of probes called significant. This procedure provides a good balance between discovery of significant spots and protection against false positives, because occurrence of the latter is held to a small proportion of the list. A balanced bootstrap analysis with 10,000 permutations was therefore applied to the tests to define the number of expected false positives by chance. Expression tables were analyzed by using hierarchical clustering, PCA, and similarity searches tools in JExpress v.2.1. ⁵ Variance normalization across the expression table was performed only prior to similarity search to identify probes with similar expression pattern and was obtained by normalizing the mean and then multiplying each vector by a scalar so that its variance becomes unity.

We compared the expression profiles of 12 salivary gland tumors and four normal tissue samples by hybridization to microarrays containing 19,200 cDNAs. Because the arrays shared a common normal tissue reference made of a pool of eight non-tumor salivary glands different from those used as test samples, the samples could be compared to identify gene expression patterns correlated with clinical features of the tumors.

We used hierarchical clustering (13) to look at the variation in gene expression among tumors. A clustering analysis was performed by using the selected 4,434 spotted cDNAs, as detailed in “Materials and Methods,” corresponding to 2,333 distinct Unigene clusters present in 75% of the arrays and with Cy5:Cy3 expression ratios that varied at least 3.0-fold in at least one array. By using the t test, we then sought to identify gene expression alterations differentiating the tumors from normal samples. t test selected 713 spotted cDNA probes (P < 0.01, average t test significant spots in bootstrap, 91.82; pFDR, 0.05) corresponding to 486 distinct Unigene clusters. The clustering based on the 486-probe set very well correlated with the tumors histotypes, with the sole exception of the mucoepidermoid carcinoma (Fig. 1). Warthin’s tumors were the most distant branch, with the non-tumor salivary glands in a separate branch from that of the pleomorphic adenoma. The single case of cell carcinoma, arisen from pleomorphic tumor, correctly clustered within the pleomorphic adenoma group. Table 1 lists differentially expressed cancer-associated genes with the most highly significant t test values. ⁶ The list of up- and down-regulated genes in most salivary gland tumors include known genes involved in neoplasia, members of signal transduction pathways, transcription factors, oncogenes, and tumor suppressor genes, cell cycle regulator, and genes that encode protein involved in cellular metabolism and adhesion. Among up-regulated transcription factors are present forkhead box O1A (FOXO1A) and PAX1.Fusion of forkhead box O1A with PAX3, a PAX1-related gene, is associated with solid tumor rhabdomyosarcoma. Overexpression of the fusion protein actives the transcription of antiapoptotic protein Bcl-XL (14), capable of suppressing cell death pathways by blocking the release of cytochrome c from mitochondria and the caspase protease cascade. By using a similarity search with Euclidean distances on the variance normalized expression table, as obtained by using the correspondent JExpress plugin, we identified genes with coherent expression patterns to those of some genes selected by the t test (Fig. 2). Interestingly, among gene subsets that display expression patterns coherent with FOXO1Athere are the PAX1, homologous to the FOX fused gene above described, a novel Bcl-2 homologue (MIL1or Bcl-rambo), and the 7b and 7c subunits of cytochrome c oxidase (COX7B, COX6C). Cytochrome c oxidase down-regulation, alongside to cytochrome c release to the cytosol, was found to be present in Bax-induced growth arrest (15). Overrepresentation of this gene subset suggests that the activation of mechanisms that block apoptosis is a major characteristic of the tumor progression in salivary glands.

One of the very few genes consistently different in every tumor is OTF1, highly expressed in non-tumor salivary glands and underexpressed in all tumors. OTF1 protein binds specifically its POU domain to octamer DNA motifs present in the promoters of several genes and regulates their expression. In addition, OTF1 participates in the cellular response to DNA damage and is involved in the regulation of stress-inducible genes (16). By using semiquantitative reverse transcription, we confirmed the expression levels of OTF1, Wee1, and EMP2(data not shown). Additionally, we identified a number of up- and down-regulated ESTs, which suggest the involvement of a large number of genes that are not yet functionally characterized in salivary gland tumors.

Defined subsets of clones that most strongly defined the division of biopsies by histotype were identified by using the t test. When we applied this test to compare each tumor to non-tumor samples, we selected 599 cDNA elements differentially expressed (P < 0.01, average t test significant spots in bootstrap, 49.95; pFDR, 0.04) in pleomorphic adenoma and 2586 elements (P < 0.01, average t test significant spots in bootstrap, 5.73; pFDR, <0.001) specifically expressed in Warthin’s tumors. The many clones associated with Warthin’s tumor is probably related to the presence of lymphoid stroma in the affected glands, because some of the selected genes are T-cell specific (i.e., NFATC3 and ITGAE, an integrin preferentially expressed on intraepithelial lymphocytes). One thousand and one hundred ninety-eight clones were differentially expressed (P < 0.001, average t test significant spots in bootstrap, 0.13; pFDR, <0.001) in pleomorphic adenoma when compared with Warthin’s tumor, of which 479 are overexpressed in Warthin’s, whereas 719 are overexpressed in pleomorphic adenoma. The specificity of clones expression is detailed in Table 1. Plag1, a gene that is often rearranged in pleomorphic tumor, is surprisingly significantly overexpressed in Warthin’s tumor (P < 0.01). In the Plag1 expression pattern group (Fig. 2), other genes present are: PLK1, progranulin, which elevated expression confers a transformed phenotype on epithelial cells, including anchorage independence in vitro and growth as tumors in nude mice (17); hypothetical protein FLJ13605 containing a tetratricopeptide repeat TPR domain, which forms scaffolds to mediate protein-protein interactions and often the assembly of multiprotein complexes and has a spatial arrangement of α-helices similar to those within 14-3-3 proteins; αL subunit of integrin leukocyte function antigen-1, a membrane glycoprotein that functions in cell-cell adhesion by heterophilic interaction with intercellular adhesion molecule 1; and polo-like serine/threonine kinase PLK1 and TXK, a member of the Tec nonreceptor subfamily of Src protein tyrosine kinases.

Finally, we selected a subset of clones that most strongly defined the division between malignant and benign lesions. By using the log likelihood test, and as log2 threshold a value of 1.5, corresponding to 3× SD from the average, we identified 16 Unigene clusters differentiating this malignant clear cell carcinoma from its related benign pleomorphic adenoma (P < 0.01, as determined by bootstrap distribution). On the contrary, this clear cell carcinoma differed from the Warthin’s tumors for a much larger number of cDNA elements (1,315), as expected for these two nonrelated tumors, whereas 133 elements were different from the non-tumor reference pool. Genes differentiating this clear cell carcinoma from its related benign tumor pleomorphic adenoma include myelin transcription factor 1-like, huntingtin, ubiquitin specific protease 9, and FLJ10355. ⁷

Because, as shown by hierarchical clustering (Fig. 1), Warthin’s tumor was very divergent from the other biopsies, we then compared the mucoepidermoid carcinoma to all of the other samples after the exclusion of Warthin’s tumor samples. By using the t test, we selected 57 cDNA probes (P < 0.01). Among relevant genes repressed in mucoepidermoid carcinoma, there is the PSG11member of the carcinoembryonic antigen family and two genes involved in apoptosis, CIDEB and ADORA1. CIDEB or cell death-inducing DFFA-like effector b, is involved in induction of apoptosis by DNA damage (18), and ADORA1, adenosine A1 receptor, is involved in induction of apoptosis by extracellular signals. Interestingly, these genes are repressed only in mucoepidermoid carcinoma. FLJ10355, which is overexpressed in the mucoepidermoid carcinomas, is also overexpressed in the other malignant tumor clear cell carcinoma and not expressed in the benign tumors. FLJ10355has some homology to Bcl9, a gene translocated in some B-cell malignancies with abnormalities of 1q21 (19) and with no recognizable protein motifs or significant homologies to any other known proteins.

In summary, the results show gene expression alterations common to all tumors and unique to each benign and malignant tumor. A larger epidemiological study is required to determine whether expression levels of the selected genes can be used as prognostic indicators or else additional genes need to be included in the analysis. Finally, the numerous ESTs of unknown function we identified could also become useful as tumor markers and represent a set of novel tumor-associated genes.