Construction of competing endogenous RNA interaction network as prognostic markers in metastatic melanoma

Patient cohort information and data processing

The melanoma patient cohort data were obtained from The Cancer Genome Atlas (TCGA), which is deposited on the GDC website (https://portal.gdc.cancer.gov/repository). Three files were downloaded containing miRNA, isoform, and clinical information, respectively. A total of 470 melanoma samples were stored in the database, including 103 in situ melanoma samples and 367 metastatic melanoma samples. The data were downloaded using the “download” command from the recommended download tool gdc-client (https://docs.gdc.cancer.gov/Data_Transfer_Tool/Users_Guide/Data_Download_and_Upload/). The clinical information and miRNA expression matrix were obtained using the R package “XML”. mRNA and lncRNA were extracted separately according to the genome annotation file. The mRNA extraction parameters were type = “gene”, gene_biotype = “protein_coding”, and gene_biotype = “protein_coding”, the lncRNA extraction parameters arguments were type = “transcript”, transcript_biotype = “lncRNA”. “3prime_overlapping_ncRNA”, “bidirectional_promoter_lncRNA”, “sense_. overlapping”, “sense_intronic” and “antisense_RNA”.

Identification of differentially expressed genes

We used the R package “edgR” to identify differential genes between metastatic and in situ melanoma samples (Robinson, McCarthy & Smyth, 2010). All q-values were corrected for statistical significance of multiple testing using FDR, with screening thresholds of |log2FC| > 1.0, FRD < 0.05, and retention of significantly different genes (Lai, 2017; McCarthy, Chen & Smyth, 2012).

Functional enrichment analysis

ClusterProfiler was used for GO/KEGG analysis (Gene Ontology/Kyoto Encyclopedia of Genes and Genomes) (Ashburner et al., 2000; Kanehisa et al., 2010; Yu et al., 2012). GO was used to describe gene functions along three aspects: biological process (BP), cellular component (CC), and molecular function (MF). KEGG was searched for pathways at the significance level of p < 0.05.

Weighted correlation network analysis

WGCNA is an algorithm used for the identification of gene co-expression networks for mRNAs or lncRNAs with different characteristics by high-throughput spectrum of expression. Neighborhood matrix assessment of weighted co-expression relationships between all datasets was performed using Pearson’s correlation analysis. We used WGCNA to analyze mRNAs and lncRNAs to obtain the mRNAs and lncRNAs with which the metastatic melanoma samples were most associated.

Construction of a ceRNA regulatory network

Some authoritative databases have provided relationship pairs such as lncRNA-miRNA and miRNA-mRNA, which can assist in the construction of a ceRNA network centered on miRNA, based on an experimental module and prediction module. MirDIP (http://ophid.utoronto.ca/mirDIP/index.jsp) was used to predict interactions between lncRNAs and miRNAs. LncBase (http://carolina.imis.athena-innovation.gr/diana_tools/web/index.php?r=lncbasev2%2Findex) databases were used to explore target mRNAs.

Construction of metastatic melanoma survival risk prognosis model

Metastatic melanoma samples with low expression of miRNA, high expression of lncRNA, and mRNA in the ceRNA network were used to perform three LASSO regression modeling filters. Factors that fit the normalized lambda value model were included in the subsequent analysis. A prognostic model was constructed using multiple COX regression analyses, and a comprehensive prognostic scoring system was obtained based on these factors. The risk score was calculated as follows: risk score = ∑βi × expRNAi, where expRNA is the expression level of RNA and β is the regression coefficient derived from the multivariate Cox regression model. Melanoma patients were classified into high-risk or low-risk groups based on the median risk score cutoff value based on this risk score formula and other factors we determined. Kaplan-Meier analysis was used to determine the overall survival between the groups and there was a difference in overall survival (OS) between the two groups. P < 0.05 was considered significant unless otherwise stated.

Cell line and cell culture

A375 (Catalog number: CRL-1619; ATCC, Manassan, VA, USA) cells were cultured in an incubator in Dulbecco’s modified Eagle’s medium (DMEM, Gibco, Grand Island, NY, USA) containing 10% fetal bovine serum (AusGeneX, Molendinar, Qld, Australia) and 1% penicillin–streptomycin (P/S; Invitrogen, Carlsbad, CA, USA) with a humidity of 5% CO₂ at 37 °C.

Fluorescence in situ hybridization (FISH)

The A375 cells were inoculated in a 12-well plate and incubated overnight in an incubator set at 37 °C with 5% CO2. The medium in the orifice plate was absorbed and washed twice with PBS, and then one ml 4% paraformaldehyde was added to each well, and fixed at room temperature for 15 min. We followed the manufacturer’s instructions for the RNA FISH kit (Catalog number: F12201/50; GenePharma, Suzhou, China). Five probes were designed and the sequences are listed in Table S1. We added 200 μl diluted DAPI working solution to each well, and stained the samples for 20 min in dark. An anti-quenching agent was added to a clean slide, and these were observed under a fluorescence microscope (630×, oil lens).

Plasmid, siRNA and primer synthesize

Plasmids of miRNA-has-miR-3662 and miRNA-shNC, siRNA of siRP11-594N15.3 and si-NC were bought from GenePharma. The overexpressed plasmids of the miR-3666 sequence are shown in Table S2. The primers of RP11-594N15.3, ZNF831, PKIA, CSF2RB, miR-3662, GAPDH, U6 and reverse universal primer were designed by Primer Premier 5 software and were purchased from the Beijing Genomics Institute (Table S3).

Expression vectors transfection

A375 cells were cultured in 10 cm petri dishes at a concentration of 7 × 10⁶/8 ml the day after transfection. Cells were stably transduced with the miR-3662 and shNC expression vectors, using Lipofectamine 2000 (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s protocol. Each plate was added to 12 μg plasmids and 24 μl Lipofectamine 2000. The medium was replaced after 6 h of transfection. Each plate was added to 80 μM siRNA for the transfection of siRNA and additional steps were performed as described above.

RNA isolation and quantitative real-time PCR (qRT-PCR)

We used TRIzol^® Reagent (Life Technologies, Carlsbad, CA, USA) to extract RNA from cells. After wiped off DNA, equal amounts (2 μg) of each RNA were separately reverse transcribed to their corresponding cDNAs using a RevertAid^™ First Strand cDNA Synthesis Kit (Thermo Fisher Scientific, Waltham, MA, USA). All the opeations are subject to the manufacturer’s instructions. CFX96 Real-Time PCR detection system (Bio-Rad, Irvine, CA, USA) was used to perform qRT-PCR with Maxima^® SYBR Green/ROX qPCR Master Mix (Thermo Fisher Scientific, Waltham, MA, USA). U6 and GAPDH were used as internal control for miR-3662 and other mRNAs respectively.

Cell counting kit-8 (CCK-8) assay

A total of 2,000 cells/well were seeded into 96-well plates, and were cultured for one, two, three, or four days. The condition of cell proliferation was assessed with a Cell Counting Kit-8 assay (C0038; Beyotime Biotechnology, China). According to manufacturer’s instruction, cells were treated with 10 μL CCK-8 (100 μL culture medium per well) solution for 2 h at 37 °C after a period of incubation. Absorbance at 450 nm was measured by using a microplate reader.

Transwell migration and invasion assays

The cells were resuspended in serum-free medium the day before experiments. For Transwell invasion assays, Transwell upper chambers (8 mm; BD Biosciences, Franklin Lakes, BJ, USA) were coated with Matrigel matrix (dilution 1:7; 356234; BD Biosciences, Franklin Lakes, BJ, USA) and 2 × 10⁵ cells were placed in. For Transwell migration assays, 2 × 10⁴ cells were seeded in the upper chamber of the Transwell without Matrigel matrix. Then, a medium containing 20% serum was added to lower chambers. Cells were incubated at 37 °C and 5% CO2 for 18 h, and then those cells remaining in upper chambers were removed using a wet cotton swab. We used 4% paraformaldehyde to fix the cells that had migrated and adhered to lower chambers for 10 min. Then hematoxylin (ZLI-9610, ZSGB-Bio) and eosin (ZLI-9613, ZSGB-Bio) were used to stain cells adhered to lower chambers for 20 min and then imaged. The number of cells was counted in ten separated high power fields with vertical cross distribution.

Identification of differential genes in metastatic and in situ melanoma samples

There were 103 samples of melanoma in situ and 367 samples of metastatic melanoma in the TCGA database (Tables S4, S5). We identified 1,413 differential lncRNAs, of which 1,059 were highly-expressed in metastatic samples and 354 in the in-situ samples. We identified 2,251 differential mRNAs, of which 1,442 were highly-expressed in the metastatic samples and 809 in the in-situ samples. We identified 81 differential miRNAs, of which 60 were highly expressed in metastatic samples and 21 in the in-situ samples (Figs. 1B–1D, Figs. S1A–S1C). The differentially expressed mRNAs in the metastatic group were enriched for GO function and KEGG pathway, immune cell activity-related terms on biological processes (BP), transmembrane transduction structures on cellular composition (CC), G protein-coupled receptor activity, passive transport and other entries on molecular function (MF) (Fig. 1E). The KEGG pathway was enriched for the PI3K-Akt signaling pathway, chemokine signaling pathway, JAK-STAT signaling pathway, cell adhesion molecules (CAMs), allograft rejection, NF-kappa-B signaling pathway, and other tumor-related pathways (Fig. 1F).

Identification of mRNA and lncRNA co-expression modules by weighted gene co-expression network analysis

We used weighted gene co-expression network analysis (WGCNA) to analyze the gene expression patterns of mRNA and lncRNA expression matrices extracted from metastatic and in-situ samples. Genes with similar expression patterns were clustered into a co-expression module (Langfelder & Horvath, 2008). We identified nine mRNA co-expression modules and 10 lncRNA co-expression modules in the two groups, and the number of genes contained in each module is shown in Table S6. The significantly related co-expression modules were selected according to the clustering of linkage levels and the number of co-expression associations. The mRNA clustering results showed that the black module and yellow module of in-situ samples, and the turquoise module of metastatic samples were identified as highly-associated co-expression modules (Figs. 2A–2C). The lncRNA clustering results showed that the pink and red module of in-situ samples, and the turquoise and brown modules of the metastatic samples were identified as highly correlated co-expression modules (Figs. 2D–2F). The corresponding interactions of the first 20 co-expressed genes in the above modules are shown in Figs. S2A–S2G. GO functional enrichment of co-expression genes in the turquoise module in mRNA showed that the module genes were associated with immune responses, and were also enriched for entries associated with immune negative regulation and cell adhesion (Fig. 2G). We removed the turquoise module in mRNA and the brown/turquoise module in lncRNA for the next step of the analysis.

Construction of ceRNA regulation network

The mRNAs from the turquoise module obtained after WGCNA analysis were intersected with the highly expressed mRNAs in the metastatic samples, and 319 mRNAs were obtained. lncRNAs from the brown/turquoise module were intersected with the highly expressed mRNAs in the metastatic samples, and 287 lncRNAs were obtained (Figs. 3A–3B). The above 319 mRNAs were subjected to Pearson’s correlation coefficient calculation with 287 lncRNAs. The results showed that 138 mRNAs had potential interactions with 89 lncRNAs when cor > 0.5 and p-value < 0.05 (Fig. S3). Using the mirDIP database to predict the target genes, the low expression miRNAs in the 21 metastatic samples were predicted to score “very high” of 3,771 mRNA target genes. The 3,771 target genes were cross-linked to 138 mRNAs by correlation analysis to obtain 11 mRNAs (CD53, CD96, CSF2RB, CXCL9, MS4A1, PIM2, PKIA, PTPN22, SH2D1A, SLA, ZNF831) (Fig. 3C). The 11 mRNAs and their corresponding six miRNAs (has-miR-346, has-miR-3662, has-miR-429, has-miR-891b, has-miR-892a, has-miR-944) were analyzed and mapped for the interaction network (Fig. 3D). These six miRNAs were used to predict 3,760 related lncRNAs in the lncBase database. The 3,760 lncRNAs were intersected with 89 lnRNAs screened in correlation analysis to obtain 17 lncRNAs (AC007386.4, AC104024.1, GTSCR1, LINC00402, LINC00861, MIAT, PRKCQ-AS1, RP11-121A8.1, RP11-164H13.1, RP11-202G28.1, RP11-203B7.2, RP11-594N15.3, RP11-598F7.3, RP11-693J15.5, RP11-768B22.2, RP5-887A10.1) (Fig. 3E). Interaction networks were mapped for the above 17 lncRNAs and six miRNAs (Fig. 3F). The final selection of 11 mRNAs and lncRNAs, the data tSNE downscaling analysis, were better able to distinguish the in-situ samples from the metastatic samples (Fig. 3G), resulting in the construction of lnRNA-miRNA-mRNA ceRNA network with the final selection of six miRNAs (Fig. 3H). The factors of our network were found to be highly expressed in metastatic samples in another melanoma cohort (GEO database; accession number: GSE65904) (Fig. 3I).

Construction of metastatic melanoma survival risk prognosis model

We clustered the 34-factors included in our ceRNA interaction network, and 17 lncRNAs, six miRNAs, and 11 mRNAs were subjected to risk prediction by LASSO regression, respectively. The regularization term controlled for the lncRNA, miRNA, and mRNA models having minimal likelihood error when λ = 7, 1, and 2, respectively. We screened lncRNAs, one miRNA, and two mRNAs for high correlation with prognosis (Figs. 4A–4F). The above 10 factors were analyzed using multiple COX regression to construct risk models. There were five factors considered to be risk factors for prognosis and the other five factors were protective factors for prognosis. lncRNA AC104024.1, MIAT, RP11-598F7.3, and miRNA hsa-mir-429 could be used as significant risk markers (P < 0.05). The concordance index of the regression model was 0.753 (Fig. 4G). In the cohort, the expressions of AC104024.1, AC007386.4, and has-mir-429 were gradually up-regulated as the risk scoring increased, while the expressions of RP11-598F7.3, PKIA, CXCL9 and MIAT were down-regulated (Fig. 4H). The KM curves of melanoma patient survival were plotted according to the expression profiles of risk and protective factors in the multivariate COX regression model. These showed that the altered expression of factors in the model was significantly correlated with the follow-up prognosis of the cohort of patients (P < 0.0001) (Fig. 4I). Four groups of ceRNA interacting triplets were obtained (RP11-594N15.3–miR3662-CSF2RB, RP11-594N15.3–miR3662-ZNF831, RP11-594N15.3–miR3662-PKIA, AC104024.1–miR346-PKIA), of which miR3662 may have potential for the treatment of metastatic melanoma (Fig. 4J).

miR-3662 down-regulates its target mRNAs and suppresses melanoma cell proliferation, migration and invasion in vitro

In order to verify whether the regulatory relationship of ceRNA was objective, we chose three groups of ceRNA interacting triplets (RP11-594N15.3–miR3662-CSF2RB, RP11-594N15.3–miR3662-ZNF831, RP11-594N15.3–miR3662-ZNF831, RP11-594N15.3–miR3662-PKIA) for experimental verification. The localization of lncRNA in the cytoplasm is the premise of its spongy function. The database “lncBase” was used to predict the interaction of lncRNA-miRNA based on a base experimental module and prediction module. The binding sites on the sequence between lncRNAs and miRNAs, and the interaction between lncRNAs and miRNAs was tested in other studies. Since miRNAs only exist in the cytoplasm, the lncRNAs predicted by the database must also theoretically exist in the cytoplasm. We conducted FISH experiment in melanoma A375 cells and selected lncRNA RP11-594N15.3 for localization to ensure that our lncRNAs conformed. Our results show that RP11-594N15.3 may be localized in the cytoplasm (Fig. 5A). Since RP11-594N15.3 is an RNA with a length of more than 3 KB, it is very difficult to construct its overexpression vector. Therefore, we decided to knock down the expression of RP11-594N15.3 with siRNA in A375 cells to observe the changes in the expression of miR-3662, CSF2RB, PKIA and ZNF831 (Fig. 5B). miR-3662 was found to be significantly overexpressed, PKIA and ZNF831 were down-regulated, and CSF2RB was not significant (Figs. 5C, 5D). Our experiments proved that the sponge function of lncRNA was objective.

In order to confirm whether the regulation of miR-3662 on target mRNAs in ceRNA regulatory triplets was objective, we transfected the miR-3662 overexpression plasmid vector and the NC plasmid vector into melanoma A375 cells, obtained A375-miR3662 and A375-NC cells, and verified the high expression of miR-3662 at the transcriptional level (Fig. 5E). The qRT-PCR of the downstream target mRNAs shows the transcriptional differences of the three targets. The transcriptional levels of PKIA and ZNF831 were down-regulated by miR-3662, but there was no significant change on CSF2RB (Fig. 5F). The overexpression of miR-3662 significantly promoted cell proliferation (Fig. 5G). We used the Transwell chamber assay to explore the functions of miR-3662 on the invasion and migration of melanoma cells. The overexpression of miR-3662 significantly suppressed cell migration through a permeable filter and invasion through Matrigel Matrix (Fig. 5H). Our results indicate that miR-3662 may repress melanoma cell proliferation, migration and invasion by down-regulating PKIA and ZNF831.