Identification of differentially expressed genes in salivary adenoid cystic carcinoma cells associated with metastasis.

INTRODUCTION: Salivary adenoid cystic carcinoma (SACC) is a frequent type of salivary gland cancer which is characterized by slow growth but high incidence of distant metastasis. We aimed to identify therapeutic targets which are associated with metastasis of SACC.
MATERIAL AND METHODS: Total RNA was isolated from a low metastatic SACC cell line (ACC-2) and a highly metastatic SACC cell line (ACC-M), which was screened from ACC-2 by combination of in vivo selection and cloning in vitro. Then the total RNA was subjected to microarray analysis. Differentially expressed genes (DEGs) were screened from ACC-M compared with ACC-2, followed by Gene Ontology function and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis. Function annotation for DEGs also was performed. A protein-protein interaction network (PPI) was constructed for DEGs.
RESULTS: A total of 1128 DEGs were identified from ACC-M cells compared with ACC-2 cells. Both up- and down-regulated DEGs were enriched in different functions in biological process (BP), cellular component (CC) and molecular function (MF). Additionally, down-regulated DEGs were mainly enriched in "Apoptosis" and "Cytokine-cytokine receptor interaction" pathways which involved IFN-α1, NTRK1 and TGF-β1. In the PPI network, PIK3CA, PTPN11 and PIK3R1 had a number of nodes greater than 10.
CONCLUSIONS: Transforming growth factor β1 might play a pivotal role during lung metastasis of SACC and be selected as a candidate target for treatment of metastatic SACC. IFNA1, NTRK1 and PIK3CA were also associated with tumor metastasis.

Introduction

Salivary adenoid cystic carcinoma (SACC) is a frequent subtype of salivary gland malignancy, which accounts for 25% of malignant tumors in the major salivary glands and 10% of all head and neck carcinomas [1]. Salivary adenoid cystic carcinoma is characterized by slow growth but high incidence of distant metastasis, particularly to the lungs. Distant metastasis occurs in 25–50% of patients and commonly in lungs but less commonly in the liver and bone [2, 3]. However, the reason for the invasiveness and metastatic dissemination of SACC remains unclear.

Although not well understood for molecular mechanisms of cancer metastasis, angiogenesis might be a possible mechanism involved [4]. Angiogenesis is the development of new blood vessels and an important process occurring in tumor growth and metastasis. The process is intricately regulated by multiple angiogenic cytokines and other factors released by tumor cells in different pathways, such as vascular endothelial growth factor (VEGF), angiopoietin-1, basic fibroblast growth factor, platelet-derived growth factor receptor, stem cell factor receptor (c-Kit), and transforming growth factor-β1 (TGF-β1) [5, 6]. As the most notable angiogenic factor, VEGF has been shown to be highly expressed in SACC, and its expression is associated with tumor size, invasion and metastasis [7-9]. Furthermore, recent studies have indicated that overexpression of inducible nitric oxide synthase and nuclear factor κB could contribute to angiogenesis by up-regulation of VEGF in many cancers [10, 11]. Additionally, c-Kit receptor, also known as CD117, was originally identified as an oncoprotein encoded by a feline sarcoma virus. Its activation would induce diverse intracellular responses such as mitogen-activated protein kinase (MAPK) and phosphoinositide-3 kinase (PI3K/Akt) pathways [12]. Thus, c-Kit regulates blood cell survival and growth control via the aforementioned two pathways. As a result, VEGF and c-Kit have been used as targets of molecule inhibitors for treatment of SACC [13, 14]. Disappointingly, the mechanism of tumorigenesis and metastasis of SACC is not clearly elucidated and therapeutic targets are rare in clinical therapy. Because conventional chemotherapy has a poor effect in the treatment of SACC, there is a great interest in determining the molecular abnormalities in SACC.

