Multi-Omics Analysis of Glioblastoma and Glioblastoma Cell Line: Molecular Insights Into the Functional Role of GPR56 and TG2 in Mesenchymal Transition

G protein-coupled receptor 56 (GPR56/ADGRG1) is an adhesion GPCR with an essential role in brain development and cancer. Elevated expression of GPR56 was observed in the clinical specimens of Glioblastoma (GBM), a highly invasive primary brain tumor. However, we found the expression to be variable across the specimens, presumably due to the intratumor heterogeneity of GBM. Therefore, we re-examined GPR56 expression in public domain spatial gene expression data and single-cell expression data for GBM, which revealed that GPR56 expression was high in cellular tumors, infiltrating tumor cells, and proliferating cells, low in microvascular proliferation and peri-necrotic areas of the tumor, especially in hypoxic mesenchymal-like cells. To gain a better understanding of the consequences of GPR56 downregulation in tumor cells and other molecular changes associated with it, we generated a sh-RNA-mediated GPR56 knockdown in the GBM cell line U373 and performed transcriptomics, proteomics, and phospho-proteomics analysis. Our analysis revealed enrichment of gene signatures, pathways, and phosphorylation of proteins potentially associated with mesenchymal (MES) transition in the tumor and concurrent increase in cell invasion and migration behavior of the GPR56 knockdown GBM cells. Interestingly, our analysis also showed elevated expression of Transglutaminase 2 (TG2) - a known interactor of GPR56, in the knockdown cells. The inverse expression of GPR56 and TG2 was also observed in intratumoral, spatial gene expression data for GBM and in GBM cell lines cultured in vitro under hypoxic conditions. Integrating all these observations, we propose a putative functional link between the inverse expression of the two proteins, the hypoxic niche and the mesenchymal status in the tumor. Hypoxia-induced downregulation of GPR56 and activation of TG2 may result in a network of molecular events that contribute to the mesenchymal transition of GBM cells, and we propose a putative model to explain this functional and regulatory relationship of the two proteins.

G protein-coupled receptor 56 (GPR56/ADGRG1) is an adhesion GPCR with an essential role in brain development and cancer. Elevated expression of GPR56 was observed in the clinical specimens of Glioblastoma (GBM), a highly invasive primary brain tumor. However, we found the expression to be variable across the specimens, presumably due to the intratumor heterogeneity of GBM. Therefore, we re-examined GPR56 expression in public domain spatial gene expression data and single-cell expression data for GBM, which revealed that GPR56 expression was high in cellular tumors, infiltrating tumor cells, and proliferating cells, low in microvascular proliferation and peri-necrotic areas of the tumor, especially in hypoxic mesenchymal-like cells. To gain a better understanding of the consequences of GPR56 downregulation in tumor cells and other molecular changes associated with it, we generated a sh-RNA-mediated GPR56 knockdown in the GBM cell line U373 and performed transcriptomics, proteomics, and phospho-proteomics analysis. Our analysis revealed enrichment of gene signatures, pathways, and phosphorylation of proteins potentially associated with mesenchymal (MES) transition in the tumor and concurrent increase in cell invasion and migration behavior of the GPR56 knockdown GBM cells. Interestingly, our analysis also showed elevated expression of Transglutaminase 2 (TG2) -a known interactor of GPR56, in the knockdown cells. The inverse expression of GPR56 and TG2 was also observed in intratumoral, spatial gene expression data for GBM and in GBM cell lines cultured in vitro under hypoxic conditions. Integrating all these observations, we propose a putative functional link between the inverse expression of the two proteins, the hypoxic niche and the mesenchymal status in the tumor. Hypoxia-induced downregulation of GPR56 and activation of TG2 may result in a network of molecular events that contribute to the mesenchymal transition of GBM cells,

