Expression of long noncoding RNAs in cancer-associated fibroblasts linked to patient survival in ovarian cancer.

Cancer-associated fibroblasts (CAFs) are the most abundant cell type in the tumor microenvironment and are responsible for producing the desmoplastic reaction that is a poor prognostic factor in ovarian cancer. Long non-coding RNAs (lncRNAs) have been shown to play important roles in cancer. However, very little is known about the role of lncRNAs in the tumor microenvironment. We aimed to identify lncRNAs expressed in ovarian CAFs that were associated with patient survival and used computational approaches to predict their function. Increased expression of 9 lncRNAs and decreased expression of 1 lncRNA in ovarian CAFs were found to be associated with poorer overall survival. A "guilt-by-association" approach was used to predict the function of these lncRNAs. In particular, MIR155HG was predicted to play a role in immune response. Further investigation revealed high MIR155HG expression to be associated with higher infiltrates of immune cell subsets. In conclusion, these data indicate expression on several lncRNAs in CAFs are associated with patient survival and are likely to play an important role in regulating CAF function.

cancer‐associated fibroblast

false discovery rate

Gene Ontology

high‐grade serous ovarian cancer

hazard ratio

Kyoto Encyclopedia of Genes and Genomes

long noncoding RNA

microRNA

protein‐protein interaction

INTRODUCTION

Epithelial ovarian cancer is the fifth leading cause of cancer death in women and the most lethal gynecological malignancy.1 High‐grade serous ovarian cancer is the most common and aggressive subtype of ovarian cancer, and despite advances in understanding the underlying genetic causes of HGSOC and improved treatment strategies, most women diagnosed with HGSOC have a poor prognosis.2 Most women are diagnosed at an advanced stage and while initial response rates to chemotherapy are high, recurrence of chemoresistant disease is a significant problem.3 Continuing to improve our understanding of the factors that influence patient prognosis and response to therapy will be beneficial to improve treatment strategies and outcomes for women diagnosed with ovarian cancer.

Until recently, most studies have focused on a greater understanding of the molecular changes present in ovarian cancer cells and how these affect tumor progression and patient outcome. However, increasingly the tumor microenvironment is gaining recognition as playing a vital part in the initiation, survival, growth, and metastasis of tumors.4 In addition, cells within the tumor microenvironment are more genetically stable than cancer cells, which potentially reduces the likelihood of continued treatment causing the accumulation of genetic changes and subsequent development of acquired resistance to therapy.5 For these reasons, the tumor microenvironment is emerging as an attractive therapeutic target to treat cancer.

In ovarian cancer, the stromal proportion present in tumors can vary from 7% to 83% of tumor tissue, and patients with a higher stromal proportion have a worse overall survival.6 Furthermore, expression profiling of HGSOC has identified a subtype that displays a “stromal expression signature”.7 Importantly, patients with this signature show higher levels of desmoplastic stroma and the poorest survival. Within the tumor microenvironment, the stroma contains multiple cell types such as endothelial cells, immune cells, and CAFs, which have all been shown to contribute to cancer progression. Cancer‐associated fibroblasts represent the most abundant cell type in the tumor stroma and are responsible for producing the desmoplastic reaction that is a poor prognostic factor in HGSOC. Additionally, CAFs have been shown to play multiple roles in ovarian cancer to promote tumor cell proliferation, migration and invasion.8, 9, 10, 11

