Identification of genes associated with matrix metalloproteinases in invasive lung adenocarcinoma.

The aim of the present study was to identify genes with similar function to that of matrix metalloproteinases (MMPs) in invasive lung adenocarcinoma (AC) and to screen the transcription factors that regulate MMPs. The gene expression dataset GSE2514, including 20 invasive lung AC samples and 19 adjacent normal lung samples, was downloaded from the Gene Expression Omnibus database. Differentially expressed genes (DEGs) were screened using the limma package in R. Genes with similar function to MMPs were identified by K-means clustering. Their correlations with MMPs were validated using Pearson correlation analysis. The expression of MMPs in lung cancer and normal tissues was evaluated by western blot analysis. Protein-protein interaction (PPI) network and transcriptional regulatory network analyses were performed with Retrieval of Interacting Genes and Database for Annotation, Visualization and Integrated Discovery, respectively. As a result, 269 DEGs were identified between invasive lung AC samples and normal lung samples, including 78 upregulated and 191 downregulated genes. Four MMPs (MMP1, MMP7, MMP9 and MMP12), which were upregulated in lung AC, were clustered into one group with other genes, including NAD(P)H quinone oxidoreductase 1, claudin 3 (CLDN3), S100 calcium-binding protein P, serine protease inhibitor Kazal type 1, collagen type XI α 1 chain, periostin and desmoplakin (DSP), following cluster analysis. Pearson correlation analysis further confirmed correlations between MMP9-CLDN3, MMP9-DSP and MMP12-DSP. PPI network analysis also indicated multiple interactions between MMPs-associated genes. Furthermore, MMPs were commonly regulated by CCAAT/enhancer binding protein α transcription factor. These findings may provide further insight into the mechanisms of MMPs in invasive lung AC.

Introduction

Lung cancer is a leading cause of cancer mortalities worldwide, with 226,160 new cases and 160,340 mortalities in 2012 (1). As a type of lung cancer, non-small cell lung carcinoma (NSCLC) consists of adenocarcinoma (AC), squamous cell carcinoma and large cell carcinoma (2). Lung AC, a heterogeneous group of tumors ranging in aggressiveness from noninvasive bronchioloalveolar carcinoma (BAC) to microinvasive tumors, mixed-type tumors and pure invasive AC, is the most common histological type of lung cancer (3,4). The first treatment for patients with early-stage lung AC is surgical resection, but the 5-year overall survival rate remains at ~80% for stage IA disease (5). Notably, patients with invasive lung AC have a lower mean 5-year survival rate compared with that of patients with BAC (59 vs. 100%, respectively) (6). Therefore, it is important to better understand the molecular mechanisms of invasive lung AC in order to develop effective preventive and curative strategies.

Cancer invasion is a highly complicated process, with the early events being the proteolytic degradation of extracellular matrix (ECM) components to provide room for infiltration, followed by migration into the adjacent tissues (7,8). Recent studies suggest that matrix metalloproteinases (MMPs), members of the matrixin subfamily of zinc metalloproteinases, are involved in the breakdown of the ECM (9,10). Over-expression of MMPs has been detected in a number of tumor types, including invasive lung AC (11). MMP2, MMP7, MMP9 and tissue inhibitor of metalloproteinase 2 are upregulated in the lung tissue of patients with primary spontaneous pneumothorax (PSP), and the imbalance of their expression may be implicated in the pathogenesis of PSP (12). The signal module and sequence variant module appearing in lung AC, in which the expression of MMP12 is upregulated but that of MMP11 is downregulated, can promote the invasion of cancer cells (13). Chemotherapy drugs can markedly inhibit the invasive ability of human lung AC cells via reducing the expression of MMP2 and MMP9 (14,15). Additionally, rosuvastatin and simvastatin can function in the treatment of lung cancer by regulating the expression of MMP2, MMP9, RAS and nuclear factor-κB-p65 (16). However, the current number of screened MMPs is limited, and further studies are required to identify genes with a similar function to that of MMPs.

In addition, it has been reported that the transcription factor signal transducer and activator of transcription 3 (STAT3) can regulate the transcription of MMP2 and promote melanoma cell invasion (17). Blocking STAT3 signaling significantly inhibits the invasion of melanoma cells (17). Storz et al (17) have demonstrated that forkhead box O3 promotes invasion and progression in numerous human solid tumors by inducing the expression of MMP9 and MMP13. However, the regulatory mechanisms of MMPs are not well understood in lung AC.

