Functional RNA-Seq Profiling of the Tumor Suppressor Gene OPCML in Ovarian Cancers: The Multifunctional, Pleiotropic Impacts of Having Three Ig Domains

1. Introduction

The IgLON family of Tumor Suppressor Genes (TSG) were first described in expression studies of brain astrocyte cell cultures [1,2]. There are five protein family members with names reflecting their initial description in neuronal growth and differentiation: Opioid Binding and Cell Adhesion Molecule (OPCML), Neurotrimin (NTM), Neuronal Growth Regulator 1 (NEGR1), Limbic System Associated Membrane Protein (LSAMP), and an unnamed gene (IgLON5). However, early functional genomics screening of gene expression profiles in tumors revealed that IgLON TSGs were important in regulating cell growth, proliferation and differentiation in a variety of tissues, including many types of tumors. In the early 2000s, OPCML was clearly shown to be silenced in the majority of ovarian cancer cases (via DNA methylation of two promoter CpG islands) [3,4,5,6,7]. Over the last two decades, much attention has been given to defining aspects of the mechanism through which OPCML controls receptor tyrosine kinase (RTK) signaling pathways through protein-protein interactions thus playing an important natural role in preventing tumor growth in many cancers [8,9,10].

In healthy cells, low-level constitutive expression of IgLON TSGs maintain an appropriate density of IgLON membrane proteins within lipid raft microdomains to modulate cell surface interactions and signal transduction pathways [10]. Because the suppression or silencing of OPCML results in the loss or absence of the protein in many cancers, this gene/protein provides an excellent target for potential replacement therapies with a site of action at the external surface of the cell membrane [10,11,12,13]. In evaluating the therapeutic potential of OPCML, it is important to keep in mind that it generally has a dimeric structure consisting of two chains each containing three immunoglobulin (Ig) domains, anchored to the cell membrane via a glycosyl-phosphatidylinositol (GPI) moiety [11]. Thus, it is likely that there are several functional activities and processes by which OPCML may impact cells through the interactions of its three Ig domains with other extracellular ligands/targets.

In this study, we assess the range of biological and cellular processes that have synergistic associations with RNA transcript expression levels of OPCML in ovarian cancer (Stage IIIC, high grade serous carcinoma). Subject cohorts and RNA-Seq data were curated by The Cancer Genome Atlas (TCGA) Research Network. Our approach uses RNA-Seq ratios for OPCML normalized to each defined gene feature to bifurcate a subject cohort into a “normal constitutive” expression ratio group (> median ratio) and a “suppressed or silenced” expression ratio group (≤ median ratio). We then test for survival (Cox Hazard Ratios) in the “normal constitutive” ratio group and mortality risk in the “suppressed or silenced” expression ratio group to identify gene sets that exhibit a potential synergy on patient outcomes, relative to OPCML gene expression. Although this is an indirect association measurement, the phenotype being tested is patient overall survival, which is a robust and well defined outcome measured as “days_to_death” after primary diagnosis.

Functional enrichment analyses of these gene sets predict potential downstream linkages between OPCML expression and cellular activities and processes. OPCML’s tumor suppressor activity is ultimately defined by the functional effects of its three Ig domains, which potentially could act along several parallel paths to suppress tumor formation. Thus, the goal of this present study is to look deeper into the multi-functional impacts of OPCML in ovarian cancer.

2. Materials and Methods

Study Data: Subject cohorts utilize RNA-Seq data generated by the The Cancer Genome Atlas (TCGA) Research Network: https://www.cancer.gov/tcga. For high-grade serous Ovarian Cancer (TCGA-OV), inclusion criteria were: Stage IIIC tumors and integer "days_to_death" value (since diagnosis) greater than 0 and less than 3,650 (10 years). These criteria resulted in an ovarian cancer cohort of 189 subjects after the removal of one individual with an extreme RNA-Seq value for OPCML transcript counts (16.2 Standard Deviation units above mean Z-score value).

Survival Analysis: Data download and initial handling used the TCGAbiolinks R package for preprocessing and GC normalization [14,15]. These values were then transformed as $lo{g}_2(x+4)$ values to more uniformly scale the RNA-Seq ratios that were calculated to normalize OPCML expression to other genes. The linear upscale by $+4$ results in a minimum working value of 2.0 when normalized RNA-Seq values are near zero and ensures that the denominator in calculating the RNA-Seq ratios was never close to zero.

