Adult-type granulosa cell tumors (aGCT) are rare ovarian sex cord tumors with few effective treatments for recurrent disease. The objective of this study was to characterize the tumor microenvironment (TME) of primary and recurrent aGCTs and to identify correlates of disease recurrence. Total RNA sequencing (RNA-seq) was performed on 24 pathologically confirmed, cryopreserved aGCT samples, including 8 primary and 16 recurrent tumors. After read alignment and quality-control filtering, DESeq2 was used to identify differentially expressed genes (DEG) between primary and recurrent tumors. Functional enrichment pathway analysis and gene set enrichment analysis was performed using “clusterProfiler” and “GSVA” R packages. TME composition was investigated through the analysis and integration of multiple published RNA-seq deconvolution algorithms. TME analysis results were externally validated using data from independent previously published RNA-seq datasets. A total of 31 DEGs were identified between primary and recurrent aGCTs. These included genes with known function in hormone signaling such as LHCGR and INSL3 (more abundant in primary tumors) and CYP19A1 (more abundant in recurrent tumors). Gene set enrichment analysis revealed that primarily immune-related and hormone-regulated gene sets expression was increased in recurrent tumors. Integrative TME analysis demonstrated statistically significant depletion of cancer-associated fibroblasts in recurrent tumors. This finding was confirmed in multiple independent datasets.

Implications:

Recurrent aGCTs exhibit alterations in hormone pathway gene expression as well as decreased infiltration of cancer-associated fibroblasts, suggesting dual roles for hormonal signaling and TME remodeling underpinning disease relapse.

Adult-type granulosa cell tumors (aGCT) are rare malignant ovarian sex cord-stromal tumors that exhibit frequent FOXL2 c.C402G (p.Cys134Trp) hotspot mutations (1, 2). Most primary aGCTs are diagnosed at an early stage and can be removed surgically with long-term survival rates exceeding 90% (3). Importantly, aGCT recurrence is difficult to predict and is nearly always incurable after relapse. Treatment of relapsed aGCTs remains a significant clinical challenge. Although knowledge about tumor-intrinsic correlates of recurrence is increasing, little is known about the tumor microenvironment (TME) drivers of aGCT recurrence.

Several published studies of primary and recurrent aGCTs have used massively parallel DNA-sequencing platforms to identify recurrent mutations in this disease (2, 4–6). These previous studies have identified enrichment of both TERT C228T promoter hotspot mutations and truncating KMT2D mutations in recurrent aGCTs (5–7). Another recent study reported that the FOXL2 homozygous genotype and chromosome instability are prevalent in recurrent aGCTs, and the homozygous FOXL2 genotype may predict a shorter time to relapse (8).

Differences in gene expression between early and advanced stages aGCT have also been investigated (9). Expression of MFAP5 was significantly higher in advanced-stage tumors whereas INSL3 and DES were significantly more abundant in early-stage tumors. At least two other prior studies have used RNA sequencing (RNA-seq) to examine differential gene expression between primary and recurrent aGCTs. One study included 10 tumor samples (6 primary and 4 recurrent) and identified 1,091 genes as differentially expressed (10), including several with known roles in estrogen signaling. A subsequent study analyzed RNA-seq data from 27 primary aGCTs, including 14 that did not recur (non-relapsed primary aGCT) and 13 that recurred (relapsed primary aGCT), and 8 recurrent (relapsed) tumors. Non-relapsed and relapsed primary aGCTs were found to have similar transcriptomic profiles whereas a comparison of primary and relapsed aGCTs samples showed three differentially expressed genes (DEG), including PLVAP (higher expression in recurrent tumors) and ASS1 and PLIN4 (higher in primary tumors; ref. 11). These prior gene expression studies have focused on tumor-intrinsic differential gene expression, and to date there has been no investigation of TME composition in this tumor type using massively parallel sequencing platforms.

In contrast, efforts to characterize aGCT TME have generally been limited by small cohort size and the use of IHC approaches that focus on a limited subset of immune cell types (12). The role of TME in aGCT development is still poorly delineated and no comprehensive comparison of TME composition among primary and recurrent tumors has been reported. The purpose of the present study was to identify tumor-intrinsic and -extrinsic differences between primary and recurrent aGCTs that may correlate with recurrence. To do this, we exploited high-quality bulk RNA-seq of cryopreserved tissue samples to perform a systematic analysis using multiple, validated TME deconvolution algorithms, bridging the knowledge gap between differential gene expression analysis and low-throughput IHC analysis of TME in primary and recurrent aGCT.

Patients and samples

This retrospective study was approved by the institutional review board at The University of Texas MD Anderson Cancer Center (protocols PA16-0891, LAB02-188, and 2022-0547) and patients provided written informed consent for tissue storage and analysis. Cryopreserved tissue samples were obtained from The University of Texas MD Anderson Cancer Center Multidisciplinary Gynecologic Cancer Tumor Bank, and 24 pathologically confirmed aGCT samples were analyzed by whole-exome sequencing (WES; ref. 6) and total RNA-seq (present study). The cohort included 8 primary aGCTs and 16 non-matched recurrent aGCTs. The somatic FOXL2 c.C402G mutation was identified in 23 of 24 tumor samples (6). Samples were identified as having truncating somatic KMT2D variants and/or focal deletions of 21q22.3 and 3p29 as described previously using the WES data (6).

Summary statistics were used to compare clinical and demographic factors by tumor category. Comparisons were made between groups using χ2, Fisher's exact, t tests, or Wilcoxon rank sum test depending on the underlying distribution of the data. Overall survival (OS) was estimated using the methods of Kaplan and Meier and compared using a log rank test. OS was defined as time from initial diagnosis to date of death or last contact. All statistical analyses of clinical data were performed using Stata/MP v17.0.

RNA-seq and data processing

