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.

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^−ΔΔCq) 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 .

Immunohistochemistry or Western blotting was used to study the protein level expression of GPR56 or TG2. in the tumor tissues or cell line. For immunohistochemistry, commercially available Tissue microarray (TMA), which included 27 GBM cases and 4 control tissue cores (BS17017a, BIOMAX, USA) or formalin-fixed paraffin-embedded (FFPE) GBM tissue sections were used. FFPE tumor tissue sections/TMA were deparaffinized, rehydrated, and endogenous peroxidases were blocked with 0.03% hydrogen peroxide. Antigen retrieval was performed using sodium citrate buffer (10 mM sodium citrate, 0.05 percent Tween 20, pH 6.0). Non-specific binding was blocked using 3% fetal bovine serum (FBS) in phosphate-buffered saline (PBS) with 0.1% Tween 20. Sections were then incubated for 2h at RT with antibodies - GPR56 (Santa- Cruz, Dallas, TX, USA, SC-9909, 1:200); TG2 (Thermo Fisher Scientific, Waltham, MA, USA, CUB7402, 1:200) followed by 1h incubation with peroxidase-labeled polymer conjugated to anti-mouse or anti-rabbit antibodies. After washing, sections were incubated with a 3,3- diaminobenzidine chromogen (DAKO, Denmark, K5007) for 2-5min and counterstained with Mayer's hematoxylin. IHC images were analyzed by a neuropathologist.

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.

Single-cell data: Neftel et al. (8) generated and analyzed single-cell 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 (MES1-like) 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 log-transformed data were z-scored using the following formula: z=X−μ where X is the expression measure of a gene in a sample, μ is the mean expression measure of a gene across all the samples, and σ is the standard deviation of expression measure of a gene across all the samples.

Human glioblastoma cell line U-373, U251, or U87-MG was cultured in high glucose Dulbecco's Modified Eagle Medium (DMEM-HG) (11995-065; Gibco, USA) and supplemented with 10% FBS (RM10409; HiMedia, USA) and with 1% Penicillin-Streptomycin (15140122-Gibco) at 37°C with 5% CO₂. U937 cells (a monocyte-like cell line) were cultured in Roswell Park Memorial Institute (RPMI)-1640 media (AL060A; HiMedia, USA) containing 20% FBS (10270106-Gibco, USA) and 1% Penicillin-Streptomycin. For the hypoxia experiment, the cells were incubated in hypoxic conditions, maintained with 1% oxygen at 37°C for 72 h.

SureSilencing shRNA plasmids targeting GPR56 and a non-targeting 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⁵ 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⁴ 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.

GPR56 control or knockdown cells were plated onto poly-L-lysine coated glass coverslips in 24 well plates at 1x10⁴ 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.

For invasion assay, we used modified two-chamber plates with an 8-μm pore size (Cat#. CLS3464, Sigma-Aldrich, USA) coated with 1 mg/ml matrigel (Cat#. E1270, Sigma-Aldrich, USA). GPR56 control and knockdown cells (10⁴ 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.

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.

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. Similar to the confidence interval, the prediction interval gives uncertainty around a single value point and allows a set of values within which a test value may fall. The linear regression model was applied to the log2 normalized CPM values. For each gene, based on the linear regression model, the log2 expression value of the gene in the test sample was predicted based on the value of the gene in the control sample at 95% prediction interval. We used the predict () function in R to do so. A gene was considered to be upregulated or downregulated if:

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.

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 500μl 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.

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⁵) 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’.

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).

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 (AC-like), 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.

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 single-sample 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, pro-inflammatory 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 U937- macrophage 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.

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₂ 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₂). 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 ).