Abstract

   Runx2 is a transcription factor involved in melanoma cell migration and
   proliferation. Here, we extended the analysis of Runt domain of Runx2
   in melanoma cells to deepen understanding of the underlying mechanisms.
   By the CRISPR/Cas9 system we generated the Runt KO melanoma cells 3G8.
   Interestingly, the proteome analysis showed a specific protein
   signature of 3G8 cells related to apoptosis and migration, and pointed
   out the involvement of Runt domain in the neoangiogenesis process.
   Among the proteins implicated in angiogenesis we identified fatty acid
   synthase, chloride intracellular channel protein-4, heat shock protein
   beta-1, Rho guanine nucleotide exchange factor 1, D-3-phosphoglycerate
   dehydrogenase, myosin-1c and caveolin-1. Upon querying the TCGA
   provisional database for melanoma, the genes related to these proteins
   were found altered in 51.36% of total patients. In addition, VEGF gene
   expression was reduced in 3G8 as compared to A375 cells; and HUVEC
   co-cultured with 3G8 cells expressed lower levels of CD105 and CD31
   neoangiogenetic markers. Furthermore, the tube formation assay revealed
   down-regulation of capillary-like structures in HUVEC co-cultured with
   3G8 in comparison to those with A375 cells. These findings provide new
   insight into Runx2 molecular details which can be crucial to possibly
   propose it as an oncotarget of melanoma.

   Subject terms: Melanoma, Protein-protein interaction networks

Introduction

   In the past half century, the incidence of metastatic melanoma has
   increased in the general population^[59]1. Western life style may
   possibly contribute to melanoma occurrence and the mortality rate
   associated to metastatic melanoma is very high. Different factors,
   endogenous or environmental, are associated to malignant transformation
   in melanoma. Among them, previous family cases of melanoma as well as
   many dysplastic nevi are strongly associated to melanoma^[60]2. In
   addition, exposure to sun or ultraviolet radiations, phenotype
   characterized by blond or red hair and pale skin, reduction of immune
   protection and chemical exposures contribute to melanoma
   onset^[61]3–[62]5.

   Several studies identified the molecular events involved in the
   malignant transformation of melanoma. Transcription regulators and
   proto-oncogenes such as serine/threonine kinase B-Raf (BRAF), mast/stem
   cell growth factor receptor c-Kit (KIT), neuroblastoma RAS viral
   oncogene homolog (NRAS), phosphatase and tensin homologue (PTEN), tumor
   suppressor p53 (P53), telomerase reverse transcriptase (TERT), and
   microphthalmia-associated transcription factor (MITF) are associated
   with metastatic melanoma. In addition, mutations in cyclin-dependent
   kinase inhibitor 2A (CDKN2A) encoding two tumour suppressor proteins,
   p16 and p14, as well as genes involved in the DNA repair mechanism or
   the melanocortin-1 receptor (MC1R) play an important role in the onset
   of melanoma^[63]2,[64]6. We recently reported new melanoma oncotargets,
   such as enolase 1, parkinsonism-associated deglycase, prostaglandin E
   synthase 3, nucleophosmin, and stathmin 1, playing a key role in
   malignant transformation^[65]7. Moreover, molecular analyses showed the
   central role of the PTEN/mutated in multiple advanced cancers 1 (MMAC1)
   in reducing cell proliferation or inducing apoptosis^[66]8, and the
   association of PTEN loss to melanoma^[67]9. Furthermore, angiogenesis
   is an important factor allowing melanoma migration and metastases.
   Mehnert et al. reported that the expression of vascular endothelial
   growth factor (VEGF) and its receptors (VEGF-R) are higher in melanomas
   compare to normal melanocytes^[68]10.

   Recently, we identified the involvement of the Runt domain of the
   Runt-related transcription factor 2 (RUNX2) gene in melanoma^[69]11.
   RUNX2 is the most important transcription factor involved in osteogenic
   differentiation. RUNX2 is also ectopically expressed in several solid
   tumors such as in thyroid, breast, pancreatic, prostate, lung, ovarian
   epithelial cancers^[70]12,[71]13 and melanoma^[72]14. We reported that
   RUNX2 promotes migration and melanoma proliferation and restricts the
   apoptosis through its Runt domain^[73]11. Importantly, the Runt domain
   downregulates P53 expression as well as single-stranded DNA-binding
   protein 1 (SSBP1), a protein that controls the stability and the
   transcriptional activity of P53^[74]15.

   Therefore, in order to identify the molecular pathways associated with
   RUNX2 in melanoma we performed a comparative shotgun proteomic analysis
   of Runt KO (3G8) and wild-type (A375) melanoma cells. We identified for
   the first time proteins associated with the Runt domain involved in
   promoting migration, proliferation and in particular in neoangiogenesis
   of melanoma cells. We then confirm the role of Runx2 in
   neoangiogenesis, analyzing the gene expression of VEGF in 3G8 and A375
   cells, and performing co-culture experiments using human umbilical vein
   endothelial cell (HUVEC) plus 3G8 or A375 cells. We observed that HUVEC
   express lower levels of CD105 and CD31 neoangiogenetic markers and have
   reduced network-like structures when co-cultured with 3G8 cells,
   confirming a role for Runx2 in neoangiogenesis.

Experimental Procedures