Gene expression profiling of several cancer types has revealed marked heterogeneity of CAFs.12, 13, 14 Therefore, distinct subtypes of CAFs could play varying roles in the tumor microenvironment and influence patient survival differently. Cancer‐associated fibroblast or stromally derived prognostic gene expression signatures have been identified in several cancer types.15, 16, 17 Furthermore, in ovarian cancer, studies have uncovered genes differentially expressed in CAFs that are predictive or prognostic biomarkers such as VCAN,18 CTGF,19 MFAP5,10 FOSB,20 EGR1, 20 and NPPB.9 However, studies investigating the DNA mutations in ovarian CAFs have concluded that somatic mutations are unlikely to contribute to gene expression changes seen in ovarian CAFs and raise the likelihood that alternative mechanisms of gene regulation occur in CAFs.21 A study published by Mitra et al showed that changes in the expression levels of miR‐31, miR‐214, and miR‐155 contribute to the reprogramming of normal fibroblasts into CAFs.22 In other cancer types, DNA methylation changes have also been shown to occur in CAFs.23, 24, 25, 26 A greater understanding of how gene expression is regulated in CAFs will help to identify new stromal biomarkers and potential therapeutic targets.

Long non‐coding RNAs represent another possible mechanism for regulating gene expression in CAFs. We have previously shown differences in lncRNA expression in ovarian CAFs compared to normal ovarian fibroblasts and that several of these lncRNAs might promote the prometastatic role of CAFs in ovarian cancer.27 Long noncoding RNAs are noncoding RNAs greater than 200 nucleotides long that do not encode protein. Once thought to be “transcriptional noise,” lncRNAs are now recognized to play crucial roles in several biological functions such as chromatin modification, transcription, and translation.28 They have also been shown to play important roles in several diseases, including cancer.29 Studies have identified lncRNAs involved in ovarian cancer that are also candidate prognostic biomarkers.30, 31 However, these studies have been restricted to whole tumor specimens or cell lines and none have examined the role of lncRNAs in the tumor microenvironment. As CAFs are known to represent a heterogeneous population of cells, with varying functional capacities that could influence patient outcome, this study aimed to investigate lncRNA expression in ovarian CAFs to determine those associated with patient outcome. We then employed a network‐based “guilt‐by‐association” approach to predict their functions. Additional analysis indicated that increased expression of one lncRNA, MIR155HG, was associated with significant increases in several immune cell subsets.

METHODS

Tissue specimens

Primary tumor specimens from 67 women diagnosed with HGSOC were obtained as previously described.19, 32 All specimens were from previously untreated HGSOC patients hospitalized at the Brigham and Women’s Hospital between 1990 and 2000. Patient specimens and corresponding clinical information were collected by written consent under protocols approved by the review board of the Brigham and Women’s Hospital Ethics Committee. All procedures were carried out in accordance with the approved guidelines and regulations. Classification was determined according to the International Federation of Gynecology and Obstetrics standards. Survival information was not available for 5 patients whose samples were excluded.

Microdissection, RNA isolation, amplification, and hybridization

Microdissection, RNA isolation, amplification, and hybridization to GeneChip Human Genome U133 Plus 2.0 Oligonucleotide arrays (Affymetrix) are previously described.19 Gene expression of endothelial cell markers (TIE‐2 and VEGFR1) and T cell markers (CD8 and CD45) were below the level of detection in our samples, indicating a lack of immune or endothelial components and enrichment for fibroblasts.19 All gene array data are available through Gene Expression Omnibus accession number GSE40595.

Statistical data preprocessing

Data preprocessing was undertaken using R Bioconducter, “affy” package. Data were normalized and background corrected using the Robust Multi‐Array Average method33 and expression values Log2 transformed. Variations across samples were assessed using the interquartile range values for each probe, and those with interquartile range less than 1 were removed for subsequent analyses. A total of 2448 probes were identified previously to be associated with lncRNAs.34 The gene symbols and titles corresponding to these probes were retrieved from GeneAnnot,35 which provides revised and improved annotations of the Affymetrix Human Genome U133 Plus 2.0 probes.

Survival analysis