These ratio values were calculated for OPCML and the 60,660 transcript features available in the TCGA data sets. For each transcript feature, the ovarian cancer cohort was bifurcated into two groups by the median ratio value: 1) a “normal constitutive” ratio group (> median ratio; higher relative OPCML expression) and 2) a “suppressed or silenced” ratio group (≤ median ratio; lower relative OPCML expression). Survival analysis of these cohort subgroups used the R packages survminer and survival to assess differences in Kaplan-Meier survival curves (Log model p-values) and Cox Hazard Ratios (HR). All HR values are calculated as the survival risk or benefit for the “normal constitutive” OPCML expression subgroup.

A challenge to working with Hazard Ratios and applying probability distributions to assess significance is that HR distributions are neither symmetric nor balanced relative to the null hypothesis (Ho) assertion that the center of the HR distribution is 1.0. The lower tail of that distribution is limited to a range of 0 to 1. The upper tail can range from 1 to +∞. We applied a simple algebraic transform in R to establish a symmetrical distribution (Eq. 1)

$$HRsym=ifelse(coxHR<1,(-1\ast {\displaystyle \frac{1}{coxHR}}+1),(coxHR-1))$$

(1)

where the symmetrical HR value (HRsym) is derived from inverting the Cox HR when it is less than 1, and then adjusting the center of the null hypothesis distribution to assert HRsym=0 when there is no difference in survival. We finally use a simple Z-score transform (Eq. 2)