To gain further insight into the molecular mechanism of tumorigenesis, bioinformatics and microarray are widely used in the study of cancers by researchers. However, microarray analysis is not well utilized in improving therapeutic outcomes for SACC. With the hope that doing so will achieve the goal of discovering an effective targeted therapy, we identified differentially expressed genes (DEGs) between a low metastatic SACC cell line (ACC-2) and a highly metastatic SACC cell line (ACC-M), which was screened from ACC-2 by combination of in vivo selection and cloning in vitro [15]. Since ACC-2 and ACC-M cells share an identical genetic background except for different metastatic behavior, it is presumed that the DEGs are metastasis-related genes, which play direct or indirect roles in the progression of metastasis. Then the DEGs were analyzed by Gene Ontology (GO) functions and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis. Then, a protein-protein interaction (PPI) network was built for DEGs. Taken together, the gene signature of metastasis could be helpful to develop novel therapeutic strategies in SACC patients.

Material and methods

RNA extraction and probe preparation

Total RNA from ACC-M and ACC-2 cells was isolated using the TRIzol method according to the manufacturer's (Invitrogen) instructions. RNA quality from each cell line was assessed by visualization of the 28S/18S ribosomal RNA ratio on 1% agarose gel. Total RNA samples were subjected to Human OneArray v6.1 (Phalanx Biotech, Taiwan, China), and all procedures were carried out according to the protocol. Briefly, 0.5 µg of RNA from two cell lines was labeled with a Cy3 fluorophore and labeled RNAs were hybridized at 37°C overnight.

Data preprocessing

The intensity of each probe was processed and normalized by the median scaling normalization method. In order to ensure that a probe was specific for one gene, we eliminated probes with multiple matching gene sequences. When several probes hybridized with transcripts from one gene, we calculated the mean values as the probe value. Normalized intensities were transformed to gene expression log2 ratios between ACC-M and ACC-2.

DEG screening

Because there was no extra replication except for one control group and one experimental group, we applied intensity alignment of probes between ACC-M and ACC-2 to identify DEGs. Genes with log fold change (FC) > 1 were considered to be significant.

GO function and KEGG pathway enrichment analysis

To identify gene functions enriched in ACC-M, we performed GO [16] function enrichment analysis for DEGs in 3 functional ontologies: biological process (BP), cellular component (CC) and molecular function (MF). KEGG pathway [17] enrichment analysis was also performed to identify significant pathways enriched in ACC-M with a platform developed by Feng-He Information Technology Co., Ltd (Shanghai, China). The p-value was calculated by hypergeometric distribution and a pathway with p < 0.01 was considered as significant.

Function annotation for DEGs

To ensure whether DEGs function in transcriptional regulation, transcription factor analysis was employed by mapping DEGs to the intersection between the TRANSFAC [18] and transcription activity term of the GO database. Combined with Tumor Suppressor Gene (TSG) [19] and Tumor Associated Gene (TAG) [20] databases, we also obtained known oncogenes and suppressor genes from identified DEGs.

Construction of PPI network

To study protein-protein association information for DEGs, the STRING database [21] was used to construct the PPI network. The selected protein pairs in PPI with an association score more than 0.9 and the number of nodes more than 3 were products of DEGs.

Results

DEG identification

To identify significant genes between ACC-M and ACC-2, DEG identification was performed. A total of 1128 DEGs were obtained including 448 up- and 680 down-regulated DEGs.

GO function enrichment analysis

To study the function changes in the process of tumor metastasis, we identified over-presented GO categories in BP, CC and MF for both up- and down-regulated DEGs. The top 5 categories for 3 type GO terms are listed in Tables I–III. From the results, up-regulated DEGs were mainly enriched in “reflex” and “synaptic transmission glycinergic” in BP, “extracellular region” and “integral to mitochondrial membrane” in CC, and “inhibitory extracellular ligand-gated ion channel activity” in MF. Down-regulated DEGs were mainly enriched in “adenylate cyclase-activating G-protein coupled receptor signaling pathway” in BP, “plasma membrane-related functions” in CC and “G-protein coupled amine receptor activity” in MF.

Table I

Gene Ontology function enrichment analysis for both up- and down-regulated differentially expressed genes in Biological Process (top 5 GO terms were listed)

GO-ID		GO-term	Clustered gene numbers	Typical gene symbol	P-value
Up-regulated	GO:0060004	Reflex	4	GJA1, GLRA1, GLRB, NPNT	0.00026
	GO:0060012	Synaptic transmission, glycinergic	2	GLRA1, GLRB	0.00042
	GO:0046903	Secretion	32	CD274, MMP7, PRKAR1A, TRPC1WNK4, WNT5A	0.00057
	GO:0043030	Regulation of macrophage activation	4	BPI, IL4, SNCA, WNT5A	0.00115
	GO:0006888	ER to Golgi vesicle-mediated transport	6	INSIG1, LMAN1, SAR1B, SEC23A, SEC23B, SPAST	0.00120