Cell cultures and co-cultures

   A375 melanoma cells (at passage 5) and normal human epidermal
   melanocytes (NHEM) were purchased from American Type Culture Collection
   (ATCC^© CTRL-1619^TM) (Rockville, MD, USA). Cell lines were cultured
   under a humidified atmosphere of 5%CO[2] and passage in growth medium:
   RMPI 1640 (Sigma-Aldrich) containing 10% FBS (fetal bovine serum)
   (Sigma-Aldrich; Merck Millipore) supplemented with antibiotics (1%
   penicillin/streptomycin) and 1% glutamine. Once 70–80% confluence was
   reached, cells were harvested using trypsin, washed and counted on a
   microscope using a Burker hemocytometer.

   Primary HUVEC (at passage 10) were kindly provided from Prof. Bazzoni
   of the University of Verona. Cells were cultured in a 75 cm^2 culture
   flask coated with Matrigel (Corning Incorporated; NY, USA) in the
   presence of M200 medium and 10% of LSGS supplement (Thermo Fisher
   scientific, Waltham, MA, USA) under a humidified atmosphere of 5%
   CO[2]. Co-cultures for CD105 gene and CD31 gene/protein expression
   analyses, as well as for net-like structure detection, were performed
   using 6 wells plate Transwell with 0.4 µm pores (Corning Incorporated;
   NY, USA). HUVEC cells were plated in the chamber, in presence of M200
   supplemented medium, while in the permeable membrane were plated A375,
   3G8 Crisped clone and NHEM cells using RPMI + 10% FBS. Each co-culture
   was performed in three replicates.

CRISPR/Cas9

   The CRISPR/Cas9 technique was applied as previously reported^[75]11. In
   brief, we used the following two specific gRNAs designed by analyzing
   the target sequence with both CHOPCHOP^[76]16,[77]17 and MIT
   ([78]http://crispr.mit.edu/) CRISPR design tools (gRNA A
   CCCATCTGGTACCTCTCCGA; gRNA B GATCGTTGAACCTTGCTACT). The above gRNAs
   were individually cloned in the PX459 V2.0 Cas9 expressing vector
   (Addgene), following the protocol described by Ran et al.^[79]18. The
   A375 cells were co-transfected and single cell clonings were obtained
   as previously reported^[80]11. To analyse the effects of CRISPR/Cas9
   editing we called the WT and CRISPR/Cas9 edited cell lines A375 WT and
   3G8, respectively. To confirm the deletion, we sequenced the PCR
   products by Sanger Sequencing. In particular, the specific Runt domain
   was amplified by PCR with the following primers: FW
   TGAAGTGGCATCACAACCCA; RV AGTCAGAGACCTACCTCGTC. The PCR products were
   purified with the FastGene® extraction kit (Nippon Genetics, Tokyo,
   Japan) and the forward PCR primer was used for Sanger sequencing
   (GenomeLab™ DTCS quick start kit and CEQ. 8000 Genetic Analysis System,
   SCIEX, San Francisco CA, USA) following manufacturer’s instructions.

Sample preparation for proteomics

   A375 wild-type (WT) and 3G8 melanoma cells were collected, washed in 1x
   phosphate-buffered saline (PBS), resuspended in RIPA buffer (50 mM
   Tris-HCl pH 8.0 with 150 mM NaCl, 1.0% Igepal CA-630, 0.5% sodium
   deoxycholate and 0.1% sodium dodecyl sulfate, Sigma-Aldrich)
   supplemented with a protease inhibitors cocktail 1x (Roche) and
   sonicated for 5–6 cycles. The lysate was gently mixed for 15 min and
   then centrifuged at 14,000 × g for 15 min at 4 °C to pellet cell
   debris. The supernatant was collected and total protein concentration
   was determined using a bicinchoninic acid protein assay
   (Sigma-Aldrich). Before mass spectrometry (MS) analysis, proteins were
   digested by trypsin following the protocol previously described^[81]19.

Mass spectrometric analysis and data processing

   Protein identification and quantification by liquid
   chromatography-tandem mass spectrometry (LC-MS/MS) analyses were
   performed as previously reported^[82]19. Briefly, LC–MS/MS analyses
   were performed using a micro-LC Eksigent Technologies (Dublin, USA)
   system interfaced with a 5600 + TripleTOF system (AB Sciex, Concord,
   Canada). The injection volume of each sample was 4.0 μL. Samples used
   to generate the Sequential Window Acquisition of All Theoretical Mass
   Spectra (SWATH-MS) spectral library were subjected to the traditional
   data-dependent acquisition (DDA), and then to cyclic data independent
   analysis (DIA) of the mass spectra, using a 25-Da window. The MS data
   were acquired with Analyst TF 1.7 (AB SCIEX, Concord, Canada). Three
   instrumental replicates for each sample were subjected to the DIA
   analysis. The MS files were searched for protein identification.
   PeakView 1.2.0.3 and Protein Pilot software v. 4.2 (AB SCIEX, Concord,
   Canada) were used to generate the peak-list. The mass spectrometry
   files were searched using Protein Pilot (AB SCIEX, Concord, Canada),
   Mascot (Matrix Science, Inc, Boston, MA).

   Samples were input in the Protein Pilot software v. 4.2 (AB SCIEX,
   Concord, Canada), which employs the Paragon algorithm, with the
   following parameters: cysteine alkylation, digestion by trypsin, no
   special factors and False Discovery Rate (FDR) at 1%. The UniProt
   Swiss‐Prot reviewed database containing human proteins (version
   2015.07.07, containing 42131 sequence entries). The Mascot search was
   performed on Mascot v. 2.4, the digestion enzyme selected was trypsin,
   with two missed cleavages and a search tolerance of 50 ppm was
   specified for the peptide mass tolerance, and 0.1 Da for the MS/MS
   tolerance. The charges of the peptides to search for were set to 2+,
   3+, and 4+, and the search was set on monoisotopic mass. The instrument
   was set to ESI‐QUAD‐TOF and the following modifications were specified
   for the search: carbamidomethyl cysteines as fixed modification and
   oxidized methionine as variable modification.

   Quantification was performed by integrating the extracted ion
   chromatogram of all the unique ions for a given peptide. SwathXtend was
   employed to build an integrated assay library with the DDA
   acquisitions, using a protein FDR threshold of 1%. Quantification was
   carried out with PeakView 2.0 and MarkerView 1.2. (ABSCIEX, Concord,
   Canada). The six peptides per protein with the highest MS1 intensity
   and six transitions per peptide were extracted from the SWATH files.
   Shared peptides were excluded as well as peptides with modifications.
   Peptides with FDR lower than 1.0% were exported in MarkerView for the
   t-test. The up- and down-regulated proteins were selected using
   p-value < 0.05 and fold change >1.5.