$$HRz={\displaystyle \frac{(HRsym-mean(HRsym)}{sd\left(HRsym\right)}}$$

(2)

where the Z-score HRsym transform (HRz) now has a mean value of 0 (null Ho “no effect” center) and a standard deviation (sd) of 1.

HRz False Discovery Rates (FDR): The high number of survival curves being assessed for each gene ratio raises concerns of high false positive rates or FDR. In this study, Monte Carlo empirical simulations were executed with the expression ratios calculated using randomized gene RNA-Seq matrices plus with randomized survival “days_to_death” values among cohort subjects. From the resulting distributions of HRz values we established the lower 0.025 quantile limits (equivalent to $\alpha$ = 0.05 for a two-tailed distribution) for HRz values that would likely be observed (with 95% confidence) due to only random sampling processes across this high-dimensional space (number of subjects * number of features). These randomized HRz thresholds of significance were measured and then applied to the real-world OPCML x$gen{e}_i$ data analysis to ensure our acceptance of significant HRz values were not in excess of a controlled false positive rate of 0.025 on either tail of the distribution (net FDR < 0.05). In this study, the vast majority of HRz values scored as significant (over 96%) were in the lower tail of the HRz distribution, i.e. survival benefit in the “normal constitutive” expression group (higher $Rati{o}_i$; HRz < 0).

In addition, we applied a screening criterion for assessing significant HRz results requiring that the Log model p-value of the survival curve for $Rati{o}_i$ stratification had to be at least 10% more significant than the Log model p-value of a corresponding survival analysis with median stratification based only on single $gen{e}_i$ RNA-Seq values. In this way, we ensured that any final decision of a functional significant $Rati{o}_i$ relationship evidenced a higher significance (lower p-value) than the sole $gen{e}_i$ expression by itself. As single gene RNA-Seq values, OPCML itself showed no significant survival differentiation in the ovarian cancer population under study.

Gene and Functional Enrichment Analyses: RNA-Seq features with significant HRz values were used to compile feature sets using custom R scripts. To better assess the functional significance of these gene sets, enrichment analyses used common curated databases: Gene Ontology (GO), Reactome (GSEA, Broad Inst., MIT) and Kyoto Encyclopedia of Genes and Genomes (KEGG). These different enrichment approaches yield different results and our approach here integrated these data streams into a functional snapshot of OPCML in Ovarian cancer [16]. The R packages clusterProfiler, enrichplot and pathfindR were used for executing analyses and visualizing results [17,18,19,20]. Most of the enriched results presented in this paper were generated by subnetwork association analysis of pathfindR using a reference protein interaction network from the BioGRID database [21]. For functional enrichment analyses, Benjamini-Hochberg FDR adjusted p-values were used at $\alpha$=0.050.

3. Results

3.1. Survival Curve HR Analyses

The assessment of RNA-Seq data from an Ovarian cancer cohort focuses on potential survival impacts of each RNA feature (n = 60,660) in independent Kaplan-Meier (KM) survival curve analyses. For each RNA feature, the cohort was first stratified by median RNA-Seq expression to establish a “baseline” contribution to survival (positive or negative) for each gene or feature. In Figure 1, example KM curves are shown for OPCML and CKS1BP3 (CDC28 Protein Kinase Regulatory Subunit 1B Pseudogene 3), neither of which evidence any significant differential survival based on their solo expression levels. In subsequent rounds of RNA-Seq ratio testing (following) these baseline p-values and hazard ratios are used to gauge any statistically significant changes when assessing RNA-Seq ratios of OPCML:$gen{e}_i$.

The survival Hazard Ratio analyses are used to identify possible synergistic associations between OPCML and other genes that show differential survival when coupled with OPCML expression. In Figure 2a, an RNA-Seq ratio of OPCML and CKS1BP3 is shown as a distribution of values among cohort subjects. Median stratification divides the cohort into two subgroups with High vs Low ratios. Although neither OPCML nor CKS1BP3 impact survival based on solo gene expression values (Fig. Figure 1), stratifying by the OPCML:CKS1BP3 ratio yields a significant survival benefit for the “High.Ratio” group (Cox HR = 0.657; model p-value = 0.0049; Fig. Figure 2b). Note, hazard ratios are always calculated as the risk ratio of the High.Ratio group relative to the Low.Ratio group. Although no functional association is yet curated for this CKS1B psuedogene, the active expression of CKS1B has been recently linked to epithelial-mesenchymal transition (EMT) in gastric cancers [22], and the regulation of EMT by OPCML via direct interaction with AXL has been described in ovarian cancers [10].

The rationale for assessing RNA-Seq ratios with survival curves is that synergistic associations between genes can exist across large interaction networks and this approach can potentially reveal linkages that may be distantly connected, e.g. from the cell membrane surface to nuclear transcription factors. Overall, Figure 2 is presented as an example of the underlying survival curve analyses used to produce the HR values that underlie our functional assessment of OPCML. In the following analyses, our goal is to execute an integrative approach across all gene features to identify patterns that reveal potential functional shifts in cellular processes. When single gene impacts become reinforced by other genes in related pathways, then a possible linkage between OPCML and specific downstream functional shifts in cellular and molecular activities becomes more apparent.

To deal with skewed and biased raw Hazard Ratio distributions, HR values were scaled equivalently (Eq. 1) and then Z-transformed (Eq. 2) to ensure normal distributions around the null hypothesis assertion of no difference in survival (HRz = 0.0). Monte Carlo iterations of randomized RNA-Seq and survival data were used to empirically determine false discovery rate thresholds. The baseline false positive rate in these randomized RNA-Seq trials averaged 0.0522. In Figure 3, the distribution of transformed Z-scores (HRz) are shown with the center of the distribution not statistically different from the null Ho assertion of no effect when HRz = 0 (t-test, mu=0, p = n.s.). In addition, the tail quantile values at 0.025 and 0.975 are equivalent and the distribution is normally distributed (Shapiro-Wilks test, Ho: Wstat = 1.0, p = n.s.).

The distribution of HRz scores (Fig. Figure 3) has two prograding tails. The FDR thresholds are empirically determined and closely approximate the theoretical 95% core distribution boundaries for a Standard Normal Z-distribution. Functionally, these boundaries define two groups of genes: a) ${\mathrm{OPCML}}^{+}$: HRz values < -2.0 identify genes with a survival benefit under constitutive OPCML expression levels (i.e., better survival in the RNA-Seq High.Ratio group [see Fig. Figure 2b]); b)${\mathrm{GENE}}_i^{+}$: HRz values > +2.0 identify genes with a survival benefit under suppressed OPCML expression levels (i.e., better survival in the RNA-Seq Low.Ratio group).