Down-regulated	GO:0007189	Adenylate cyclase-activating G-protein coupled receptor signaling pathway	9	ADCY1, ADCY3, ADRB1, ADRB2, OPRM1, VIP	2.08E-05
	GO:0051046	Regulation of secretion	31	IL1A, TGF-β1, CORIN, FAM3D, TNFSF15,TRAF6	7.08E-05
	GO:0002526	Acute inflammatory response	13	APOL2, CD163, IL1A, IL20RB, VCAM1, VTN	7.23E-05
	GO:0019233	Sensory perception of pain	10	NTRK1, OPRD1, ADRB1,CHRNB2, P2RX4, VIP	8.45E-05
	GO:0051049	Regulation of transport	55	IL1A, TGF-β1, BCL2, BEST1, PIK3R1, TRIB3	0.00021

Table III

Gene Ontology function enrichment analysis for both up- and down-regulated differentially expressed genes in Molecular Function (top 5 GO terms were listed)

GO-ID		GO-term	Clustered gene numbers	Typical gene symbol	P-value
Up-regulated	GO:0005237	Inhibitory extracellular ligand-gated ion channel activity	3	GABRB2, GLRA1, GLRB	8.71E-05
	GO:0016933	Extracellular-glycine-gated ion channel activity	2	GLRA1, GLRB	0.00253
	GO:0016934	Extracellular-glycine-gated chloride channel activity	2	GLRA1, GLRB	0.00253
	GO:0015276	Ligand-gated ion channel activity	9	CNGA3, GABRA1, GLRB, P2RX2, TRPC1, TRPC5	0.00256
	GO:0022834	Ligand-gated channel activity	9	CNGA3, GABRA1, GLRB, P2RX2, TRPC1, TRPC5	0.00256
Down-regulated	GO:0008227	G-protein coupled amine receptor activity	6	ADRB1, ADRB2, HTR1A, HTR2A, TAAR2, TAAR6	0.00227
	GO:0004939	β-Adrenergic receptor activity	2	ADRB1, ADRB2	0.00298
	GO:0004867	Serine-type endopeptidase inhibitor activity	9	CPAMD8, OVOS, SERPINA1, SERPINA5, SPINK9, SPOCK1	0.00327
	GO:0004252	Serine-type endopeptidase activity	12	CORIN, F11, PRSS54, ST14, MST1L, PCSK5	0.00328
	GO:0051537	2 iron, 2 sulfur cluster binding	4	AIFM3, CISD3, GLRX2, NDUFS1	0.00329

KEGG pathway enrichment analysis

To gain further insight into the function of DEGs, we used the platform developed by Feng-He Information Technology Co., Ltd (Shanghai, China) to identify the significant pathways. With the selected criteria we finally obtained 9 pathways (Table IV). The up-regulated DEGs were mainly enriched in the “protein processing in endoplasmic reticulum” pathway, while down-regulated DEGs were significantly enriched in “Calcium signaling pathway” and “Olfactory transduction”. From the pathway information in the KEGG database, we found that just “Apoptosis” and “Cytokine-cytokine receptor interaction” pathways were associated with tumorigenesis.

Table IV

Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis for both up- and down-regulated differentially expressed genes