Total RNA was extracted from cryopreserved tissue of aGCTs samples using the Qiagen RNEasy Mini Tissue kit. For RNA isolation frozen tissue samples were cryosectioned and a section adjacent to the portion submitted for RNA-seq underwent hematoxylin and eosin staining to assess tumor cell content. These adjacent sections were reviewed by a gynecologic pathologist (R.R. Broaddus) and only tumor samples with an estimated tumor cell content >60% without significant necrosis were included in the present study. Libraries were prepared from cDNA using the NuGEN Ovation Ultralow Library System V2. Paired-end bulk RNA-seq was performed on the high-throughput sequencing Illumina HiSeq 2000 platform. RNA-seq reads were aligned to the hg19/GRCh37 reference human genome using the STAR software (version 2.6.0b, RRID:SCR_004463; ref. 13) with default parameters. The annotation was obtained in GTF format from Ensembl database (https://useast.ensembl.org/index.html), version Homo_sapiens.GRCh37.75.gtf. Gene expression values were quantified as the number of reads that align to each gene (read count). Raw count data of each gene were obtained using HTSeq (version 0.11.1; ref. 14). RNA-seq metrics are presented in Supplementary Table S1. RNA-seq data are available through the EGA database (EGAS00001006478).

Differential gene expression analysis

A read count matrix was used to identify significantly DEGs. Normalization and differential expression analysis were performed using R package “DESeq2” (v.1.32.0, RRID:SCR_015687; ref. 15) with default settings except for parameter “minReplicatesForReplace” that was set to “Inf.” Minimal pre-filtering rule was applied to count matrix to keep only rows that have at least 10 reads total (removing low expression genes that have no counts, or only a single count across all samples). We performed multiple comparisons adjustment using the Benjamini–Hochberg approach to acquire adjusted P values. DEGs with the adjusted P value of <0.05 and absolute fold change >2 were considered as statistically significant. The heatmap of DEGs was constructed using “pheatmap” R package (16). To visualize our results, Z-scores were computed on a gene-by-gene basis by subtracting the mean and then dividing by the standard deviation using scale = “row” pheatmap option.

Functional enrichment analysis

To analyze the functions and the potential pathways of DEGs, Gene Ontology (GO) term (17) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway (18) enrichment analyses were performed using the “clusterProfiler” R package (v.4.0.5, RRID:SCR_016884; ref. 19). Over-representation analysis (ORA) was performed with enrichGO and enrichKEGG clusterProfiler functions.

Because DEGs are the only input for ORA, this method may overlook the genes that show non-significant changes in expression. ORA does not account for genes with small changes in expression that might be also biologically relevant. To overcome these challenges and to solve the problem of undetectable, minor changes in gene expression Gene Set Enrichment Analysis (GSEA) method was applied (20). This method considers sets of genes simultaneously ranking all genes based on difference in expression between studied groups what allows to detect situations where all genes in a predefined gene set change in a small but coordinated way. GSEA was conducted according to the method described (20), using GO terms (21), KEGG pathways and Molecular Signatures Database (MSigDB; ref. 22) gene sets (H: hallmark; C3: regulatory target; C6: oncogenic signatures; C7: immunologic signature). To perform GSEA, we used gseGO, gseKEGG, and GSEA functions of the “clusterProfiler” R package. In our study, we performed pre-ranked GSEA where all genes were ordered according to Fold Change enrichment calculated by DESeq2 (recurrent vs. primary tumors). A ranked gene list with corresponding fold change values was then provided to the “clusterProfiler” R package. The parameters were set as following: eps = 0, minGSSize = 10, maxGSSize = 500, pvalueCutoff = 0.05. One thousand permutations were performed to obtain randomized enrichment score and calculate the normalized enrichment score. The P value associated with each gene set was adjusted by using the Benjamini–Hochberg multiple testing procedures; the gene sets with an adjusted P value of <0.05 were considered as significant.

Next, we conducted single sample GSEA (ssGSEA) using the “ssGSEA” method of “GSVA” R package (RRID:SCR_021058; ref. 23) and gene sets from MSigDB. As input for the GSVA algorithm, a gene expression matrix in the form of RNA-seq counts and a database of gene sets were provided. Differences between primary and recurrent groups were found using scores calculated by “ssGSEA” and applying the Wilcoxon test.

TME analysis

Several validated cell-type estimation methods were used to estimate levels of noncancerous cells in the TME. To analyze the tumor-infiltrated cell fractions, we used CIBERSORTx (RRID:SCR_016955; ref.24), EPIC (25), MCP-counter (26), quanTIseq (RRID:SCR_022993; ref. 27), and xCell (28) methods. To perform deconvolution and quantify the proportions of immune cells in aGCT samples with CIBERSORTx (24), gene expression profiles were uploaded to the CIBERSORTx web server (https://cibersortx.stanford.edu/), Leukocyte gene signature matrix (LM22) was used as signature matrix (29), with this LM22 matrix CIBERSORTx algorithm allows specific discrimination of 22 human hematopoietic cell phenotypes. CIBERSORTx deconvolution algorithm uses a set of reference gene expression values (a signature with 547 genes) and, based on those values, infers gene expression profiles and provides an estimation of the abundances of member cell types in a mixed cell population, including in data from bulk tumor samples (24). CIBERSORTx can be run in two modes, in our study “Absolute mode” was used because this mode allows comparison between samples, output absolute levels are represented as scores of arbitrary units that reflects the absolute proportion of each cell type.

To perform TME analysis with EPIC (25) method, we used the web interface (http://epic.gfellerlab.org), “tumor-infiltrating cells” option was used as reference profile. EPIC includes RNA-seq–based gene expression reference profiles from immune cells and other nonmalignant cell types found in tumors: B cells, cancer-associated fibroblasts (CAF), T cells CD4+, T cells CD8+, endothelial cells, macrophages, and NK cells. In addition to these cell types, EPIC returns the fraction of “otherCells” what represents the remaining cells that mainly correspond to cancer cells when applying EPIC on tumor samples. EPIC output values can be interpreted as a cell fraction relative to all cells in sample (the sum of all cell fractions is always equal to 1 per sample).

The quanTIseq method quantifies the fractions of 10 immune cell types from bulk RNA-seq data, this method was validated in blood and tumor samples using simulated, flow cytometry, and IHC data (27). Along with EPIC, quanTIseq can estimate the absolute fractions of immune cells referred to the whole-cell mixture, but in addition unlike EPIC, quanTIseq also considers immune cells relevant for cancer immunology like regulatory T cells (Treg) cells, dendritic cells, and classically (M1) and alternatively (M2) activated macrophages. In our study, quanTIseq analysis (the software available at http://icbi.at/quantiseq) was performed in “preproc” mode where FASTQ files of raw RNA-seq reads from tumor samples were used as input (option--tumor = TRUE, other parameters were used in default settings). In this mode, RNA-seq data processing was performed by quanTIseq itself, including adapter sequences removal, reads quality control, gene counts and TPM (transcripts per million) generation, and annotation.

MCP-counter (26) algorithm was used as an alternative method to EPIC for estimation of not only immune cells abundance (this method allows to quantify eight immune cell types) but also of two stromal cell populations: CAFs and endothelial cells. MCP-counter analysis was performed using the “MCPcounter” R package (26).

The xCell (28) method infers 64 immune and stromal cell types and provides scores, comparable between samples. For running the xCell method, http://xcell.ucsf.edu/ web server was used. For xCell and all methods above (except quanTIseq) TPM-transformed matrix generated by quanTIseq was used as input.

The difference in the noncancerous cell fractions between primary and recurrent aGCTs was calculated using the Wilcoxon test, cell fractions with P < 0.05 considered as significantly different between conditions. For all TME computational deconvolution method results, we performed Pearson correlation analysis on the fractions of noncancerous cells infiltrating studied samples.

Independent validation datasets

For the independent external validation, we obtained previously published publicly available pre-processed bulk RNA-seq data for 10 freshly frozen aGCT tissue samples (6 primary and 4 recurrent) from Mendeley Data (https://data.mendeley.com/datasets/vfvp9gcpzd/draft?a=989a06cc-2a6e-4909-a589-9d150c308b90; ref. 10). Available transcript abundances were converted into a gene-level TPM matrix using the “tximport” R package (RRID:SCR_016752; ref. 30). We also obtained estimated transcript abundances for bulk RNA-seq data generated for 34 fresh-frozen aGCTs by Pilsworth and colleagues (31) from the authors through personal communication. All gene expression values expressed as TPM were provided to TME deconvolution tools. Differences in library preparation and sequencing platforms across study precluded data integration at the level of transcript abundance estimates. However, for each separate dataset the fold change as the ratio between recurrent and primary scores/fractions for each cell type and each method was calculated and compared between datasets to evaluate the consistency of the results from different studies.

Data availability

RNA-seq data are available through the EGA database (EGAS00001006478).

Differential expression analysis: recurrent versus primary tumors

Differential gene expression analysis comparing primary with recurrent aGCTs identified a total of 31 DEGs (25 protein-coding genes and 6 long non-coding RNA transcripts) for which expression was significantly different between two studied groups (adjusted P < 0.05, fold change > 2). Eighteen transcripts showed increased expression and 13 genes had decreased expression in the recurrent tumors when compared with the primary tumors. Differential gene expression analysis results are presented as an expression heatmap of unsupervised hierarchical clustering of all 31 DEGs along with clinicopathologic annotations (Fig. 1A) and as a volcano plot of DEGs between recurrent and primary aGCT samples (Fig. 1B). The clinical summary for studied cohort is presented in Table 1, the detailed description for each patient presented as a Supplementary Data (Supplementary Table S2).

Figure 1.

Differential gene expression analysis comparing primary with recurrent aGCTs identified a total of 31 DEGs, including genes involved in ovarian hormone signaling. A, Heatmap of unsupervised hierarchical clustering of 31 DEGs. Rows represent genes, and columns represent aGCT samples. The variance stabilizing transformation was applied for expression data. Colors of heatmap display gene expression Z-scores that are computed on a gene-by-gene basis: Increased expression (red) and decreased expression (blue). Asterisk identifies sample without FOXL2 c.C402G mutation. B, Volcano plot of DEGs between recurrent and primary aGCT samples. The x-axis indicates the fold change (Log2 scaled), the y-axis shows the adjusted P value (Log10 scaled). Each point represents a different gene, the red/blue color of the points categorize the genes with increased/decreased expression in recurrent tumors (DESeq2, adjusted P < 0.05 is considered as statistically significant). The genes with the highest fold change and genes involved in ovarian sex cord hormonal signaling are highlighted in bold and broadly discussed in the article. C and D, Violin plots showing the expression levels of key genes with increased expression in primary (C) and recurrent (D) tumors. aGCT, adult-type granulosa cell tumor; DEGs, differentially expressed genes.

Figure 1.

Differential gene expression analysis comparing primary with recurrent aGCTs identified a total of 31 DEGs, including genes involved in ovarian hormone signaling. A, Heatmap of unsupervised hierarchical clustering of 31 DEGs. Rows represent genes, and columns represent aGCT samples. The variance stabilizing transformation was applied for expression data. Colors of heatmap display gene expression Z-scores that are computed on a gene-by-gene basis: Increased expression (red) and decreased expression (blue). Asterisk identifies sample without FOXL2 c.C402G mutation. B, Volcano plot of DEGs between recurrent and primary aGCT samples. The x-axis indicates the fold change (Log2 scaled), the y-axis shows the adjusted P value (Log10 scaled). Each point represents a different gene, the red/blue color of the points categorize the genes with increased/decreased expression in recurrent tumors (DESeq2, adjusted P < 0.05 is considered as statistically significant). The genes with the highest fold change and genes involved in ovarian sex cord hormonal signaling are highlighted in bold and broadly discussed in the article. C and D, Violin plots showing the expression levels of key genes with increased expression in primary (C) and recurrent (D) tumors. aGCT, adult-type granulosa cell tumor; DEGs, differentially expressed genes.

Close modal
Table 1.

Characteristics and clinical features of patients with primary and recurrent aGCT.

Primary (n = 8)Recurrent (n = 16)
CharacteristicN (%)N (%)P
Age at diagnosis     0.009 
N 16  
 Mean (SD) 58.38 (8.65) 44.75 (11.99)  
 Median (min–max) 62.00 (45.00–70.00) 43.50 (32.00–76.00)  
Tumor size at diagnosis (mm)     0.386 
N 11  
 Mean (SD) 163.50 (124.38) 104.55 (37.98)  
 Median (min–max) 145.00 (45.00–400.00) 100.00 (30.00–160.00)  
Tumor size at sample collection (mm)     0.209 
 N 16  
 Mean (SD) 163.50 (124.38) 97.69 (58.43)  
 Median (min–max) 145.00 (45.00–400.00) 92.50 (40.00–260.00)  
Menopause status at diagnosis     0.024 
 Post-menopausea 87.50 28.57  
 Pre-menopause 12.50 10 71.43  
 Unknown    
Initial stage at diagnosis     0.618 
 Advanced stage 25.00 15.38  
 Early stage 75.00 11 84.62  
 Unknown    
Anatomical tumor site     <0.001 
 Abdomen pelvis 0.00 10 62.50  
 Large Bowel 0.00 12.50  
 Liver 0.00 6.25  
 Ovary 100.00 12.50  
 Spleen 0.00 6.25  
Primary (n = 8)Recurrent (n = 16)
CharacteristicN (%)N (%)P
Age at diagnosis     0.009 
N 16  
 Mean (SD) 58.38 (8.65) 44.75 (11.99)  
 Median (min–max) 62.00 (45.00–70.00) 43.50 (32.00–76.00)  
Tumor size at diagnosis (mm)     0.386 
N 11  
 Mean (SD) 163.50 (124.38) 104.55 (37.98)  
 Median (min–max) 145.00 (45.00–400.00) 100.00 (30.00–160.00)  
Tumor size at sample collection (mm)     0.209 
 N 16  
 Mean (SD) 163.50 (124.38) 97.69 (58.43)  
 Median (min–max) 145.00 (45.00–400.00) 92.50 (40.00–260.00)  
Menopause status at diagnosis     0.024 
 Post-menopausea 87.50 28.57  
 Pre-menopause 12.50 10 71.43  
 Unknown    
Initial stage at diagnosis     0.618 
 Advanced stage 25.00 15.38  
 Early stage 75.00 11 84.62  
 Unknown    
Anatomical tumor site     <0.001 
 Abdomen pelvis 0.00 10 62.50  
 Large Bowel 0.00 12.50  
 Liver 0.00 6.25  
 Ovary 100.00 12.50  
 Spleen 0.00 6.25  

Note: Bold text indicates P values of <0.05, which were considered as significant.

aOne caused by surgery.

Among the genes with significantly increased expression in primary tumors were insulin-like peptide 3 (INSL3) and luteinizing hormone/choriogonadotropin receptor (LHCGR; Fig. 1C). Genes with increased expression in recurrent tumors included neural EGF-like 2 (NELL2), IL1 receptor like 1 (IL1RL1), and aromatase (CYP19A1; Fig. 1D). The full list of 31 DEGs is shown in Supplementary Table S3. To have a larger overview of the transcriptomic profile between primary and recurrent tumors, a heatmap of unsupervised hierarchical clustering of the top 1,000 most variable genes is presented in Supplementary Fig. S1. Comparative gene expression plots for DEGs identified in previous studies (10, 11) are presented in Supplementary Fig. S2. Because previous exome-based analysis of copy-number alterations identified recurrent focal deletions of 21q22.3 and 3p29 in this cohort of 24 aGCT samples (9/24 tumors, 38% and 7/24 tumors, 29%, respectively; ref. 6), we determined the positions of DEGs. None of the identified DEGs were in 3p29 and 21q22.3 focal deletion regions.

Functional annotation and pathway enrichment analysis

To better understand the underlying biological processes modulated by these 31 DEGs, we performed an ORA with respect to the GO terms (21) and KEGG pathways for identified DEGs. No significantly enriched (adjusted P < 0.05) GO terms and KEGG pathways were identified for set of 31 DEGs through ORA, likely due to small gene set input.

Next, GSEA was performed to compare primary and recurrent aGCT samples using sets of genes grouped by categories. This analysis yielded several gene sets, most of which had genes with increased expression in recurrent tumors. The top 10 GSEA GO terms for biological processes, molecular functions, and cellular components are presented as a dot plot (Supplementary Fig. S3A). Regulation of hormone levels (GO:0010817), humoral immune response (GO:0006959), and steroid metabolic process (GO:0008202) were the top three enriched biological processes. Blood microparticle (GO:0072562) was the most significantly enriched cellular component.

GSEA performed with gene sets from KEGG pathways identified significantly enriched gene sets in recurrent tumors. The top three KEGG pathways were Complement and coagulation cascades (hsa04610, adjusted P = 2.14×10−18), retinol metabolism (hsa00830, adjusted P = 1.82×10−13), and drug metabolism–cytochrome P450 (hsa00982, adjusted P = 1.82×10−13; Supplementary Fig. S3B; GSEA plots with the top 5 gene sets also presented in Supplementary Fig. S4).

In addition, calculation of enrichment scores on single sample basis in ssGSEA performed with “GSVA” R package and consequent primary versus recurrent groups comparison were performed. After applying the ssGSEA method for all available MsigDB gene sets, comparison between primary and recurrent groups revealed statistically significant difference only in one gene set MIR758_5P from collection “C3: regulatory target gene sets,” which contain genes predicted to be targets of microRNA hsa-miR-758-5p. Knowledge of microRNA hsa-miR-758-5p and its role in cancer progression are limited. This miRNA was differentially expressed in metastatic compared with non-metastatic seminomatous germ cell tumors (32). However, the mechanistic connection between expression of this microRNA and aGCT recurrence is still unclear. Further investigation of microRNA hsa-miR-758-5p and its target genes in the context of aGCT pathology is needed.

Analysis of tumor microenvironment

Tumor-infiltrating non-cancerous cells can shape the TME and are closely linked to metastasis potential, progression, and recurrences. We therefore performed an integrated analysis of TME composition, estimated using an integrated analysis of several validated tools. Although these computational methods have been comprehensively benchmarked (33), each of them can identify different cell subtypes with varying degrees of accuracy. To circumvent this limitation, we therefore estimated the fractions of noncancerous cells in aGCT samples using several available computational algorithms and compared their results.

Using EPIC (25) results, we identified statistically significant difference between primary and recurrent tumors in fractions of endothelial cells (P = 0.011) and CAFs (P = 0.045). There was no statistically significant difference between fractions of “otherCells” that mainly correspond to cancer cells between the two conditions. EPIC results are presented as stacked bar plot showing the proportion of cells infiltrating 24 aGCT samples and boxplot showing the fraction of noncancerous cell type identified by the EPIC method in primary and recurrent tumor samples (Fig. 2A and B). The correlation analysis based on EPIC results showed that the fractions of CAFs and endothelial cells were positively correlated (Fig. 2C). Fractions of T cells CD4+ were also positively correlated with endothelial cells.

Figure 2.

EPIC method results display noncancerous cell fractions in aGCT samples and lower content of stromal components in recurrent tumors in comparison with primary tumors. A, Stacked bar plot showing the proportion of noncancerous cells infiltrating 24 aGCT samples. Column names of the plot are the aGCT sample ID. Asterisk identifies sample without FOXL2 c.C402G mutation. Sample order is same as in Fig. 1A. B, Boxplot showing the fraction of each cell type identified by the EPIC method in primary and recurrent tumor samples. The difference in the noncancerous cell fractions between these two groups was calculated using the Wilcoxon test; *, P < 0.05. C, Heatmap showing the correlation between noncancerous cells fractions. The shade of each color circle represents the corresponding correlation value between two cell types, numbers in each circle, indicating the Pearson correlation coefficient. aGCT, adult-type granulosa cell tumor; CAFs, cancer-associated fibroblasts.

Figure 2.

EPIC method results display noncancerous cell fractions in aGCT samples and lower content of stromal components in recurrent tumors in comparison with primary tumors. A, Stacked bar plot showing the proportion of noncancerous cells infiltrating 24 aGCT samples. Column names of the plot are the aGCT sample ID. Asterisk identifies sample without FOXL2 c.C402G mutation. Sample order is same as in Fig. 1A. B, Boxplot showing the fraction of each cell type identified by the EPIC method in primary and recurrent tumor samples. The difference in the noncancerous cell fractions between these two groups was calculated using the Wilcoxon test; *, P < 0.05. C, Heatmap showing the correlation between noncancerous cells fractions. The shade of each color circle represents the corresponding correlation value between two cell types, numbers in each circle, indicating the Pearson correlation coefficient. aGCT, adult-type granulosa cell tumor; CAFs, cancer-associated fibroblasts.

Close modal

Applying another deconvolution-based approach CIBERSORTx (24) in ”Absolute mode,” we found statistically significant difference in infiltration levels of resting CD4+ memory T cells, they had significantly higher fractions in recurrent tumors (P = 0.001, Supplementary Fig. S5A).

Using xCell algorithm (28), which allows to identify 64 immune and non-immune cell types, we found that the proportions of endothelial cells, hematopoietic and mesenchymal stem cells, were significantly lower in recurrent tumors than those in primary tumors (Supplementary Fig. S5B). Using results provided by xCell, correlation analysis reveals a cluster of several positively correlated cells (Supplementary Fig. S6A). CIBERSORTx method results showed positive correlation between neutrophils and resting NK cells (Supplementary Fig. S6B).

With the estimation of quanTIseq (27), the infiltration levels of B cells, M1 and M2 macrophages, monocytes, neutrophils, NK cells, non-regulatory CD4+ T cells, CD8+ T cells, Treg cells, and myeloid dendritic cells in two groups of aGCTs were retrieved. According to quanTIseq results the infiltration levels of neutrophils were significantly higher in recurrent tumors than those in the primary tumors (P = 0.022; (Supplementary Fig. S7A and S7B). Dendritic cell infiltration levels were lower in recurrent tumors and identified also as significantly different between two conditions (P = 0.009). There was no significant difference of other immune cell types between the two groups. According to quanTIseq results, macrophages M1 positively correlated with neutrophils, NK cells positively correlated with Tregs, T cells CD8+ positively correlated with monocytes (Supplementary Fig. S7C).

MCP-counter (26) showed no statistically significant difference between primary and recurrent samples in infiltration levels of considered groups of cells (Supplementary Fig. S8A); however, there was a trend of decreasing CAFs and endothelial cells (P = 0.12 and P = 0.14, respectively) and increasing neutrophils (P = 0.11) abundance in recurrent tumors. Correlation between CAFs and endothelial cells was seen with MCP-counter results (P = 0.005). Levels of neutrophils correlated with monocytes (P = 6.08×10−4), as well as cytotoxic lymphocytes correlated with T cells CD8+ (P = 8.47×10−6; Supplementary Fig. S8B).

Because the TME gene set scores were unique to each deconvolution method, we compared and evaluated the consistency of the outputs across methods by calculating the fold change as the ratio between recurrent and primary scores/fractions for each cell type and each method. Various noncancerous cell estimates associated with primary or recurrent tumors are presented in Fig. 3. Tumor infiltration analysis by different tools showed consistent results for CAFs, endothelial cells, neutrophils, and macrophages; moreover, applying several different methods, fractions of these cell types were significantly different between primary and recurrent conditions.

Figure 3.

Various noncancerous cells estimates preferentially associated with primary tumors (blue) or recurrent tumors (red) showing consistent infiltration results for neutrophils, macrophages M1, endothelial cells, and cancer-associated fibroblasts. Tumor microenvironment infiltration results based on CIBERSORTx (absolute mode), EPIC, MPC-counter, quanTIseq (“preproc” mode; FASTA files with raw reads used), and xCell predictions. Only cell types identified by at least two methods are presented. Only direct measurements were considered. Benchmark recommended (Sturm et al, 2019) cell-type/algorithm estimates highlighted in gray. The difference in the noncancerous cell fractions between primary and recurrent groups was calculated using the Wilcoxon test.

Figure 3.

Various noncancerous cells estimates preferentially associated with primary tumors (blue) or recurrent tumors (red) showing consistent infiltration results for neutrophils, macrophages M1, endothelial cells, and cancer-associated fibroblasts. Tumor microenvironment infiltration results based on CIBERSORTx (absolute mode), EPIC, MPC-counter, quanTIseq (“preproc” mode; FASTA files with raw reads used), and xCell predictions. Only cell types identified by at least two methods are presented. Only direct measurements were considered. Benchmark recommended (Sturm et al, 2019) cell-type/algorithm estimates highlighted in gray. The difference in the noncancerous cell fractions between primary and recurrent groups was calculated using the Wilcoxon test.

Close modal

To explore the question of whether tumor sample purity is correlated with CAF abundance estimates, we investigated correlations between the fraction of CAFs and tumor purity using two orthogonal sequencing-based metrics: Raw allelic frequency of FOXL2 c.C402G mutation (variant allele frequency, VAF) and tumor purity as estimated by Sequenza based on whole-exome local copy number and SNP B allele frequency (ref. 34; Supplementary Fig. S9A–S9F). The raw VAF of the heterozygous FOXL2 c.C402G mutation would be presumed to be 0.5 in a completely pure aGCT tumor sample. Across the three TME deconvolution algorithms that detect CAFs, we did not identify statistically significant correlations between CAFs fraction and either method for estimating tumor purity (P > 0.05 for all comparisons), except between xCell CAFs estimations and VAF of FOXL2 c.C402G mutation (P = 0.0025). For FOXL2 c.C402G VAF correlation estimates, the sample without a detectable FOXL2 c.C402G mutation (GCT008) was excluded along with a sample where the VAF was 0.9 (GCT001) as this likely reflects FOXL2 copy-number changes that confound a simple relationship between FOXL2 c.C402G VAF and tumor purity. Correlation between CYP19A1 expression and tumor purity estimated by either method also was not found (Supplementary Fig. S9G and S9H).

We next validated our results using two independent datasets from Haltia and colleagues (10) and Pilsworth and colleagues (31) studies, which provided bulk RNA-seq data for fresh-frozen primary and recurrent aGCT samples. Notably, a high degree of consistency was seen among cross-study estimates of the fold change in CAF between primary and recurrent aGCTs with a mean fold change estimate of 0.38 (range, 0.28–0.57) for EPIC (Fig. 4). Similar results were also obtained for cross-study estimate of CAF fold change using MCP-Counter, with a mean fold change estimate of 0.46 (range, 0.23–0.71; Supplementary Fig. S10).

Figure 4.

Fold change by noncancerous cell-type estimates between recurrent and primary tumors calculated by EPIC demonstrates depletion of cancer-associated fibroblasts fractions in recurrent tumors in three independent aGCT cohorts. Black lines indicate range of observed fold-change estimates for recurrent versus primary tumors. Area of data points is proportional to cohort size.

Figure 4.

Fold change by noncancerous cell-type estimates between recurrent and primary tumors calculated by EPIC demonstrates depletion of cancer-associated fibroblasts fractions in recurrent tumors in three independent aGCT cohorts. Black lines indicate range of observed fold-change estimates for recurrent versus primary tumors. Area of data points is proportional to cohort size.

Close modal

The purpose of this study was to leverage RNA-seq data to perform a robust and comprehensive comparative analysis of TME composition in primary and recurrent aGCT. Comparison of the gene expression profiles revealed a relatively small number of genes that differ significantly between the primary and recurrent tumors. This relative homogeneity across conditions may reflect the shared etiology with respect to cell type of origin and the presence of FOXL2 c.C402G mutation, which can be considered as a shared aGCT-initiating event. This homogeneity is consistent with observations in another comparative study of primary and recurrent aGCTs (11).

DEGs: higher expression in primary tumors

Of the 13 genes with higher expression in primary aGCTs, INSL3 is notable with around 64-fold lower expression in recurrent tumors (Fig. 1). Of note, a previous study comparing gene expression between stages 1 and 3 aGCTs determined that expression of INSL3 was approximately 70-fold higher in stage 1 aGCT in comparison with stage 3 (9). INSL3 is a member of the relaxin family of neo-hormones that are involved in reproduction-related processes in mammals (35). INSL3 is largely produced by the theca cells, and it is thought to be essential for androstenedione synthesis, which is the major steroid precursor for the granulosa cells to produce estrogens. LHCGR is expressed in theca, interstitial, differentiated granulosa, and luteal cells in the normal ovary (36) and in our study also showed significant 9-fold higher expression in primary compared with recurrent tumors. Because both INSL3 and LHCGR genes have been implicated in ovarian development, this decreased expression in recurrent tumors may reflect pathways of de-differentiation in which normal hormonal regulatory pathways and biosynthetic functions become altered during tumor progression. Alternatively, samples from primary aGCTs were all derived from ovarian sites (Table 1) and may include some portion of normal ovarian sex cords, which would not be expected to be found in recurrent tumors from metastatic sites. Thus, certain genes may seem enriched using bulk RNA-seq due to differences in the admixture of normal and cancerous cells.

aGCTs can produce estradiol; this hormone could play a role in this cancer through its nuclear receptors, that is, ERα (encoded by ESR1) and ERβ (encoded by ESR2). Estradiol mediates its effects through ERα-dependent genomic mechanisms and ERβ/ERα-dependent extra-nuclear mechanisms. Although it was previously shown that estradiol can promote progression of GCTs, with a clear dependence on ERα (37), our analysis did not identify a significant difference in ESR1 expression level between primary and recurrent aGCTs (Supplementary Fig. S2B). Moreover, we did not observe statistically significant differences in expression levels of genes involved in estrogen signaling that showed higher expression in primary tumors compared with recurrent tumors in Haltia and colleagues’ study (ref. 10; Supplementary Fig. S2A). Those differences in observations may be due to the differences between the cohorts and cohort size (our cohort is almost 2.5 times bigger: 10 Haltia and colleagues’ samples used for RNA-seq vs. 24 samples in current study) as well as algorithms underlying methods that were used for RNA-seq data processing. Moreover, in differential expression analysis as a threshold for statistical significance Haltia and colleagues (10) used an adjusted P value of <0.25. In our work, we used more strict criteria and a value of 0.05. With threshold 0.25 in our work CREB3L1 that was identified as DEG in Haltia and colleagues’ study also could be considered as a gene with statistically significant differences in expression levels between primary and recurrent tumor groups (adjusted P = 0.16). In accordance with this prior work, we found that ESR2 was the most highly expressed of the estrogen receptors (Supplementary Fig. S2B).

In primary aGCTs, we found significantly higher expression of long noncoding RNA LINC01116 that promotes the development of breast cancer by binding directly to miR-145, resulting in the upregulation of ESR1, a target gene of miR-145 (38). Interestingly, it has been shown that overexpressed LINC01116 promotes epithelial ovarian cancer progression via increasing proliferation and migration, and inhibiting cell apoptosis (39).

The expression levels of genes ASS1 and PLIN4, which were increased in primary aGCTs in a prior study (11), were not significantly different in our study (Supplementary Fig. S2C). These differences can be due to the differences between the studied cohorts or formalin fixation transcriptomic effects in archival tissue samples. Unlike the previous Andersson's study (11) where archival aGCT FFPE samples were used, in our work we had fresh-frozen tissue. This is a major difference—FFPE artifacts can significantly confound RNA-seq analysis. We should emphasize that our samples were all cryopreserved and not formalin fixed.

DEGs: higher expression in recurrent tumors

Of the 18 genes with high expression in recurrent tumors, NELL2, GDF6, and IL1RL1 show the greatest increase (more than 16-fold). NELL2 is involved in the onset of female puberty by regulating the release of gonadotropin-releasing hormone, facilitating the maintenance of the normal reproductive cycle in mammals (40). Information about the role of NELL2 in ovarian cancer and specifically in aGCTs is limited. However, studies in non–small cell lung cancer (41) and Ewing sarcoma (42) provide some evidence that NELL2 may be involved in cell proliferation regulation and metastasis. GDF6 is involved in follicular development and it was previously shown that increases in GDF6 expression were associated with increased metastasis and decreased survival rate in patients with melanoma (43).

IL1RL1 (also known as ST2) is a member of the IL1 cytokine family, with IL-33 as its ligand. Studies performed by Tong and colleagues (44) have shown that both IL33 and IL1RL1 were highly upregulated in epithelial ovarian cancer tumors compared with normal ovary and ovarian benign tumors, and the expression levels were further increased in tumor tissues at the metastatic site. In our study, only IL1RL1 expression level was significantly elevated in recurrent tumors in comparison with primary tumors, we have not found statistically significant difference in expression levels of IL33 between considered groups. Elevated levels of IL1RL1 soluble form (referred to as soluble ST2 or sST2) have also been detected and associated with poor prognosis in hepatocellular carcinoma (45) and breast cancer (46), suggesting that IL1RL1 is likely to be involved in the development of a variety of cancers.

In a previous study of aGCTs, CYP19A1 (also known as aromatase) expression was highly variable among primary and recurrent tumors with no specific trend (10), but in our work, CYP19A1 showed statistically significant higher expression in recurrent than in primary tumors. This is in contrast with other genes involved in normal ovarian sex cord hormonal signaling (INSL3 and LHCGR) in which higher expression was seen in primary tumors. It has previously been shown that CYP19A1 expression can change in response to hormonal modulation, such as with systemic aromatase inhibitor treatment, and that germline SNPs can also influence this effect (47). We speculate that a complex interplay between genotype and exposure to hormonal therapies may account for the aggregate increase in CYP19A1 expression we observed in recurrent versus primary tumors. The treatment history of patients in this cohort was heterogeneous, as is common for recurrent aGCTs where standard-of-care approaches are limited, and thus a complete exploration of this phenomenon awaits a future, focused analysis of a larger cohort.

Difference in TME between primary and recurrent aGCTs

Tumors are not only composed of cancer malignant cells but are embedded in a complex network of stromal cells and immune cells. It also comprises the extracellular matrix, which consists of chemokines, inflammatory cytokines, and other secreted molecules. Tumor-infiltrating noncancerous cells can shape the TME and are closely linked to tumor initiation, growth, metastasis potential, and recurrences. In recent years, the TME has been reported to play a vital role in the tumorigenesis of ovarian cancer (48). However, the role of TME in aGCT development is poorly delineated.

This tumor infiltration analysis using multiple robust TME deconvolution tools showed consistent trends for myeloid-derived cells (neutrophils and macrophages), endothelial cells, and CAFs (Fig. 3). Fractions of neutrophils and macrophages were higher in recurrent tumors compared with primary tumors, whereas fractions of CAFs and endothelial cells were lower. Of these observations, statistically significant differences in CAF abundance were confirmed using additional independent datasets.

Macrophages and neutrophils are immune cells that in recent years have attracted increasing attention because of their cancer-promoting effects (49, 50). In tumors stimulated by different factors, monocytes can differentiate into two main phenotypes: Antitumor M1 macrophages or pro-tumor M2 macrophages (51, 52). In our study, different deconvolution tools were able to identify specific subtypes of macrophages, demonstrating that fractions of macrophages are increased in recurrent tumors. Although in our 24 samples cohort no one method showed statistically significant differences between fractions of macrophages and macrophages M1 between primary and recurrent conditions, xCell, quanTIseq, and CIBERSORTx showed the tendency in the same direction for macrophages M1. We believe that macrophages (both M1 and M2) should be the subject of further investigation because the increase of macrophages in recurrent aGCTs was observed not only in the current work but also in the study made on 14 patients with aGCT (FFPE tumor samples: 6 primary and 8 recurrent tumors) presented in the AACR abstract (53). In line with our results, in the mentioned work an integrated multiplexing approach for the immunoprofiling of the TME of aGCT showed an increase in macrophage density in recurrent tumors compared with primary tumors. However, an M2 phenotype and not an M1 was identified.

Neutrophils also contribute to different cancer types (50), including ovarian cancer metastasis (54). In an ovarian cancer meta-analysis study consisting of 12 studies and 3,854 patients, elevated pretreatment neutrophil-to-lymphocyte ratio has been proposed to serve as a predicative factor of poor prognosis for patients with ovarian cancer (55). In our study, neutrophil fraction was significantly higher in recurrent tumors in comparison with primary tumors according to the quanTIseq tool.

In contrast with macrophages and neutrophils, CAFs and endothelial cells showed significantly lower infiltration level of recurrent tumors in comparison with primary tumors. CAFs are a group of heterogeneous cells and are the most prevalent stromal components of the TME (56). In different ovarian cancer types, high abundance of CAFs was associated with poor prognosis, and previously it was shown that CAFs promote tumor progression via various mechanisms. CAFs cells contribute to ovarian cancer initiation and development (57), cancer growth (58), and resistance to treatment (59). Specific subpopulations of CAFs accumulation in mesenchymal high-grade serous ovarian cancers are associated with an immunosuppressive environment (60). Interestingly, in contrast with most of the previous studies on different ovarian cancer types, our results showed that the high abundance of CAFs was associated with primary tumors not with recurrent tumors, a finding that we found to be consistent across two additional independent datasets.

It is conceivable that a lower number of CAFs in the recurrent tumor could be a result of a different tumor/stroma ratio with purer tumors on recurrence; however, we would like to emphasize that in the present study we used only tumor samples with an estimated tumor cell content >60% without significant necrosis. Also, we have not observed a significant correlation between CAFs and multiple sequencing-based tumor purity metrics (Supplementary Fig. S9A–S9F).

Because CAFs consist of heterogeneous cell populations, using computational approaches we could detect a mixture consisting of a broad range of CAFs that express different markers. It is crucial to note that these markers are not specific and many CAFs may not express all these markers at the same time, reflecting the high degree of heterogeneity and possible opposite CAFs functions in the context of TME. On the basis of our computational deconvolution results, we cannot distinguish specific CAF subtypes responsible for the difference between primary and recurrent aGCT groups but despite most of the previous studies providing evidence of CAFs pro-tumor effect, there are also reports about the antitumor role of some CAFs subtypes in different cancer types. For example, a diminished CD68+ CAF subset contributes to the poor prognosis of patients with oral squamous cell carcinoma (61). A study focusing on Shh-deficient signaling demonstrated a protective role of αSMA+ myofibroblasts and indicated that depletion of CAFs induces immunosuppression and accelerates pancreas cancer with reduced survival (62). Similar depletion of antitumor CAFs observed in our study can lead to TME remodeling and may correlate with aGCTs recurrence. Considering the above, further investigation of the functions of CAFs and their different subtypes in the context of aGCTs TME and the roles they may conceivably play in immune surveillance seems to be an interesting and promising area for future investigation.

As CAFs, endothelial cells and specifically microvascular endothelial cells also showed significantly lower fractions in recurrent tumors in comparison with primary tumors according to most of the used methods. The reason for the decrease in the level of these cells in the recurrent tumors is not obvious. Clinical relevance of microvascular endothelial cells in aGCTs previously has not been reported. In colorectal cancer, abundance of microvascular endothelial cells is associated with worse survival in the cohort that did not receive chemotherapy, whereas patients with microvascular endothelial cell abundant colorectal cancer better respond to chemotherapy (63). In our study, only the xCell method was able to identify microvascular endothelial cells. A reduced fraction of microvascular endothelial cells in recurrent aGCTs may be related to worse response to chemotherapy. This observation warrants further study.

It also should be noted that some of the cell types estimated by the xCell algorithm (28) were not expected to be in aGCT samples such as osteoblasts, hepatocytes, melanocytes, etc. (Supplementary Fig. S5B). One explanation of such findings could be that the xCell algorithm is based on composite marker genes and results must be interpreted within the context of the tumor specimen or tissue type being studied.

On the basis of the results above, we conclude that immune response and related processes play a major role in aGCT progression and malignancy; however, the specific contribution of noncancerous cells to aGCT progression is still not clear. Comparing different methods and considering benchmark studies, we can see that it is not easy to detect fractions of myeloid cells only by computational approach, but our findings suggest that myeloid cells fractions are higher in recurrent tumors. The role of stromal components in aGCT recurrency is still unclear but our results and results of external validation clearly show that CAFs (or at least some CAFs populations) are less abundant in recurrent tumors.

Strengths and limitations of this study

A significant strength of this study is the utilization of a relatively large cohort of well-annotated aGCT samples collected as fresh-frozen tissue. Not only does this facilitate comparative analysis of a rare cancer type but also eliminates transcriptomic artifacts associated with sequencing formalin-fixed, paraffin-embedded samples. Another significant strength of this study is the integrative analysis of TME composition using multiple computational methods, which has been omprehensively benchmarked and showed reliability for bulk RNA-seq data deconvolution (33). When these methods are concordant, it suggests consistent and robust differences in the TME between primary and recurrent tumors.

However, the results obtained by computational methods for macrophages, neutrophils, and CAFs should be interpreted with caution. Some macrophages can share features of both M1 and M2 phenotypes and co-express M1/M2 markers (64). The identification of macrophages based solely on classical M1/M2 expression markers can lead to misinterpretation of the true role of macrophages in recurrent aGCTs, and additional validation should be applied. Neutrophils in TME also can be polarized to antitumor N1 or pro-tumor N2 phenotype (65). This can complicate the interpretation of the data obtained for neutrophils as well. In addition, some tools indicated that the fractions of macrophages and neutrophils were relatively low, and the results obtained for these cell types in different validation cohorts are controversial.

CAFs are heterogeneous cell populations. Therefore, it is essential to classify CAFs by combining different markers to identify their biological characteristics. In our study, we only report that CAFs fractions identified by computational methods are lower in recurrent tumors, but we cannot make any robust conclusions about specific subtypes of CAFs because the methods we used do not provide those predictions. Another limitation is differences among patients in this cohort regarding the number and type of prior treatment lines.

This is fundamentally a correlative study and our results do not preclude other clinical or pathologic variables from exhibiting co-variance with either recurrence or specific changes in the TME (Table 1). Indeed, it is possible that differences in age in diagnosis we observe between patients with primary and recurrent tumors in our cohort could be related to recurrence risk through a mechanism that is mediated by changes in the tumor microenvironment. Further exploration of clinical and molecular co-variates will be a focus of future studies.

Finally, a limitation of our study is the use of non-matched aGCT samples. Matched primary-recurrent samples would facilitate a more direct comparison of gene expression dynamics over the growth and progression of a single tumor, and such a study will be an important future area of investigation. Similarly, technological advancement in single-cell RNA-seq makes this an attractive platform for future investigation of aGCT gene expression heterogeneity and TME. However, these approaches typically require dissociation of fresh tissue before cryopreservation and thus were not feasible for use on the banked tissue samples available for this study. Methods such as imaging mass cytometry or CosMX Spatial Molecular Imaging as well allow the identification of distinct cell types and unique TME internal structure. Finally, all data in this study are observational and correlative in nature. Future studies are needed to determine mechanistic links between TME composition, tumor differential gene expression, and propensity for recurrence.

Conclusions

The TME composition is significantly different between primary and recurrent aGCTs, with recurrent tumors having increased myeloid fraction (neutrophils and macrophages) in conjunction with reduced prevalence of stromal cells such as CAFs and endothelial cells. Rigorous differential gene expression analysis identifies genes involved in ovarian hormone signaling (INSL3, CYP19A1, and LHCGR) as significantly altered in recurrences, suggesting that changes in hormonal signaling pathways as well as TME remodeling may underpin aGCT recurrence.

B.M. Fellman reports grants from NIH MDACC support grant during the conduct of the study. R.R. Broaddus reports grants from National Cancer Institute during the conduct of the study. D.M. Gershenson reports grants from Novartis outside the submitted work; Genentech consulting, member of Protocol Steering Committee; non-compensated Aadi advisory board member, one-time personal fee; Onconova advisory board member, one-time personal fee; Verastem advisory board member, one-time personal fee Springworks advisory board member, one-time personal fee; Elsevier: royalty for educational publication; UpToDate: royalty for educational publication; chair, International Consortium for Low-Grade Serous Ovarian Cancer; equity interest: Johnson and Johnson, Bristol Myers Squibb, and Procter and Gamble. R. Hillman reports grants from Sumitomo Pharma Co., Ltd. during the conduct of the study. No disclosures were reported by the other authors.

E. Khlebus: Conceptualization, data curation, software, formal analysis, validation, investigation, visualization, methodology, writing–original draft, project administration. V.K. Vuttaradhi: Investigation. T. Welte: Investigation. N. Khurana: Investigation. J. Celestino: Investigation. H.C. Beird: Investigation, writing–review and editing. C. Gumbs: Investigation. L. Little: Investigation. A.F. Legarreta: Data curation. B.M. Fellman: Formal analysis. T. Nguyen: Investigation. B. Lawson: Investigation. S. Ferri-Borgogno: Resources, validation, writing–review and editing. S.C. Mok: Resources, validation, writing–review and editing. R.R. Broaddus: Investigation, writing–review and editing. D.M. Gershenson: Investigation, writing–review and editing. P.A. Futreal: Funding acquisition. R.T. Hillman: Conceptualization, resources, supervision, funding acquisition, methodology, writing–original draft, project administration.

This work was supported by grants from the Cancer Prevention and Research Institute of Texas (CPRIT) to R.T. Hillman (RR200045) and to P.A. Futreal (R1205). R.T. Hillman is supported by the Jennifer “Jenny” Song Fund for Granulosa Cell Tumor Research. The University of Texas MD Anderson Cancer Center Multidisciplinary Gynecologic Cancer Tumor Bank is supported in part by NIH P50CA83639 SPORE in Ovarian Cancer. This research was in part supported by the National Institutes of Health through M.D. Anderson's Cancer Center Support Grant CA016672. We thank Xingzhi Song and Jianhua Zhang for the help with RNA-sequencing data analysis.

The publication costs of this article were defrayed in part by the payment of publication fees. Therefore, and solely to indicate this fact, this article is hereby marked “advertisement” in accordance with 18 USC section 1734.

Note: Supplementary data for this article are available at Molecular Cancer Research Online (http://mcr.aacrjournals.org/).

1.
Shah
SP
,
Köbel
M
,
Senz
J
,
Morin
RD
,
Clarke
BA
,
Wiegand
KC
, et al
.
Mutation of FOXL2 in granulosa-cell tumors of the ovary
.
N Engl J Med
2009
;
360
:
2719
29
.
2.
Da Cruz Paula
A
,
da Silva
EM
,
Segura
SE
,
Pareja
F
,
Bi
R
,
Selenica
P
, et al
.
Genomic profiling of primary and recurrent adult granulosa cell tumors of the ovary
.
Mod Pathol
2020
;
33
:
1606
17
.
3.
Lauszus
FF
,
Petersen
AC
,
Greisen
J
,
Jakobsen
A
.
Granulosa cell tumor of the ovary: a population-based study of 37 women with stage I disease
.
Gynecol Oncol
2001
;
81
:
456
60
.
4.
Alexiadis
M
,
Rowley
SM
,
Chu
S
,
Leung
DTH
,
Stewart
CJR
,
Amarasinghe
KC
, et al
.
Mutational landscape of ovarian adult granulosa cell tumors from whole-exome and targeted TERT promoter sequencing
.
Mol Cancer Res
2019
;
17
:
177
85
.
5.
Pilsworth
JA
,
Cochrane
DR
,
Xia
Z
,
Aubert
G
,
Färkkilä
AEM
,
Horlings
HM
, et al
.
TERT promoter mutation in adult granulosa cell tumor of the ovary
.
Mod Pathol
2018
;
31
:
1107
15
.
6.
Hillman
RT
,
Celestino
J
,
Terranova
C
,
Beird
HC
,
Gumbs
C
,
Little
L
, et al
.
KMT2D/MLL2 inactivation is associated with recurrence in adult-type granulosa cell tumors of the ovary
.
Nat Commun
2018
;
9
:
1
7
.
7.
Hillman
RT
,
Lin
DI
,
Lawson
B
,
Gershenson
DM
.
Prevalence of predictive biomarkers in a large cohort of molecularly defined adult-type ovarian granulosa cell tumors
.
Gynecol Oncol
2021
;
39
:
5567
.
8.
Kraus
F
,
Dremaux
J
,
Altakfi
W
,
Goux
M
,
Pontois
L
,
Sevestre
H
, et al
.
FOXL2 homozygous genotype and chromosome instability are associated with recurrence in adult granulosa cell tumors of the ovary
.
Oncotarget
2020
;
11
:
419
.
9.
Alexiadis
M
,
Chu
S
,
Leung
D
,
Gould
JA
,
Jobling
T
,
Fuller
PJ
.
Transcriptomic analysis of stage 1 versus advanced adult granulosa cell tumors
.
Oncotarget
2016
;
7
:
14207
19
.
10.
Haltia
UM
,
Pihlajoki
M
,
Andersson
N
,
Mäkinen
L
,
Tapper
J
,
Cervera
A
, et al
.
Functional profiling of FSH and estradiol in ovarian granulosa cell tumors
.
J Endocr Soc
2020
;
4
:
bvaa034
.
11.
Andersson
N
,
Haltia
UM
,
Färkkilä
A
,
Wong
SC
,
Eloranta
K
,
Wilson
DB
, et al
.
Analysis of non-relapsed and relapsed adult type granulosa cell tumors suggests stable transcriptomes during tumor progression
.
Curr Issues Mol Biol
2022
;
44
:
686
98
.
12.
Pierini
S
,
Tanyi
JL
,
Simpkins
F
,
George
E
,
Uribe-Herranz
M
,
Drapkin
R
, et al
.
Ovarian granulosa cell tumor characterization identifies FOXL2 as an immunotherapeutic target
.
JCI insight
2020
;
5
:
e136773
.
13.
Dobin
A
,
Davis
CA
,
Schlesinger
F
,
Drenkow
J
,
Zaleski
C
,
Jha
S
, et al
.
STAR: ultrafast universal RNA-seq aligner
.
Bioinformatics
2013
;
29
:
15
21
.
14.
Anders
S
,
Pyl
PT
,
Huber
W
.
HTSeq—a Python framework to work with high-throughput sequencing data
.
Bioinformatics
2015
;
31
:
166
9
.
15.
Love
MI
,
Huber
W
,
Anders
S
.
Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2
.
Genome Biol
2014
;
15
:
1
21
.
16.
Kolde
R
.
pheatmap: pretty heatmaps
.
2019
.
Available from:
https://CRAN.R-project.org/package=pheatmap.
17.
The Gene Ontology Consortium
.
The gene ontology resource: 20 years and still going strong
.
Nucleic Acids Res
2019
;
47
:
D330
8
.
18.
Kanehisa
M
,
Goto
S
.
KEGG: kyoto encyclopedia of genes and genomes
.
Nucleic Acids Res
2000
;
28
:
27
30
.
19.
Wu
T
,
Hu
E
,
Xu
S
,
Chen
M
,
Guo
P
,
Dai
Z
, et al
.
clusterProfiler 4.0: a universal enrichment tool for interpreting omics data
.
Innov J
2021
;
2
:
100141
.
20.
Subramanian
A
,
Tamayo
P
,
Mootha
VK
,
Mukherjee
S
,
Ebert
BL
,
Gillette
MA
, et al
.
Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles
.
Proc Natl Acad Sci
2005
;
102
:
15545
50
.
21.
Ashburner
M
,
Ball
CA
,
Blake
JA
,
Botstein
D
,
Butler
H
,
Cherry
JM
, et al
.
Gene ontology: tool for the unification of biology
.
Nat Genet
2000
;
25
:
25
29
.
22.
Liberzon
A
,
Birger
C
,
Thorvaldsdóttir
H
,
Ghandi
M
,
Mesirov
JP
,
Tamayo
P
.
The molecular signatures database hallmark gene set collection
.
Cell Syst
2015
;
1
:
417
25
.
23.
Hänzelmann
S
,
Castelo
R
,
Guinney
J
.
GSVA: gene set variation analysis for microarray and RNA-seq data
.
BMC Bioinf
2013
;
14
:
1
15
.
24.
Newman
AM
,
Steen
CB
,
Liu
CL
,
Gentles
AJ
,
Chaudhuri
AA
,
Scherer
F
, et al
.
Determining cell type abundance and expression from bulk tissues with digital cytometry
.
Nat Biotechnol
2019
;
37
:
773
82
.
25.
Racle
J
,
de Jonge
K
,
Baumgaertner
P
,
Speiser
DE
,
Gfeller
D
.
Simultaneous enumeration of cancer and immune cell types from bulk tumor gene expression data
.
Elife
2017
;
6
:
e26476
.
26.
Becht
E
,
Giraldo
NA
,
Lacroix
L
,
Buttard
B
,
Elarouci
N
,
Petitprez
F
, et al
.
Estimating the population abundance of tissue-infiltrating immune and stromal cell populations using gene expression
.
Genome Biol
2016
;
17
:
1
20
.
27.
Finotello
F
,
Mayer
C
,
Plattner
C
,
Laschober
G
,
Rieder
D
,
Hackl
H
, et al
.
Molecular and pharmacological modulators of the tumor immune contexture revealed by deconvolution of RNA-seq data
.
Genome Med
2019
;
11
:
1
20
.
28.
Aran
D
,
Hu
Z
,
Butte
AJ
.
xCell: digitally portraying the tissue cellular heterogeneity landscape
.
Genome Biol
2017
;
18
:
1
14
.
29.
Newman
AM
,
Liu
CL
,
Green
MR
,
Gentles
AJ
,
Feng
W
,
Xu
Y
, et al
.
Robust enumeration of cell subsets from tissue expression profiles
.
Nat Methods
2015
;
12
:
453
7
.
30.
Soneson
C
,
Love
MI
,
Robinson
MD
.
Differential analyses for RNA-seq: transcript-level estimates improve gene-level inferences
.
F1000Res.
2015
;
4
:
1521
.
31.
Pilsworth
JA
,
Todeschini
A-L
,
Neilson
SJ
,
Cochrane
DR
,
Lai
D
,
Anttonen
M
, et al
.
FOXL2 in adult-type granulosa cell tumour of the ovary: oncogene or tumour suppressor gene?
J Pathol
2021
;
255
:
225
31
.
32.
Ernst
S
,
Heinzelmann
J
,
Bohle
RM
,
Weber
G
,
Stöckle
M
,
Junker
K
, et al
.
The metastatic potential of seminomatous germ cell tumours is associated with a specific microRNA pattern
.
Andrology
2020
;
8
:
1687
98
.
33.
Sturm
G
,
Finotello
F
,
Petitprez
F
,
Zhang
JD
,
Baumbach
J
,
Fridman
WH
, et al
.
Comprehensive evaluation of transcriptome-based cell-type quantification methods for immuno-oncology
.
Bioinformatics
2019
;
35
:
i436
i445
.
34.
Favero
F
,
Joshi
T
,
Marquard
AM
,
Birkbak
NJ
,
Krzystanek
M
,
Li
Q
, et al
.
Sequenza: allele-specific copy number and mutation profiles from tumor sequencing data
.
Ann Oncol
2015
;
26
:
64
70
.
35.
Ivell
R
,
Anand-Ivell
R
.
Insulin-like peptide 3 (INSL3) is a major regulator of female reproductive physiology
.
Hum Reprod Update
2018
;
24
:
639
51
.
36.
Ascoli
M
,
Fanelli
F
,
Segaloff
DL
.
The lutropin/choriogonadotropin receptor, a 2002 perspective
.
Endocr Rev
2002
;
23
:
141
74
.
37.
Cluzet
V
,
Devillers
MM
,
Petit
F
,
Pierre
A
,
Giton
F
,
Airaud
E
, et al
.
Estradiol promotes cell survival and induces Greb1 expression in granulosa cell tumors of the ovary through an ERα-dependent mechanism
.
J Pathol
2022
;
256
:
335
48
.
38.
Hu
HB
,
Chen
Q
,
Ding
SQ
.
LncRNA LINC01116 competes with miR-145 for the regulation of ESR1 expression in breast cancer
.
Eur Rev Med Pharmacol Sci
2018
;
22
:
1987
93
.
39.
Fang
YN
,
Huang
ZL
,
Li
H
,
Tan
WB
,
Zhang
QG
,
Wang
L
, et al
.
LINC01116 promotes the progression of epithelial ovarian cancer via regulating cell apoptosis
.
Eur Rev Med Pharmacol Sci
2018
;
22
:
5127
33
.
40.
Ryu
BJ
,
Kim
HR
,
Jeong
JK
,
Lee
BJ
.
Regulation of the female rat estrous cycle by a neural cell-specific epidermal growth factor–like repeat domain containing protein, NELL2
.
Mol Cells
2011
;
32
:
203
7
.
41.
Wang
Y
,
Li
M
,
Zhang
L
,
Chen
Y
,
Zhang
S
.
m6A demethylase FTO induces NELL2 expression by inhibiting E2F1 m6A modification leading to metastasis of non–small cell lung cancer
.
Mol Ther Oncolytics
2021
;
21
:
367
76
.
42.
Jayabal
P
,
Zhou
F
,
Lei
X
,
Ma
X
,
Blackman
B
,
Weintraub
ST
, et al
.
NELL2–cdc42 signaling regulates BAF complexes and Ewing sarcoma cell growth
.
Cell Rep
2021
;
36
:
109254
.
43.
Venkatesan
AM
,
Vyas
R
,
Gramann
AK
,
Dresser
K
,
Gujja
S
,
Bhatnagar
S
, et al
.
Ligand-activated BMP signaling inhibits cell differentiation and death to promote melanoma
.
J Clin Invest
2018
;
128
:
294
308
.
44.
Tong
X
,
Barbour
M
,
Hou
K
,
Gao
C
,
Cao
S
,
Zheng
J
, et al
.
Interleukin-33 predicts poor prognosis and promotes ovarian cancer cell growth and metastasis through regulating ERK and JNK signaling pathways
.
Mol Oncol
2016
;
10
:
113
25
.
45.
Bergis
D
,
Kassis
V
,
Ranglack
A
,
Koeberle
V
,
Piiper
A
,
Kronenberger
B
, et al
.
High serum levels of the interleukin-33 receptor soluble ST2 as a negative prognostic factor in hepatocellular carcinoma
.
Transl Oncol
2013
;
6
:
311
8
.
46.
Lu
DP
,
Zhou
XY
,
Yao
LT
,
Liu
CG
,
Ma
W
,
Jin
F
, et al
.
Serum soluble ST2 is associated with ER-positive breast cancer
.
BMC Cancer
2014
;
14
:
1
8
.
47.
Wang
L
,
Ellsworth
KA
,
Moon
I
,
Pelleymounter
LL
,
Eckloff
BW
,
Martin
YN
, et al
.
Functional genetic polymorphisms in the aromatase gene CYP19 vary the response of breast cancer patients to neoadjuvant therapy with aromatase inhibitors
.
Cancer Res
2010
;
70
:
319
28
.
48.
Yang
Y
,
Yang
Y
,
Yang
J
,
Zhao
X
,
Wei
X
.
Tumor microenvironment in ovarian cancer: function and therapeutic strategy
.
Front Cell Dev Biol
2020
;
11
:
758
.
49.
Yin
M
,
Shen
J
,
Yu
S
,
Fei
J
,
Zhu
X
,
Zhao
J
, et al
.
Tumor-associated macrophages (TAMs): a critical activator in ovarian cancer metastasis
.
Onco Targets Ther
2019
;
12
:
8687
.
50.
Xiong
S
,
Dong
L
,
Cheng
L
.
Neutrophils in cancer carcinogenesis and metastasis
.
J Hematol Oncol
2021
;
14
:
1
17
.
51.
Grivennikov
SI
,
Greten
FR
,
Karin
M
.
Immunity, inflammation, and cancer
.
Cell.
2010
;
140
:
883
99
.
52.
Gupta
V
,
Yull
F
,
Khabele
D
.
Bipolar tumor-associated macrophages in ovarian cancer as targets for therapy
.
Cancers
2018
;
10
:
366
.
53.
Juncker-Jensen
A
,
Hilliard
TS
,
Stavrou
N
,
Nagy
M
,
Parnell
E
,
Kuo
J
, et al
.
Abstract LB-C18: an integrated multiplexing approach for the immunoprofiling of the tumor microenvironment of ovarian granulosa cell tumors
.
Mol Cancer Ther
2019
;
18
:
LB
C18
.
54.
Lee
W
,
Ko
SY
,
Mohamed
MS
,
Kenny
HA
,
Lengyel
E
,
Naora
H
.
Neutrophils facilitate ovarian cancer premetastatic niche formation in the omentum
.
J Exp Med
2019
;
216
:
176
94
.
55.
Huang
QT
,
Zhou
L
,
Zeng
WJ
,
Ma
QQ
,
Wang
W
,
Zhong
M
, et al
.
Prognostic significance of neutrophil-to-lymphocyte ratio in ovarian cancer: a systematic review and meta-analysis of observational studies
.
Cell Physiol Biochem
2017
;
41
:
2411
8
.
56.
Wang
Z
,
Yang
Q
,
Tan
Y
,
Tang
Y
,
Ye
J
,
Yuan
B
, et al
.
Cancer-associated fibroblasts suppress cancer development: the other side of the coin
.
Front Cell Dev Biol
2021
;
9
:
146
.
57.
Schauer
IG
,
Sood
AK
,
Mok
S
,
Liu
J
.
Cancer-associated fibroblasts and their putative role in potentiating the initiation and development of epithelial ovarian cancer
.
Neoplasia
2011
;
13
:
393
405
.
58.
Ko
SY
,
Barengo
N
,
Ladanyi
A
,
Lee
JS
,
Marini
F
,
Lengyel
E
, et al
.
HOXA9 promotes ovarian cancer growth by stimulating cancer-associated fibroblasts
.
J Clin Invest
2012
;
122
:
3603
17
.
59.
Zhang
F
,
Cui
JY
,
Gao
HF
,
Yu
H
,
Gao
FF
,
Chen
JL
, et al
.
Cancer-associated fibroblasts induce epithelial-mesenchymal transition and cisplatin resistance in ovarian cancer via CXCL12/CXCR4 axis
.
Future Oncol
2020
;
16
:
2619
33
.
60.
Givel
AM
,
Kieffer
Y
,
Scholer-Dahirel
A
,
Sirven
P
,
Cardon
M
,
Pelon
F
, et al
.
miR200-regulated CXCL12β promotes fibroblast heterogeneity and immunosuppression in ovarian cancers
.
Nat Commun
2018
;
9
:
1
20
.
61.
Zhao
X
,
Ding
L
,
Lu
Z
,
Huang
X
,
Jing
Y
,
Yang
Y
, et al
.
Diminished CD68+ cancer-associated fibroblast subset induces regulatory T-cell (Treg) infiltration and predicts poor prognosis of oral squamous cell carcinoma patients
.
Am J Pathol
2020
;
190
:
886
99
.
62.
Özdemir
BC
,
Pentcheva-Hoang
T
,
Carstens
JL
,
Zheng
X
,
Wu
CC
,
Simpson
TR
, et al
.
Depletion of carcinoma-associated fibroblasts and fibrosis induces immunosuppression and accelerates pancreas cancer with reduced survival
.
Cancer Cell
2014
;
25
:
719
34
.
63.
Oshi
M
,
Huyser
MR
,
Le
L
,
Tokumaru
Y
,
Yan
L
,
Matsuyama
R
, et al
.
Abundance of microvascular endothelial cells is associated with response to chemotherapy and prognosis in colorectal cancer
.
Cancers
2021
;
13
:
1477
.
64.
Singhal
S
,
Stadanlick
J
,
Annunziata
MJ
,
Rao
AS
,
Bhojnagarwala
PS
,
O'Brien
S
, et al
.
Human tumor-associated monocytes/macrophages and their regulation of T-cell responses in early-stage lung cancer
.
Sci Transl Med
2019
;
11
:
eaat1500
.
65.
Masucci
MT
,
Minopoli
M
,
Carriero
MV
.
Tumor associated neutrophils. Their role in tumorigenesis, metastasis, prognosis and therapy
.
Front Oncol
2019
;
9
:
1146
.
This open access article is distributed under the Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International (CC BY-NC-ND 4.0) license.