The expression levels of all lncRNAs across samples were separated into low vs high expression using a fuzzy clustering algorithm, wherein each data point belongs to a cluster to some degree that is specified by a membership degree.36 The memberships are nonnegative, and for a fixed sample, they sum to 1. For each lncRNA, the fuzzy clustering algorithm was set to identify 2 clusters of samples representing higher vs lower expression levels. Samples belonging to either of clusters with the membership degree greater than 0.7 were included, ie, the remaining samples were considered as “undetermined” and excluded from subsequent analyses, resulting in 2 distinct groups of low vs high expression for each lncRNA. Fuzzy clustering was undertaken using R “cluster” package.

Kaplan‐Meier analysis and the log‐rank test were used to assess the association between the expression level of each lncRNA in CAFs and the patients’ overall survival. The prognostic value of each lncRNA’s expression levels as well as debulking status and chemotherapy response (sensitive vs the rest) was determined with univariate Cox proportional hazard modelling, and those significantly related to survival were incorporated into a multivariate analysis. In order to overcome the multicollinearity among lncRNA expression, principle component regression analysis was performed.37 Accordingly, instead of lncRNA expression profiles, the corresponding principal components were used as covariates in the multivariate Cox regression analysis. The Wald test was used to assess the statistical significance of the Cox models (α = 0.05). All survival analyses were carried out using R “survival” package; tied event times were handled by Breslow’s approximation.38

Functional prediction of lncRNAs associated with patient survival

Potential functions of prognostic CAF‐expressed lncRNAs were predicted using a network‐based “guilt‐by‐association” approach as follows.

Construction of coexpressed “interactome”

A coexpression network was first constructed where nodes are the identified lncRNAs and all protein‐coding genes and edges represent significant correlations, ie, |Pearson’s correlation coefficient| > 0.7, correlation adjusted P‐value < 10E‐6. The gene coexpression network was then mapped on a cellular interactome comprising protein‐protein and gene regulatory interactions. Experimentally validated human PPIs detected in more than 2 experiments were combined with highly ranked predicted PPIs (ie, FDR > 60%) as predicted by kotlyar et al39 to form a comprehensive PPI database. An experimentally derived gene regulatory network was secured from ORTI,40 a comprehensive repository of mammalian transcriptional interactions. The resulting network is a coexpressed interactome comprising coexpressed protein‐protein or gene regulatory interactions plus coexpressed lncRNA‐gene associations.

Network modules were identified using “community” detection algorithms where communities are groups of nodes with dense connections internally and sparser external connections.41 Communities represent transcripts that are more likely to be involved in distinct similar biological processes and thus can be used to assign functions to lncRNAs associated with them. All network analyses were performed using R “igraph” package. Different community detection algorithms were tried out, eg, Louvain,42 greedy,43 infomap,44 and walktrap.45 The Louvain algorithm found clusters with a relatively higher modularity score and thus was used to report the results.

Functional enrichment analysis

Network modules containing at least 1 lncRNA underwent GO and pathway enrichment analysis using the R “enrichR” package which implements Fisher’s exact test and FDR adjustment on a wide range of gene set libraries.46 We used biological processes (GO_Biological_Process_2017b, comprising 10 125 GO terms on 13 247 genes) and KEGG pathways (KEGG_2016, comprising 293 pathways on 7010 genes) as EnrichR datasets to predict putative functions of lncRNAs. For ease of visualization, GO terms enriched by each module were summarized into representative subsets of the terms using REViGO.47

Validation of MIR155HG in stromal‐enriched whole tumor specimens

To determine whether MIR155HG was prognostic in whole tumor samples, we used the same cohort as was used for the CIBERSORT analysis. Tothill et al7 have previously clustered these samples into 6 molecular subtypes (C1‐C6) using k‐means clustering and identified that C1 is enriched by genes associated with stromal cell types, enabling us to validate MIR155HG prognostication in an independent cohort.

We followed our pipeline to preprocess GSE9899 raw data and to compare survival differences between patients with high and low MIR155HG expression using Kaplan‐Meier analysis and the log‐rank test in each of the subgroups. Average MIR155HG gene expression was also compared between the subgroups and differences were compared using the nonparametric Wilcoxon test to account for nonnormality of MIR155HG mean expression across samples within each subgroup.