KEGG-ID		KEGG-term	Clustered gene numbers	Typical gene symbol	P-value
Up-regulated	KEGG:4141	Protein processing in endoplasmic reticulum	11	HSPA2, HSPA5, HSPA8, LMAN1, MAN1A1, PDIA4	0.00063
Down-regulated	KEGG:4020	Calcium signaling pathway	15	ADCY1, ADCY3, ATP2B1, GRIN2C, HTR2A, P2RX1	0.00047
	KEGG:4740	Olfactory transduction	24	ADCY3, OR1A2, OR1L8, OR2A14, OR2G3, OR2V2	0.00120
	KEGG:4210	Apoptosis	9	TGF-β1, IFNA1, NTRK1, PIK3R1, PIK3R2, TNFSF10	0.00167
	KEGG:4080	Neuroactive ligand-receptor interaction	18	ADRB1, ADRB2, CHRNB2, F2, GRID1, GRIN2C	0.00243
	KEGG:4910	Insulin signaling pathway	11	G6PC2, IRS1, PCK2, PIK3R1, PIK3R2, SHC2	0.00434
	KEGG:4540	Gap junction	8	ADCY1, ADCY3, ADRB1, HTR2A, PDGFC, PLCB	0.00761
	KEGG:5414	Dilated cardiomyopathy	8	TGF-β1, PRKACG, ADCY1, ADCY3, ADRB1, TNNT2	0.00761
	KEGG:4060	Cytokine-cytokine receptor interaction	16	TGF-β1, IFNA1, NTRK1, CCL16, CCR9, TNFSF10	0.00998

Function annotation of DEGs

Combined with genes in Tumor Suppressor Gene (TSG) and Tumor Associated Gene (TAG) databases, we screened known tumor suppressor genes and oncogenes in both up- and down-regulated DEGs. In 25 DEGs enriched in tumor-related “Apoptosis” and “Cytokine-cytokine receptor interaction” pathways, we found that IFN-α1 (interferon α1) was a tumor suppressor gene but neurotrophic tyrosine kinase, receptor, type 1 (NTRK1) was an oncogene. Just transforming growth factor β1 (TGF-β1) was associated with lung diseases but not with transcriptional regulatory function.

Protein-protein interaction network

After PPI network construction with the criteria of an association score higher than 0.9 and the number of nodes more than 3, we selected 32 differentially expressed proteins in the network (Figure 1). Moreover, there were 3 DEGs with the number of nodes higher than 10 including 2 up-regulated proteins – phosphatidylinositol-4,5-bisphosphate 3-kinase, catalytic subunit α (PIK3CA) and protein tyrosine phosphatase, non-receptor type 11 (PTPN11) – and one down-regulated protein, phosphoinositide-3-kinase, regulatory subunit 1α (PIK3R1).

Figure 1

Protein-protein interaction network for differentially expressed genes. Red color represents up-regulated genes, and green color represents down-regulated genes. Diamonds represent genes with the number of nodes more than 10. Square represents oncogene which associates with “Apoptosis” and “Cytokine-cytokine receptor interaction” pathways

Discussion

Salivary adenoid cystic carcinoma is a common malignant tumor that arises from the secretory epithelial cells of salivary glands and has a unique potential for cell invasion and distant migration, particularly to the lungs. To date, there are very few prospective studies of molecular biomarkers which can be used to predict prognosis of SACC or therapeutic response. In this study, DEGs were screened from a highly metastatic cell line (ACC-M) compared with a low metastatic cell line (ACC-2), followed by functional annotation and PPI network construction for DEGs. DEGs were significantly enriched in “Apoptosis” and “Cytokine-cytokine receptor interaction” pathways which were associated with tumors in the KEGG database, which involved 3 significant DEGs: NTRK1, IFNA1 and TGF-β1. Moreover, PIK3CA, PIK3R1 and PTPN11 were also important for SACC in the PPI network.

Transforming growth factor β1 is a multifunctional polypeptide which plays vital roles in apoptosis, cell proliferation and epithelial-mesenchymal transition [22-26]. Previous studies have suggested that TGF-β1 signaling could promote invasion and metastasis of breast cancer and glioma [27, 28]. In addition, it has been proved that overexpression of active TGF-β1 in vivo accelerates lung metastases of transgenic mammary tumors with an effect on tumor size or tumor cell proliferation [29]. Moreover, it was found that TGF-β1 expression was significantly increased in human primary SACC samples with metastasis [30]. Further study showed that TGF-β1 could promote migration and invasion of SACC via TGF-β1/Smad signaling and induce epithelial-mesenchymal transition in normal stromal cells or epithelial cells of SACC [31]. On the other hand, TGF-β1 is a pluripotent cytokine with dual roles and is also considered as a tumor suppressor in tumorigenesis [32]. Inconsistent with the previous study, our microarray test showed that TGF-β1 was down-regulated in a highly metastatic SACC cell line. From the results of pathway enrichment analysis, TGF-β1 was found to be enriched in the “Cytokine-cytokine receptor interaction” pathway. Interaction between the cytokine and cytokine receptor leads to multiple biological responses including prevention of tumor cell apoptosis and cell survival [33]. As a result, we speculate that a decrease of the cytokine TGF-β1, as a tumor suppressor, regulates lung metastasis of SACC through the “cytokine-cytokine receptor interaction” pathway.

