Identification of potential biomarkers of lung adenocarcinoma brain metastases via microarray analysis of cDNA expression profiles.

Brain metastases originating from lung adenocarcinoma (LAD) occur frequently. The aim of the current study was to assess potential biomarkers for the prognosis of lung adenocarcinoma brain metastasis (LAD-BM) through the analysis of gene expression microarrays. The current study downloaded two gene expression datasets, GSE14108 and GSE10245, from the Gene Expression Omnibus database. From GSE14108 and GSE10245, 19 LAD-BM samples and 40 primary LAD samples were selected for analysis. To identify the differentially expressed genes (DEGs), the current study compared the two sample groups, using the limma R package. Subsequently, pathway enrichment analysis was conducted using the Cluster Profiler R package, and the construction of the protein-protein interaction (PPI) network was executed utilizing the Search Tool for the Retrieval of Interacting Genes database. The microRNA-target network was built using the TargetScore R package. Then, these networks were established and visualized using Cytoscape software. An array of 463 DEGs was identified in the LAD-BM samples, including 256 upregulated and 207 downregulated genes. Based on functional term enrichment analysis using the Gene Ontology database and signaling pathway enrichment analysis using the Kyoto Encyclopedia of Genes and Genomes database, it was identified that the overlapping DEGs were primarily involved in chemokine-associated signal transduction, which may mediate lung cancer cell metastasis to the brain. Chemokine ligand 2, lysozyme, matrix metalloproteinase-2 (MMP-2), lysyl oxidase (LOX) and granzyme B were identified as potential biomarkers according to a topological analysis of the PPI networks. Two notable nodes, MMP-2 and LOX, appeared in the PPI network and were key points in the microRNA-target network, as they were regulated by hsa-let-7d. Many DEGs and microRNAs were regarded as prognostic biomarkers for lung adenocarcinoma metastasis in the current study. These DEGs were primarily associated with chemokine-mediated signaling pathways. In addition, MMP-2 and LOX were predicted to be targets of hsa-let-7d.

Introduction

Brain metastasis (BM) is a common complication of malignant cancers in adulthood (1) and serves as a major cause of cancer-associated cases of mortality (2,3). Lung cancer, particularly non-small cell lung cancer (NSCLC) (4), is the most common primary cancer that metastasizes to the brain. The percent of patients with NSCLC that suffer from BM is ~10% at the onset of illness (5), with another 25–40% demonstrating advancing BM during the disease course (6). Although the precise occurrence is unknown, it is predicted that 35–50% of all patients with NSCLC are afflicted with BM, based on various studies (4,7). Lung adenocarcinoma (LAD) is a common NSCLC with a frequent incidence of developing brain metastases and accounts for over 50% of all NSCLC brain metastases (1). Once lung cancer cells have spread to the brain, the patient has a very poor prognosis. The median overall survival of untreated and treated brain metastases patients is 4–11 weeks and 4–15 months, respectively (8,9). Prior investigations have established that there are two mechanisms underlying the involvement of LAD metastasis to the brain, including activation of the Wnt/T-cell factor signal pathway (with homeobox B9 and lymphoid enhancer binding factor 1 expression markedly associated with metastatic ability) and the upregulation of plasminogen activator (PA) inhibitory serpins (10,11). Recently, Singh et al (12) revealed that Sparc/osteonectin, cwcv, and kazal-like domains proteoglycan 1 and twist family bHLH transcription factor 2 are necessary moderators of brain metastasis-initiating cells in LAD and play a key role in the metastatic process in the brain.

MicroRNAs (miRNAs or miRs) are non-coding single stranded nucleotides between 19–24 nucleotides in length (13), which may be important factors in the regulation of tumor invasion and metastasis, including that of LAD, by regulating gene expression (14). Previous studies have identified that miR-95, miR-378 and miR-145 may regulate certain key pathways in the process of lung cancer cell metastasis to the brain (15–17). However, there are few reports of miRNAs associated with LAD metastasis.

A previous study has identified differentially expressed genes (DEGs) between primary LAD and metastatic brain tumors by cDNA microarray (18). The results of this previous study suggest that those DEGs, including genes coding for cytoskeletal proteins, cellular antigens and plasma membrane proteins, may serve important roles in cell-cell interactions. However, tumor cell metastasis to the brain is a complicated multi-step process and other important molecules that may be involved, and how they regulate the mechanisms of metastasis, remain largely unknown. Therefore, the current study analyzed the GSE14108 and GSE10245 microarray datasets to identify DEGs in primary LADs and BM. Furthermore, the current study performed enrichment analysis and protein-protein interaction (PPI) network construction. Additionally, the hub genes were combined with the TargetScan database to construct a miRNA-target regulatory network. These comprehensive bioinformatics methods provided an opportunity to identify effective biomarkers for the prognosis of LAD metastasis and to better illuminate the underlying molecular pathogenesis of BM.