CIBERSORT analysis to determine immune cell infiltrates

CIBERSORT is an analytical tool designed to accurately estimate the immune cell subsets present in whole tumor samples from their gene expression profiles.48 We used the default LM22 signature matrix consisting of 547 genes that accurately distinguish 22 mature human hematopoietic populations and activation states, including 7 T cell types, naïve and memory B cells, plasma cells, natural killer cells, and myeloid subsets.

Using CIBERSORT, we quantified immune cell infiltrates from a cohort of 83 previously characterized high‐grade ovarian tumors7 described as having a stromal expression signature and an increased density of fibroblasts (C1) with data available from Gene Expression Omnibus, accession number GSE9899. Differences were compared between tumors with high and low MIR155HG expression. Data are expressed as absolute fractions of each immune cell type and differences between groups were measured by t test, with significance determined by a P value of less than .05.

RESULTS

Long noncoding RNA expression levels in ovarian CAFs associated with patient prognosis

Characteristics of patients and tumor samples included in this study are shown in Table 1. Microdissected CAF samples were quality controlled for CAF composition using 14 gene markers of CAFs and fibroblasts previously reported.49 To this end, differential expression analysis was carried out comparing CAFs vs matched microdissected epithelial tumors; a total of 7161 genes were identified to be significantly upregulated in CAFs (adjusted P value < .05 using moderated t test with FDR correction) covering 86% (12/14) of CAF markers, which indicates significant enrichment for CAF composition (ie, P value = 1.02e‐5, Fisher’s exact test). Additionally, expression profiles of marker genes in each individual CAF sample were assessed and visualized as shown in Figure S1, clearly illustrated the high expression level of CAF markers across all samples.

Table 1

Clinical characteristics of patients with ovarian cancer and tumor samples

Characteristic	n = 62	Description
Age at diagnosis, years (mean ± SD)	60.94 ± 12.37	—
Stage (III/IV), grade	55/7, 3	—
Debulking (optimal/suboptimal)	49/13	Optimal debulking corresponds to <1 cm residual tumor
Site, histological type	Ovary, serous	—
Chemoresponse (R/S/R‐S/Ref)	18/24/7/4	R, resistant (recurred < 6 months); S, sensitive (recurred > 6 months); R‐S, resistant‐sensitive (recurred at 6 months); Ref, refractory (never responded)

Kaplan‐Meier survival analysis indicated increased or decreased expression levels of 10 lncRNAs in ovarian CAFs were associated with patients’ overall survival (P‐value < .05). The symbols, chromosomal locations, and titles of these 10 lncRNAs are listed in Table 2. Kaplan‐Meier plots of each lncRNA identified distinct survival trends between samples in groups of high vs low expression (Figure 1). These groups were identified using a fuzzy clustering algorithm that was set to categorize samples into 2 clusters representing high vs low expression levels. The box plots beside each Kaplan‐Meier plot show how samples were clustered into 2 groups of high and low expression. Fuzzy clustering assigns membership grades to each sample indicating the degree to which it belongs to each cluster. The middle gray box in each plot represents “uncategorized” and “removed” samples that did not strongly belong to either group (membership degree < 0.7). We also undertook an identical analysis on the gene expression data using microdissected epithelial tumor cells from matched patient samples. Other than CRNDE, none of the identified lncRNAs was significantly associated with patients’ overall survival (Figure 2), confirming the CAF‐specific prognostic utility of the identified lncRNAs. Additionally, CAF expression patterns of the lncRNAs were not significantly correlated with their expression in matched epithelial samples (Figure 2), corroborating CAF distinct regulatory mechanisms.

Table 2

Long noncoding RNAs (lncRNAs) associated with overall survival in ovarian cancer patients based on Kaplan‐Meier analysis (P value < .05)