NTRK1, as a high affinity receptor for nerve growth factor (NGF), is a member of the neurotrophin receptor family. After binding with NGF, NTRK1 undergoes dimerization and autophosphorylation of tyrosine residues [34]. Then apoptosis, proliferation and differentiation-related proteins including phosphatidylinositol 3-kinase, Rac, Ras and mitogen activated protein kinase, respectively, are activated [35]. These proteins can stimulate apoptosis of tumor cells. In our study, NTRK1 was found to be down-regulated and enriched in the pathway “Apoptosis”. As a result, decreased NTRK1 may contribute to anti-apoptosis of cancer cells and tumorigenesis. Additionally, NTRK2 has been found to be implicated in the pathogenesis of lung cancers [36]. Thus, NTRK2 may be associated with lung metastasis of SACC.

INF-α1, a member of the interferon family, is known as a tumor suppressor gene. A previous study showed that the cytokine mediates its apoptotic effects by inducing expression of a tumor necrosis factor (TNF)-related apoptosis-inducing ligand [37]. Thus, down-regulated INF-α1 may accelerate tumor development. PIK3CA encodes the p110α catalytic subunit of phosphatidylinositol-3-kinase (PI3K), which has been shown to be expressed in cancers, such as lung and ovarian cancer [38, 39]. Previous studies have indicated that PIK3CA is regulated by transcription factor p53, which is involved in apoptosis and cell cycle arrest [40]. Moreover, inactivation of p53 leads to up-regulation of PIK3CA and appears to be an early step in ovarian carcinogenesis [39]. Furthermore, PIK3CA has been used as a candidate target for inhibition of metastatic tumor growth in bone-metastatic tumors [41]. Consequently, we infer that PIK3CA might be associated with lung metastasis in SACC.

In conclusion, we have identified several DEGs with change of expression in SACC. Transforming growth factor-β1 might be a unique molecule in lung metastasis and could be selected as a candidate target for clinical therapy. NTRK1, IFNA1 and PIK3CA are also associated with tumor development of SACC. Nevertheless, a limitation has to be mentioned. There is no biological replicate in the study. The analysis is merely a qualitative experiment, but it can partly elucidate molecular mechanisms of lung metastasis of SACC. Further experiments are needed, enrolling more biological replicates and verifying the results of bioinformatics.

Table II

Gene Ontology function enrichment analysis for both up- and down-regulated differentially expressed genes in Cellular Component (top 5 GO terms were listed)

GO-ID		GO-term	Clustered gene numbers	Typical gene symbol	P-value
Up-regulated	GO:0005576	Extracellular region	70	CCL20, CD274, IGSF10, MMP7, NGFR, WNT5A	0.00040
	GO:0032592	Integral to mitochondrial membrane	4	CPT1A, ETFDH, MFF, TMEM70	0.00528
	GO:0044421	Extracellular region part	38	CCL20, CD274, IL4, MMP26, MMP7, WNT5A	0.00539
	GO:0031301	Integral to organelle membrane	11	A4GALT, CPT1A, HSPA5, INSIG1, QSOX1, TMEM70	0.00573
	GO:0005593	FACIT collagen	2	COL14A1, COL16A1	0.00655
Down-regulated	GO:0031226	Intrinsic to plasma membrane	62	ADAM29, ADCY3, CCR9, CD163, NTRK1, TPO	0.00048
	GO:0005887	Integral to plasma membrane	59	ADAM29, ADCY3CCR9, CD163, NTRK1, TACR2	0.00089
	GO:0044459	Plasma membrane part	84	ADAM29, ADCY3, CCR9, CD163, NTRK1, OPRD1	0.00143
	GO:0044304	Main axon	6	ADORA2A, MYO1D, NCMAP, ROBO1, SCN8A, SPOCK1	0.00253
	GO:0005576	Extracellular region	90	CCL16, CD163, IL1A, IL22, TEK, TGF-β1	0.00299