These two groups, ${\mathrm{OPCML}}^{+}$ and ${\mathrm{GENE}}_i^{+}$, are fundamentally different, yet connected to OPCML via RNA-Seq ratios. We highlight these groups here because they will be a focal point of subsequent analyses with green vs. orange color codes to separate functional results. The number of significant HRz scores for the ${\mathrm{OPCML}}^{+}$ group by far outnumber the ${\mathrm{GENE}}_i^{+}$ group count (2311 vs. 537; Fig. Figure 3). A fundamental interpretation of this result is simply that survival benefits are more frequently evident when OPCML expression is at “normal constitutive” levels. When OPCML is “suppressed or silenced”, there are other genes (${\mathrm{GENE}}_i^{+}$ group) that may potentially compensate for the loss of OPCML expression or parallel pathways that reinforce the impacts of OPCML but on downstream components of RTK signaling pathways.

3.2. Functional Enrichment Analyses

The collection of all significant gene features identified by survival RNA-Seq ratio HRz scores (Fig. Figure 3) were used for functional enrichment analyses. These were executed using the pathfindR package to identify active sub-network structures based on a defined protein interaction network (PIN), the results of which were then used for enrichment analyses [18]. For this project, four PIN database sources were used in serial combination with five pathway gene mapping databases (total independent enrichment analyses = 20). In Figure 4a, the count of enriched functions identified in each combination of PIN x pathway are shown. The PIN analysis sets were fairly equivalent in count results. For subsequent analyses of functions, results from all five protein interaction networks were compiled together, filtering out multiple instances of the same enriched function by keeping the one with the highest reported p-value. This collection of unique entries resulted in a data table with a total of 242 significantly enriched identified functions, pathways and processes.

In Figure 4b, differential composition of identified genes within enrichment functions is shown as a distribution using the ratio of ${\mathrm{OPCML}}^{+}$ gene counts to ${\mathrm{GENE}}_i^{+}$ counts in each function. The median of the distribution is 4.0 (i.e., $4$ genes to $1$) and is equivalent to the ratio of their total counts which is 4.3 ($2311$ genes to $537$). Thus the median position is close to the expected null position of a distribution determined solely by random assortment based solely on the proportional representation of the two groups of genes. Using the median distribution value (orange line in Fig. Figure 4b) functions were classified into two similar groups: 4x or greater representation of component ${\mathrm{OPCML}}^{+}$ genes, and < 4x ${\mathrm{OPCML}}^{+}$ genes, which translates into a greater representation of ${\mathrm{GENE}}_i^{+}$ genes. This classification scheme was used to assess whether functions highly enriched with ${\mathrm{OPCML}}^{+}$ genes were distinct in terms of cellular processes from other functions identified as significant in the analysis.

In Figure 5, hierarchical clustering of functions by their enriched gene composition does not show any visually distinct segregation of clusters based on the ${\mathrm{OPCML}}^{+}$and ${\mathrm{GENE}}_i^{+}$representation in the final gene sets. Here, 20 top-level clusters are identified (A-T) and three color classification groups are highlighted as examples (a,b,c). In general, the distribution of functions having high proportions of ${\mathrm{OPCML}}^{+}$genes are well mixed with those with more ${\mathrm{GENE}}_i^{+}$genes. The ${\mathrm{OPCML}}^{+}$functions do not separate into unified large (> 10) clusters. The overall topology of the tree shows that there is a diversity of enriched functions with many having very limited overlaps in gene composition.

Three example clusters are indicated in Figure 5. In cluster a, there are 10 functions ($6$and $4$) that have a high degree of gene overlap among them. The functions in this cluster group are described in Table 1. This cluster involves Toll-like Receptor (TLR) pathways. These entries are ranked by p-values and the top 6 have a >4x majority of ${\mathrm{OPCML}}^{+}$genes. The 4 functions with more ${\mathrm{GENE}}_i^{+}$composition form a discrete subgroup, but have sufficient overlap in enriched genes to be highly related in terms of pathway and process activities.

In cluster b (Fig. Figure 5), there are 31 functions evenly split on the ${\mathrm{OPCML}}^{+}$and ${\mathrm{GENE}}_i^{+}$ classifications. This cluster is focused on PI3K/AKT signaling in cancer and related IGF1 and FGFR1 signaling (Table 1). We know that OPCML drives degradation of RTK signaling receptors [13] in cancers. This large cluster aligns with how “constitutive normal” expression of OPCML contributes to increased survival in patients. It is interesting to note that when ranked by p-values in Table , the top ten functions are all in the >4x ${\mathrm{OPCML}}^{+}$gene count class. Again suggestive of this proximal cellular function by OPCML.