The present study aimed to identify differentially expressed genes (DEGs), particularly MMPs-associated genes, using gene expression profile data of invasive lung AC and adjacent normal samples collected from a publicly available database. Besides, the transcriptional factors regulating MMPs-associated genes were also investigated.

Materials and methods

Ethical statement

The present study was approved by the Ethics Committees of the Beijing Shijitan Hospital (Beijing, China) and the Chest Hospital Affiliated to Shanghai Jiaotong University (Shanghai, China), and it was performed in accordance with the ethical standards (18). In addition, written informed consent was obtained from all patients, prior to enrollment in the present study.

Microarray data

The gene expression dataset GSE2514, which is based on two platforms [GPL81 (MG_U74Av2) Affymetrix Murine Genome U74A Version 2 Array and GPL8300 (HG_U95Av2) Affymetrix Human Genome U95 Version 2 Array], was downloaded from Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/) database (19). Microarray data obtained with the platform GPL8300 (HG_U95Av2; Affymetrix Human Genome U95 Version 2 Array; Affymetrix, Inc., Santa Clara, CA, USA) was downloaded from the GEO database (19) in our study. It contained 20 lung AC samples and 19 adjacent lung samples (~1 cm away from the tumor site), which were obtained from 5 male and 5 female patients undergoing lobectomy (9) and wedge resection (1). Of these patients, 9 had a history of tobacco smoking. The ages of the patients ranged from 45 to 73 years. Their tumors were all invasive lung AC. The majority of tumors were low-to-intermediate grade and low stage, although 2 stage III tumors were included in the analysis. Probe annotation files were also acquired.

Preprocessing and differential analysis