Symbol	Chromosome	Title
CRNDE	chr16q12.2	Colorectal neoplasia differentially expressed (nonprotein coding)
DANCR	chr4q12	Differentiation antagonizing non‐protein coding RNA
LOC642852	chr21q22.3	Uncharacterized LOC642852
MALAT1	chr11q13.1	Metastasis associated lung adenocarcinoma transcript 1 (nonprotein coding)
MEG3	chr14q32	Maternally expressed 3 (nonprotein coding)
MGC2752	chr19q13.43	Uncharacterized LOC100653267
MIR155HG	chr21q21.3	MIR155 host gene (nonprotein coding)
NEAT1	chr11q13.1	Nuclear paraspeckle assembly transcript 1 (nonprotein coding)
TP73‐AS1	chr1p36.32	TP73 antisense RNA 1 (nonprotein coding)
XIST	chrXq13.2	X (inactive)‐specific transcript (nonprotein coding)

Figure 1

Differentially expressed long noncoding RNAs (lncRNAs) in cancer‐associated fibroblasts associated with significant differences in overall survival among ovarian cancer patients. Higher expression of 9 lncRNAs was associated with shorter survival, whereas increased expression of MIR155HG was associated with longer survival as depicted in the Kaplan‐Meier curves. Box plots show results from the fuzzy clustering algorithm that separates high and low expression into 2 distinct groups

Figure 2

Expression of long noncoding RNAs (lncRNAs) that are prognostic in cancer‐associated fibroblasts (CAFs) are not prognostic in matched tumor epithelium of ovarian cancer patients. With the exception of CRNDE, expression of the lncRNAs in tumor epithelium were not associated with differences in patient survival, as depicted in the Kaplan‐Meier curves. Expression of lncRNAs was not correlated between microdissected CAF samples and matched tumor epithelium

The statistics of univariate cox proportional hazards analysis are shown in Figure 3A. Expression profiles of CRNDE, MALAT1, MEG3, TP73‐AS1, and XIST, as well as chemoresponse and tumor debulking, were statistically significant predictors of mortality in univariate analysis. The identified lncRNAs showed mutually significant correlations with each other (Figure 3B), but not with chemoresponse or tumor debulking (Figure S2). To adjust for existing collinearity among lncRNAs, the first principal component of these 10 lncRNAs (capturing 98% of variations) as well as clinical characteristics (ie, chemoresponse and debulking) were used as regressors in multivariate analysis. The first principal component of the lncRNAs (P value = .000116, HR = 0.74) and chemoresponse (P value = .000168, HR = 0.22) were significant predictors of survival in multivariate cox analysis. Debulking status approached, but did not reach statistical significance (P value = .067754, HR = 1.91).

Figure 3

Univariate Cox proportional hazards and correlation analyses. A, Expression levels of CRNDE, MALAT1, MEG3, TP73‐AS1, and XIST (highlighted) were significant predictors of mortality among patients with ovarian cancer as well as debulking status and chemoresponse. B, The majority of lncRNAs were significantly positively correlated with each other (blue shaded boxes, darker blue represents stronger correlations). Nonsignificant correlations are depicted by small gray circles; significant correlations are depicted by larger black circles as outlines in the correlation significance scale

Prognostic lncRNAs in ovarian CAFs enriched for pathways known to be involved in CAF function

“Guilt‐by‐association” assigns putative functions to coding/noncoding transcripts based on genes coexpressed with them. It relies on the idea that genes with similar expression patterns across multiple samples are more likely to be coregulated, share similar functions, or are involved in similar biological processes.50 Coexpression network analyses have been previously used to predict functions of lncRNAs.51, 52, 53