Materials and methods

Data resources

The raw data from the GSE14108 and GSE10245 microarray datasets were downloaded from the Gene Expression Omnibus (GEO) (www.ncbi.nlm.nih.gov/geo/) database. From the datasets, 19 LAD-BM samples (GSE14108) (19) and 40 LAD samples (GSE10245) (20) were used in the current study. The platform used for the detection of these microarray data was the GPL570 Human Genome U133 Plus 2.0 Array (Affymetrix; Thermo Fisher Scientific, Inc., Waltham, MA, USA).

Pre-treatment and DEG analysis

Raw microarray data were pre-processed using the Affy package (version 1.48.0; http://bioconductor.org/packages/release/bioc/html/affy.html); robust multi-array average (RMA) was used for background adjustment and the linear scaling method was applied for normalization. Then, the standardized data were summarized based on the Perfect Match-Mismatch difference model (21). The intensity of gene expression was calculated from the hybridization signal of the probe-set (containing multiple probes). When multiple probes corresponded to the same gene, the average value of probe expression was considered the gene expression value. Using the limma package of R software (version 3.26.9; http://bioconductor.org/packages/release/bioc/html/limma.html), DEGs were identified by comparing the LAD-BM samples with the LAD samples (22). All DEGs were selected according to the following criteria: |log2 (fold-change)| ≥1.5 and a false discovery rate of adjusted P<0.001.

Functional enrichment analysis of DEGs

To identify the potential biological processes that were affected, the Bioconductor package ‘Cluster Profiler’ of R software (version 3.2.2; http://www.bioconductor.org/packages/release/bioc/html/clusterProfiler.html) (23,24) was used to classify the enriched Gene Ontology (GO) terms. Information in the Kyoto Encyclopedia of Genes and Genomes (KEGG) (http://www.genome.jp/kegg/pathway.html) database was used for the pathway enrichment analysis of DEGs (25). P<0.05 was considered to indicate a statistically significant selection of GO terms and KEGG pathways.

PPI network analysis of the DEGs

The PPI networks of the DEGs were constructed utilizing the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING; string-db.org) database (26) and visualized using Cytoscape software (version 3.6.1; http://www.cytoscape.org/). From 787,896 pairs of human protein interactions containing 16,730 genes, DEG-containing interactions were obtained. STRING (https://string-db.org) utilized a combined score (0–1) (27) to assess reliability. Spearman's correlation coefficient was implemented to assess the edge scores, which evaluated the probability of two co-expressed gene pairs in the current study. Each protein was regarded as a node in the network, and the degree of a node was regarded as the number of interactions with other nodes. Hub genes were nodes with ≥50 degrees.

Establishment of the miRNA-target regulatory network

The TargetScore (version 1.12.0) R package (28) was used to predict biological miRNA-gene interactions based on the TargetScan database (http://www.targetscan.org). The current study combined ‘targetscan context+score’ with ‘probabilities of conserved targeting’ of the hub gene to calculate the gene ‘targetscore’ value, utilizing a variational Bayesian-Gaussian mixture model (28). The miRNAs were identified as miRNA-gene interactions according to the criterion ‘targetscore’ >0.4. The miRNA-gene regulatory network was subsequently visualized using Cytoscape software (version 3.6.1).

Results

Identification of DEGs in LAD-BM

Overall, the current study identified 463 DEGs (256 upregulated and 207 downregulated) between the LAD-BM and LAD samples. A volcano plot (Fig. 1) demonstrated the distribution of DEGs, and a heat map is presented in Fig. 2.

Figure 1.

Volcano plot of DEGs. The blue dots represent 207 downregulated DEGs, the red dots represent 256 upregulated DEGs, and the green dots represent non-DEGs. DEG, differentially expressed gene.

Figure 2.

Heat map of gene expression in brain metastases and lung adenocarcinoma samples. The expression level of each gene was normalized and the relative value to the median among 59 samples is demonstrated by color. Red represents high expression and green indicates low expression.

GO terms and KEGG pathways of DEGs

In total, 65 GO terms and 23 KEGG pathways were enriched in DEGs according to the criteria P<0.05. The 20 most significantly enriched GO terms are presented in Table I. The majority of the enriched GO terms contained ‘chemokine-mediated signaling pathway’, ‘collagen catabolic process’ and ‘blood vessel morphogenesis’. Certain enriched GO terms also included ‘angiogenesis’, ‘extracellular matrix organization’ and ‘extracellular structure organization’. The enriched KEGG pathways of the DEGs are listed in Table II. The majority of the enriched KEGG pathways, including ‘cytokine-cytokine receptor interaction’ and ‘chemokine signaling pathway’ were directly involved in the process of lung cancer metastasis to the brain.

Table I.

Top 20 most significantly enriched GO terms for biological processes of differentially expressed genes in brain metastases samples compared with lung adenocarcinoma.

GO ID	GO name	DEG number of genes involved	P-value
GO:0070098	Chemokine-mediated Signaling pathway	13	3.65×10−9
GO:0030574	Collagen catabolic process	12	2.77×10−8
GO:0048514	Blood vessel morphogenesis	30	2.80×10−8
GO:0006959	Humoral immune response	18	6.75×10−8
GO:0044236	Multicellular organismal metabolic process	15	6.94×10−8
GO:0044243	Multicellular organismal catabolic process	12	7.55×10−8
GO:0002685	Regulation of leukocyte migration	15	1.12×10−7
GO:0001525	Angiogenesis	26	1.35×10−7
GO:0030198	Extracellular matrix organization	25	1.45×10−7
GO:0043062	Extracellular structure organization	25	1.52×10−7
GO:0070661	Leukocyte proliferation	20	3.00×10−7
GO:0002688	Regulation of leukocyte chemotaxis	12	3.02×10−7
GO:0030595	Leukocyte chemotaxis	16	3.15×10−7
GO:0050920	Regulation of chemotaxis	15	3.27×10−7
GO:0060326	Cell chemotaxis	18	4.57×10−7
GO:0032963	Collagen metabolic process	13	4.88×10−7
GO:0048247	Lymphocyte chemotaxis	9	5.20×10−7
GO:0050795	Regulation of behavior	17	6.23×10−7
GO:0044259	multicellular organismal macromolecule metabolic process	13	7.89×10−7
GO:0002548	Monocyte chemotaxis	9	1.39×10−6

GO, gene ontology; ID, identifier; DEG, differentially expressed gene.

Table II.

The 10 most significantly enriched signaling pathways of differentially expressed genes.

KEGG pathway no.	Signaling pathway	DEGs involved	DEG number of genes involved	P-value
hsa04060	Cytokine-cytokine receptor interaction	PDGFRA, CCL2, IL7R, CCL18, CXCR4, HGF, CXCL12 and others	23	1.99×1008
hsa04145	Phagosome	CTSS, HLA-DRA, FCGR2B, HLA-DPA1, HLA-DQB1, NOS1, COMP and others	13	3.13×1005
hsa05323	Rheumatoid arthritis	IL17A, CCL2, TNFSF13B, HLA-DRA, HLA-DPA1, CXCL12, MMP1 and others	10	2.55×1005
hsa05150	Staphylococcus aureus infection	C1R, C1S, FCGR2A, FPR3, HLA-DRA, FCGR2B, HLA-DPA1, HLA-DQB1	8	2.69×1005
hsa04672	Intestinal immune network for IgA production	CXCR4, TNFSF13B, HLA-DRA, HLA-DPA1, CXCL12, TNFRSF17, HLA-DQB1	7	7.61×1005
hsa05320	Autoimmune thyroid disease	IFNA14, IFNA7, HLA-DRA, HLA-DPA1, HLA-DQB1, GZMB	6	0.001016
hsa05164	Influenza A	CCL2, SOCS3, CASP1, HLA-DRA, HLA-DPA1, CCL5, TNFSF10, HLA-DQB1 and others	11	0.001679
hsa04062	Chemokine signaling pathway	CCL2, GNG2, CCL18, CXCR4, CXCL12, CCL5, CXCL9, CCL19, CCL8, CXCL13, CXCL11	11	0.002944
hsa04940	Type I diabetes mellitus	GAD2, HLA-DRA, HLA-DPA1, HLA-DQB1, GZMB	5	0.002649
hsa05152	Tuberculosis	IFNA14, IFNA7, FCGR2A, CTSS, HLA-DRA, FCGR2B, HLA-DPA1 and others	10	0.005802

KEGG, Kyoto Encyclopedia of Genes and Genomes; DEG, differentially expressed gene; IgA, Immunoglobulin A.

PPI network of DEGs

Using the STRING database, a PPI network was constructed and is presented in Fig. 3. It included 942 pairs and 317 nodes. A total of five hub genes were identified, including CCL2 (degree=60), LYZ (degree=60), MMP-2, (degree=58), LOX (degree=53) and GZMB (degree=50).

Figure 3.

Protein-protein interaction network of differentially expressed genes. Red nodes represent the hub genes (≥50 degrees).

miRNA-target gene regulatory network

The current study focused on the five hub genes, and further assessed their miRNA-target associations. The miRNA-target regulatory network was based on interactions in the TargetScan database. A miRNA-gene-regulated network was established, including 23 connections and five hub genes, which are presented in Fig. 4. MMP-2 and LOX were significant targets that were identified by and interacted with many miRNAs. MMP-2 was predicted to be the target of the following six miRNAs: hsa-miR-17, hsa-miR-449b, hsa-miR-135a, hsa-miR-530f, hsa-miR-1285 and hsa-let-7d. LOX was predicted to be the target of the following 12 miRNAs: hsa-miR-181a, hsa-miR-155, hsa-miR-26a, hsa-miR-148b, hsa-miR-530f, hsa-miR-182, hsa-miR-200c, hsa-miR-429, hsa-miR-200b, hsa-let-7d, hsa-miR-30b and hsa-miR-27a.

Figure 4.

MicroRNA-gene regulatory network of five hub genes in the protein-protein interaction network.

Discussion

In the current study, the significantly enriched GO terms of DEGs contained ‘chemokine-mediated signaling pathway’, ‘collagen catabolic process’ and ‘blood vessel morphogenesis’. As demonstrated by Table II, the enriched KEGG pathways were tumor-associated biological processes, including ‘cytokine-cytokine receptor interaction’ and ‘chemokine signaling pathway’. The overlapping DEGs, including CCL2, CXCL12 and CXCR4, which were mainly enriched in the GO terms and KEGG pathways, were primarily involved in chemokine-associated signal transduction. Numerous studies have previously demonstrated that chemokines serve a pivotal role in cancer invasion, angiogenesis and metastases (29–31). Emerging studies have revealed that the CXCR4/CXCL12 signaling axis is involved in tumor cell migration and angiogenesis, facilitating the formation of BM (32–34). Hartmann et al (35) identified that CXCR4 cooperates with integrins to mediate adhesion and chemo-resistance in small cell lung cancer cells. CXCR4 and CXCL12, which are overexpressed in brain metastases, have been identified to be correlated with brain-specific metastasis and poor prognosis in NSCLC patients (36). Furthermore, the proteolysis of interstitial collagen has been recognized as a contributing factor that participates in tumor cell invasion and metastasis (37). Collagen is one of the chief structural proteins of the extracellular matrix (ECM) and provides a major obstacle impeding cancer cell migration. Thus, it has been speculated that collagen catabolism serves an important role in facilitating the spread and invasion of cancer cells to host organs (38). As demonstrated in Table I, the second most enriched GO term was ‘collagen catabolic process’. Indeed, collagen degradation can be observed histologically at the periphery of certain aggressive tumors (39). Thus, the remodeling of the ECM, which is caused by the degradation of collagen, serves an important role in tumor invasion and migration.

The PPI networks based on topology analysis revealed five potential key genes, including CCL2, LYZ, MMP2, LOX and GZMB, which may play important roles in lung cancer BM. CCL2 encodes the chemokine (C-C motif) ligand 2, which interplays with its receptor, C-C chemokine receptor type 2 (CCR2), and subsequently promotes tumor progression and metastasis caused by the Snail-induced epithelial-to-mesenchymal transition (EMT) (40). Previous research has reported that the release of CCL2 is induced by activating protease-activated receptor 2 through matrix metalloproteinase-1 expression, which initiates the activation of the thromboxane A2 receptor in LAD cells (41). Furthermore, activation of the thromboxane A2 receptor also increases the expression of CCL2 and then recruits macrophages, thereby stimulating human LAD cell invasion (42). In addition, the CCL2/CCR2 axis cooperates with interleukin-6 to enhance EMT by activating signal transducer and activator of transcription 3-Twist signaling, and then boosts lung cancer progression and metastasis (43).

The matrix metalloprotease (MMP) family includes a series of functionally similar enzyme molecules that degrade protein substrates based on a highly conserved mechanism (44). MMP-2, also known as gelatinase A, has the ability to degrade matrix proteins including gelatin, type IV and type V collagens and elastin (45). MMP-2 is involved in tumor metastasis as it degrades vascular basement membranes (46). Recently, a study demonstrated that the epigenetic activation of MMP-2 induced by the interaction of megakaryocytic leukemia 1 with histone methyltransferase SET1 can also facilitate the migration and invasion of ovarian tumor cells (47).

The LOX gene encodes lysyloxidase, which is a secretory copper-dependent amine oxidase (48). Previous studies have identified that the overexpression of LOX can also promote tumor cell proliferation, invasion and metastasis in various cancer types (49–51). The downregulation of LOX markedly enhances the expression of E-cadherin and decreases the expression of vimentin (49), increasing the tendency of tumor cells to metastasize. Furthermore, Wilgus et al (52) confirmed that the high expression of LOX is associated with invasion and poor prognosis in patients with LAD.

The most well-known function of lysozyme, encoded by the LYZ gene, is anti-infection (53). Previous research has revealed that the expression of LYZ is associated with an unfavorable prognosis and may be a potential prognostic factor in male breast cancer (54). However, the relationship between LYZ and lung cancer is still not clear.

Additionally, granzyme B is a cytotoxic T-lymphocyte-associated serine esterase. The inhibition of granzyme B and interferon-γ activates the transforming growth factor-β/Smad pathway and then inhibits T-cell mediated cancer clearance in vivo, emphasizing the role of the perforin/granzyme pathway in cancer clearance (55). Interestingly, in specific settings, granzyme B and perforin are associated with the regulatory T cell-mediated inhibition of cancer clearance in vivo (56). However, Smyth et al (57) discovered that granzyme A and B are not required for cytotoxic T cell- and natural killer cell-induced cancer rejection, including spontaneous and experimental cancers. Therefore, more convincing studies are required to confirm the function of granzymes in cancer immune surveillance and rejection.

A miRNA-target regulatory network was constructed, and many potential miRNAs were identified. MMP-2 and LOX were the leading targets identified by and interacting with multiple miRNAs. MMP-2 and LOX were predicted to be targets of hsa-let-7d and hsa-miR-530f. So far, the function of hsa-miR-530f is unclear. Hsa-let-7d reportedly inhibits cancer pathogenesis (58). Furthermore, let-7d may inhibit cancer cell migration, invasion and metastasis, by directly targeting PBX3, COL3A1 and CCL7 (58,59). Meanwhile, a recent study reported that let-7d may suppress proliferation and invasion of trophoblast cells, by targeting MMP-2 (60). However, to the best of our knowledge, the targeting of MMP-2 by hsa-let-7d has not yet been reported in lung cancer. In the current study, MMP-2 and LOX were predicted to be targets of hsa-let-7d, implying that hsa-let-7d may target the two genes during lung cancer metastasis.

The current study also had specific limitations. Certainly, the current results may be more meaningful if matched cases were used to validate the expression of selected genes. Clinically, the most common primary origin of BM is the lung (61), and LAD is usually characterized by the early development of BM, even though certain primary tumors may be observed without symptoms (62). Unfortunately, the current study did not have the opportunity to obtain enough BM and matched primary tumors to analyze the expression of selected genes in tumor patients. Studies aimed at identifying candidate genes that lead to tumor metastasis are difficult to perform, particularly using human tumor samples. This is due to the unpredictable time point of tumor BM and the difficulty in obtaining tumor tissue samples from patients who have never received chemotherapy or radiotherapy, which is required to ensure that the analysis and interpretation of the results is not confused with drug or radiation-induced changes in gene expression (63). Furthermore, for many years, patients with BM have not generally been considered candidates for surgical intervention, and neurosurgeons have been reluctant to surgically treat these patients due to the wide dissemination of cancer and limited survival (64). Prospective or multi-center cooperation may be necessary to collect a well-defined BM group and matched primary tumor cases for further comparative analysis of candidate genes (65). Comprehensive bioinformatics analyses of gene expression profiles based on public databases (e.g. GEO and TCGA) have yielded important insights into potential prognostic biomarkers for tumors (66–68). The current study re-analyzed microarray data from 19 LAD-BM samples and 40 primary LAD samples in the GEO database. All DEGs in LAD-BM compared with the control group were identified via a bioinformatics-based method. Furthermore, the current study performed GO term and pathway enrichment analyses, and PPI network construction. The current study also combined the DEG data with information on miRNAs in the TargetScan database to predict miRNA-target interactions. These analyses aided the identification of key genes associated with BM, including CCL2, LYZ, MMP2, LOX and GZMB, and the essential miRNA, hsa-let-7d. Through these comprehensive bioinformatical methods, the current study may contribute to understanding the molecular mechanism underlying BM, facilitating the identification of potential gene targets for the diagnosis and treatment of patients with BM.

In conclusion, many DEGs and miRNAs were regarded as potential biomarkers for the prognosis of LAD metastasis in the present study. These DEGs were mainly associated with chemokine-mediated signaling pathways. In addition, MMP-2 and LOX were predicted to be targets of hsa-let-7d. However, these predictive results require additional experiments to confirm their function.