Bioinformatics analysis to identify the critical genes, microRNAs and long noncoding RNAs in melanoma.

Melanoma, which is usually induced by ultraviolet light exposure and the following DNA damage, is the most dangerous skin cancer. The purpose of the present study was to screen key molecules involved in melanoma.Microarray data of E-MTAB-1862 were downloaded from the ArrayExpress database, which included 21 primary melanoma samples and 11 benign nevus samples. In addition, the RNASeq version 2 and microRNA (miRNA) sequencing data of cutaneous melanoma were downloaded from The Cancer Genome Atlas database. After identifying the differentially expressed genes (DEGs) using Limma package, enrichment analysis and protein-protein interaction (PPI) network analysis were performed separately for them using DAVID software and Cytoscape software. In addition, survival analysis and regulatory network analysis were further performed by log-rank test and Cytoscape software, respectively. Moreover, real-time reverse transcription polymerase chain reaction (RT-PCR) was performed to further verify the expression patterns of several selected DEGs.A total of 382 DEGs were identified in primary melanoma samples, including 206 upregulated genes and 176 downregulated genes. Functional enrichment analysis showed that COL17A1 was enriched in epidermis development. In the PPI network, CXCL8 (degree = 29) and STAT1 (degree = 28) had higher degrees and could interact with each other. Survival analysis showed that 21 DEGs, 55 long noncoding RNAs (lncRNAs) and 32 miRNAs were found to be associated with prognosis. Furthermore, several regulatory relationships were found in the lncRNA-gene regulatory network (such as RP11-361L15.4 targeting COL17A1) and the miRNA-gene regulatory network (such as hsa-miR-375 targeting CCL27 and hsa-miR-375 targeting insulin-like growth factor 1 receptor [IGF1R]). Real-time RT-PCR results showed that the overall direction of differential expression was consistent except COL17A1.CXCL8 interacted with STAT1, CCL27, and IGF1R targeted by hsa-miR-375, and COL17A1 targeted by RP11-361L15.4 might function in the development and progression of melanoma, which should be verified by more detailed experiments.

Introduction