Figure 4A depicts a schematic view of a network‐based “guilt‐by‐association” approach followed in this work. We first constructed a coexpression network whose nodes include all protein‐coding genes as well as 10 lncRNAs identified by the survival analysis; edges represent coexpression relationships (|Pearson’s correlation coefficient| >0.7). This network held 6791 nodes and 5 557 325 edges and shows a relatively low degree of modularity (0.18) where 91% of nodes fall into 3 gigantic clusters. A high degree of modularity, however, has often been reported in biological networks.54 We mapped the coexpression network on a cellular interactome comprising PPIs and gene regulatory interactions to derive a coexpressed interactome. The corresponding network held 6791 nodes and 43 545 edges whose modularity was improved to 0.47. Overall, 17 clusters (excluding singletons, ie, clusters including only 1 member) were identified; 3 clusters contained the identified lncRNAs. DANCR, LOC642852, MALAT1, MEG3, MGC2752, TP73‐AS1, and XIST were coclustered in a relatively large module of 1711 genes. MIR155HG fell into a separate cluster (size = 242 genes), which is in concordance with the correlation pattern of lncRNAs visualized in Figure 2B CRNDE was clustered with only 2 protein‐coding genes, which is not sufficient for enrichment analysis, and NEAT1 was identified as a singleton. NEAT1 does not show sufficiently high correlation with any gene, indicating that the correlation‐based guilt‐by‐association analysis cannot reveal its function and complementary analyses are required. We therefore undertook functional enrichment analysis on the 2 former clusters comprising 8 of 10 lncRNAs of interest and identified cellular processes and pathways overrepresented (adjusted P value < .05) by the corresponding 2 subnetworks as predicted by functions of the constituent lncRNAs. Figure 4B shows the clusters and summarizes representative enriched GO terms. A list of all enriched GO terms as well as overrepresented pathways and cluster composition are available in Document S1. The cluster containing DANCR, LOC642852, MALAT1, MEG3, MGC2752, TP73‐AS1, and XIST showed an enrichment for pathways primarily involved in metabolic processes as well as autophagy and cilium assembly. The cluster containing MIR155HG was enriched for GO terms associated with the immune system, particularly pathways associated with T cell activation, antigen processing and presentation, leukocyte migration, and activation of an immune response. This cluster was also enriched for pathways involved in ECM organization, cell death, metabolic processes, and cytokine signaling. The KEGG pathways related to infectious diseases, immune diseases, and the immune system were also enriched, suggesting a role for MIR155HG in regulating the immune microenvironment.

Figure 4

Functional enrichment analysis. A, Workflow for prediction of long noncoding RNA (lncRNA) function using the proposed network‐based guilt‐by‐association approach that incorporates known protein‐protein and gene regulatory interactions to derive a coexpressed interactome. Coexpressed communities were then identified and modules containing at least 1 lncRNA were subject to functional enrichment analysis. B, Two clusters were identified containing at least 1 lncRNA. For each cluster, representative Gene Ontology terms are listed in the tables and are grouped by a broad functional classification. Overall, the node size is proportional to degree of nodes and nodes are colored blue to red by log2 of degree. Nodes are labeled by the gene/lncRNA name if degree > 50

Validation of MIR155HG expression and prognostication in stromal‐enriched whole tumor specimens

We were interested to know whether MIR155HG would remain prognostic using gene expression data from whole tumor samples. Therefore, we examined its prognostic value using an independent dataset.7 In this dataset, samples were classified into 6 subtypes (C1‐C6) based on gene expression. As shown in Figure 5A, MIR155HG was only prognostic in the C1 subtype (log‐rank test P value = .0199). This subtype is described as having a stromal expression signature and an increased density of myofibroblasts. Therefore, our observation that MIR155HG is only prognostic in the Cluster 1 subgroup likely reflects the higher contribution of CAFs to the observed gene expression in this subtype.

Figure 5

Analysis of MIR155HG in whole tumor specimens. A, High expression of MIR155HG was able to predict prolonged survival in the C1 subtype of ovarian tumors that have been shown to have a stromal gene expression signature, but was not prognostic in the other subtypes. B, C1 and C2 subsets showed significantly higher MIR155HG expression than the other subtypes