Western blotting

   Western blotting was performed as described in our former study^[83]20.
   Briefly, protein samples from two biological replicates were diluted
   1:1 in Laemmli’s sample buffer (62.5 mM Tris-HCl, pH 6.8, 25% glycerol,
   2% SDS, 0.01% Bromophenol Blue), heated for 5 min at 90 °C, and
   separated by sodium dodecyl sulfate−polyacrylamide gel electrophoresis
   (SDS-PAGE), followed by transfer onto polyvinylidene difluoride (PVDF)
   membranes. Amido Black staining (Sigma-Aldrich) was used to confirm
   equal protein loading in different lanes. Then the primary and
   secondary antibodies were incubated with PVDF membranes (Supplemental
   Table [84]1), and chemiluminescent signal acquisitions were done with
   the ChemiDoc MP Imager, a CCD imager, using the Image Lab 5.2.1
   software (Biorad).

Bioinformatics analysis

   Identified proteins were annotated according to the information
   associated to their main functions available under gene ontology (GO)
   categories using the Database for Annotation, Visualization and
   Integrated Discovery (DAVID) (v6.8)
   ([85]http://david.abcc.ncifcrf.gov/)^[86]21. Functional annotation
   chart report was used to identify GO biological processes, molecular
   function and cellular component. Gene enrichment analysis was performed
   by setting the threshold of EASE Score, a modified Fisher Exact
   P-Value, to 0.1 and p-value smaller than 0.01, while the threshold of
   minimum gene counts belonging to an annotation term was set to 2.

   Moreover, Ingenuity Pathway Analysis (IPA, Ingenuity Systems, Redwood
   City, CA) was performed, as previously described^[87]22, in order to
   identify the perturbed canonical pathways associated with deregulated
   proteins of 3G8 cells. Briefly, proteins associated with canonical
   pathways were estimated as significant using the Fisher’s exact test
   (p-value ≤ 0.01). In order to further explain the observed expression
   changes in 3G8 cells, the IPA Upstream Regulator analysis was performed
   to identify a putative cascade of upstream regulators. The upstream
   regulators were assumed as valid effectors of gene/protein expression
   if the corresponding p-value obtained by Fisher’s exact test was ≤0.01.
   Activation z-score algorithm was used to allow for prediction whether
   an upstream regulator is activated (z ≥ 2) or inactivated (z ≤ −2)
   based on the direction of expressional change of the associated genes.

   Finally, protein-protein interaction network analysis was carried out
   using STRING platform ([88]http://string-db.org)^[89]23 by setting the
   specie under investigation (Homo sapiens) and medium confidence level
   (score 0.4). We retrieved known interactions experimentally determined
   (pink edge) and on curated database (light blue edge), excluding all
   other prediction methods implemented in STRING (such as co-expression
   and text-mining). Additional white nodes and network depth were kept to
   the minimum value^[90]1, to exclude as many false positive interactions
   as possible. Genes were subjected to k-means clustering for five
   clusters and disconnected nodes were excluded.

   The Cancer Genome Atlas (TCGA) of skin cutaneous melanoma dataset that
   consisted of the maximum number of sequenced cases (471 patients) was
   interrogated through the cBio Cancer Genomics Portal
   ([91]http://www.cbioportal.org/) for alterations in our query list
   consisting in proteins associated to angiogenesis, i.e. fatty acid
   synthase (FASN), chloride intracellular channel protein-4 (CLIC4), heat
   shock protein beta-1 (HSPB1), Rho guanine nucleotide exchange factor 1
   (ARHGEF1), D-3-phosphoglycerate dehydrogenase (PHGDH), myosin-1c
   (MYO1C), and caveolin-1 (CAV1) for missense, truncating and in-frame
   mutations, amplification and deletions, mRNA up- and downregulation and
   protein up- and down regulation by reverse phase protein array (RPPA)
   assay. A survival analysis was also carried out comparing the cases
   with alterations in all, or in each, of the proteins involved in
   angiogenesis to the ones without alterations to look at the median
   overall survival.

Real time RT-PCR

   Total RNA extraction and reverse transcription were performed as
   previously reported^[92]11. PCRs were performed in a total volume of
   25 µl using 20 ng of cDNA from each sample. For endoglin (CD105,
   Hs00923996_m1), platelet endothelial cell adhesion molecule (CD31,
   Hs1065279_m1), and VEGF-A (Hs00900055_m1) expression analysis, we used
   TaqMan® probes (Life Technologies LTD, CA, USA) and TaqMan® Universal
   Master Mix (Life Technologies LTD, CA, USA). Gene expression was
   normalized to housekeeping Beta-2 Microglobulin (B2M, Hs99999901_s1)
   gene, calculating the relative fold expression differences. TaqMan SDS
   analysis software was used to analyse the Ct values. Three independent
   experiments with three replicates for each sample were performed.

Immunofluorescence

   After 7 days of co-culture, cells were fixed and processed according to
   the manufacturer’s protocols. As primary antibody we used CD31 (Dako,
   Glostrup, Denmark), diluted 1:20 in Antibody Dilution Buffer and
   incubated 2 hours at room temperature. Slides were then incubated with
   goat mouse fluorescein conjugated (Cat. AP124F, Millipore) secondary
   antibody for 2 hours at room temperature. Nuclear staining was
   performed by ProLong™ Gold Antifade Mountant with DAPI. Images were
   recorded using a EVOS Cell Imaging Systems Leica (Thermo Fisher
   scientific, Waltham, MA, USA) at 4X and 10X.

Net-like formation assay

   Six-wells Transwell plates, prepared as previously described, were
   evaluated after 7 days. The membrane was removed and HUVEC cells,
   plated in the chamber, were photographed under polarized phase at 10x
   magnification through a Leica DMi1 Microscope (Leica, Wetzlar,
   Germany). Completely enclosed spaces were counted to quantify net-like
   structures^[93]24.

Capillary-like tube formation on matrigel-3D matrix

   Capillaries formation was evaluated in Matrigel (Corning, NY, US)
   according to manufacturer’s instructions. Briefly, 100 µl of Matrigel
   (5 mg/mL) were aliquoted in each 24 well plate and incubated at 37 °C
   for 30 minutes. After trypsinization 4.80 × 10^5/mL of HUVEC were
   suspended in 270 μL Matrigel matrix solution (5 mg/mL) and incubated at
   37 °C for 30 to 45 min. After 12 hours 24 well transwell inserts
   (Sarsted, Nümbrecht Germany) were added to the HUVEC seeded plate and
   incubated 37 °C. After 24 hours of co-culture incubation under
   normoxia, tubes formation was observed using DMi1 microscope (Leica,
   Wetzlar, Germany) at 5x and 10x magnifications.

GFP V2a kit tube formation assay

   Tube formation was also evaluated by a Matrigel-green fluorescent
   protein (GFP) fluorescence coupling system. 24 well plate was coated
   using 5 mg/mL of Matrigel (Corning, NY, US) and incubated at 37 °C for
   30 minutes. HUVEC-GFP cells were plated according to the datasheet
   instructions and after 24 hours 24 well transwell inserts (Sarsted,
   Nümbrecht Germany) were co-cultured for 24 hours. Tubes formation was
   observed under Evos FL auto microscope (Leica, Wetzlar, Germany) using
   a GFP channel at 4x and 10x magnifications and the tubular structures
   formed in the matrigel were counted in 5 random fields.

Experimental design and statistical rationale

   The SWATH-MS-based proteomics analysis was performed on a total number
   of two biological replicates (obtained by seeding and growing the cells
   in separate times) for both 2G8 cells and control A375 cells. Three
   instrumental replicates for each sample were then subjected to the DIA
   analysis. Co-culture experiments with HUVEC cells, as well as RT PCR
   analysis, were performed in three independent experiments with three
   replicates for each sample.

   We used the Student Paired t-test to compare the variation of a
   variable between two groups, while one-way analysis of variance (ANOVA)
   was used for multiple comparisons. Differences were considered
   statistically significant with p < 0.05. Statistical analyses were
   performed using SPSS for Windows, version 16.0 (SPSS Inc, Chicago, IL,
   USA).

Results

Deletion of the Runt domain in A375 melanoma cells

   After applying the CRISPR/Cas9 procedure to A375 cells, we performed a
   single cloning and expansion. Remarkably we obtained only 11 clones.
   Therefore, we analysed the sequences of the PCR products obtained by
   amplifying the Runt domain region (Fig. [94]1A).

Figure 1.

   [95]Figure 1
   [96]Open in a new tab

   Deletion of the Runt domain in A375 melanoma cells. (A) PCR products.
   11 clones have been generated by applying the CRISPR/Cas9 technique.
   The PCR products, analysed by agarose gel electrophoresis, have been
   obtained by amplifying the Runt domain region. (B) Sanger sequencing
   chromatogram. The image shows the comparison between the WT and the 3G8
   sequences with the deletion in genomic context.

   As the sequence in 10 of 11 clones showed none or little changes, such
   as in-frame insertion or deletion (Supplementary Fig. [97]1) compared
   to A375, we chose to investigate the clone 3G8 in which a deletion
   corresponding to 114 bp was observed (Fig. [98]1B). Before further
   experiments, we confirmed by immunoblotting in 3G8 cells the expression
   of the remaining protein at comparable level of the Runx2 full length
   (Fig. [99]2I), as we previously published^[100]11.

Figure 2.

   [101]Figure 2
   [102]Open in a new tab

   Runt domain-dependent effects on the proteome of melanoma cells. (I)
   Immunoblot analysis of Runx2 and of some deregulated proteins
   (indicated by the corresponding gene name) in A375 and 3G8 cells. The
   primary antibodies, along with the sources and working dilutions are
   listed in Supplementary Table [103]1. The blots were cropped to focus
   upon the specific proteins indicated. Amido black-stained PVDF membrane
   is also presented to show equal loading of protein. (II) Functional
   annotation of 3G8 deregulated proteins according to DAVID. The enriched
   cellular component (a), biological process (b), and molecular function
   (c) terms are reported.

Quantitative analysis of Runt domain-dependent effects on the proteome of
melanoma cells

   To assess the global effects of Runt domain deletion from RUNX2, a
   comparative shotgun proteomic analysis of A375 and 3G8 melanoma cells
   was carried out on two biological replicates. We identified a total of
   1809 proteins in the first biological replicate, and 1671 in the second
   one, with 1445 identified proteins shared in common. In total 1169
   protein groups were quantified in the two biological replicates with a
   peptide confidence cut-off of 99% (FDR < 1%); among these 112 were
   deregulated. In particular, 41 proteins were found to be up-regulated
   in 3G8 melanoma cells compared to controls, while 71 were found
   down-regulated (p < 0.05, fold-change cut-off > 1.5). The lists of:
   “identified peptides and their area”, “identified proteins and their
   normalized area”, “quantified and differentially expressed proteins”
   are presented in Supplementary Tables [104]2, [105]3, [106]4
   respectively.

   To validate proteomic data, we performed immunoblot analysis to
   determine the relative amount of up and down-regulated candidates in
   3G8 cells compared to A375 cells (Fig. [107]2I).

   These proteins were selected on the basis of several factors, such as
   fold of variation, novelty and relevance. The results of immunoblot
   experiments were largely consistent with those observed from the
   quantitative proteome analysis.

Modulation of biological processes and pathways induced by inhibition of
Runx2 DNA-binding function in melanoma cells

   To reveal signalling pathways and interaction networks perturbed by
   inhibition of Runx2 DNA-binding function, we performed a bioinformatics
   analysis on the proteins whose expression significantly changed between
   A375 and 3G8 melanoma cells.

   Cellular component analysis annotated more than half of differentially
   expressed 3G8 proteins to the extracellular exosome (51,8%), followed
   by a smaller percentage of annotated protein to cytosol localization
   (49,1%) (Fig. [108]2IIa). In addition, the GO enrichment analysis of
   biological process suggested that the deregulated 3G8 proteins have a
   role in oxidation-reduction (10%), cell-cell adhesion (5,4%) and
   nucleosome assembly (4,5%) (Fig. [109]2IIb). According to the molecular
   function analysis, most of the deregulated proteins in 3G8 cells
   exhibited “binding” activity, in particular protein (64,3%) and poly(A)
   RNA (25%) binding activity (Fig. [110]2IIc).

   The complete list of GO terms for cellular component, biological
   process and molecular function analysis, together with p-value and
   genes, is reported in Supplemental Table [111]5.

   To reveal functional enrichment and deregulated pathways within the
   proteins that exhibited a changed abundance in 3G8 cells we performed
   an IPA analysis (Supplemental Table [112]6). Overall, IPA software
   revealed that proteins modulated in 3G8 cells show a significant link
   with many different pathways. The twenty top canonical pathways
   (p-value ≤ 0.001) are shown in Fig. [113]3. Among them, the “Granzyme A
   signalling” is the most statistically significant one. Notably, an
   overall increase in the activity of “Rho GDI signalling” and “actin
   cytoskeleton pathways” was detected (orange histograms), in addition to
   an overall decrease in the activity of “sirtuin”, “Rho GTPases” and
   “protein kinase A” signalling pathways (blue histograms). Moreover, the
   IPA software allowed us to examine also the potential upstream
   regulators associated with the above-described proteomic profile.
   Interestingly, the most significant upstream regulators include the
   transforming growth factor beta 1 (TGFB1; p = 7.65E-09, z = 0.148),
   tretinoin (p = 6.24E-07, z = 1.437) and mir-122 (p = 5.98E-06,
   z = 1.044), as well as the ligand-dependent nuclear estrogen receptor
   (ESR1; p = 8.30E-05, z = 0.289), the transcription regulators N-myc
   (MYCN; p = 8.86E-05, z = −1.39) and tumor protein p53 (TP53;
   p = 1.12E-04, z = 0.218). In addition, the transcription regulator
   retinoblastoma (RB1; p = 3.40E-03, z = −2.05) is predicted to be
   significantly inhibited. A complete list of transcriptional upstream
   regulators with significant p-values (p ≤ 0.01) can be found in
   Supplemental Table [114]6. To determine the key proteins in the
   function network, the STRING online tool was used to analyse the
   protein-protein interaction networks. The analysis resulted in a
   network (showed in Fig. [115]4) with an average local clustering
   coefficient (i.e. indication of the embeddedness of single nodes) of
   0.374, and an average node degree (i.e. number of interactions, at the
   score threshold, that a protein has on the average in the network) of
   0.9. The low average local clustering coefficient suggests the presence
   of low connected neighborhoods in the network, as indicated by the
   presence of small groups of interacting proteins (2 or 3 proteins). The
   results showed that 5 clusters were determined by the K-means method,
   represented by different colored nodes in Fig. [116]4, with a maximum
   and minimum size of 42 and 3 genes respectively. In particular, the
   green cluster has 4 members and was found to be enriched in histone
   subunits, the yellow cluster has 5 members 3 of which are components of
   the cytoskeleton, while the light-blue cluster has 3 members which play
   a role in ubiquitin pathway

Figure 3.

   [117]Figure 3
   [118]Open in a new tab

   Identification of deregulated canonical pathways in 3G8 cells based on
   IPA. The top 20 significantly altered canonical pathways associated
   with deregulated proteins of 3G8 cells. The y-axis corresponds to the
   −log of the P-value (Fisher’s exact test) and the ratio (orange points)
   the number of genes in a given pathway that meet cut off criteria,
   divided by the total number of genes that map to that pathway. The
   orange and blue colored bars indicate predicted pathway activation, and
   predicted inhibition, respectively (z-score). White bars are those with
   a z-score at or very close to 0. Gray bars indicate pathways where no
   prediction can currently be made.

Figure 4.

   [119]Figure 4
   [120]Open in a new tab

   Protein network of modulated proteins in 3G8 cells by using STRING.
   Schematic view of known and predicted protein-protein interactions
   according to the STRING database (v. 10). The minimum required
   interaction score was set as the medium confidence (0.400) and the
   disconnected nodes in the network were hidden. Pink edges represent
   interactions experimentally determined, while light blue edges
   represent known interaction from curated database. The network map was
   clustered in 5 clusters using the K-means method

Runt domain increases angiogenesis ability in melanoma cells

   In order to prove the involvement of Runt domain in inducing
   neoangiogenesis, as suggested by the proteomic analysis, we analysed
   the expression of VEGF gene in A375 and 3G8 melanoma cells.
   Interestingly, as shown in Fig. [121]5A we found a downregulation of
   VEGF expression in 3G8 cells (p < 0.01) as compared to A375 cells.
   Therefore, we investigated the expression of CD105 and CD31, the
   markers of tumour neoangiogenesis, in HUVEC co-cultured with NHEM, A375
   or 3G8 melanoma cells. The data obtained showed that HUVEC cells
   express lower levels of CD105 and CD31 mRNAs and reduced expression of
   CD31 protein when co-cultured with 3G8 as compared to A375 cells
   (Fig. [122]5B).

Figure 5.

   [123]Figure 5
   [124]Open in a new tab

   Runt domain promotes angiogenesis in melanoma cells. (A) Real-Time PCR
   analysis shows decreased gene expression of VEGF in 3G8 melanoma cells
   as compared to A375 control cells. (B) Real-time PCR (left) and
   immunofluorescence (right) respectively show significant
   down-regulation of CD31 and CD105 mRNAs and reduced CD31 protein
   expression in HUVEC cells co-cultured with 3G8 cells compared to that
   co-cultured with A375 cells. (C) In the transwell co-culture setup
   HUVEC cells were cultured as a monolayer in the chamber and NHEM, A375,
   or 3G8 cells were plated in the upper well. Quantification (left) and
   representative images (right) of the tube formation assay are shown;
   magnification 10x. (D) Representative phase contrast and
   epifluorescence images of capillary-like structure formation in HUVEC
   (upper panels) and HUVEC-GFP (bottom panels) co-cultured with A375 or
   3G8 cells on matrigel-3D matrix; magnification 40x (left). The branch
   points of the capillary-like tubes were counted and the quantitative
   data are presented (right). Statistical significance was calculated
   based on three independent experiments (*p < 0.05; **p < 0.01, p-values
   between depicted groups). NHEM = normal melanocytes.

   Moreover, further confirming the role of Runt domain in
   neoangiogenesis, we detected by transwell co-culture setup a reduced
   number of network-like structures in HUVEC co-cultured with 3G8
   (p < 0.01, Fig. [125]5C) as compared to A375 cells. The pro-angiogenic
   effect of Runt domain was also assessed by capillary-like tube
   formation assays in the matrigel-3D matrix. We observed that HUVEC
   co-cultured with 3G8 cells had a significant reduction in the capillary
   forming features with few connecting structures when compared with A375
   cells (p < 0.05, Fig. [126]5D). All these findings support the role of
   the Runt domain in promoting the neoangiogenesis.

TCGA analysis indicates that angiogenic proteins modulated in 3G8 cells are
altered in melanoma patients

   Upon querying the TCGA provisional database for skin cutaneous
   melanoma, consisting of 471 patients, our results indicate that our
   target list of proteins associated to angiogenesis (that depend on Runt
   domain) was found to be altered in 51.36% of the total patients
   (Fig. [127]6A). In particular, FASN was found altered in 17.95% of
   patients, while CLIC4, HSPB1, ARHGEF1, PHGDH, MYO1C, CAV1 were found
   altered in 6.05%, 9.6%, 2.92%, 9.6%, 7.72%, 15.87% of patients
   respectively. Using TCGA cohort survival data in univariate analysis,
   alterations in all these genes were not associated with reduced patient
   overall survival (Fig. [128]6B), while there was significantly lower
   overall survival among melanoma patients with alterations of MYO1C gene
   (Logrank test p-value: 0.0019) (Fig. [129]6C).

Figure 6.

   [130]Figure 6
   [131]Open in a new tab

   cBioPortal query of the TCGA provisional database for skin cutaneous
   melanoma dataset for alterations in genes of the list of targets
   implicated in angiogenesis (i.e. FASN, CLIC4, HSPB1, ARHGEF1, PHGDH,
   MYO1C, CAV1) and correlation with overall survival. (A) Combined
   alterations of 7 genes for our assay targets were seen in 242 cases
   (51.36%) of melanoma in the TCGA dataset. (B) Patients with alterations
   in these genes had no a statistically significant lower survival when
   compared to cases with no alterations, while (C) patients with
   alterations in MYO1C gene had a statistically significant lower
   survival when compared to cases with no alteration.

Discussion

   Our previous study showed that Runt domain KO melanoma cells are
   characterized by reduced proliferation, epithelial-mesenchymal
   transition and metastasis, suggesting that the Runt domain of Runx2 may
   represent an oncotarget in melanoma^[132]11. In the present study, with
   the aim of deepening the understanding of Runt domain-dependent effects
   on melanoma cells biology, we carried out a proteomic analysis of A375
   and del-Runt (3G8) melanoma whole-cell extracts to investigate the
   intracellular molecular mechanisms characterizing the deletion of the
   DNA-binding domain from RUNX2.

   The data obtained revealed the regulation of some key pathways in 3G8
   cells. One of the interesting observations was that the most
   significant deregulated pathways in cells with the deletion of Runt
   domain of Runx2 is the “Granzyme A signalling”. The identified proteins
   of 3G8 cells involved in this signalling are all downregulated histone
   subunits (Supplemental Table [133]6) which are typical targets of
   granzyme A degradation^[134]25. Granzyme A is a selective protease that
   triggers a form of cell death characterized by all of the morphological
   features associated with apoptosis (membrane blebbing, chromatin
   condensation and nuclear fragmentation) without activating the caspase
   cascade. Accordingly, in our previous study, we observed an increased
   apoptosis in Runt domain KO cells; therefore, the downregulation of
   this protease corroborates the negative role of the Runt domain in
   inducing apoptosis in melanoma cells. Interestingly, it has been
   reported that Runx3 regulates the expression of granzyme B^[135]26, and
   that RUNX2^fl/fl mice have defects in granzyme B expression^[136]27. To
   the best of our knowledge, our finding represents the first indication
   that a link between Runx2 and Granzyme A may also exist. Consistent
   with our previous study^[137]11 the other data obtained by canonical
   pathway enrichment analysis confirmed the effect of the Runt domain of
   Runx2 on cell migration, metastasis and proliferation of tumour cell.
   Indeed, among the top enriched and induced pathways we found the
   “signalling of RhoGDI” a guanine nucleotide dissociation inhibitor
   (GDI) specific for the “Rho family of small GTPases”, a signalling
   instead detected as inhibited in 3G8 cells. Emerging evidences suggests
   that RhoGDI (activated in 3G8 cells) plays dual opposite roles in
   tumour progression^[138]28, while Rho GTPases (inhibited in 3G8 cells)
   act as promoters of cell migration^[139]29. Also the “sirtuin
   signalling pathways” which promotes melanoma metastases^[140]30, and
   the “protein kinase A signalling” implicated in the initiation and
   progression of tumours^[141]31, resulted to be inhibited pathways in
   3G8 cells. In particular, among the 3G8 cells deregulated proteins that
   we identified by MS and immunoblotting analyses, there are the
   pro-apoptotic BH3-interacting domain death agonist (BID, +1.84 fold
   change), as well as Ras GTPase-activating-like protein (IQGAP1, −1.54
   fold change) which is involved in melanoma cell migration^[142]32, and
   annexin-VI (ANXA6, +1.93 fold change) which acts as a multifunctional
   scaffold in cell motility^[143]33 and represents a tumour suppressor in
   melanoma^[144]34.

   Interesting data were also suggested by the upstream regulator analysis
   performed in silico by IPA (Supplemental Table [145]6). This analysis
   indicated that the most significant upstream regulators is TGF-β1, of
   which Runx2 represents one of the major targets^[146]35. We found
   interesting that the modulation of some of the TGF-β1 targets
   identified in this study, such as FASN, ANPEP, CLIC4, HSPB1, ARHGEF2,
   PHGDH, MYO1C and CAV1, are involved in processes that counteract
   angiogenesis. Interestingly, the involvement of Runx2 in regulation of
   angiogenesis and vasculogenic factors (such as VEGF) has already been
   demonstrated during bone development^[147]36 as well as in
   cancer^[148]37–[149]41, but never in melanoma.

   Notably, TGF-β1 markedly decreases the expression of fatty acid
   synthase^[150]42 which we found downregulated in 3G8 cells (−1.95 fold
   change). This enzyme besides playing a pivotal role in lipid
   metabolism, is also implicated in the induction of
   angiogenesis^[151]43, in particular in cancer cells^[152]44.
   Accordingly, FASN inhibitors play a key role in antiangiogenic
   treatments^[153]45; in particular, they reduce the tumour cell-mediated
   formation of HUVEC capillary-like structures in melanoma^[154]46.
   Accordingly, we found a reduced number of network-like structures in
   HUVEC cocultured with 3G8 compared to cocultures with WT melanoma cells
   (Fig. [155]5C). We also assed HUVEC capillary-like tube formation when
   in 3D co-culture with 3G8 or A375 cells using both a Matrigel (which
   provides similar conditions to that observed in the tumour
   microenvironment) as well as a Matrigel-GFP fluorescence coupling
   system (which allowed direct visual evaluation of pro-angiogenic effect
   of melanoma cells via comparing fluorescent tubular structures formed
   by GFP-expressing HUVECs, while obviating noise effect resulting from
   any nonspecific cell attachment to HUVEC targets). As depicted in
   Fig. [156]5D the formation of capillary-like structures on 3D-matrix
   was significantly reduced in HUVEC cells co-cultured with 3G8 as
   compared to A375 cells. Among the downregulated TGF-β1 targets that we
   identified in 3G8 cells there is also aminopeptidase N (ANPEP/CD13,
   −1.59 fold change) a ubiquitously expressed membrane peptidase. ANPEP
   is involved in melanoma angiogenesis and in melanoma cell
   invasion^[157]47. In particular, the inhibition of ANPEP in melanoma
   cells correlates with the anti-tumour angiogenesis effect of
   bestatin^[158]48. ANPEP is also an important mediator of resistance to
   inhibition of BRAF, one of the most aggressive oncogenes found in
   melanoma which modulates angiogenesis^[159]49. The proteomic analysis
   also suggested the downregulation of chloride intracellular channel
   protein-4 (CLIC4, −1.82 fold change), a p53 and TGF-β regulated protein
   which promotes angiogenesis supporting acidification of vacuoles along
   the intracellular tubulogenic pathway^[160]50 and induces tubular
   morphogenesis^[161]51. Emerging evidences indicate that targeting CLIC4
   could represents a strategy to diminish some of the tumour enhancing
   effects of the cancer stroma^[162]52 and could further suppress the
   invasion in melanoma cells^[163]53. We also showed the downregulation
   in 3G8 cells of heat shock protein beta-1 (HSPB1, −1.55 fold change), a
   molecular chaperone highly expressed in many cancers^[164]54. It has
   been reported that by combining HSPB1 silencing and BRAF inhibition
   melanoma cells are fully committed to death^[165]55 and, accordingly,
   reduced level of HSPB1 correlates with a less aggressive phenotype and
   improved survival in patients with melanoma^[166]56. Interestingly, the
   downregulation of HSPB1 is associated to reduced endothelial cell
   proliferation and decreased secretion of key mediator of angiogenesis
   (such as VEGF-A, VEGF-C, and basic fibroblast growth factor)^[167]56.
   Accordingly, we observed a lower gene expression of VEGF-A in KO Runt
   cells (3G8) as compared to A375 melanoma cells (Fig. [168]5A). We also
   detected a reduced number of network-like structures, as well as
   reduced expression of CD105 and CD31 endothelial markers in HUVEC cells
   co-cultured with 3G8 as compared to that co-cultured with A375 melanoma
   cells (Fig. [169]5B,C). CD105 and CD31 are angiogenesis-associated
   molecules and they are considered direct markers of the levels of
   neoangiogenesis. In addition, CD105 and CD31 expression indicates the
   density of intra-tumorous vessels^[170]57 and it is considered
   important prognostic factor in cancer^[171]57,[172]58, Therefore, our
   results strongly suggest that Runt domain promotes the neoangiogenesis
   in melanoma.

   We also showed in 3G8 cells downregulation of Rho guanine nucleotide
   exchange factor 1 (ARHGEF1, −1.50 fold change) and 2 (ARHGEF2, −2.70
   fold change) which are fine-tuners of angiogenesis and vascular
   function^[173]59; D-3-phosphoglycerate dehydrogenase (PHGDH, −1.89 fold
   change), a key enzyme of the serine synthesis pathway whose knockdown
   inhibits angiogenesis in vitro and in vivo (as a result of decreased
   nucleotide and heme synthesis)^[174]60; as well as of unconventional
   myosin-Ic (MYO1C, −1.81 fold change) whose depletion resulted in an
   anti-angiogenic effect due to reduced cell-surface and total cellular
   level of VEGFR2 the predominant receptor in angiogenic
   signalling^[175]61. Finally, we also found upregulation of Caveolin-1
   (CAV1, +2.11 fold change) as a result of the deletion of the Runt
   domain of Runx-2 in 3G8 cells. CAV1 has a dual role in
   melanoma^[176]62: in advanced stages its expression is associated with
   enhanced metastatic potential, while, at earlier stages, it functions
   as a tumour suppressor^[177]63. Interestingly, some researchers showed
   a relationship between its induction and reduced angiogenesis^[178]64.
   Of note, we also checked cBioPortal and found that these proteins
   related to angiogenesis (i.e. FASN, CLIC4, HSPB1, ARHGEF1, PHGDH,
   MYO1C, CAV1) are altered in a total of 242 melanoma patients out of a
   total of 471 (51.36%). The most common alteration found in patients for
   all these targets is the mRNA transcriptional upregulation of their
   encoding genes, hence the fact that the deletions of Runt domain leads
   to a downregulation at protein level for all of them (with the
   exception of CAV1 which is upregulated) seems particularly interesting
   in confirming the key role of Runx2 in the development of melanoma.
   Moreover, by using cBio Cancer Genomics Portal, the overall survival
   rates of patients with melanoma were compared between tumour samples
   with or without alterations in each and in all these proteins. However,
   the results indicate that our list of proteins involved in
   angiogenesis, whose expression is connected to Runt domain, do not
   constitute a “signature” to predict survival of patients with melanoma.

   Our approach provides a comprehensive view of molecular details that
   depend on Runx2 DNA binding function in melanoma cells. It should
   however be emphasized that the findings here reported are related to
   the consequences of deleting the Runt domain of Runx2, and not to Runx2
   full knockout. Thus, it is possible that some alterations are not
   directly caused by lack of the Runt domain, but for instance, they may
   involve the remaining non-Runt domains of the protein which may still
   interact with regulatory proteins^[179]65. Indeed, the mutant Runx2
   protein in 3G8 cells could still alter gene expression through indirect
   mechanisms, which are still relevant in understanding the role that
   Runx2 plays in melanoma. In conclusion, these findings support the idea
   that the Runt domain of Runx2 is involved in increasing the risk of
   tumour proliferation/migration and angiogenesis, and may be an ideal
   oncotarget in melanoma cancer cells.

Supplementary information

   [180]Supplementary Info (Table 1, Fig 1)^ (407.3KB, pdf)
   [181]41598_2019_44552_MOESM2_ESM.xlsx^ (68.6MB, xlsx)

   Supplemental_Table_2 list of peptides area and of identified peptides
   [182]41598_2019_44552_MOESM3_ESM.xlsx^ (2MB, xlsx)

   Supplemental_Table_3 list of normalized protein area and of identified
   proteins.xlsx
   [183]41598_2019_44552_MOESM4_ESM.xls^ (408.5KB, xls)

   Supplemental_Table_4 List of differentially expressed proteins.xls
   [184]41598_2019_44552_MOESM5_ESM.xls^ (41KB, xls)

   Supplemental_Table_5 GO enrichment analysis.xls
   [185]Supplemental_Table_6 IPA analysis.xls^ (411.5KB, xls)

Acknowledgements