The raw array data in the CEL file (Affymetrix, Inc.) were transformed into recognizable gene expression data using the robust multi-array average algorithm from the affy package in R (20). Upon normalization, probes were mapped to genes according to the annotation files. The levels of probes corresponding to one gene were averaged and used as the final gene expression value. DEGs were screened using the limma package (http://www.bioconductor.org/packages/release/bioc/html/limma.html) (21) in R based on the cut-offs of P<0.05 and |log fold-change| >1.

Cluster analysis

Cluster analysis of DEGs was performed with Cluster 3.0 (http://bonsai.ims.u-tokyo.ac.jp/mdehoon/software/cluster) using the K-means clustering algorithm (22), which was conducted on data with K (the number of clusters)=5.

Pearson correlation analysis

The method of Pearson correlation analysis (23) in the cor function in R (https://www.r-project.org/) was used to perform the correlation analysis between the expression of DEGs and MMPs. A correlation coefficient of >0.8 was used as the cut-off criterion.

Western blot analysis

A total of 6 lung AC and matched adjacent lung tissue biopsy samples were obtained in June 2014 from Shanghai Lung Cancer Center, Chest Hospital Affiliated to Shanghai Jiaotong University (Shanghai, China). Tissues were washed with ice-cold PBS and lysed in ice-cold radioimmunoprecipitation assay lysis buffer (Sangon Biotech Co., Ltd., Shanghai, China). Extracted proteins were quantified using a BCA Protein Assay kit (Sangon Biotech Co., Ltd., Shanghai, China) and separated with SDS-PAGE with a 10% separating gel and a 5% stacking gel. Subsequently, proteins were transferred onto polyvinylidene difluoride membranes (EMD Millipore, Billerica, MA, USA). Then, membranes were blocked in 5% nonfat dry milk and probed with primary rabbit antibodies against MMPs (anti-MMP1, cat. no. sc-58377; anti-MMP7, cat. no. sc-80205; anti-MMP9, cat. no. sc-13520; and anti-MMP12, cat. no. sc-8839; all diluted to 1:500; Santa Cruz Biotechnology, Inc., Dallas, TX, USA) overnight at 4°C. Horseradish peroxidase-conjugated anti-rabbit immunoglobulin G (cat. no. 111-035-003; 1:5,000; Jackson ImmunoResearch Laboratories, Inc., West Grove, PA, USA) served as a secondary antibody and was incubated for 1 h at 4°C. The anti-β-actin antibody (cat. no. 8227; 1:10,000; Abcam, Cambridge, CA, USA) was used as the control. The position of protein bands was developed with ECL chemoluminescence kit (Merck Millipore) and visualized under a ChemiDoc MP imaging system (Bio-Rad Laboratories, Inc., Hercules, CA, USA).

Functional and pathway enrichment analysis

Gene Ontology (GO; http://www.geneontology.org/) (24) and Kyoto Encyclopedia of Genes and Genomes (KEGG; http://www.genome.jp/kegg/) (25) pathway enrichment analyses were performed for the MMPs-associated genes using the Database for Annotation, Visualization and Integrated Discovery (DAVID) online tool (26) to reveal altered biological functions in lung AC. P<0.05 and a false discovery rate (FDR) adjusted by the Benjamini and Hochberg method (27) of <0.01 were set as the thresholds.

Selection of disease-associated genes

Diseases associated with the DEGs were identified by gene set enrichment analysis against the Genetic Association Database Disease Class (http://geneticassociationdb.nih.gov/) using the annotation server DAVID (P<0.05) (26).

Construction of protein-protein interaction (PPI) network

PPI network analysis for the DEGs was carried out with the Search Tool for the Retrieval of Interacting Genes (STRING; http://string-db.org/), and the PPI network involving MMPs-associated genes was selected. Interaction associations with a combined score of ≥0.4 were retained and the PPI network was visualized by Cytoscape software (version 2.8; http://www.cytoscape.org) (28).

Transcriptional regulatory network analysis

Transcriptional regulatory network analysis was performed with DAVID (26) for the group of DEGs including MMPs, and the transcriptional regulatory network was visualized by Cytoscape software (version 2.8; http://www.cytoscape.org) (28).

Results

DEGs analysis

Upon normalization of the raw data (Fig. 1), the DEGs were screened using the limma package in R. As a result, 269 DEGs were identified between lung AC and adjacent normal lung samples, including 78 upregulated and 191 downregulated genes.

Figure 1.

Box plot for normalized gene expression data. Horizontal axis and vertical axis separately represent sample names and expression values. The blue box and the black line represent expression values and their median, respectively. The medians (black lines) are almost all at the same level, indicating a good performance of normalization.

MMPs-associated genes screening

MMP1, MMP7, MMP9 and MMP12 were clustered into one group with a number of other DEGs (Table I), which were upregulated in lung AC compared with their expression in normal lung tissues. These results were further confirmed by Pearson correlation analysis and western blotting. The present findings also demonstrated that MMP9 expression exhibited a significant correlation with DSP (correlation coefficient=0.8309295) and CLDN3 (correlation coefficient=0.8058015) expression. Additionally, MMP12 expression was significantly correlated with DSP expression (correlation coefficient=0.8127249). The expression of MMP1, MMP7, MMP9 and MMP12 was also observed to be upregulated in lung AC compared with that in adjacent lung tissues (Fig. 2).

Table I.

List of genes clustered with MMPs.

ID	Log (fold-change)	Mean expression	t	P-value	False discovery rate
RAMP1	1.803939312	6.214583105	−10.555338710	4.33×10−13	3.50×10−11
FAP	1.561195027	6.500446054	−9.026036586	3.63×10−11	1.48×10−9
CLDN3	1.625063261	9.253029938	−7.959600860	9.43×10−10	2.41×10−8
CP	1.533166764	7.200364641	−7.654471654	2.45×10−9	5.46×10−8
ZWINT	1.534376695	6.188206312	−7.484374926	4.18×10−9	8.53×10−8
MDK	1.627179122	9.344798091	−7.324035903	6.94×10−9	1.32×10−7
MMP9	1.708966359	7.597559232	−6.771606537	4.04×10−8	5.82×10−7
COL3A1	1.643158431	9.138801925	−6.532240647	8.72×10−8	1.11×10−6
COL11A1	2.387669841	5.767529739	−6.412529728	1.28×10−7	1.54×10−6
TOX3	1.770307511	4.159327973	−6.212832065	2.44×10−7	2.68×10−6
GABBR1 UBD	1.619335627	7.242024866	−6.016224500	4.60×10−7	4.64×10−6
MMP12	1.644563956	5.138055675	−5.757801490	1.06×10−6	9.61×10−6
POSTN	1.502940090	7.019690170	−5.698113750	1.28×10−6	1.13×10−5
SPINK1	2.517220534	7.945875667	−5.613361864	1.69×10−6	1.44×10−5
NQO1	1.527576974	7.820504924	−5.374097053	3.64×10−6	2.77×10−5
SLC7A5	1.517121297	7.685388227	−5.301930776	4.59×10−6	3.37×10−5
MMP7	2.087445737	6.987333841	−5.300686125	4.61×10−6	3.38×10−5
MMP1	1.553980336	4.289661258	−5.235072252	5.68×10−6	4.05×10−5
DSP	1.669305242	7.306405073	−5.162024033	7.18×10−6	4.94×10−5
CEACAM5	2.090196480	7.403316921	−4.560737891	4.80×10−5	2.43×10−4
S100P	1.513031206	7.826887416	−3.965427875	2.97×10−4	1.12×10−3

t, (the sample mean-the expected value for the population mean)/the standard error of the mean. MMP, matrix metalloproteinase; ID, identity; RAMP1, receptor activity modifying protein 1; FAP, familial adenomatous polyposis; CLND3, claudin 3; CP, ceruloplasmin; ZWINT, ZW10 interacting kinetochore protein; MDK, midkine; COL3A1, collagen type III α 1; COL11A1, collagen type XI α I; TOX3, TOX high mobility group box family member 3; UBD, ubiquitin D; POSTN, periostin; SPINK1, serine peptidase inhibitor, Kazal type 1; NQO1, NAD(P)H quinone oxidoreductase 1; SLC7A5, solute carrier family 7 member 5; DSP, desmoplakin; CEACAM5, carcinoembryonic antigen-related cell adhesion molecule 5; S100P, S100 calcium-binding protein P.

Figure 2.

Expression of MMPs (MMP1, MMP7, MMP9 and MMP12) in lung AC tissues and adjacent NTs evaluated by western blotting. AC, adenocarcinoma; NT, normal tissue; MMP, matrix metalloproteinase.

Functional enrichment analysis was applied for the above group of DEGs with the DAVID online tool. As shown in Table II, collagen metabolism (P=4.37×10−6) and multicellular organism macromolecule metabolic processes (P=5.98×10−6) were enriched for these DEGs.

Table II.

GO terms over-represented for the MMPs-associated genes.

		Gene
Category	Term	N	%	P-value	Genes	False discovery rate
GO_BP	GO:0032963 collagen metabolic process	4	20	4.37×10−6	MMP9, COL3A1, MMP7 and MMP1	0.005547327
GO_BP	GO:0044259 multicellular organism macromolecule metabolic process	4	20	5.98×10−6	MMP9, COL3A1, MMP7 and MMP1	0.007594962
GO_CC	GO:0005578 proteinaceous ECM	7	35	2.30×10−6	MMP9, COL3A1, MMP7, POSTN, COL11A1, MMP12 and MMP1	0.002198433
GO_CC	GO:0031012 ECM	7	35	3.56×10−6	MMP9, COL3A1, MMP7, POSTN, COL11A1, MMP12 and MMP1	0.003399869

Go, Gene Ontology; BP, biological process; CC, cellular component; ECM, extracellular matrix; MMP, matrix metalloproteinase; COL3A1; collagen, type III α 1; POSTN, periostin; COL11A1, collagen type XI α I.

Disease-relevant genes selection

Relevant diseases to the MMPs-associated DEGs were retrieved with DAVID. Under the Genetic Association DB Disease Class, MMP9, NAD(P)H quinone oxidoreductase 1 (NQO1), MMP12 and MMP1 were linked with lung AC, while MMP9, MMP12 and MMP1 were associated with lung function (Table III).

Table III.

Diseases associated with MMPs-associated genes.

Term	P-value	Genes
Nasopharyngeal cancer	8.72×10−6	MMP9, MMP7, NQO1, MMP12 and MMP1
Abdominal aortic aneurysm	1.18×10−4	MMP9, COL3A1, MMP12 and MMP1
Brain cancer	6.15×10−4	MMP9, MMP7, NQO1 and MMP1
Gastric ulcer	6.37×10−4	MMP9, MMP7 and MMP1
Subarachnoid hemorrhage	1.68×10−3	MMP9, MMP12 and MMP1
Lung function	3.48×10−3	MMP9, MMP12 and MMP1
Ovarian cancer	3.78×10−3	MMP9, MMP7, NQO1 and MMP1
Bladder cancer	5.85×10−3	MMP9, NQO1, MMP12 and MMP1
Colorectal cancer	6.11×10−3	MMP9, MMP7, NQO1, MDK and MMP1
Coronary artery luminal dimensions	7.91×10−3	MMP7 and MMP12
Osseointegrated implant failure	7.91×10−3	MMP9 and MMP1
Alzheimer's disease dementia, vascular	1.18×10−2	MMP9 and MMP1
Chronic obstructive pulmonary disease	1.84×10−2	MMP9, MMP12 and MMP1
Rheumatoid arthritis	1.85×10−2	MMP9, MMP7, MMP12 and MMP1
Uterine leiomyoma	1.97×10−2	MMP9 and MMP1
Adenomyosis endometriosis	1.97×10−2	MMP9 and MMP7
Breast cancer	2.09×10−2	MMP9, POSTN, NQO1, MMP12 and MMP1
Cervical artery dissection, spontaneous	2.36×10−2	MMP9 and COL3A1
Aneurysm	2.36×10−2	MMP9 and MMP12
Lung cancer	2.51×10−2	MMP9, NQO1, MMP12 and MMP1
Coronary artery disease	3.07×10−2	MMP9, COL3A1 and MMP12
H. pylori infection stomach cancer	3.52×10−2	MMP9 and MMP7
Left ventricular remodeling	3.90×10−2	MMP9 and MMP1

MMP, matrix metalloproteinase; MDK, midkine; COL3A1, collagen type III α 1; POSTN, periostin; NQO1, NAD(P)H quinone oxidoreductase 1.

PPI network analysis

A PPI network was constructed for the MMPs-associated DEGs using STRING (Fig. 3). The results revealed that interaction associations existed among MMP9, NQO1 and MMP7. In addition, the other genes clustered into the MMPs group also interacted with each other.

Figure 3.

Protein-protein interaction network for MMPs-associated genes. MPPs are marked in circles and other genes in squares. Lung cancer-associated genes are presented as white and other genes as yellow. MMP, matrix metalloproteinase; GABBR1, γ-aminobutyric acid type B receptor subunit 1; UBD, ubiquitin D; FAP, familial adenomatous polyposis; COL11A1, collagen type XI α I; COL3A1, collagen type III α 1; POSTN, periostin.

Transcriptional regulatory network analysis indicated that 9 genes, including MMP1, MMP9 and NQO1, were commonly regulated by the CCAAT/enhancer binding protein (C/EBP) α (CEBPA) transcription factor (Fig. 4).

Figure 4.

Transcriptional regulatory network containing MMPs. The transcription factor CEBPA is marked as white, while lung cancer-associated genes are presented as green and other genes as yellow. The dashed arrows indicate regulatory associations. MMP, matrix metalloproteinase; FAP, familial adenomatous polyposis; NQO1, NAD(P)H quinone oxidoreductase 1; POSTN, periostin; COL11A1, collagen type XI α I; COL3A1, collagen type III α 1; TOX3, TOX high mobility group box family member 3; CLDN3, claudin 3; CEBPA, CCAAT/enhancer binding protein (C/EBP) α.

Discussion

In the present study, cluster analysis demonstrated that 4 MMPs, including MMP1, MMP7, MMP9 and MMP12, were upregulated in lung AC samples, and may be involved in ECM metabolic processes, thus promoting cancer invasion. These results were validated by western blotting and were in accordance with previous studies (29–32). MMP1 is the most highly expressed interstitial collagenase for degrading fibrillar collagens (33). Overexpression of MMP1 has been associated with tumor invasion and metastasis by modulating the polarization of T helper (Th)1/Th2 inflammatory responses (29). Downregulation of MMP-7 mediated by antisense oligonucleotide changes the ultrastructure of lung AC A549 cells, leading to decreased microvilli, endoplasmic reticulum dilation, swelling of mitochondria and formation of apoptotic bodies, which eventually inhibits invasion in lung AC A549 cells (30). Similarly, the expression levels of MMP9 and MMP12 were observed to be higher in NSCLC than in normal samples (31). Upregulation of MMP12 and MMP9 may be one of the mechanisms to promote lung cancer cell invasion (32).

In addition, the present study identified several DEGs that were clustered into one group with MMPs, thus indicating the same function in cancer cell invasion. A number of these genes have been linked to lung cancer in the following studies. NQO1 is a member of the NQO family, and altered expression of this protein has been reported in numerous tumors (34). Mutations in this gene contribute to susceptibility to various forms of cancer, including lung cancer (35,36). S100 calcium-binding protein P (S100P) is a member of the S100 family of proteins containing two EF-hand calcium-binding motifs, which are involved in the regulation of cell cycle progression and differentiation (37). Bartling et al (31) reported that S100P expression is mainly increased in AC, and that S100P upregulation is detected in early rather than in advanced tumor stages. Overexpression of S100P may lead to cancer cell invasion by changing the expression levels of several cytoskeletal proteins (38). Serine peptidase inhibitor, Kazal type 1 (SPINK1) is a secreted serine protease inhibitor (39). Overexpression of SPINK1 is associated with aggressiveness of prostate cancer (40). A previous study by Soon et al (34) revealed that SPINK1 is an invasion factor associated with prognosis in breast cancer patients. However, its role in lung cancer remains unknown. The present study suggested that SPINK1 may also be involved in lung AC invasion. Ubiquitin D (UBD, also known as F-associated transcript 10) is a member of the ubiquitin-like modifier family, and appears to be upregulated in hepatocellular carcinoma, as well as in gastrointestinal and gynecological cancers (41). Increased cytoplasmic UBD is significantly associated with depth of colon cancer invasion (42). Few studies have investigated the roles of UBD in lung AC invasion, including the present study. Collagen type XI α I (COL11A1) has also been demonstrated to be overexpressed in NSCLC, and has been identified as a potential invasion-associated gene in cancer (43). However, its mechanism of function is not understood yet. The present study revealed that COL11A1 could interact with the osteoblast-specific factor periostin (POSTN), although this requires further experimental validation. In addition, the associations between MMPs and various genes (including MMP9-CLDN3, MMP9-DSP and MMP12-DSP) were also further confirmed by Pearson correlation analysis. CLDN3, a component of tight junctions, has been reported to be upregulated in NSCLC, and may be important in invasion (44). Agarwal et al (39) reported that CLDN3-mediated increased invasion may be accomplished through the activation of MMP2. In the present study, CLDN3 was significantly correlated with MMP9 via Pearson correlation analysis, indicating a potential link between CLDN3 and MMP9 in lung AC cell invasion. DSP acts as a tumor suppressor via inhibiting the Wnt/β-catenin signaling pathway in NSCLC (45). Thus, the expression patterns of DSP and MMPs are theoretically opposite, which has been demonstrated in the human epithelial carcinoma cell line A431 during the epithelial-mesenchymal transition, with upregulated MMPs and downregulated DSP levels (46).

According to the present transcriptional regulatory network analysis, MMPs-associated genes were regulated by CEBPA. CEBPA is a basic leucine zipper domain transcription factor that not only can bind to certain promoters and enhancers as a homodimer, but can also form heterodimers with the associated proteins CEBPB and CEBPG (47). It has been reported that CEBPB is an important mediator in the activation of MMP genes (including MMP1, MMP3 and MMP10) in A549 lung carcinoma cells stimulated with the inflammatory cytokine interleukin-1β (48). Therefore, it could be speculated that CEBPA may also serve an important role in upregulating the expression of MMP1 and MPP9, thus affecting the invasion of lung cancer. Although the present results indicated that CEBPA could modulate the expression of fibroblast activation protein α, CLDN3, collagen type III α 1, COL11A1, TOX high mobility group box family member 3, NQO1 and POSTN, no experimental evidence could be obtained.

In conclusion, the present study identified several genes such as SPINK1 and UBD that may serve important roles in lung AC invasion with MMPs (including MMP1, MMP7, MMP9 and MMP12). Several MMPs-associated genes were observed to regulated by the CEBPA transcription factor. These findings may provide various underlying targets for prevention of lung AC invasion. However, further experimental investigations or studies on other datasets are required to validate the present observations.