INTRODUCTION
Glioblastoma (GBM) is one of the most aggressive forms of primary brain tumors, with a poor prognosis of 12-15 months and just 3-5% of survival over five years (1). The dismal prognosis is attributed mainly to the complex inter-and intratumor heterogeneity, contributing to drug resistance (2,3). Large-scale transcriptomic studies to characterize the heterogeneity of GBM led to the identification of molecular subtypes of GBM, namely proneural (PN), classical (CL), and mesenchymal (MES) (3,4). The PN and MES subtypes have been more consistently identified, with PN relating to a more desirable outcome and MES to poor survival (5,6). Recent work on GBM heterogeneity reports that molecular features of histologically defined anatomical areas of the tumor (leading edge (LE), infiltrating tumor (IT), microvascular proliferation (MVP), cellular tumor (CT), and pseudo-palisading cells around necrosis (PAN) are highly distinct and conserved regardless of whether they are derived from the same or different patients (7). Also, a recent analysis based on single-cell RNA sequencing revealed that GBM cells exhibit four cellular states (neuralprogenitor like (NPC-like), oligodendrocyte progenitor like (OPC-like), astrocyte like (AC-like), and mesenchymal-like (MES-like) and a single malignant cell can generate all four states (8,9). Thus, such inter-/intra tumoral heterogeneity and cellular plasticity may have significant implications for treatment and likely outcomes.
Adhesion G protein-coupled receptors (aGPCRs) are the second largest family of G protein-coupled receptors (GPCRs) involved in many cellular functions and serve as drug targets for various clinical conditions (10,11). aGPCRs are characterized by an unusually long extracellular domain that interacts with other cells as well as mediates interactions between cells and the extracellular matrix (ECM) (12), seven transmembrane domains, and an intracellular C-terminus. In recent years, GPR56 has emerged as an important aGPCR with a multitude of functions in health and disease. GPR56 was originally discovered for its role in genetic disorder bilateral frontoparietal polymicrogyria (BFPP) (13,14), resulting in severe cortical damage malformation, thus implying a role in brain development. Further, GPR56 is also implicated in neuronal myelination and myelin repair (15,16), Oligodendrocyte precursor cell (OPC) proliferation (17), radial axonal sorting by Schwann cells (16), and others. Altered expression of GPR56 has been reported in multiple cancers, including glioblastoma, melanoma, breast cancer, and colon cancer (18)(19)(20)(21)(22). Its expression has been found to influence cell adhesion, migration, and proliferation in a variety of cancer cell types as well as epithelial-mesenchymal transition (EMT) and radio resistance (19,(23)(24)(25)(26)(27). Whereas these reports clearly imply GPR56 role in cancer, it seems to function as an oncogenic factor in some (19,28) and as a tumor suppressor in others (26,29).
Based on the earlier studies in brain development and melanoma, collagen III and TG2, respectively, have been identified as key interactors of GPR56 (13,30). While GPR56 interacts with collagen III during brain development, its interaction with TG2 plays a role in both developing brain (OPC proliferation) and cancer. TG2 was discovered to be the ECM-ligand of GPR56 in melanoma tissues where it was shown that the binding of TG2 to the receptor initiated internalization and lysosomal degradation of TG2, resulting in decreased fibronectin deposition in the ECM (13,15,30). The interaction of TG2-GPR56 suggests an inhibitory role of GPR56 in the progression of melanoma (27). TG2 is a multifunctional protein, present intracellularly, on the cell surface as well as in the ECM. TG2 crosslinks fibronectin and other adhesion proteins and, in some instances modulates their interaction with integrins on the cell surface with a role in cell migration, growth, and differentiation (31). TG2 can drive cancer stem cell survival and tumor formation either via ECM modifications that activate stromal cells or by modulation of integrin signaling inhibiting the hippo pathway in tumor cells (32). In GBM, high levels of TG2 expression are associated with aggressive MES phenotype, and a recent study suggested that TG2 inhibition prevents differentiation into mesenchymal subtype (33). GPR56 is upregulated in GBM and other astrocytomas, as we and others have reported (20,34,35). However, we observed that GPR56 expression is heterogeneous across GBM specimens and within the tumor, which was also supported by the observation of Moreno et al., who showed that while GPR56 is highly expressed in the PN subtype of GBM, it is downregulated in the MES subtype (26). Given its role in the developmental processes and implications in cancer, it forms an attractive target to study. Thus it is important to delineate the mechanistic basis of these varied influences of the receptor and the pathways involved in normal development and cancer.
In this study, we tried to understand the heterogeneous expression and functional significance of GPR56 in GBM using public domain data and experimental multi-omics data, including transcriptomics, proteomics, and phosphoproteomics, of GPR56 knockdown U373 GBM cells. Our study suggests a role of GPR56 and TG2 interaction in PN to MES transition in GBM and provides insights on the putative molecular events involved, rendering this interaction a potential therapeutic target.

Sample Collection and Processing
Glioblastoma samples were obtained from patients at Mazumdar Shaw Medical Center, Bengaluru, India, during surgery, with informed consent. All procedures were carried out in compliance with the recommended protocols and with the approval of the Institutional Ethics Committee, NHH/MEC-CL-2015/384(A). The samples were snap-frozen and formalin-fixed, and only tumors histopathologically classified as GBM as per WHO guidelines (2016) were used. Brain tissue specimens obtained from temporal lobe epilepsy surgeries were used as experimental controls (Supplementary Table 1), were obtained from the brain bank at the National Institute of Mental Health and Neurosciences (NIMHANS), Bangalore, India.

qRT-PCR
Expression of individual genes at the RNA level in cells or tissue specimens was assessed using qRT-PCR. For this purpose, total RNA was isolated from the cells or tissues using Macherey Nagel Total RNA and Protein Isolation Kit (Cat#: 740933.50, Macherey-Nagel, Germany). The extracted RNA was then converted into cDNA using a High-Capacity cDNA Reverse Transcription Kit (Cat#: 4368814, Applied Biosystems, Lithuania) as per the manufacturer's instructions. qRT-PCR was carried out using the KAPA Biosystem's SYBR Fast qPCR universal master mix (2X) (Cat#. KK4602, KAPA Biosystems, MA, USA) and the Roche Light Cycler 480 Real-time PCR system (Roche Diagnostics, Germany). The Second derivative maximum (2 −DDCq ) method was used to normalize the qPCR results (36), and relative change in expression was evaluated using the geometric mean of the selected reference genes (18S and RAB7A). Primers used for the qRT-PCR reaction are given in Table 1.
Western blot was performed on whole-cell protein extracts, transferred to 0.45 µ m pore size PVDF membrane. The membrane was blocked for 1h in a blocking buffer containing 3% BSA. The primary antibodies were diluted in the blocking buffer as indicated below and incubated with the membrane overnight. The blot was then washed three times with TBST (Tris-buffered saline and Tween 20%) and incubated with respective secondary antibodies in the blocking buffer. This was followed by three washes with TBST, followed by two washes with TBS, and developed using an ECL kit. Beta-actin and GAPDH were used as loading control of TG2 and GPR56, respectively. The details of the antibodies used in the analysis are given in Table 2 below.

Public Domain Data Retrieval and Analysis
Single-cell data: Neftel et al. (8) generated and analyzed singlecell RNA-seq data from 28 pediatric and adult glioblastoma tumors to identify four major neoplastic cell types defined by six gene modules; 1. Mesenchymal -Hypoxia independent (MES1like) and hypoxia dependent (MES2-like) mesenchymal related gene sets, 2. Astrocytic -astrocytic (AC-like) marker gene set, 3. Oligodendroglial progenitor-like (OPC-like) lineage marker gene set and 4. Neural -stem and progenitor cell gene set (NPC1-like and NPC2-like) as well as two cell cycling modules, namely G1S and G2M (CC). The processed data, which is available in Transcripts per Million (TPM) format at the Broad Institute Single-Cell Portal (GEO -SCP393; https://www.broadinstitute. org), was downloaded and used to study the expression profiles of GPR56 across GBM meta-modules. Ivy Glioblastoma Atlas Project (Ivy GAP): In the Ivy GAP analysis (7) (https://glioblastoma.alleninstitute.org), laser capture microdissection of the GBM tumor areas followed by RNA-seq analysis was used to generate molecular signatures of cells present in five major anatomic features of GBM visible by H and E staining, namely, the LE, IT, MVP, NE, and PAN The processed Ivy GAP quantitation data was downloaded from the website and log2 transformed. Where required, the logtransformed data were z-scored using the following formula: z=X−m where X is the expression measure of a gene in a sample, m is the mean expression measure of a gene across all the samples, and s is the standard deviation of expression measure of a gene across all the samples.

Generating GPR56 shRNA Mediated Knockdown Cells
SureSilencing shRNA plasmids targeting GPR56 and a nontargeting control were purchased from SABiosciences (Qiagen) and contained the following targeting sequences: shGPR56.2(21 bp): GAACCGACATGCTGGGAGATT and shNON-CODING (21 bp): GGAATCTCATTCGATGCATAC. Plasmids were linearized with NaeI (NEB) for 3h at 37°C, the enzyme inactivated by incubation at 65°C for 20 mins. The linearized DNA was purified using phenol-chloroform extraction, followed by DNA precipitation with isopropanol. Reverse transfection of 8x10 5 U373 cells/well was performed using 0.4 µg linearized shRNA plasmid DNA with Fugene6 (Roche) in Opti-MEM using a 24 well plate. Cells were incubated for 24h and then reseeded in a six-well plate at 5x10 4 cells/well in DMEM with 10% FBS containing 200µg/ml hygromycin (Selection medium). The selection medium was initially removed daily until only a few viable cells remained in each well. Once colonies were established, cloning cylinders (Sigma-Aldrich, USA) were placed around single colonies on the drained plate and sealed with 1% (w/v) low-melting agarose (Gibco, USA). Cells were then treated with accutase, removed, and seeded into a fresh 24 well plate supplemented with the selection medium.

Confocal Microscopy
GPR56 control or knockdown cells were plated onto poly-Llysine coated glass coverslips in 24 well plates at 1x10 4 cells/well in DMEM with 10% FBS for 72 h. Cells were washed with PBS, fixed with 4% paraformaldehyde in PBS for 5 mins to conserve epitopes. Cells were washed three times with PBS, blocked with 1% BSA in PBS, for 30 mins prior to incubation with 1µg/ml anti GPR56 antibody (R&D Systems) for 2h. Cells were washed three times with PBS, incubated for 1h with 4µg/ml donkey anti-sheep Alexa Fluor 568 (Invitrogen, USA) in 1% BSA/PBS, and washed three times with PBS for 5 min. Nuclei were counterstained with a DAPI-containing vectashield mounting medium (Vector Laboratories, USA). Cells were analyzed using a Leica SP5 confocal microscope (Leica Microsystems, Germany) using a 63x oil immersion.

Invasion Assay
For invasion assay, we used modified two-chamber plates with an 8-mm pore size (Cat#. CLS3464, Sigma-Aldrich, USA) coated with 1 mg/ml matrigel (Cat#. E1270, Sigma-Aldrich, USA). GPR56 control and knockdown cells (10 4 cells) were added in serum-free media onto the top chamber. In the lower chamber, complete media was used as a chemoattractant. After incubating for 24 hours, the non-invading upper chamber cells were removed by gentle wiping with a cotton-tipped swab. Invading cells on the lower surface of the membrane were fixed with 4% formaldehyde in PBS for 10min, permeabilized with methanol for 20min, stained with 0.4% crystal violet (Cat#. V5265, Sigma-Aldrich, USA) for 15min, counted, and photographed at a magnification of 10X. The fold increase in invasion was calculated by dividing the total number of invading cells from the GPR56 knockdown cells group by those from the GPR56 control group.

Wound Healing Assay
GPR56 control and GPR56 knockdown cells were cultured in a 24 well plate until they attained 80% confluence state; serum-starved for 24 hours and wounded using disposable 200ul pipette tip. Cell migration to the cell-free area was observed and imaged every 24h for 2 days at 10X magnification using a phase-contrast microscope. The migration distance was measured using the integrated Carl Zeiss software (Zeiss, Germany) as follows: The wound area (in pixels) was measured by taking the average of the areas (n=3) delineated by the wound boundaries. The percentage of wound closure is calculated using the formula (A0-An/A0) *100, where A0 is the average area at the 0 hours and An is the average area at the nth hour. Results were derived from three independent experiments, each performed in triplicates.

RNA-Seq Analysis
Total RNA isolated as described above was quantified and assessed for quality using Agilent BioAnalyzer to ensure that all samples had an RNA integrity number (RIN) of 7 or more. The poly-A enriched RNA library's construction was performed according to the manufacturer's protocol using the NEB Ultra RNA-seq Library Prep kit protocol (NEB, USA). All libraries were quantitated using Qubit High Sensitivity Assay (Invitrogen, USA), and RNA-seq was carried out using Illumina HiSeqX by a 151 bp paired-end sequencing as per the manufacturer's instructions to generate 30M, 100bp paired-end reads. Using the STAR (37) algorithm with default parameters, the sequenced transcriptome was aligned to the hg19 reference genome. Gene expression was quantified using the ENSEMBL reference with bedtools (38). For count data normalization, the Count per Million (CPM) method was applied using the following formula: CPM = (count/sum (count)*1000000). The CPM data was log2 transformed after adding a pseudo count of 1 for further analysis. A two-fold change in expression was considered to identify differentially expressed genes. Since we used single sample-pair for analysis, we used a linear regression model and Prediction Interval approach to identify differentially expressed genes and to assess the significance. The z-score for each gene was calculated using the following formula: Pathway analysis: To know the significance of differentially expressed genes we categorized them by GO-molecular function, and biological pathways using the Protein ANalysis THrough Evolutionary Relationships (PANTHER v13.1, pantherdb.org) (39,40).
Single-sample gene set enrichment analysis (ssGSEA): ssGSEA was run for PN and MES sub-type gene signatures reported by Verhaak et al. (3) to develop an enrichment score. The low/high GPR56-associated signature was defined by differentially expressed genes between GPR56 knockdown and control cells.

Sample Preparation and iTRAQ Labeling
GPR56 control and knockdown cells were grown to 70% confluence, starved in serum-free medium for 12 h, and then lysed in cell lysis buffer (2% SDS in 50 mM triethyl ammonium bicarbonate (TEABC) with sonication. Protein concentration was estimated using the BCA method (Pierce; Waltham, MA, USA). 200 µg of protein from the GPR56 control or knockdown cells were reduced using 5mM DTT at 60°C for 20 minutes. Subsequently, alkylation was carried out using 15mM iodoacetamide for 15 minutes in the dark. After reduction and alkylation, the proteins were precipitated with ice-cold acetone and incubation at −20°C overnight. For enzyme digestion of the proteins, sequencing grade trypsin reconstituted in 50 mM Triethyl ammonium bicarbonate was added to the dried protein in the ratio of 1:20 (trypsin: protein in µg), and the digestion was carried out at 37°C for overnight. Trypsin digested peptides were then subjected to iTRAQ labeling using an iTRAQ 8-plex kit (AB SCIEX Pte Ltd., USA) as per the manufacturer's instructions. Labeling tag details are as follows: GPR56 control sample batch1 with 114, GPR56 knockdown batch1 sample with 115, GPR56 control sample batch 2 with 118, and GPR56 knockdown batch 2 with 119 tags. Reactions were quenched with 10 mM glycine. All the four labeled samples were combined, desalted using C18 StageTip and vacuum dried.

Phosphopeptide enrichment
For total proteome analysis, 10 percent (40 µg) of the pool was used and subjected to MS analysis in duplicates, as described below. Remaining pool -90 percent (360 µg) was used to enrich phosphopeptides using a metal affinity-based Phosphopeptide Enrichment Kit (Pierce, Thermo Fisher Scientific, USA). Briefly, the dried peptides were dissolved in 150 ml of binding buffer. TiO 2 beads were washed twice with washing buffer, and a total of 300 µg of tryptic peptide solution was incubated with an appropriate amount (tryptic peptide: TiO 2 = 1:1, w/w) of TiO 2 beads by end-over-end rotation at room temperature for 30 min. The phosphopeptide-bound beads were collected by brief centrifugation, washed twice with 500 ml washing buffer, and transferred to a C18 StageTip (Thermo Fisher Scientific) placed on the top of a 1.5-ml centrifuge tube and was centrifuged to remove the wash buffer, and phosphopeptides were collected from the resin with elution buffer. The eluents were dried and stored at −80°C until further LC-MS/MS analysis.

Mass Spectrometry (MS)
The tryptic peptides or the phosphopeptide fraction were subjected to LC-MS/MS analysis on Orbitrap Fusion Tribrid mass spectrometer interfaced with Easy nano-LC II (Thermo Scientific, Bremen, Germany) and were analyzed in duplicates. The peptides were first loaded on a preanalytical column (2cmx75µm, Magic C18 Aq) (Michrom Bioresources, Inc.) using solvent A (0.1% formic acid). Peptides were then resolved on an analytical column (50cm x75µm, Magic C18 Aq) using a gradient of 5-38% of solvent B (95% acetonitrile, 0.1% formic acid) at a flow rate of 280 nL/min for 120 min. The data-dependent acquisition of MS spectra in the range of 400-1600 m/z was carried out using Orbitrap mass analyzer with a mass resolution of 120,000 at MS level and 60,000 at MS/MS level; higher energy collision dissociation was selected for fragmentation with 37% normalized collision energy. The automatic gain control was set to 2x10 6 ions for full MS and 1x10 6 ions for MS/MS. Internal calibration was executed using lock mass from ambient air (m/z 445.1200025).

Data Analysis
Protein identifications and quantifications of differentially expressed proteins were carried out as follows. The MS data was analyzed using Proteome Discoverer (PD; Thermo Fisher Scientific, Version 2.2). MS/MS search was carried out using the SEQUEST search engine against the NCBI RefSeq database version 89 (containing 425211 entries). Search parameters included trypsin as an enzyme with 1 missed cleavage allowed; precursor and fragment mass tolerance were set to 20 ppm and 0.1 Da, respectively; methionine oxidation was set as a dynamic modification while methylthio modification at cysteine and iTRAQ modification at N-terminus of the peptide was set as static modifications. Signal to noise ratio applied was 1.5 or more. Peptide identifications were obtained by setting a target FDR threshold of 1% at the peptide level, using a decoy database. Protein abundance values obtained from the PD output, are based on the ratios of relative intensities of the iTRAQ reporter ions from the control and knockdown cells, released during MS/MS fragmentation of each peptide. The intensities were checked and found to be conformed to less than 40% coefficient of variation.
The abundance values of proteins or the phosphopeptides for phosphoproteome experiments of the control and knockdown were normalized by dividing the abundance values of the proteins or the phosphopeptides (in phosphoproteome analysis) in each column by the column mean. The normalized values were then log2 transformed, and fold change in abundance for proteins were calculated. A 1.5 fold change in expression was considered to identify differentially expressed proteins and phosphopeptides. We used a single sample pair in duplicates for our study. Thus, as in the RNA analysis, to identify and add statistical strength to the differentially abundant proteins and phosphopeptides, the Prediction Interval (PI) model was applied to the log2 normalized abundance values of proteins or of phosphopeptides (for the phosphoproteome data). We analyzed the two replicates separately. Prediction interval (PI) was calculated from the linear regression model applied using the log2 values from the control sample, and the log 2 abundance values of the proteins or phosphopeptides in the test sample were predicted based on the values in the control sample at a 95% prediction interval (upper and lower PI) to do so.
For determining the upregulated or downregulated protein, we used the predict () function in R as in the transcriptomics data: protein up= protein test > upper PI and protein test − protein control > 0:58 protein down= protein test < lower PI and protein test − protein control < −0:58 The z-score for each protein was calculated using the following formula: For further analysis and interpretation, proteins with ≥2 unique peptides or at least 2 peptide-spectrum matches (PSM) were considered for single peptide identifications. Single peptide based identifications with 1PSM per peptide were included only if supported with concordant transcript-level evidence. All phosphopeptides considered in the subsequent analysis were with ≥2 PSM with less than 40% coefficient of variation.

Cytokine Array
Cytokine profiling was performed using the Proteome Profiler TM Array Human (XL) Cytokine Array kit (ARY022B-R&D Systems). Arrays were incubated overnight, at 4°C, with 500ml of the conditioned media (CM) from GPR56 control or knockdown U373 cells, and the arrays were processed as per the manufacturer's instructions. The array images were analyzed using ImageJ to determine signal intensities. The mean pixel density of the cytokine/chemokine spots of GPR56 knockdown CM were expressed relative to GPR56 control CM values to determine the differential expression.

Macrophage Infiltration Assay
The conditioned medium of GPR56 control or knockdown U373 cells cultured in parallel at 70% confluency in DMEM/10% FBS was collected. U937 cells were cultured in RPMI1640 medium and induced by 5nM PMA (phorbol 12-myristate-13-acetate) (P1585, Sigma-Aldrich) for 48 hours as per the protocol described earlier (41,42) to differentiate them into macrophages. To check the level of macrophage infiltration, a standard transwell Boyden chamber invasion assay was performed (43), where U937 derived macrophages (n= 4x10 5 ) seeded in the top chamber (with 10% FBS) migrated towards secreted chemoattractants in the conditioned media from GPR56 KD cells containing 10% FBS) in the bottom chamber. Cell migration across the membrane was evaluated as described for 'Invasion assay'.

Network Analysis
Differentially altered gene and protein entities from the top five deregulated pathways, altered cytokines in the CM of the GPR56 knockdown cells, and the proteins with differentially altered phosphopeptides detected in our analysis were merged and used to examine protein-protein interactions. STRING database version 11 (http://string-db.org) was applied to construct the network. All the entities were uploaded into the STRING tool, and the default interactions were extracted with a confidence score >0.4. Then, the network was visualized using Cytoscape software (3.8.0).

GPR56 Expression in GBM Is Heterogeneous With Lowered Expression in Hypoxic Regions
Early studies have shown that GPR56, an aGPCR, is overexpressed in GBM (20,34). However, when we reassessed the mRNA expression of GPR56 by qRT-PCR in multiple GBM biopsy specimens (n=28) and at the protein level using immunohistochemistry (IHC) on a commercial tissue microarray (n=27), we found expression of GPR56 to be heterogeneous across specimens (Figures 1A-C and Supplementary Tables 1, 2). Moreno et al. recently reported similar findings and identified that GPR56 expression varies between subtypes of GBM (26).
To understand the heterogeneous expression of GPR56 further, we examined the expression pattern across different tumor compartments, using the spatial gene expression data from Ivy GAP resource (7). Figure 1D shows the abundance of GPR56 transcripts across various tumor areas. As seen, GPR56 was found to be distinctly upregulated in the cellular and infiltrating tumor area (CT and IT), while it was downregulated in microvascular proliferating (MVP) and pseudopalisiding cells around necrosis (PAN) areas. The IHC analysis also supported that GPR56 expression was low in the PAN (Red arrows) MVP regions and high in CT areas (green arrows). In addition, we examined the public domain GBM single-cell gene expression dataset (8) to understand its expression in different neoplastic cell types better. The single-cell neoplastic data identified four major groups of neoplastic cells, identified by six gene-based modules, i.e., hypoxia-independent (MES1-like) and hypoxia-dependent mesenchymal (MES2-like), astrocytic (AClike), oligodendroglial lineage (OPC-like), and neural-stem and progenitor type (NPC1-like and NPC2-like) as well as two cell cycling (CC) modules (G1S and G2M -CC) ( Figure 1E). We observed that in neoplastic cells, GPR56 expression was high in cycling (G1S) and in AC-like, OPC-like, and NPC-like cells but was lost in mesenchymal-like cells (MES1, MES2), particularly in hypoxia dependent MES2-like cells.

Multi-Omics Analysis of GPR56 Knockdown U373 GBM Cells Reveal Cellular and Molecular Changes Consistent With the Mesenchymal Transition
We generated GPR56 knockdown U373 GBM cells to do transcriptomic, proteomic, and phosphoproteomic analysis in order to gain insight into the molecular changes associated with GPR56 function. GPR56 expression levels were assessed following isolation of various single cell clones from shRNA GPR56 transfected cells using Western blotting. Figure 2A demonstrates that there was a significant reduction in GPR56 levels in clones 3 and 4, as assessed by loss of signal intensity for full-length GPR56 at 72 kDa and the processed N-GPR56 at 60 kDa. We have used shRNA GPR56 knockdown, clone 3 (GPR56-KD) for the rest of the analysis. Confocal analysis of control (GPR56-NC) and GPR56-KD cells revealed punctate staining for GPR56 in U373 control cells ( Figure 2B, top panel) and loss of staining in GPR56-KD (bottom panel). This was further verified at the mRNA level using qRT-PCR ( Figure 2C).
We next tested if this alteration in GPR56 expression affected chemotactic motility in invasion potential of GPR56-KD cells using a matrigel invasion chamber and compared invasion to the respective control cells. Figure 2D shows a significant increase in invasive behavior of the GPR56-KD cells. Similarly, wound closure was accelerated in GPR56 knockdown cells compared to control cells using a scratch wound assay. At 24 and 48h, about 69.3% ± 3.2%, and 94.06% ± 10.9% of the wound got closed respectively in the GPR56 knockdown cells compared to 46.1% ± 11.7% and 70.6% ± 4.8%, respectively during the same period in the GPR56 control cells ( Figure 2E). Additionally, functional assays were also performed with GPR56 shRNA clone 4, and the results further support the observations made with clone 3 (Supplementary Figure 1). These data indicate that GPR56 silencing plays a role in the invasive and migratory properties of the GBM cells.
To gain insights into molecular changes associated with GPR56 knockdown response, we carried out transcriptome, proteome, and phosphoproteome analysis of GPR56 control and knockdown cells. Figure 3 shows differential expression of RNA and proteins and proteins with differential phosphorylation in these cells. The transcriptome analysis of GPR56 knockdown U373 cells showed 1010 differentially expressed transcripts with 697 upregulated and 313 downregulated ( Figure 3A and Supplementary Table S3A). With the protein samples obtained from the same GPR56 knockdown cells, we used a liquid chromatography-based tandem mass spectrometry (LC-MS/MS) approach to identify differentially altered proteins and phosphoproteins. We could access 2221 proteins in the proteomic analysis (< 1% FDR at peptide level); 46 were differentially expressed, out of which 28 were upregulated, and 18 were downregulated. ( Figure 3B and Supplementary Table  S3B). The total coverage of proteins observed was relatively low (presumably due to the LC-MS/MS analysis carried out without pre fractionation of the tryptic peptides before mass spectrometry, as discussed under Methods). However, for all differentially expressed proteins identified, we observed that there was a positive correlation between protein and transcript abundances (average r=0.45; Figure 3D).
The phosphoproteomic analysis ( Figure 3C) identified 4471 phosphopeptides mapping to 1997 proteins; about half of these were new identifications, not detected in the global proteomic analysis. Out of 1997, 70 proteins had differentially altered phosphopeptides, with 36 showing elevated levels and 34 lower phosphorylation levels as compared to the control (Supplementary Table S3C). The phosphorylation was distributed across all types of phosphorylation sites (Serine, threonine, and tyrosine). Phosphopeptide differentials could originate either due to a change in the protein abundance per se or due to phosphorylation events as such. For this purpose, we examined the correlation between corresponding protein abundances (using global protein analysis) and the abundances of the phosphopeptides mapping to the respective proteins (average r=0.4; Figure 3E). We observed that the abundance of some phosphopeptides matched with corresponding proteins with an overall change in abundance. In contrast, 30 phosphopeptides were differentially altered, while the overall abundance of the corresponding protein itself remained unchanged, implying differential phosphorylation of these sites, which potentially leads to functional alterations of the proteins (Supplementary Table S4).
For functional annotations of differential expressions observed, we used RNA-seq data (differentially expressed genes; DEGs) and the PANTHER database (Version 15.0) on account of higher coverage and a larger number of differential entities observed. In the molecular function group, the DEGs were enriched mainly in DNA binding, catalytic activity, molecular regulators, and transcriptional regulator activity ( Figure 4A). Most DEGs are catalogued within the biological processes mapped to cellular processes, biological regulation, metabolic processes, and signaling. The top signaling pathways enriched are inflammation mediated by chemokine and cytokine signaling pathway, CCKR signaling pathway, angiogenesis, WNT, integrin pathway ( Figure 4B). The gene entities mapping to these pathways are given in the Table (Figure 4C). These enriched pathways are known to be associated with the mesenchymal transition or status in general (44). We have indeed detected several genes associated with EMT in general (45), such as TWIST1, MMP2, ITGB1, FN1, TG2, etc., from the transcriptome data, with details provided in the Supplementary  Table S5.
Based on the transcriptome data on GBM generated by the TCGA group, Verhaak et al. (3) defined four subtypes of GBM with distinct gene signatures, namely PN, N, CL, and MES subtypes, with the PN and the MES subtypes being more commonly observed. Given the enrichment of pathways associated with mesenchymal transition, we performed singlesample gene set enrichment analysis (SSGSEA) using these signatures of the molecular subgroups of GBM. We observed that differential gene expression data from GPR56 knockdown cells were enriched positively for genes associated with the MES subtype and negatively to the PN subtype ( Figure 4D, Supplementary Table S6). This observation is consistent with an earlier report with GBM cells (26). Further, we looked at the differentially phosphorylated protein group (Supplementary Table S4; n=30 proteins) in GPR56 silenced cells, and we found that it included proteins that have been reported to be involved in epithelial-mesenchymal transition (EMT), in general ( Table 3). In GBM, proneural-mesenchymal transition (PMT) is considered the equivalent of epithelial-mesenchymal transition associated with other aggressive cancers. Taken together, these results are consistent with an association of loss of GPR56 with the mesenchymal transition.
Tumor cells associated with mesenchymal transition are believed to promote tumor development by releasing growth factors and cytokines (62). Considering the mesenchymal-like shift observed upon loss of GPR56, we sought to examine the secretome of GPR56 knockdown cells. For this purpose, the conditioned medium (CM) of GPR56 control and knockdown cells were analyzed using cytokine arrays. Our analysis showed increased secretion of tumor invasion markers, proinflammatory cytokines, pro-angiogenic factors in the CM of GPR56 knockdown cells ( Figure 5A). Further, one of the effects of cytokines released by mesenchymal cells is their ability to recruit immune cells (63). To determine the role of GPR56 expression by GBM cells in recruiting immune cells, we performed an invasion assay using PMA-induced U937macrophage cells in the Boyden top compartment with the CM of GPR56 knockdown cells as the chemotactic medium in the bottom chamber and compared invasion infiltration of macrophage cells to the CM of control cells. We observed the conditioned medium from the GPR56 knockdown cells enhanced infiltration of macrophages ( Figure 5B). Thus, both cytokine array data and macrophage infiltration assay are consistent with the enrichment of soluble pro-inflammatory mediators in the CM of GPR56 knockdown cells.

GPR56 and TG2 Interaction May Have a Determinant Role in the Mesenchymal Transition of GBM
Interestingly, our multi-omics findings revealed that upon loss of GPR56, TG2 a known ligand of GPR56 in cancer was upregulated ( Figures 6A, B), suggesting coordinated expression of the two. The inverse expression level of GPR56 and TG2 was also observed in the Ivy Gap GBM RNA-seq data set ( Figure 6C). We validated this inverse expression using IHC in GBM specimens ( Figures 1E and 6D), supporting the earlier finding of an 'antagonistic' relationship between the two proteins in melanoma (29). Based on the finding that GPR56 expression was significantly lower in MVP and PAN areas (more distinct in MVP; Figure 6C), which are usually hypoxic (7), we queried plausible association between hypoxic state and expression of GPR56 and TG2. In solid tumors generally, hypoxia triggers molecular events leading to epithelial to mesenchymal transition (EMT), resulting in metastasis (64). In GBM, EMT is not prevalent as observed in other tumors, and hypoxia can cause similar molecular changes leading to the transition of the less aggressive proneural GBM to highly aggressive, chemo and radiation-resistant mesenchymal GBM, referred to as proneural to mesenchymal transition (PMT) (44). Literature survey indeed revealed that the GPR56 gene was found to be a target of hypoxia-inducible factors, and hypoxia downregulates GPR56 in some cancer cells (65,66). Thus, we induced hypoxia (1% O 2 for 72h) in three GBM cell lines in vitro and analyzed the expression of GPR56 and TG2 ( Figure 6E). In line with earlier reports and our observation in Ivy Gap, we observed lowered levels of GPR56 and elevation of TG2 expression in all GBM cell lines cultured under hypoxic conditions compared to normoxia (21% O 2 ). Similarly, loss of GPR56 is associated with the mesenchymal transition (26), while inhibition of TG2 has been reported to reverse mesenchymal differentiation of glioma stem cells (67). However, the present analysis does not provide any evidence for the direct role of hypoxia-inducing factors (HIF) in the co-regulation of GPR56 and TG2 expression and needs to be verified with direct investigations using HIF knockout cells. Nevertheless, together these observations permit us to propose a link between tumor hypoxia, reciprocal expression of GPR56 and TG2, and mesenchymal transition.
To gain further insights and assess protein interaction partners, we used differentially expressed gene and protein entities from the top five deregulated pathways (n = 55), all proteins with differentially altered phosphopeptides (n = 70), and altered cytokines (n = 16), from knockdown cells and constructed protein-protein interaction network using STRING (Supplementary Table S7). With a total of 138 altered entities, the STRING output captured 155 interactions (p-value: 1.0e-16), which were then visualized in Cytoscape. Our network indicated that GPR56-TG2 interaction is a core component from among all the entities belonging to pathways and processes that are aligned with the mesenchymal transition, as discussed above. Together, all our observations discussed above support a regulatory link of the GPR56 and TG2 interaction in the mesenchymal transition in GBM (Figure 7).

DISCUSSION
In this study, beginning with the heterogeneous expression of GPR56 in GBM specimens, we examined its expression in histologically distinct areas of GBM tumor tissue (7) as well as in single-cell expression data from GBM (8). Further, upon silencing of GPR56 in GBM cells, we performed a detailed functional and multi-omics molecular analysis of the GPR56 knockdown GBM cells to understand the consequential molecular changes and thereby the role of GPR56 in GBM.
The main finding of our analysis is: 1, GPR56 expression in GBM is heterogeneous (1A and B), supported by Ivy GAP spatial gene expression data showing higher expression in the cellular tumor and infiltrating tumor region and lower in the MVP and PAN regions. (Figures 1D, E). The single-cell expression data of GBM (8) further indicates low expression of GPR56 in the hypoxia-dependent mesenchymal cell types (MES1, MES2) ( Figure 1E). 2. GPR56 silencing in GBM cell line U373 increased their invasive, migratory behavior ( Figures 2B, C) as well as resulted in changes in gene, protein, and phosphoprotein expression, mapping to pathways and processes which are consistent with the mesenchymal state (Figures 3 and 4; Table 3). 3, The GPR56 knockdown cells also revealed an altered pattern of chemotactic and pro-inflammatory cytokine release which may mediate enhanced recruitment of immune cells (Figures 5A, B). 4, Interestingly, we observed the expression of TG2 to be reciprocal to that of GPR56, in the knockdown cells as well as in spatial tumor areas ( Figure 6). Similarly, hypoxic conditions reduced GPR56 expression and increased TG2 expression in GBM cells in-vitro ( Figure 6E), consistent with an 'antagonistic' relationship between the two proteins earlier reported in melanoma (29). Thus these observations point towards the appearance of mesenchymal state on suppression of GPR56 expression, and the interdependent reciprocal expression of GPR56 and TG2 observed may be a key feature in this process. We have tried to integrate all this data along with literature information in order to tease out molecular insights into the possible interplay of GPR56 and TG2 in the transitioning of cells towards the mesenchymal state.
The mesenchymal transition involves complex and highly coordinated molecular changes leading to altered cell adhesion and migration behavior of the cells, which are governed by a number of master regulators (MRs) including intracellular transcription factors such as STAT, Gli, ZAB induced by TNFa mediated nuclear factor kappa B (NF-kB) activation (44). In GBM, GPR56 is shown to inhibit the NF-kB pathway (26), while TG2 has been reported to activate it (68)(69)(70). On the other hand, TG2 expression itself is activated by NF-kB, thus creating an auto feedback loop (71). Given the regulatory loop between TG2 and NF-kB and the apparent role of GPR56 in limiting TG2 levels, GPR56 may be an important regulator of mesenchymal transition. Hypoxic conditions not only suppress GPR56 expression but may be linked to other molecular changes associated with the transition of the cell towards the mesenchymal phenotype. Thus, integrating our observations and the published information, we propose a putative model explaining the interplay of GPR56, TG2, and other interactors that may form an essential molecular network underlying MES transition in GBM (Figure 8). 3 | List of proteins detected with differential phosphorylation in GPR56-KD cells and known to be associated with epithelial to mesenchymal transition in general.

No
Gene Symbol Expression in GPR56 KD cells Function in EMT PMID Publication Phosphoproteomics data was analyzed as described under methods to identify proteins with altered phosphopeptides and not their overall abundance, indicating the change in their phosphorylation status in the GPR56-KD cells. They were mapped to the epithelial-mesenchymal transition process as per published literature, and the proteins are listed in the Table. In support of the scheme shown in Figure 8, in the GBM proneural cells, GPR56 is present at higher levels and may negatively regulate the NF-kB pathway-driven mesenchymal transition, as reported by Moreno et al. (26). Furthermore, it has been reported that PN glioma stem cells, on overexpression of TG2, upregulate mesenchymal MRs and mesenchymal markers (33). The low level of TG2 in the PN state is thus consistent with the "suppression" of the mesenchymal state as an additive mechanism. Thus, the inverse expression of GPR56 and TG2 through internalization of TG2 by GPR56 and its degradation, as reported in melanoma (29), may be an important aspect even in GBM. TG2 protein exists extracellularly, on the cell surface as well as intracellularly. Further, TG2 can assume two distinct mutually exclusive conformationsopen or closed forms regulated by binding of allosteric modulators Ca 2+ and GTP. Conformation-sensitive FRET experiments in live cells have shown that TG2 assumes a 'closed' conformation when localized in the perinuclear area and an 'open' conformation when localized near the plasma membrane. This may suggest that localized conformational transition (open form) occurs, presumably linked to Ca 2+ channel activation, provoking a shift in intracellular Ca 2+ /GTP concentration (72). The 'closed' signaling-active GTP-bound conformation of TG2 can drive intracellular signaling, whereas the 'open' transamidaseactive conformation of TG2 can enhance cancer cell survival after externalization and through altering ECM assemblies (32). Indeed, stromal TG2 (ECM) can promote tumor growth, and this may be suppressed by GPR56 expression on tumor cells as has been reported in melanoma cells whereby GPR56-TG2 binding results in internalization of the latter, followed by degradation through the endosomal pathway as discussed above (29). Thus, both transcriptional suppression and protein degradation may contribute to the negative regulation of TG2 in PN cells.
Hypoxia may trigger several molecular changes through the induction of hypoxia-inducible factors, possibly as an adaptive response. These include 1. Deregulation of Ca 2+ homeostasis which may be facilitated by changes in Ca 2+ signaling (influx) elements that have been reported to accompany GBM progression (73). 2. Release of inflammatory cytokines and their activation leading to molecular changes consistent with the mesenchymal state of the cells. For example Moreno et al. (26) showed that cytokine TNFa promoted a decrease in GPR56 expression, which may relieve the inhibition of the NF-kB pathway-mediated mesenchymal transition in GBM. Cytokines secreted in the microenvironment and known to be involved in EMT in general were also found to be upregulated in the transcriptome data and may induce TG2 expression (74,75). The elevated TG2 may lead to intracellular activation of the NF-kB pathway as an additive mechanism and induction of other MES-regulators. 3. Ca2+ mediated conformational change of TG2 accompanied by its extracellular translocation may enhance the extracellular TG2 pool, which may also be contributed to by endothelial cells and immune cells -a major source of secreted TG2. Together this may result in remodeling of ECM through increased crosslinking of ECM proteins by TG2 and interaction with cell surface integrins and activation of integrin signaling. 4. Intracellular TG2 in its GTP-bound 'closed' form (31) may also directly interact with integrins and modulate integrin signaling to bring about changes in migration-related downstream events. All these events may contribute to and have a role in mesenchymal transition in GBM. In summary, our study reveals some new insights on the regulatory role of GPR56 and TG2 in the mesenchymal transition of GBM. TG2 has been documented in the literature as a natural ligand of GPR56 (30). Our findings suggest that GPR56 may play a role in regulating the dynamics of TG2 levels and activity in tumor cells during the mesenchymal transition. Their reciprocal expression may be the dominant determinant of the ECM architecture on the one hand, and intracellularly, NF-kB mediated activation of pro-mesenchymal signaling cascades on the other. Elevation of TG2 levels concurrent with the loss of GPR56 is most likely due to transcriptional activation of the TG2 gene as well as its prolonged persistence in ECM due to the absence of GPR56mediated endosomal pathway and degradation. Since TG2 is already looked at as a potential therapeutic target in cancer (33,76) as its expression promotes chemo/radioresistance and invasive functions by inducing mesenchymal transition, a more in-depth understanding of the interplay of these two molecules in the therapeutic perspective, would be important to investigate. However, there is no clarity on the regulatory factors or mechanisms that underlie suppression of GPR56 expression in the mesenchymal state. Although TNFa alpha seems to promote downregulation of GPR56 in the mesenchymal GBM cells (26), the exact mechanism is not clearly understood. Association with the hypoxic condition and involvement of hypoxia-inducing factors (HIF 1 and 2) are strong indicators emanating from our analysis, but at this point, they are more of a probability to be experimentally confirmed. There could also be other regulatory factors involved. For example, miRNA-10a has been reported to be associated with the mesenchymal state, and temozolomide resistance in GBM cells (77), and interestingly, miRNA10a -  Table S3C), and the cytokines released ( Figure 5A) were used to construct the network using the STRING tool. Green indicates EMT-related molecules; pink, Phosphoproteins Blue, both (Phospho and EMT related) Grey, others. Direct key node interactors of GPR56 and TG2 are highlighted in Red.
GPR56 target interaction are revealed by the miRWalk target prediction tools, raising the possibility of its role in the regulation of GPR56 expression. Thus, GPR56 and its regulation present an attractive subject for future investigations.

CONCLUSION
Expression of GPR56 in Glioblastoma (GBM) was found to be heterogeneous, with the heterogeneity arising from expression of the receptor in spatially different tumor tissue and cell types (high in proliferating cells and low in hypoxic mesenchymal cells). On the basis of cellular assays and multi-omics analysis of GPR56 silenced U373 GBM cells, we infer that GPR56 plays a vital role in GBM cell invasion, migration, and mesenchymal transition. Furthermore, analysis of GPR56 silenced cells and spatial gene expression data of GBM tumors also revealed the expression of Transglutaminase 2 (TG2 -a known interactor of GPR56) to be inversely correlated with that of GPR56. GBM cell lines cultured under hypoxic conditions further supported the reciprocal regulation of the expression of the two genes. Integrating these observations, we have proposed a putative mechanistic link between the inverse expression of GPR56 and TG2, the hypoxic niche, and the regulation of mesenchymal transition in GBM.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: The RNA-seq data used in this study are available at the Gene Expression Omnibus (NCBI-GEO) database (http://www.ncbi.nlm.nih.gov/gds) under the accession ID GSE192874. The mass spectrometry proteomics data used in this study have been deposited to the ProteomeXchange Consortium (http://proteomecentral.proteomexchange.org) via the PRIDE (http://www.ebi.ac.uk/pride) partner repository with the dataset identifier PXD031569.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by Institutional (MSCTR) Ethics Committee, NHH/

AUTHOR CONTRIBUTIONS
RG performed and coordinated various experiments, integrated and analyzed data and wrote the manuscript. PS carried out transcriptomic and bioinformatics analysis. DN assisted in IHC and WB. AE and GS helped with proteomic and mass spectrometry experiments. KC, performed surgeries on GBM patients and provided patient samples and patient data. AL evaluated histology and immunohistochemistry sections. KV participated in data analysis and interpretations. LB generated stable knockdown cells and performed their initial characterization. DA and VK supervised the development of knockdown cells by LB, helped in data interpretation, critically reviewed and edited the manuscript. NS supervised public domain and experimental transcriptomics and bioinformatics analysis, jointly worked with RS in data interpretation and overall development of the study. RS conceptualized and designed and supervised the overall study, total data analysis, integration of the results and wrote the manuscript. All authors have read and approved the final version of the manuscript.