We also compared MIR155HG expression levels between the different subtypes. MIR155HG expression was significantly higher in the C1 and C2 subtypes (Figure 5B). Both these subtypes have previously been associated with higher levels of infiltrating CD3+ T cells, with the C1 subtype showing high levels of stromal CD3+ T cells and the C2 subtype showing a high level of intratumoral CD3+ T cells.7 The higher MIR155HG expression seen in these subtypes supports our functional prediction analysis and the CIBERSORT results showing higher CD3+ T cells in tumors with high MIR155HG expression.

High MIR155HG expression associated with increase in immune cell subsets in stromal‐enriched whole tumor specimens

Based on the enrichment analysis and the longer survival seen in patients with high MIR155HG expression, we hypothesized that MIR155HG is associated with differences in immune cell infiltrates within the tumor. In order to investigate this further, we used CIBERSORT48 to examine the immune infiltrates present in whole tumor specimens obtained from a cohort of 285 ovarian cancer patients,7 separated by their MIR155HG expression. As shown in Figure 6, tumors with high MIR155HG expression had significantly higher numbers of plasma cells, CD8+ T cells, CD4+ memory activated T cells, follicular helper T cells, gamma delta T cells, M1 macrophages, and eosinophils.

Figure 6

CIBERSORT immune cell infiltrate analysis. Graph displays the absolute fraction (%) ± SEM of different immune cell subsets within ovarian tumors, separated into MIR155HG low and high expression groups. *Significant difference (P value < .05)

DISCUSSION

Cells within the microenvironment of solid tumors are not passive bystanders in tumor progression and metastasis but play an active and essential role. In many tumor types, including ovarian, CAFs are known to influence tumor cell behavior by increasing tumor cell survival, proliferation, migration, and invasion.8, 9, 10, 11 Cancer‐associated fibroblasts also interact with the other cell types present in the tumor microenvironment to promote angiogenesis and help tumor cells evade immune destruction.13, 55 Given these findings, a greater understanding of the molecular features of CAFs and their potential role in the clinical behavior of tumors is essential when designing new therapies that target CAFs.

We recently reported that several lncRNAs are differentially expressed in ovarian CAFs compared to normal ovarian fibroblasts and that several of these lncRNAs contribute to the prometastatic phenotype of CAFs.27 In the current study, we investigated whether differences in lncRNA expression in CAFs influence patient outcome in HGSOC. Given the importance of CAFs in ovarian cancer, there is a rationale for exploring the molecular aberrations present in CAFs as these could provide valuable prognostic information. We identified 10 lncRNAs with variable expression in ovarian CAFs that were associated with overall survival in HGSOC. In addition, an expression signature based on these 10 lncRNAs was an independent predictor of patient survival. Several of these lncRNAs are already known to play a role in either ovarian cancer,56, 57, 58 or other cancer types59; however, these previous studies have only examined lncRNAs in whole tumor specimens, therefore, it is not clear whether it is expression in the tumor cells or the microenvironment that is associated with patient survival. By analyzing expression data from microdissected CAFs and matched tumor cells, we were able to show that for 9 out of 10 of our lncRNAs, differential expression in CAFs specifically, and not tumor cells, was associated with patient survival. This suggests for the first time that these lncRNAs could play an important role in CAFs and the tumor microenvironment.