In cluster c (Fig. Figure 5), there are 8 functions. In contrast to clusters a and b, 7 of the 8 are in the >4x ${\mathrm{OPCML}}^{+}$gene count class, and the cluster break is at a high level in the tree (at branch P) and thus represents a distinctive functional grouping. In Table 1, this cluster focuses on external membrane cAMP signaling and connected pathways. Of note is the function “Calcitonin-like ligand receptors” (CLRs; Reactome 419812) with a xFold-enrichment of 16.2x, the highest in this study. CLRs have potent tumor micro-environment impacts in terms of stimulating both angiogenesis and cell proliferation [note: need to add ref]. The tumor suppressor activity of OPCML, particularly in its tumor micro-environment location, is clearly antagonistic to the action of CLRs in tumors.

Overall, hierarchical clustering of enriched functions by gene content is an efficient way to organize and visualize functional groupings within a large set of results (n=242). However, the challenge becomes how to work through all the cluster groupings to find key pathways and processes that are relevant to our understanding of the breadth of cellular impacts OPCML may have in suppressing tumor proliferation in ovarian cancer. The short summaries of the three groups selected from Figure 5 (Table ) are obviously important impacts OPCML has on ovarian tumors, but these were identified because they are structurally significant to the topology of the dendrogram. We need a quantitative strategy to winnow through all 242 enriched functions to distill the greatest relevance of these pathways for OPCML’s range of potential actions in ovarian cancer.

To guide identifying functions with high relevance to Ovarian Cancer, we utilized the PubMed database (US National Institute of Health, National Center for Biotechnology Information (NCBI)) to count the number of research publications over a 25-year period (1998 to 2023) relating our gene sets to ovarian cancer. For every single $gen{e}_i$ found in our enriched functions, a PubMed search was executed via API to count all articles where that $gen{e}_i$ appeared in combination with the subject “ovarian cancer”. In Figure 6, total $gen{e}_i$ PubMed counts are shown for both ${\mathrm{OPCML}}^{+}$and ${\mathrm{GENE}}_i^{+}$groups. A primary assumption is that the genes with the highest number of publications in the field of ovarian cancer are likely to be those that are “known” to be relevant to the progression, severity and treatment of these tumors. First, although both gene groups are well represented in the ovarian cancer literature, the ${\mathrm{OPCML}}^{+}$group has 2.4x more publication counts among the top 30 genes than the ${\mathrm{GENE}}_i^{+}$group (2,766 vs. 1,164).

Second, the top six genes, each with more than 100 publications (KIT, EGF, STAR, GAL, FGF2, VIP), have established roles in ovarian cancer. These genes have been studied and published because of their known relevance to understanding this cancer. Third, these gene-level PubMed counts can be combined into an Ovarian Cancer relevancy score for each enriched function as a roll-up score (geometric average) from the constituent genes found in an enriched function. This scoring approach, calculated for all enriched functions identified, allows for a triage ranking of functions relative to their potential involvement with onset and progression of ovarian cancer.

Our “known” Ovarian Cancer Relevancy score (OVca) was used to sort enriched functions from Figure 5 and collapse the complexity and overlaps of the dendrogram relationships into a ranked order for functions with >4x majority of ${\mathrm{OPCML}}^{+}$genes and those with less than a 4x majority (i.e., more ${\mathrm{GENE}}_i^{+}$component genes). Within any function cluster with more than 4 members, selection was limited to the top 4 functions by OVca score. Table 2 presents top 40 ranked functions. Because of the volume of ovarian cancer research published involving constituent genes of these ranked functions, the relevance of these pathways to ovarian cancer is high by default. In Figure 6, KIT, EGF and FGF2 are in the top 5 genes by publication count. In Table , functions with the highest OVca scores are those where KIT, EGF and FGF2 were identified as constituent genes. For most molecular oncologists, there are likely no surprises at the top of this table. However, the strength of this table is that all results are based upon a statistical linkage between OPCML:$gen{e}_i$ RNA-Seq ratios and survival benefits observed in patients with “normal constitutive” levels of OPCML expression. OPCML is linked to all the member genes of these functions in some fashion. Thus, cellular level impacts of OPCML are far reaching through many mechanisms of tumor progression and disease severity.