Melanoma is the most aggressive type of skin cancer arising from the malignant transformation of melanocytes.[ Malignant melanoma has a strong tendency to spread to other parts of the body and causes serious illness and death.[ In 2016, there was an estimation of 76,380 new diagnoses of cutaneous melanoma and 10,130 deaths related to melanoma in the United States.[ Despite significant breakthroughs and advances in early diagnosis and prevention as well as targeted therapies, the prognosis of melanoma remains unoptimistic. In addition, evidence shows that the 5-year survival rate decreases to 63% in patients with melanoma with regional metastases and it is only 16% in patients with distant metastases.[ Therefore, identifying molecular pathogenic basis for melanoma and novel treatment strategies remains an active area of research.

Recently, many researches have focused on the roles of mRNAs, microRNAs (miRNAs), and long noncoding RNAs (lncRNAs) in the tumorigenesis of melanoma. For instance, neural precursor cell expressed, developmentally downregulated 9 , which is significantly overexpressed in human metastatic melanoma, promotes invasion ability in vitro and metastasis ability in vivo of melanocytes.[ Recently, Zeng et al[ demonstrated that trans-activating response region DNA-binding protein (TDP-43) was a novel oncogene in melanoma and it can regulate melanoma proliferation and metastasis. MiR-211 inhibits expressions of the entire melastatin locus insulin-like growth factor type 2 receptor, transforming growth factor β receptor type II, and nuclear factor–activated T cell 5 , suggesting that miR-211 suppresses melanoma invasion in the progression of human melanoma.[ A study reported that miR-7-5p could inhibit melanoma cell proliferation and metastasis in part through its direct suppression of nuclear factor kappa B (NF-κB) subunit RelA expression and NF-κB signaling.[ The lncRNA SPRY4 intronic transcript 1 (SPRY4-IT1) is overexpressed in melanoma cells, accumulates in cell cytoplasm, and affects cell dynamics, indicating that SPRY4-IT1 may be related to the molecular mechanisms of human melanoma.[ The lncRNA BRAF–activated noncoding RNA is reported to overexpress in human melanoma, and be implicated in the proliferation of melanoma cell in vivo and in vitro through mediating the activation of mitogen-activated protein kinase (MAPK) pathway.[ Overexpressed lncRNA growth arrest-specific transcript 5 can reduce the invasion and migration of melanoma SK-Mel-110 cells through inhibiting the expression and activity of matrix metalloproteinase 2 to some degree.[ However, a number of melanoma-associated cellular molecules remain unexploited, and the precise molecular mechanisms of melanoma are not yet fully understood.

Bioinformatic data mining of gene expression microarray data provides a useful tool for revealing many previously unrecognized significant genes implicated in the pathogenesis of diseases.[ For instance, a recent study analyzed previously published melanoma gene expression profiles from available database and identified several lncRNAs and mRNAs that may be associated with melanoma tumorigenesis and metastasis.[ In 2016, Eriksson et al[ carried out analysis of gene expression profiles to identify the common changes in metastatic and nonmetastatic primary melanomas and benign nevi, finding that collagen triple helix repeat containing-1 regulated by B-Raf proto-oncogene, serine/threonine kinase (BRAF), nuclear factor of activated T cells-c2, and transforming growth factor β played an essential role in migration and invasion of melanoma cells and could be used as therapeutic target for metastatic melanomas. In order to identify more genes, miRNAs and lncRNAs associated with melanoma, the microarray data deposited by Eriksson et al, as well as the RNASeq version 2 (RNASeqV2) and miRNA sequencing (miRNASeq) data of cutaneous melanoma from The Cancer Genome Atlas (TCGA) database were downloaded and reanalyzed. We screened the differentially expressed genes (DEGs) between primary melanoma samples and benign nevus samples. Afterwards, enrichment analysis and protein-protein interaction (PPI) network analysis were conducted for the DEGs. In addition, survival analysis and regulatory network construction were further carried out to select the key molecules in melanoma.

Materials and methods

Data source

The flow diagram of bioinformatics analysis is shown in Figure 1. Microarray data under the accession number of E-MTAB-1862 were extracted from the European Bioinformatics Institute ArrayExpress database (http://www.ebi.ac.uk/arrayexpress), which were sequenced on the platform of A-AFFY-44—Affymetrix GeneChip Human Genome U133 Plus 2.0 (HG-U133_Plus_2). The E-MTAB-1862 dataset included 21 primary melanoma samples (8 females and 13 males, mean Breslow thickness = 6.8) and 11 benign nevus samples (6 females and 5 males) isolated from skin. In addition, the RNASeqV2 and miRNASeq data of cutaneous melanoma were downloaded from TCGA database. This study just reanalyzed the microarray data downloaded from public database and performed bioinformatics analysis. The authors declare that no experiments were performed on humans or animals for this investigation. Thus, ethics approval or consent to participate was not applicable.

Figure 1

The flow diagram of bioinformatics analysis. DEG = differentially expressed gene, lncRNA = long noncoding RNA, PPI = protein-protein interaction, RNASeqV2 = RNASeq version 2, TCGA = The Cancer Genome Atlas.

Data preprocessing and DEGs screening

After E-MTAB-1862 was downloaded, the raw data were converted to expression matrix using the Affy package[ in R. Combined with the annotation information from the platform, probes were converted into gene symbols. When multiple probes mapped to the same gene symbol, their mean value was considered as the gene expression value. Through aligning to the GENCODE lncRNA data,[ the RNASeqV2 data were further annotated.

The DEGs between primary melanoma samples and benign nevus samples were screened by the linear models for microarray data (Limma) package[ in R. The adjusted P value (all named false discovery rate, FDR) <.05 and |log2fold change (FC)| >1 were used as the cut-off criteria.

Functional and pathway enrichment analysis

Gene Ontology (GO, http://www.geneontology.org/) describes genes and gene products from molecular function (MF), cellular component (CC), and biological process (BP) aspects.[ Kyoto Encyclopedia of Genes and Genomes (http://www.genome.jp/kegg/), which integrates information of systemic functional, genomic, and chemical aspects, is known for pathway mapping.[ The DEGs were performed with functional and pathway enrichment analyses using DAVID software (http://david.abcc.ncifcrf.gov/),[ with the FDR <0.05.

PPI network construction

The Search Tool for the Retrieval of Interacting Genes (STRING, http://string-db.org/) database contains a large amount of known and predicted interaction information.[ Cytoscape software (http://www.cytoscape.org/) is developed for constructing a unified framework via integrating biomolecular interaction networks and gene expression data.[ The interactions among the proteins encoded by the DEGs were predicted by the STRING database,[ and then the PPI network was visualized using the Cytoscape software.[ The proteins involved in more interactions were key nodes in the network.

Survival analysis

The median value (each variable value in the statistical population was sorted in the order of size to form a series, and the variable value in middle position was the median value) of each lncRNA, miRNA, and mRNA in tumor samples was calculated. Then, the samples with expression higher than the median value were divided into high expression group, and the samples with expression lower than the median value were divided into low expression group. Using the KMsurv package[ and the survival package[ in R, the Kaplan-Meier curve was constructed for the high and low expression groups. Afterwards, the significant difference in prognosis between the 2 groups was obtained by log-rank test,[ with P value <.05 as the threshold.

Predicting the target genes of the risk lncRNAs and miRNAs

The target genes of lncRNAs were predicted based on coexpression relationships, and the genes with FDR <0.05 were obtained as the potential targets of the lncRNAs. What's more, miRanda (http://www.mocrorna.org),[ MirTarget2 (http://mirdb.org/miRDB),[ PicTar (http://pictar.mdc-berlin.de/),[ PITA (http://genie.weizmann.ac.il/pubs/mir07/mir07_data.html),[ TargetScan (http://www.targetscan.org/),[ miRecords (http://www.mirecords.umn.edu),[ and Mirwalk (http://mirwalk.uni-hd.de/)[ databases were used to search the target genes of prognostic risk miRNAs. The genes being found in 1 or more databases were considered as the targets of prognostic risk miRNAs. Finally, the regulatory network was constructed using Cytoscape software.[

Quantitative real-time polymerase chain reaction validation

In an attempt to validate differential expression of genes obtained by analyzing the microarray data, we selected a total of 5 genes, 2 upregulated and 3 downregulated, and quantified their expression levels by real-time polymerase chain reaction (PCR). Real-time PCR was performed on 5 primary melanoma samples and 10 benign nevus samples. Total RNA from the sample tissue was isolated using Trizol (TAKARA, Japan) protocol. First-strand cDNA was generated from 0.5 μg total RNA using PrimerScript RT Master Mix (Takara Biotechnology Co, Ltd, Dalian, China) according to the manufacturer's instructions and quantitative PCR was performed using SYBR Premix Ex Taq (Perfect Real-Time) Kit (Takara Biotechnology, Dalian, PR China) according to the manufacturer's protocol. PCRs were performed in triplicate in a total reaction volume of 20 μL, including 10 μL SYBR Premix Ex Taq (2×), 1 μL of PCR Forward Primer (10 μM), 1 μL of PCR Reverse Primer (10 μM), and 8 μL of cDNA (diluted in double-distilled water). The quantitative real-time PCR reaction was set at an initial step of 3 minutes at 50°C and 3 minutes at 95°C; and 95°C (10 seconds) and 60°C (30 seconds) in a total of 40 cycles. All experiments were done in triplicate and all samples were normalized to glyceraldehyde-3-phosphate dehydrogenase (GAPDH). Primers used for amplification were as follows: CXCL8 sense, 5′- CCTGATTTCTGCAGCTCTGTG-3′ and CXCL8 antisense, 5′-CCAGACAGAGCTCTC.

TTCCAT-3′; STAT1 sense, 5′-GGTTAACGTTCGCACTCTGTG-3′ and STAT1 antisense, 5′-GCTGCTGAAGTTCGTACCAC-3′; COL17A1 sense, 5′-GCAGCAGCGGCTACATA

AAC-3′ and COL17A1 antisense, 5′-AGAAGAGTTGCCACTGGAGC-3′; insulin-like growth factor 1 receptor (IGF1R) sense, 5′-CTCTGGCCGACGAGTGGAG-3′ and IGF1R antisense, 5′-CTCGGTAATGACCGT

GAGCTT-3′; CCL27 sense, 5′- AGGCTGAGCAACATGAAGGG-3′ and CCL27 antisense, 5′- TGCTGGGTGGCAGTAGGAAT-3′; GAPDH sense, 5′- CAGTGCCAGCCTCGTCTCAT-3′ and GAPDH antisense, 5′-AGGGGCCATCCACAGTCTTC-3′. Statistical analysis of the PCR data was conducted using the relative expression software tool algorithm, which uses a pairwise-fixed reallocation and randomization test to determine significance (28). Data are expressed as the mean ± SEM from 3 independent experiments. Statistical analysis was carried out with SPSS 22.0 software (SPSS Inc, Chicago, IL) using a Student t test. A value of P < .05 was considered as statistically significant.

Results

Data preprocessing and DEGs analysis

The expression values of 20,545 genes were obtained from the 32 samples included in E-MTAB-1862. After further annotating the RNASeqV2 data downloaded from TCGA database, the expression values of 465 lncRNAs were also acquired from 448 samples. Among the 448 samples, 154 samples had annotation information of death time and were used for the subsequent survival analysis.

Under the thresholds of |log2FC >1 and FDR <0.05, total of 382 DEGs (including 206 upregulated genes and 176 downregulated genes) were screened in primary melanoma samples compared with benign nevus samples. Therefore, the number of upregulated genes was more than that of downregulated genes.

The top 5 GO and pathway terms enriched for the upregulated genes are listed in Table 1, including immune response (GO_BP, P = 1.70E-14), extracellular region part (GO_CC, P = 1.74E-07), chemokine activity (GO_MF, P = 1.09E-05), and toll-like receptor signaling pathway (pathway, P = 7.68E-05). In addition, the GO terms significantly enriched for the downregulated genes including epidermis development (GO_BP; P = 8.98E-07; which involved collagen, type XVII, alpha 1, COL17A1), integral to plasma membrane (GO_CC, P = 4.19E-04) and sequence-specific DNA binding (GO_MF, P = 5.64E-03) (Table 2). For the downregulated genes, there were no significantly enriched pathways. The enriched GO_BP terms for the upregulated genes and the downregulated genes separately are shown in Figures 2 and 3.

Table 1

The top 5 functional and pathway terms enriched for the upregulated genes.

Table 2

The top 5 GO terms enriched for the downregulated genes.

Figure 2

The Gene Ontology_biological process (GO_BP) terms enriched for the upregulated genes.

Figure 3

The Gene Ontology_biological process (GO_BP) terms enriched for the downregulated genes.

A total of 226 nodes (including 137 upregulated genes and 89 downregulated genes) and 598 interactions were involved in the PPI network constructed for the DEGs (Fig. 4). Especially, chemokine (C-X-C motif) ligand 8 (CXCL8, degree = 29) and signal transducers and activators of transcription 1 (STAT1, degree = 28) were key nodes in the PPI network for they had higher degrees, and they could interact with each other.

Figure 4

The protein-protein interaction network constructed for the differentially expressed genes.

Survival analysis was performed for the DEGs, lncRNAs, and miRNAs in melanoma samples. As a result, a total of 21 DEGs, 55 lncRNAs, and 32 miRNAs were found to be associated with prognosis.

Regulatory network construction

Based on coexpression relationships, the target genes of all prognostic risk lncRNAs were predicted (such as RP11–361L15.4 targeting COL17A1) and the lncRNA-gene regulatory network is shown in Figure 5. Meanwhile, 78 target genes were obtained for 6 (hsa-miR-100, hsa-miR-10a, hsa-miR-183, hsa-miR-184, hsa-miR-22, and hsa-miR-375) of the prognostic risk miRNAs (such as hsa-miR-375 targeting chemokine (C-C motif) ligand 27, CCL27; and hsa-miR-375 targeting IGF1R). In addition, the miRNA gene regulatory network is shown in Figure 6.

Figure 5

The regulatory network of prognostic risk long noncoding RNAs (lncRNAs) and their target genes. Purple diamonds represent lncRNAs. Red and green circles stand for upregulated and downregulated genes, respectively. Dotted lines indicate negative correlation between lncRNAs and target genes.

Figure 6

The regulatory network of prognostic risk microRNAs (miRNAs) and their target genes. Purple rectangles represent miRNAs. Red and green circles stand for upregulated and downregulated genes, respectively.

Validation of several DEGs by real-time RT-PCR

Aiming to further verify the expression patterns of selected DEGs, real-time PCR, which allows quantitative analysis of mRNA expression, was applied. A total of 5 primary melanoma samples and 10 benign nevus samples (control) were obtained from 15 additional persons who were not part of the microarray analysis. As shown in Figure 7, the overall direction of differential expression was consistent. Among 5 genes assayed, 2 genes (STAT1 and IGF1R) yielded concordance with the results of microarray analysis (STAT1 was significantly upregulated and IGF1R was significantly downregulated). The other 2 genes (CXCL8 and CCL27) showed no significant changes and COL17A1 showed a slight change in the opposite direction, which may be affected by the small sample size.

Figure 7

Real-time reverse transcription polymerase chain reaction (RT-PCR) validation of several identified differentially expressed genes. Control indicates benign nevus samples. Melanoma indicates primary melanoma samples. ∗, P < .05; ∗∗, P < .01.

Discussion

In the current study, a total of 382 DEGs were identified in primary melanoma samples, including 206 upregulated genes and 176 downregulated genes. In the PPI network, CXCL8 (degree = 29) and STAT1 (degree = 28) had higher degrees and could interact with each other. Previous study demonstrates that CXCL8 expression acts in mediating multiple cellular phenotypes related to growth and metastasis of melanoma.[CXCL8 can cause cell proliferation and angiogenesis through both receptors (chemokine [C-X-C motif] receptor 1, CXCR1; and chemokine [C-X-C motif] receptor 2, CXCR2) and mediate invasive and migration of melanoma cells via CXCR2; hence, blocking the CXCL8 signaling axis may be a promising therapeutic method for metastatic melanoma.[ The STAT1/STAT3 balance in tumor cells and host lymphocytes is significantly regulated by interferonα2b, and the ratio of pSTAT1/pSTAT3 in tumor cells at baseline may be a promising predictor of therapeutic effect in patients with cutaneous melanoma.[STAT1 is speculated to be involved in the progression of late-stage melanoma; thus, interference with the STAT1 pathway may have curative effects for patients with late-stage melanoma.[ Real-time PCR in the present study had further confirmed that STAT1 was overexpressed in melanoma. Collectively, we suggested that CXCL8 and STAT1 might function in the development and progression of melanoma via interacting with each other. More experimental investigations are required to verify the role of CXCL8 and STAT1 and their interaction.

Survival analysis showed that 21 DEGs, 55 lncRNAs, and 32 miRNAs were found to be associated with prognosis. Furthermore, several regulatory relationships were found in the miRNA gene regulatory network (such as hsa-miR-375 targeting CCL27 and hsa-miR-375 targeting IGF1R). Ectopic expression of miR-375 inhibits motility, proliferation, and invasion, and causes shape changes of melanoma cells, indicating that miR-375 may be implicated in the development and progression of melanomas.[CCL27 and chemokine (C-C motif) receptor 10 (CCR10) in human melanomas may help tumor cells to grow, invade tissue, evade immune response, and spread to lymph nodes.[ Villanueva et al[ found that IGF1R/phosphatidylinositol 3-kinase (PI3K) signaling is boosted in resistant melanomas, and combined treatment with MAPK and IGF1R/PI3K inhibitors can overcome acquired resistance to BRAF.[ Uveal melanomas express MET, epidermal growth factor receptor, hepatocyte growth factor, and IGF-1R, in addition, epidermal growth factor and IGF1 may have correlation with the pathogenesis of uveal melanoma.[ In human melanoma, IGF1 is found to be a novel growth factor for autocrine-driven proliferation in vitro, and IGF1-IGF1R autocrine pathway may be a therapeutic target for the disease.[ Therefore, CCL27 and IGF1R targeted by hsa-miR-375 might also associate with the mechanisms of melanoma. The validation of real-time PCR strongly supported the downregulation of IGF1R; the expression of CCL27 and the role of hsa-miR-375 in melanoma should be further verified by multiple approaches.

Functional enrichment analysis showed that COL17A1 was enriched in epidermis development. Fibrillar collagen can act as a powerful barrier to viral distribution and matrix-modifying treatments significantly promote the treatment response of human melanoma.[ Van Kempen et al[ firstly report that the proangiogenic effect of type I collagen in a de novo tumor model keeps ahead of the development of invasive melanomas. Collagen types XVII, which is a type of epithelial adhesion protein, has correlation with the epidermal basement membrane and plays a key functional role.[ Moreover, collagen XVII endodomain accumulated in melanocytic tumors is related to malignant transformation and regulates antibody-induced melanoma apoptosis.[ The majority of evidence has demonstrated that microarrays are invaluable discovery tools with acceptable reliability for genome-wide gene expression screening.[ However, microarray data alone often lack the degree of specificity, which was needed to make accurate biological conclusions.[ Genes found to be differentially expressed by microarray analysis should be validated with other well-established technologies for analyzing gene expression, such as real-time reverse transcription (RT)-PCR. Dallas et al[ reported a 13% to 16% nonconcordance between quantitative reverse transcription PCR (qRT-PCR) and microarray data, which may emphasize the continuing requirement for caution in validation of gene expression data. In this study, the qRT-PCR results showed a slight change in the opposite direction of COL17A1. However, considering the role of collagen XVII demonstrated by the previous study,[ we still suggested the significant role of COL17A1 in melanoma. On the contrary, in the lncRNA-gene regulatory network, COL17A1 was targeted by the lncRNA RP11–361L15.4, indicating that COL17A1 targeted by RP11-361L15.4 might function in the pathogenesis of melanoma. More experiments should be performed to explore the expression of COL17A1 and the association of RP11-361L15.4 and COL17A1.

In conclusion, we have conducted an in-depth bioinformatics analysis and identified a total of 382 DEGs in primary melanoma samples. CXCL8 and STAT1 might function in the development and progression of melanoma via interacting with each other. CCL27 and IGF1R targeted by hsa-miR-375 might also associate with the mechanisms of melanoma. Furthermore, COL17A1 targeted by RP11-361L15.4 might function in the pathogenesis of melanoma. A subset of selected genes was validated by real-time RT-PCR and several genes were confirmed their differential expression in small size of melanoma samples. More in vivo or in vitro experiments with larger sample size were required to further explore the exact roles of the identified genes, miRNAs and lncRNAs in melanoma, which may represent potential tumor markers or drug targets.