Even though several of the lncRNAs identified in this study have previously been studied in ovarian cancer or other cancers, their function in CAFs is not clear. In addition, LOC642852 and MGC2752 are not well characterized. In order to elucidate the potential functions of these lncRNAs in CAFs we used a network‐based guilt‐by‐association approach. The majority of lncRNAs clustered together, suggesting they play similar roles in CAFs. This cluster was associated with pathways involved in metabolism, autophagy, and cilium assembly. Cancer‐associated fibroblasts are already known to play a role in the metabolic reprogramming of the tumor microenvironment in order to favor cancer growth and metastasis.60 Through alterations in their metabolic activity, CAFs take on a catabolic phenotype that is then able to provide nutrients to anabolic cancer cells.61 In ovarian cancer, CAFs have been shown to have altered metabolism compared to normal ovarian fibroblasts and targeting this altered metabolism resulted in tumor regression.62 An essential part of tumor metabolism, the process of autophagy, is activated in CAFs as a potential mechanism to allow CAFs to provide metabolic products to feed cancer cells, and high levels of autophagy in the tumor microenvironment has been associated with cancer progression.63 In ovarian cancer, autophagy could protect ovarian CAFs against oxidative stress.64 Supporting our findings, both MEG3 and MALAT1 have been shown to induce autophagy in ovarian cancer.65, 66 Cilium assembly pathways were also enriched in this cluster. Primary cilia are important for signaling between stromal cells and adjacent tumor cells and autophagy promotes cilia formation.67 Enrichment of pathways associated with metabolism, autophagy, and cilium assembly suggests that lncRNAs belonging to this cluster might be important in creating a metabolic environment conducive to ovarian tumor growth. Higher levels of these lncRNAs could be indicative of more metabolically active and aggressive tumors, resulting in worse patient survival. However, further studies are required to validate the functional roles of these lncRNAs in ovarian CAF metabolism and autophagy.

The other cluster identified in this study contained the lncRNA associated with longer survival, MIR155HG. MIR155HG was originally identified as a proto‐oncogene in B‐cell lymphomas68 and is known to regulate many immune and inflammatory processes.69 The role of MIR155HG has not been well‐studied in cancer; however, recent studies have shown increased expression to be associated with worse survival in glioma and pancreatic adenocarcinoma patients, but improved survival in colorectal cancer patients.70, 71, 72, 73 In tumor cells, MIR155HG appears to have an oncogenic role and is associated with increased cell growth and decreased apoptosis.73, 74 However, MIR155HG could play a different role in CAFs. Our functional prediction analysis showed the cluster containing MIR155HG was highly enriched for immune‐related pathways involved in T cell activation and activation of an immune response. Given patients with higher expression of MIR155HG in CAFs survived for longer, MIR155HG could be an important component in the interaction between CAFs and immune cells, and might be promoting or permitting an antitumor immune response. The CIBERSORT analysis indicated high MIR155HG expression to be associated with increased immune cell subsets previously shown to be associated with improved survival,75 further experiments manipulating MIR155HG expression in CAFs are required to determine whether or not MIR155HG is directly involved in the induction of an antitumor immune response. Interestingly, MIR155HG’s associated miRNA, miR‐155, is a well‐known regulator of immunity and its expression has been shown in preclinical mouse models to be essential for mounting an antitumor immune response.76, 77 This finding could support a role for MIR155HG in antitumor immunity as many miRNA host genes have been shown to have similar functions to their associated miRNA,78 however, this requires further investigation.

The concept that lncRNAs are involved in the regulation of the immune system and more specifically in the regulation of tumor immunity is relatively recent and therefore their importance in the tumor microenvironment and their potential clinical utility is not yet well known. Given that high MIR155HG expression in CAFs is associated with increased patient survival as well as increased infiltration of antitumor immune cell subsets, MIR155HG expression could represent a useful biomarker to predict response to immunotherapy.

In summary, we have identified that variable expression of several lncRNAs in ovarian CAFs is linked to patient survival. Functional prediction using computational models highlights several potential ways these lncRNAs are regulating the ovarian tumor microenvironment to create an environment for tumors to grow and metastasize and evade immune destruction. Given the crucial role of the tumor microenvironment in cancer initiation and progression, continuing to understand the complexity of CAFs is essential to identify novel biomarkers and improved ways to therapeutically target the tumor microenvironment.

CONFLICT OF INTEREST

The authors have no conflict of interest.

FigS1‐S2

Click here for additional data file.

Supinfo

Click here for additional data file.