Whereas Table 2 generally reinforces what is already know about OPCML and Ovarian Cancer, a parallel ranking of functions with <4x majority of ${\mathrm{OPCML}}^{+}$genes (more ${\mathrm{GENE}}_i^{+}$genes present) provides a second perspective in Table 3. First, the publication counts for ${\mathrm{GENE}}_i^{+}$genes related to ovarian cancer are lower and consequently the top row OVca score values are less than in Table 2. Despite these lower OVca scores, p-values and xfold-enrichment values are significant. Second, ${\mathrm{GENE}}_i^{+}$group genes have a significant survival benefit when OPCML expression is low and these genes are expressed at a level above the median value in the patient cohort. These pathways in Table 3 are parallel to the impacts of OPCML such that when OPCML is silenced or suppressed, the expression of these pathways may ameliorate or compensate for the OPCML deficit and confer a survival benefit to patients.

4. Discussion

4.1. OPCML & IgLON Family Protein Structure

Much like other transmembrane proteins that function as receptors, IgLON family proteins (OPCML, NEGR1, LSAMP, NTM and IGLON5) have extracellular extensions of multiple Ig domains that arose from gene duplication. This multi-modular characteristic affords versatility to each IgLON family member in terms of specific protein-protein interactions each one orchestrates. All IgLONs possess three (3) sequential Ig domains (Fig. Figure 7) which is likely a stereospecific condition for IgLONs to engage in surface interactions with other juxtaposed protein extracellular domains (e.g, some RTKs possess two extracellular Ig domains). In addition, this configuration also enables IgLONs to form dimers via the first Ig domains (D1). Moreover, this further affords the possibility of alternate orientations of IgLONs participating in cell surface interactions either on the same cell (cis) or potential interactions with other neighboring cells (intercellular/trans).

The OPCML monomer has three Ig binding domains (Fig. Figure 7). The crystal structure was solved by Birtley et al. [11] where an excellent structural description of OPCML can be found. In Figure 7 known mutations in OPCML associated with different cancers are highlighted by the amino acid alpha carbons. Although very rare in occurrence, there were 119 clinically reported OPCML mutations presented by Birtley et al. [11]. More significantly, Birtley et al. investigated the impact of point mutations at 8 specific clinically relevant amino acid substitution sites (across all three Ig domains) and demonstrated that the reported clinical mutations at these sites in different tumors inactivated the tumor suppressor activity of OPCML. Thus, each of the three Ig domains appear to be important to the overall functionality of the OPCML protein.

The genesis of this current work, was the recognition that OPCML’s tumor suppressor activity is mechanistically a collection of functional effects defined by three Ig domains that in combination act along parallel paths to suppress tumor formation. The functional effects of OPCML are distributed across all three Ig domains independently, and thus the archetypal and pleiotropic tumor suppressor action of OPCML results from the integration of impacts from each of these three Ig domains. Thus, Birtley et al. [11] demonstrated the clinical relevance and multifunctionality of the three Ig domains and their extraordinary potency of integrated functions to inhibit tumor formation in ovarian cancer cells. The goal of this present study is to look deeper into the multi-functional impacts of OPCML in cancers through its multiple Ig domains.

4.2. RTKs, Ligands & Signaling

One of the principal functions of OPCML as a tumor suppressor protein is to inactivate receptor tyrosine kinases (RTKs) by dephosphorylation of the activated RTK and subsequently stimulating the degradation of that receptor tyrosine kinase (RTK). In ovarian cancers, of that receptor tyrosine kinases (RTK). In ovarian cancers, OPCML drives inactivation and degradation of the RTKs HER2, HER4, FGFR1, FGFR3, EPHA2, and AXL [23], and the specificity of this network is underscored by identification of a group of RTKs that are unaffected by OPCML. There may, however, be more RTKs impacted by OPCML to be described with further research. AXL is a good example of this inactivation/degradation pathway. When AXL is complexed with its ligand GAS6 in ovarian cancer cells, OPCML chaperones the AXL-GAS6 complex into external membrane lipid-raft domains where the complex becomes dephosphorylated through apposition to a raft restricted phosphatase and then targeted for ubiquitin-dependent degradation [10,13]. In the last decade, there has been a growing recognition of the important role that AXL and other RTKs play in ovarian cancers and their role in signaling through PI3K/AKT and MEK/ERK signaling pathways. By driving the dephosphorylation of RTKs within lipid raft domains, OPCML prevents low avidity aggregation of RTK species and their internal downstream signaling machinery that typically allows a high probability of the initiation and maintenance of signaling, and in this way inhibiting sustained ERK signaling, which is crucial for the spectrum of pleiotropic functions exhibited by RTK driven cancers. Thus, we know that OPCML plays a critical role in regulating cancer progression through its interaction with RTKs in the external membrane surface of ovarian cancer cells.

In this work, we demonstrate that in addition to this primary effect on RTKs, OPCML can influence a broader range of cellular activities that reinforce the potency of its Tumor Suppression effects. Beyond OPCML’s chaperoning effects on RTKs into cell-membrane lipid rafts, the analyses reported here regarding functional pathway enrichment (Fig. Figure 5 and Tables 2&3) emphasize deeper secondary association impacts of OPCML with other cellular processes that also impact tumor growth.

In Figure 8, a summary set of gene relationships and functional associations are compiled to illustrate the transcriptomic linkages between OPCML and genes that control broad functional aspects of tumor growth. These analyses are all based on RNA-Seq transcriptome data relative to overall survival in ovarian cancer patients. Here we see a downstream tumor suppressor role for OPCML beyond the direct interaction with membrane RTKs previously reported, and confirmed/validated orthogonally by previously reported pathways. The RNA-Seq ratio analyses reveal significant linkages with genes that span a range of tumor growth activities. Overall, the proteins shown in Fig. 8 are organized by cellular location and general function groupings: 1) known RTK signaling ligands, 2) extracellular matrix proteins, 3) other extracellular ligands, 4) other membrane receptors and proteins, 5) a few intracellular proteins involved with signaling and 6) two transcription factors and one nuclear protein. Note that the circle symbols in Fig. 8 with an asterisk indicate ${\mathrm{OPCML:GENE}}_i^{+}$ratios that have a direct and significant survival hazard ratio (HRz) for patients.

Highlighted in Figure 8 are the RTKs KIT, ROS1 and TEK that are significantly associated with OPCML RNA-Seq ratios in this study. These associations have not as yet been reported in any direct experimental studies. However, patient survival HRz scores reveal a new and significant RNA-Seq level linkage between OPCML and KIT, while TEK and ROS1 are prominent in several of the enriched function pathway results.

A key observation in Figure 8 is that there is a diverse range of transcriptomic links between OPCML and other proteins across mainly the extracellular matrix and cell membrane compartments. In this RNA-Seq analysis, we see far broader and more downstream associations at a transcriptional level (c.f. direct protein:protein interaction studies), indicative of the multifunctional “downstream-ripple” effect of OPCML expression. The enumeration of functionally enriched pathways (Figure 5 and Table 2 & Table 3) demonstrates a wide diversity of functions at granular resolution. In contrast, the level of proteins that impact tumor growth processes in Figure 8 reveal how many of these granular pathways align with OPCML’s general cellular activities as a tumor suppressor gene.

Additionally, many of the proteins in Figure 8 impact various immune system functions. Related to OPCML’s prominent external membrane presence, there appear to be many secondary associations connecting OPCML to immune status and activities within the tumor microenvironment. A recent study of NEGR1 colorectal cancers (CRC, mouse model system) showed that a decrease in NEGR1 expression resulted in CRC progression through the PI3K/AKT pathway in addition to a strong immune system interaction resulting in M2 macrophage polarization [24]. This role of immune system modulation may indicate a new role for OPCML specifically, and IgLON family proteins in general, that we are just beginning to describe and understand. Any degree of functional impact on immune modulation via OPCML could have relevant implications for cancer immunotherapies.

4.3. OPCML: Tumor Suppressor “Sphere of Influence”

Given that 49 out of 58 proteins (84%) in Figure 8 are operative externally to the cell or within the cell membrane, we conclude that a cell expressing OPCML could be described as establishing an OPCML "sphere of influence" extending from the extracellular matrix to the internal side of the membrane proximal to the PI3K/AKT and MEK/ERK signal transduction pathways. There are a total of 29 extracellular/membrane functional proteins with linkage to OPCML in contrast to only six intracellular proteins. The impact of OPCML expression in ovarian cancer cells supports strong relationships to external, microenvironment facing signal networks, less so downstream cytoplasm signal networks.

Alongside the regulation of classic outside-in signaling via cell surface receptors, extracellular proteins can impact three important processes in tumor growth: interactions between tumor cells, interactions between tumor and immune system cells, and modulation tumor microenvironments. In Fig. 8, we have identified 15 ligands and 8 proteins in the extracellular space that are important players in these processes. The ligands are mostly chemokines and cytokines (POMC, SHH, IL12A, IL4, CCL11, CXCL6, IL5, CXCL1, LEP and TNFSF11) that influence immune system activity and regulate tissue inflammatory responses. OPCML’s spatial relationship with the extracellular compartment with its three Ig domains allows for the opportunity to interact with other proteins and small molecules in this space. As a result, we see transcriptional linkages between some of these ligands and OPCML. Other proteins present here are involved in matrix remodeling, which can alter tumor invasiveness and thus metastasis (MMP13, S100A8, S100A9 and ELN).

In the membrane compartment there are also associations with G protein-coupled receptors (GPCRs: PTGIR, OPRM1, ADRA1A, HTR2A, HTR2B, FPR2 and NMUR2). We interpret these linkages as potential feedback loops between GPCRs and OPCML to influence OPCML’s cell-surface adhesion molecule role. In addition to these GPCRs, other significant signaling membrane proteins that are impacted by OPCML are FGR, CD36 and NGFR.

Overall, Figure 8 illustrates a potential outward radiating effect of OPCML from a cell expressing the protein to impact the surrounding extracellular matrix. It is possible in such a "sphere of influence" model that the OPCML membrane complement in one cell would regulate oncogenic RTK signals directly impacting that cell, as well as shaping the external mircoenvironments adjoining that cell in such a way that would further reduce other oncogenic processes impacting it. Thus by extension, any cell actively expressing OPCML could be shaping the microenvironment of its neighboring cells. The potential exists for an OPCML tumor suppression "sphere of influence" in which the composition of the extracellular matrix is determined by the collective OPCML titers of adjoining cells, whether these are cancer cells, immune cells or supportive stromal cells. Thus, the range of impacts on tumor growth and potential influence over neighboring cells, establish OPCML as a critical, multifunctional, pleiotropic tumor suppressor gene.

What our analyses show here are the transcriptional linkages that exist across many genes in diverse pathways that demonstrate significant associations with OPCML’s impact on mortality risk (overall survival) in ovarian cancer. This transcriptional association network is obviously quite different from previously demonstrated direct protein/protein interaction networks of physically interacting molecules for OPCML, although as previously stated it does also represent an orthogonal approach for validation of these previously published observations. The RNA-Seq data reveal concordant transcriptional events where subsets of the dozens of transcription factors known to be active in ovarian cancer can shift mRNA production of many genes simultaneously. It is interesting to note that the vast majority of genes shown in Figure 8 are in the ${\mathrm{OPCML}}^{+}$group vs. the ${\mathrm{GENE}}_i^{+}$group. We detect a reduced mortality risk when OPCML expression levels are in a normal, constitutive expression range and the expression of many of the other signaling pathway genes are below the cohort median. This suggests that the seesaw balance of down vs. up-regulation of cell signaling in cancers extends across many parallel functional pathways.

4.4. Conclusions

In this study, we have pursued a broader assessment of the potential functional impacts of OPCML and its component three Ig domains in ovarian cancer. Our bioinformatic approach of using RNA-Seq ratios, survival curves and enriched pathway analyses provides a novel perspective on the scope of cross-talk associations and cellular impacts that OPCML can influence. The net effect of OPCML is more extant than its proximal role in regulating the turnover of membrane RTKs. These normal activities span a breadth of pathways and processes, all of which contribute to the efficacy of OPCML’s tumor suppressor function.

These results highlight the role of OPCML’s regulation of RTK signaling pathways and their downstream signaling (PI3K/AKT and MEK/ERK) whilst also identifying three new potential RTK linkages to KIT, TEK and ROS1 in ovarian cancers. In addition, we show here that multiple known signaling receptor ligands are also linked in terms of gene expression to OPCML. Functional pathway enrichment analyses indicate a broad array of secondary linkages or cross-talk impacts of OPCML that translate to survival benefits for patients when OPCML gene expression is near normal constitutive levels. This information shows the broad reach of tumor suppressor activities engaged by this one gene/protein. Overall, OPCML is a multi-functional, pleiotropic, tumor suppressor that can inhibit tumor growth via multiple, secondary pathways and processes in a way that significantly impacts cancer patient survival outcomes consistent with its functional biology.