Systematic Growth Factor Profiling Platform for 3D Tumor Models Reveals Estradiol-Responsive Cellular Mechanisms of Immunotherapy Resistance

Systematic Growth Factor Profiling Platform for 3D Tumor Models Reveals Estradiol-Responsive Cellular Mechanisms of Immunotherapy Resistance | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results Systematic Growth Factor Profiling Platform for 3D Tumor Models Reveals Estradiol-Responsive Cellular Mechanisms of Immunotherapy Resistance Kwanghwan Lee , Minsung Kim , Si On Lim , Dong-Ju Shin , Yun Shin , Jung-Joo Choi , Maria Lee , Hyun Ju Kang , Jeong-Won Lee , Jin-Ku Lee doi: https://doi.org/10.1101/2025.09.04.674374 Kwanghwan Lee 1 Genomic Medicine Institute, Medical Research Center, Seoul National University , Seoul 03080, Republic of Korea 2 Department of Anatomy and Cell Biology, Seoul National University College of Medicine , Seoul 03080, Republic of Korea Find this author on Google Scholar Find this author on PubMed Search for this author on this site Minsung Kim 3 Department of Biomedical Sciences, Seoul National University , Seoul 03080, Republic of Korea Find this author on Google Scholar Find this author on PubMed Search for this author on this site Si On Lim 3 Department of Biomedical Sciences, Seoul National University , Seoul 03080, Republic of Korea Find this author on Google Scholar Find this author on PubMed Search for this author on this site Dong-Ju Shin 3 Department of Biomedical Sciences, Seoul National University , Seoul 03080, Republic of Korea Find this author on Google Scholar Find this author on PubMed Search for this author on this site Yun Shin 3 Department of Biomedical Sciences, Seoul National University , Seoul 03080, Republic of Korea Find this author on Google Scholar Find this author on PubMed Search for this author on this site Jung-Joo Choi 4 Department of Obstetrics and Gynecology, Gynecologic Cancer Center, Samsung Medical Center, Sungkyunkwan University School of Medicine , Seoul 06351, Republic of Korea Find this author on Google Scholar Find this author on PubMed Search for this author on this site Maria Lee 5 Department of Obstetrics and Gynecology, Seoul National University College of Medicine, Seoul National University Hospital , Seoul 03080, Republic of Korea Find this author on Google Scholar Find this author on PubMed Search for this author on this site For correspondence: marialee{at}snu.ac.kr hj0601.kang{at}snu.ac.kr garden.lee{at}samsung.com jinkulee{at}snu.ac.kr Hyun Ju Kang 1 Genomic Medicine Institute, Medical Research Center, Seoul National University , Seoul 03080, Republic of Korea 2 Department of Anatomy and Cell Biology, Seoul National University College of Medicine , Seoul 03080, Republic of Korea Find this author on Google Scholar Find this author on PubMed Search for this author on this site For correspondence: marialee{at}snu.ac.kr hj0601.kang{at}snu.ac.kr garden.lee{at}samsung.com jinkulee{at}snu.ac.kr Jeong-Won Lee 4 Department of Obstetrics and Gynecology, Gynecologic Cancer Center, Samsung Medical Center, Sungkyunkwan University School of Medicine , Seoul 06351, Republic of Korea 6 Samsung Advanced Institute for Health Sciences and Technology, Sungkyunkwan University School of Medicine, Seoul 06351, Republic of Korea Find this author on Google Scholar Find this author on PubMed Search for this author on this site For correspondence: marialee{at}snu.ac.kr hj0601.kang{at}snu.ac.kr garden.lee{at}samsung.com jinkulee{at}snu.ac.kr Jin-Ku Lee 1 Genomic Medicine Institute, Medical Research Center, Seoul National University , Seoul 03080, Republic of Korea 2 Department of Anatomy and Cell Biology, Seoul National University College of Medicine , Seoul 03080, Republic of Korea 3 Department of Biomedical Sciences, Seoul National University , Seoul 03080, Republic of Korea Find this author on Google Scholar Find this author on PubMed Search for this author on this site For correspondence: marialee{at}snu.ac.kr hj0601.kang{at}snu.ac.kr garden.lee{at}samsung.com jinkulee{at}snu.ac.kr Abstract Full Text Info/History Metrics Preview PDF Abstract Current organoid culture systems face critical limitations: standardized growth factor formulations fail to capture patient-specific signaling requirements, and single-cell-type approaches overlook tumor-stromal interactions essential for understanding immunotherapy resistance. To address these challenges, we developed an innovative biofabrication platform that systematically integrates patient-derived three-dimensional (3D) cultures with comprehensive growth factor profiling across 128 combinations. Through systematic screening of 23 ovarian cancer patient samples and single-cell RNA sequencing, we identified two estradiol-responsive cellular populations that coordinate immunosuppression: a malignant cell fraction (MAL.PDCD5) and a cancer-associated fibroblast (CAF) fraction (FB.TNFSF10). MAL.PDCD5 cells suppress immune infiltration by downregulating antigen presentation, while FB.TNFSF10 cells promote immunosuppression through enhanced TGF-β signaling. Spatial transcriptomic analysis revealed striking mutual exclusivity between FB.TNFSF10 cells and T/NK cells, providing evidence of active immune cell exclusion. Most significantly, FB.TNFSF10 abundance emerged as a robust predictor of immune checkpoint inhibitor therapy resistance across multiple cancer cohorts, independent of conventional biomarkers. This biofabrication platform provides a scalable framework for drug screening and biomarker discovery, with immediate applications in precision medicine for patient stratification and combination therapy development. 1. Introduction The biofabrication of physiologically relevant tumor models represents a critical advancement in cancer research, enabling the recapitulation of complex cellular interactions that govern treatment responses [ 1 , 2 ]. While traditional 2D cell culture systems have provided valuable insights into cancer biology, they fail to capture the three-dimensional architecture and cellular heterogeneity characteristic of native tumor microenvironments (TMEs) [ 3 , 4 ]. Recent advances in biofabrication technologies, including organoid culture systems, bioprinting, and microfluidic platforms, have opened new possibilities for engineering tumor models that more accurately represent in vivo conditions [ 5 – 8 ]. However, existing 3D culture protocols face significant limitations. Most rely on empirically derived, standardized media formulations that fail to capture patient-specific growth factor requirements, using one-size-fits-all combinations without systematic optimization [ 9 , 10 ]. This approach potentially misses critical signaling dependencies that vary between patients and tumor subtypes [ 11 , 12 ]. Additionally, the complexity of the TME extends beyond cancer cells to include various stromal components, immune cells, and signaling molecules that collectively influence tumor progression and therapeutic responses [ 13 , 14 ]. Cancer-associated fibroblasts (CAFs) have emerged as critical orchestrators of TME dynamics [ 15 ], exhibiting remarkable plasticity and adopting distinct functional states that can either promote or suppress antitumor immunity [ 16 ]. Yet current biofabrication approaches have not fully captured the regulatory mechanisms governing CAF heterogeneity, particularly in hormone-responsive contexts [ 17 ]. To address these limitations, we developed an innovative biofabrication platform that systematically integrates patient-derived organoid cultures with comprehensive growth factor profiling across 128 combinations. This approach overcomes empirical media limitations and provides unprecedented insights into patient-specific signaling requirements. We pioneered the integration of single-cell RNA sequencing with 3D culture systems to investigate estradiol-mediated immune regulation in cancer, uniquely combining patient-derived cultures, systematic growth factor screening, and molecular profiling to identify specific cellular populations that drive immunosuppression and therapy resistance. This biofabrication strategy provides a scalable framework for modeling hormone-responsive tumor microenvironments and offers new opportunities for drug screening and biomarker discovery in cancer immunotherapy. 2. Materials and methods 2.1. Reagents and tools table Basal media composition and growth factors View this table: View inline View popup 2.2. Biofabrication of patient-derived tumor models Tissue Collection and Cell Culture After obtaining informed consent, tumor specimens and clinical records (Supplementary Table 1) were collected from patients undergoing surgery at Samsung Medical Center (SMC) and Seoul National University Hospital (SNUH) under approved Institutional Review Board protocols (IRB file numbers SMC 2019-07-163-019 and H-2011-127-1173). Tissue samples were snap-frozen in liquid nitrogen for genomic analysis or processed for organoid generation following established protocols [ 18 ]. Fresh tissue samples were mechanically dissociated and enzymatically digested with collagenase (Sigma-Aldrich) at 37°C under agitation for up to 30 min. When necessary, red blood cells (RBCs) were eliminated using RBC lysis buffer (Roche). Patient-derived CAFs [ 19 ] were established from the tissue dissociates of samples O010, O019, O032, and O044 by seeding 5 × 10 5 cells per well in the assay media (Supplementary Table 2). 2.3. Systemic growth factor (GF) profiling GFs (Supplementary Table 3) were dispensed into 384-well plates (761001; Nest) using an I.DOT Liquid Handler (Dispendix) with 40 µL of basal medium (M0, Supplementary Table 2) per well. Cells suspended in 50% Matrigel (Corning) were printed onto 384-pillar plates (Cellvitro, MBD) in a temperature-controlled chamber (4-8°C) to prevent premature Matrigel solidification. After 30 min of gelation at 37°C, the pillar plate was combined with the GF assay plate and cultured at 37°C with 5% CO 2 for 7 days. Cell viability was assessed using CellTiter-Glo 3D reagent (Promega) 36, and growth rates were calculated as the ratio of day 7 to day 0 luminescence signals. 2.4. Data analysis The cell proliferation measurements for the 128 GF combinations were performed in quadruplicate. Quality control was implemented using two criteria: (i) outlier removal using the 1-standard deviation method, excluding values outside (mean ± 1 standard deviation) of four-replicate data, and (ii) filtering of combinations with a coefficient of variation of ≥30% after outlier removal. Growth rates were calculated relative to day 0 measurements and converted to z-scores for normalization. Hierarchical clustering of samples based on growth scores was performed using the ComplexHeatMap R package. Detailed reagent information is provided in the Supplementary Table 11. 2.5. DNA/RNA isolation and next-generation sequencing analysis DNA and RNA were extracted from 17 tissue samples using an AllPrep DNA/RNA Mini Kit (Qiagen) or a DNeasy Blood & Tissue Kit and RNeasy Mini Kit (Qiagen), according to the manufacturer’s instructions. Whole-exome sequencing libraries were prepared using an Enzymatic Fragmentation Kit v1.0, Mechanical Fragmentation Kit (Twist Bioscience), or SureSelectXT Low Input Target Enrichment Kit (Agilent). RNA-seq libraries were prepared using the Illumina RNA Prep with Enrichment (Illumina) or TruSeq RNA Access (Illumina) library preparation kits. For RNA-seq library preparation for cell line experiments, a TruSeq Stranded mRNA Sample Prep Kit (Illumina) was used. 2.6. Variant calling from whole-exome sequencing Germline and somatic variant calling were conducted using the nf-core/Sarek 3.1.2 pipeline [ 20 ]. Briefly, raw sequencing data in FASTQ files were aligned to the human GRCh38 reference genome using BWA-MEM [ 21 ]. Following alignment, duplicate reads were marked using GATK Mark Duplicate Spark. The base quality scores were recalibrated using GATK BaseRecalibrator and GATK ApplyBQSR. Somatic variants, including single-nucleotide variants and small insertions/deletions, were identified using Mutect2 [ 22 ]. Germline variants were identified using Strelka2 [ 23 ]. Single-nucleotide variants and insertions/deletions were annotated using the Ensembl Variant Effect Predictor [ 24 ]. Finally, MAF files were generated using VCF2MAF and oncoplots were created using the R package maftools. 2.7. Bulk RNA-seq RiboDetector was used to remove rRNA sequences from the FASTQ files [ 25 ]. Subsequently, the STAR aligner [ 26 ] was used to map the RNA-seq reads to the human reference genome (GRCh38). Read counts were processed to obtain the expected counts using the RSEM package [ 27 ]. Batch effects arising from the two RNA-seq platforms were removed using ComBat-Seq in the R package sva [ 28 ]. Finally, the normalized expression values (TPM) were calculated from the expected batch effect removed counts. 2.8. GF responsiveness genomic interaction analysis The GF responsiveness for each GF–sample pair was assessed as the difference between the mean growth scores in the presence and absence of GF: responsiveness (ΔG) = mean growth score of the sample (with GF) − mean growth score of the sample (without GF). The significance of the interaction between GF responsiveness and mutation signatures was assessed using the Mann–Whitney U test, and that of the interaction between GF responsiveness and gene expression was measured using the Spearman correlation coefficient. 2.9. TCGA cohort analysis We used the TCGA-OV (Ovarian Serous Cystadenocarcinoma; n = 421) dataset for analysis. Gene expression data (STAR-Counts) were downloaded using the R package TCGAbiolinks and clinical data were obtained from cBioPortal. To predict the C1/C2 subtypes of TCGA patients with ovarian cancer, we performed k-means clustering (k = 2) using the Python package scikit-learn with 150 DEGs from our 18 samples as input. The gene expression values of the TCGA samples were normalized using variance-stabilizing transformation in the R package DESeq2. For the two clusters resulting from k-means clustering, we mapped the cluster enriched in ECM and developmental functions to C1, and the cluster enriched in immune functions to C2. Functional enrichment of DEGs in TCGA patients was performed using the gProfiler [ 29 ]. The epigenetically regulated RNA expression-based (EREG.EXPss) score was developed by Malta et al. [ 30 ]. was used as the stemness score for the TCGA samples. The C1 and C2 scores were calculated using GSVA for the C1 and C2 upregulated DEGs. 2.10. Analysis of GPER signaling in human ovarian cancer cell lines Four human ovarian cancer cell lines were used: A2780, A2780/CP20, PA-1, and TOG-21G cells (Supplementary Table 11). The cells were cultured in a cell culture medium (Supplementary Table 2). The authenticity of each cell line was verified using the PowerPlex Fusion 5C® System (Promega). The absence of mycoplasma was confirmed using the MycoAlert PLUS Mycoplasma Detection Kit (Lonza), following the manufacturer’s instructions. For GPER signaling analysis, one day after seeding, the cells were treated with the GPER inhibitor G15 (Selleckchem) at a concentration of 50 μM for 12 h. Detailed information on the reagents used is provided in the Supplementary Table 11. 2.11. Cell deconvolution analysis Three cell-type ESs (ImmunoScore, StromaScore, and MicroenvironmentScore) were calculated based on TPM values using the xCell web service ( http://xCell.ucsf.edu/ ) [ 31 ]. CIBERSORTx was used for the cell-type deconvolution of TCGA samples [ 32 ]. First, we created a signature matrix of cell subtypes in our scRNA-seq data using the ‘Create Signature Matrix’ mode. As an input, we used the normalized gene expression values of 200 down sampled cells for each subtype. Next, we calculated the cell fractions in each sample using the ‘Impute Cell Fractions’ mode with the TPM values of the TCGA samples. 2.12. Subtype analysis of TCGA patients The consensus molecular subtypes of TCGA patients were determined using the R package consensusOV, which uses a random forest classifier trained on consistently subtyped tumors from multiple datasets [ 33 ]. We used the TPM value of each gene as the input and set the other options as defaults. Three tumor immune phenotypes were obtained from Desbois et al. [ 34 ]. 2.13. Gene Set Enrichment Analysis (GSEA) in TCGA R package GSVA was used to calculate the ESs of the GSEA gene sets. The input was the raw count value of each gene, and the Poisson function was used as the kernel function (kcdf = Poisson) with the other options set as the default. The curated gene sets for the four reactome pathways were downloaded from MSigDB [ 35 ]: GPER1 signaling (R-HSA-9634597), ESR-mediated signaling (R-HSA-8939211), Wnt signaling (R-HSA-195721), and Wnt signaling in cancer (R-HSA-4791275). 2.14. Single-cell preparation and scRNA-seq Cryopreserved patient-derived organoids were thawed and cultured in a 50% Matrigel dome in a medium with or without estradiol for 3 days. Sequencing data were processed using CellRanger and demultiplexed based on major voting in souporcell [ 36 ], vireoSNP [ 37 ], and demuxlet [ 38 ]. The resulting barcode matrix of 24,508 cells was processed using Seurat 5.0.0 [ 39 ], which is a toolkit for scRNA-seq data analysis. Low-quality cells were excluded based on the following criteria: (i) the number of expressed genes < 250; (ii) number of expressed genes ≥ 6000; (iii) RNA count ≤ 500; (iv) mitochondrial transcripts (indicative of apoptosis) ≥ 20%; and (v) log10 (genes per unique molecular identifier) ≤ 0.8. Unique molecular identifier count data were normalized using log transformation. In total, 20 and 309 single cells were retained, respectively (Supplementary Table 7). The top 2000 highly variable genes were selected to aggregate samples into a merged dataset. The merged cell-by-gene matrix was processed using SCTransform (v2) [ 40 ] in the Seurat package, which included functions for the normalization, regression, and identification of variable genes. Principal component analysis was applied to the processed matrix. Next, we applied harmony integration [ 41 ], which effectively reduced technical batch effects while preserving biological variation. The RunHarmony function returns a Seurat object updated with the corrected harmony coordinates. Subsequently, the main cell clusters were identified using Seurat’s FindClusters function with the following options: reduction = harmony and resolution = 0.35, and visualized using Uniform Manifold Approximation and Projection. 2.15. PRIMUS cell subtype analysis We applied PRIMUS [ 42 ] to identify the phenotypic cell subtypes of malignant cells, fibroblasts, and macrophages, while accounting for patient-specific components and technical noise. To determine the optimal number of clusters, we ran PRIMUS for k values ranging from 1 to 15, each with five random initial parameter sets. Based on the Bayesian information criterion, we selected k = 3, 6, and 2 for the malignant cells, fibroblasts, and macrophages, respectively. Fibroblasts expressing both epithelial and fibroblast markers were excluded from this study. Representative gene sets were identified using the methods reported by Zhang et al. [ 43 ]. through the following steps. (i) Marker genes for each cluster were identified using differential expression analysis by comparing query clusters with other clusters. Genes with a log2(fold change) ≥ 0.05 and FDR < 0.05 were selected from each comparison. (ii) Using CS-core [ 44 ], a co-expression matrix was calculated to identify significantly co-expressed genes (FDR < 0.05). (iii) From the co-expression matrix, co-expression networks were extracted using the WGCNA algorithm [ 45 ], and the component genes of each co-expression network were used as representative gene sets. (iv) Finally, a gene set overrepresentation analysis was conducted for the remaining genes in each representative gene set. 2.16. Single-cell GSEA The ES gene set for each cell was calculated using the R package escape. C1 and C2 scores were determined as the percentile rank difference between the ESs of the C1- and C2-upregulated DEGs. The gene set ‘Oxidative Phosphorylation’ was obtained from the cancer hallmark gene sets. The gene set ‘Antigen Presentation’ and TGF-β CAF markers were sourced from Hornburg et al. [ 46 ]. The CAF.S1 markers were downloaded from Kieffer et al. [ 47 ]. 2.17. Cell-cell interaction analysis To identify and visualize cell-cell interactions among cell subtypes and identify major signaling differences between estradiol-treated and untreated cells, the R package CellChat was used according to the developer’s vignette ( https://github.com/sqjin/CellChat ) [ 48 ]. We split the integrated scRNA-seq data into estradiol-treated and untreated groups and ran CellChat independently. Next, we conducted a comparative analysis of the two CellChat datasets shown in the vignette ‘Comparison analysis of multiple datasets using CellChat.’ 2.18. Survival analysis of TCGA-OV patients Survival information for TCGA patients was obtained from the ‘Curated survival data’ file in the TCGA Ovarian Cancer dataset of UCSC Xena ( https://xenabrowser.net/datapages/ ), as reported by Liu et al. [ 49 ]. Survival and Survminer R packages were used to create Kaplan– Meier curves, and the significance of survival differences was determined using the log-rank test. 2.19. Quantitative analysis of FB.TNFSF10 fractions and pathway enrichment To quantify the FB.TNFSF10 cell fractions, bulk RNA-seq data were deconvoluted using CIBERSORTx, with a cell signature matrix derived from our scRNA-seq data. The enrichment of FB.TNFSF10 cell markers was assessed using ssGSEA with the GSVA package. Three marker gene sets were used to ensure robust ssGSEA results. The first set comprised markers identified by differential expression analysis, using the Wilcoxon rank-sum test in the Seurat package. The second set included markers derived from co-expression network analysis, as shown in Supplementary Figure 13B . Third set of integrated markers from both approaches. included markers derived from the co-expression network analysis. Third set of integrated markers from both approaches. 2.20. Pancancer analysis of TCGA patients Transcriptomic data from 11,096 samples representing 33 cancer types were obtained from the TCGA PancanAtlas database ( https://gdc.cancer.gov/about-data/publications/pancanatlas ). These batch-corrected gene expression data were used to perform GSEA of TGFβ-associated Cancer-Associated Fibroblast (TGFβ.CAF) and CAF.S1 marker signatures. In addition, the pre-computed xCell scores for the TCGA pancancer samples were retrieved from the xCell web portal ( http://xCell.ucsf.edu/ ). 2.21. Xenium data analysis We obtained preprocessed 5k Xenium data of an ovarian cancer sample from the 10x Genomics database ( https://www.10xgenomics.com/datasets/xenium-prime-ffpe-human-ovarian-cancer ). For the fibroblast subtype prediction, cells categorized as “tumor-associated” and “stromal-associated fibroblasts” were selected. We used the DeepCC tool [ 50 ] and trained a deep-learning model using the fibroblast transcriptome from the scRNA-seq dataset. To analyze the correlation between fibroblast subtypes and T/NK cells, we divided the Xenium image into a 7×8 grid of tiles. Fibroblast fractions were calculated from the total fibroblast count within each tile, whereas T/NK cell fractions were calculated from the total cell count per tile to prevent bias from variation in fibroblast abundance. The ssGSEA analysis for the TGFβ CAF gene set was performed using the ’escape R package’ with parameters ’method = UCell’ and ’maxRank = 5000’. For visualization of Xenium cells, we used the scanpy and squidpy Python packages. 2.22. ICI-treated patients’ data collection and preprocessing We analyzed bulk RNA sequencing (RNA-seq) data from tumor samples across 16 cohorts with ICI treatment: Kim et al. [ 51 ]; Huang et al. [ 52 ]; VanAllen et al. [ 53 ]; Riaz et al. [ 54 ]; Chen et al. [ 55 ]; Liu et al. [ 56 ]; Hugo et al. [ 57 ]; Gide et al. [ 58 ]; Prat et al. [ 59 ]; Lee et al. [ 60 ]; Cho et al. [ 61 ]; Jung et al. [ 62 ]; Hwang et al. [ 63 ]; Miao et al. [ 64 ]; Phillips et al. [ 65 ]; and IMvigor210 [ 66 ]. For eight datasets (Gide, Hugo, Jung, Kim, Phillips, Lee, Riaz, and Cho), raw fastq files were obtained from patient samples and aligned to the GRCh38 human reference genome using STAR (version 2.7.9a). Read quantification was performed using the RSEM pipeline (version 1.3.3). For the IMvigor210 dataset, raw read counts were directly downloaded. Gene expression levels were calculated by normalizing the counts using the trimmed means of M-value (TMM) normalization implemented in the edgeR R package (version 3.32.1) [ 67 , 68 ]. For the remaining datasets (Huang, VanAllen, Liu, Chen, Hwang, and Miao), we used processed gene expression levels provided by Lee et al. Processed expression data from the Prat dataset were retrieved from the Gene Expression Omnibus (GEO) using the GEOquery R package (version 2.58.0). For drug response classification, patients were categorized based on the RECIST criteria: those with a partial or complete response were classified as responders (R), whereas those with stable or progressive disease were classified as non-responders (NR). The Van Allen cohort was an exception, in which the response classification followed the criteria established by Lee et al . 3. Results 3.1. Integrative genomics and growth factor-omics analysis Cancer exhibits significant molecular heterogeneity among patients, resulting in diverse GF responses that can affect treatment outcomes. Understanding these different responses is crucial for developing personalized therapeutic strategies. However, systematic approaches to characterize patient-specific GF responsiveness are limited. To address this gap, we designed a comprehensive study to identify the GFs that influence cancer growth and their associated genomic patterns. Our approach uniquely combines growth factor screening with genomic profiling to establish correlations between cancer cell responsiveness to GF combinations and underlying molecular features. This integrative strategy aimed to uncover previously unrecognized patterns of GF responsiveness and their molecular determinants ( Figure 1A ). Download figure Open in new tab Figure 1. Schematic overview and growth factor-omics analysis A. The study aims, samples, experimental scheme for growth factor-omic data acquisition, and genomics analysis are described. B. Ranking of cell growth results for 128 combinations of seven GFs. The bar plot shows the mean cell growth compared that with on day 0. The heatmap shows the combinations of the seven GFs, with gray indicating the absence of the corresponding GF. We collected tissue samples from 31 patients (EOC) who underwent comprehensive growth factor-associated molecular profiling using multiple analytical approaches. The experimental design included parallel analyses of bulk RNA sequencing of 21 frozen specimens, and tissue dissociation of 27 samples. From the dissociated tissues, we successfully established GF profiles for 18 specimens and isolated CAFs from four samples (O010, O019, O032, and O044; Supplementary Table 1). To investigate the effect of hormonal signaling, which is increasingly recognized as a critical modulator of tumor behavior, we performed single-cell RNA sequencing (scRNA-seq) on four dissociated tissue samples under estradiol-positive and estradiol-negative conditions. This approach allowed us to examine the cellular heterogeneity and hormone-responsive populations at an unprecedented resolution. Analysis of the successfully processed samples (n = 23) revealed a predominance of high-grade serous carcinomas (n = 19), whereas the remaining samples comprised single cases of clear cell, mucinous, endometrioid, and low-grade serous subtypes (Supplementary Table 1). For GF screening, we developed a systematic approach using a supplemented basal medium (Supplementary Table 2) containing seven key GFs in 128 combinations (Supplementary Table 3). This comprehensive screening strategy enabled the identification of complex GF interactions that may have been missed in simpler experimental designs. Analysis of EOC-derived cell growth responses revealed a median growth rate of 5.93, with a 3.45-fold difference between the highest and lowest rates, highlighting the substantial variability in GF responsiveness among patients ( Figures 1B and 2A , Supplementary Table 4). Download figure Open in new tab Figure 2. Differential cell growth induced by growth factors (GFs). A. Growth of epithelial ovarian cancer (EOC)-derived cells obtained from 18 patients after 7 days of incubation in medium containing growth factors (GFs) in 128 combinations. Box plots span the first to third quartiles and whiskers represent the 1.5 interquartile range. The central line indicates the median value. B. Cell growth efficiency of the top three growth factor-omics conditions (WREsNo (top), WEsEF (middle), and WREs (bottom)) compared to day 0 and basal medium (gray, M0). W, Wnt; R, R-spondin; Es, oestradiol; No, noggin; E, EGF; F, FGF10. C. Representative images of tumor growth under WRE sNo conditions from three independent experiments. Top: Day 0 (D0). Middle: basal medium (M0) on day 7. Bottom: WREsNo on Day 7. Scale bars, 200 μm. D. Effects of seven GFs on EOC cell growth. Cell growth was normalized for each sample using the z-score. Box plots span the first to third quartiles and whiskers represent the 1.5 interquartile range. The central line indicates the median value. The significance of the difference between the GF-positive (red line) and GF-negative (gray line) conditions was assessed using the Mann–Whitney U test. We identified specific GF combinations that consistently promoted robust cell growth across samples. Combinations containing Wnt + R-spondin + estradiol + Noggin (WREsNo), Wnt + estradiol + EGF + FGF10 (WEsEF), and Wnt + R-spondin + estradiol (WREs) induced the highest cell growth rates ( Figures 2A-C ). Statistical analysis revealed that WREs significantly increased growth rate ( P < 0.001), whereas N-acetyl-L-cysteine (NAC) reduced growth ( P = 0.029). FGF10, EGF, and Noggin showed no significant independent impact on growth, suggesting complex interactions between GFs rather than simple additive effects ( Figure 2D ). These findings provide insights into the patient-specific patterns of GF responsiveness and establish a foundation for understanding how these patterns relate to underlying genomic features and potential therapeutic responses. 3.2. Molecular associations with GF responsiveness Understanding the molecular basis of differential GF responses is essential for predicting treatment outcomes and for developing personalized therapeutic strategies. Although previous studies have identified individual genetic markers associated with treatment responses, a comprehensive analysis of how molecular features are related to GF responsiveness is lacking. To address these knowledge gaps, we performed hierarchical clustering analysis of standardized growth rates across our patient cohort ( Figure 3A ) [ 69 ]. This unbiased approach revealed two clusters, C1 and C2, suggesting differences in GF response patterns. Detailed analysis showed that C1 samples exhibited significantly higher growth rates in response to estradiol and Wnt signaling, whereas C2 samples showed higher responsiveness to R-spondin ( Figure 3B ). Other GFs, including NAC, FGF10, EGF, and Noggin, exhibited no significant association with either cluster, indicating specificity in the molecular determinants of the GF response ( Supplementary Figure 1A ). Download figure Open in new tab Figure 3. Hierarchical clustering of GF responsiveness reveals two major clusters. A. Hierarchical clustering of normalized cell growth in 18 ovarian cancer samples revealed two clusters: C1 and C2. B. Comparison of the growth scores of samples in clusters C1 (red line) and C2 (green line) cultured in the presence or absence of the indicated GFs. Box plots span from the first to third quartiles, and the whiskers represent the 1.5 interquartile range. The central line indicates the median value. P -values were calculated using the Mann–Whitney U test. C. Correlation of growth scores between the indicated GFs. Correlations were analyzed using PCC. D. Responsiveness of four HGSOC molecular subtypes (DIF, PRO, IMR, and MES) to the indicated GFs. The heatmap indicates the amount of change in average normalized cell growth. Significance was assessed using the Mann–Whitney U test. E, F, G. Bland–Altman (MA) plot of the correlation between gene expression and responsiveness to estradiol (E), Wnt (F), and R-spondin (G). The x-axis represents the mean gene expression across the samples and the y-axis shows the Spearman correlation coefficient (SCC). The dot size indicates the significance of the SCC. H, I, J. GSEA results of positive (red) and negative (blue) SCC between gene expression and responsiveness to estradiol (H), Wnt (I), and R-spondin (J). Gene Ontology: Biological Process gene sets were used in GSEA. Further investigation revealed a significant inverse correlation between R-spondin and estradiol responsiveness (Pearson correlation coefficient = –0.75, P = 1.9 × 10⁻³; Figure 3C and Supplementary Figure 1B ). This unexpected finding suggests that these pathways may have antagonistic effects, providing new insights into potential combination treatment strategies. The strong correlation also indicated these factors were the main determinants distinguishing C1 and C2 clusters. To understand the molecular basis of these differences, we performed a differential gene expression (DEG) analysis, which identified 50 representative genes, with 25 genes specific to each cluster ( Supplementary Figure 2 and Supplementary Table 5). This analysis revealed distinct molecular signatures associated with GF responsiveness. Estradiol responsiveness was significantly higher in the proliferative (PRO) and mesenchymal (MES) subtypes, whereas the differentiated (DIF) subtype showed lower estradiol and higher R-spondin responsiveness. Immunoreactive (IMR) and MES subtypes are predominantly R-spondin responsive ( Figure 3D and Supplementary Figure 3 ) [ 70 , 71 ]. Elevated estradiol responsiveness was associated with upregulation of HDCA2 and RPL23 , downregulation of TMSB4X and HSPB1 , and enrichment of cytoplasmic translation and cell adhesion gene sets, whereas immune-related gene sets were depleted ( Figures 3E and 3H ). Enhanced Wnt responsiveness correlated with increased FAT1 and TPP1 expression, decreased HSPD1 and CD9 expression, and the enrichment of angiogenesis-related gene sets ( Figures 3F and 3I ). The heightened R-spondin response resulted in increased ARF5 and UBE2H expression, reduced KRT8 and COX5B expression, and enrichment of axoneme assembly and lymphocyte chemotaxis genes ( Figures 3G and 3J ). Integrated analysis of the mutation profiles and GF responsiveness revealed significant correlations between specific genetic alterations and GF sensitivity (Supplementary Figures 4 and 5). 3.3. Characteristics of C1/C2-type samples in TCGA cohort Although our initial analysis identified distinct GF response patterns and associated molecular signatures in our patient cohort, validation using a larger independent dataset is crucial to establish the broader relevance of these findings. To validate the molecular characteristics of the C1/C2 clusters and assess their associated clinical features, we classified the samples from the TCGA-OV cohort (n = 421) into C1 or C2 using k-means clustering based on expression patterns. We identified 193 C1-type and 228 C2-type patients with no significant cluster-specific mutation patterns ( Supplementary Figure 6 ). Analyses of demographic and clinical characteristics revealed several significant associations. C1-type patients exhibited significantly higher age at diagnosis and aneuploidy scores than C2-type patients ( P = 5.4 × 10⁻¹⁰ and P = 1.8 × 10⁻¹², respectively). However, tumor mutation burden, ESR1 expression, winter hypoxia score, MSI score, and mutation counts did not differ significantly between the groups ( Supplementary Figure 7A ). Overall and progression-free survival rates were similar ( Supplementary Figure 7B ). DEG analysis revealed that 2,604 and 486 genes were upregulated in the C1 and C2 treatments, respectively ( Figure 4A and Supplementary Table 6). C1 upregulated genes were enriched in ’collagen-containing ECM’ (FDR = 1.22 × 10⁻⁵), whereas C2 genes were enriched in ’immune response’ (FDR = 8.59 × 10⁻¹⁹, Figure 4B ). The stemness score was significantly higher in C1 ( P = 0.016; Figure 4C ). C1 showed a higher StromaScore ( P = 0.0013), whereas C2 exhibited higher ImmunoScore and MicroenvironmentScore ( P = 7.4 × 10⁻⁹ and P = 4.9 × 10⁻⁷, respectively; Figure 4D ). Download figure Open in new tab Figure 4. C1/C2 signature characteristics of the TCGA-OV cohort. A. Top 50 DEGs between the C1- and C2-type samples in TCGA-OV cohort. The heatmap shows the z-score-normalized gene expression of the top 25 upregulated and downregulated genes. B. GO enrichment results of DEGs between the C1-(red) and C2-type (green) samples in the TCGA-OV cohort. C. The stemness scores of the C1-type (red line) and C2-type (green line) samples were calculated using the EREG.EXPss scores. Box plots span from the first to third quartiles, and the whiskers represent the 1.5 interquartile range. The central line indicates the median value. Statistical significance was calculated using the Mann–Whitney U test. D. Cell-type ESs of C1-type (red line) and C2-type (green line) samples were calculated using xCell: ImmuneScore (composite score of immune cell types), StromaScore (composite score of stromal cell types), and MicroenvironmentScore (composite scores of immune and stromal cell types). Box plots span from the first to third quartiles, and the whiskers represent the 1.5 interquartile range. P -values were calculated using the Mann–Whitney U test. E. Enrichment results of C1-type (red line) and C2-type (green line) samples for four Reactome pathways: GPER1 signaling (R-HSA-9634597), ESR-mediated signaling (R-HSA-8939211), WNT signaling (R-HSA-195721), and WNT signaling in cancer (R-HSA-4791275). Box plots span from the first to third quartiles, and the whiskers represent the 1.5 interquartile range. The central line indicates the median value. P -values were calculated using the Mann–Whitney U test. F. Stacked bar plot displaying the distribution of HGSOC consensus subtypes in C1- and C2-type samples. Statistical significance was calculated using the chi-square test. G. ESs of immune phenotypes (tumor immune phenotype [TIP]: desert, excluded, or infiltrated) of TCGA patients in the C1-type (left) and C2-type (right) samples. Box plots span from the first to third quartiles, and the whiskers represent the 1.5 interquartile range. The central line indicates the median value. P -values were calculated using the Mann–Whitney U test. Gene set variation analysis revealed significantly higher enrichment scores for G-protein-coupled estrogen receptor (GPER) signaling and Wnt-related signaling in C1 ( P = 0.04, P = 0.0057, and P = 7.1 × 10⁻⁵, respectively; Figure 4E ). ESR-mediated signaling showed no significant difference ( P = 0.98), suggesting that GPER may be a key mediator of estradiol responsiveness [ 72 – 74 ]. C2 samples showed higher fractions of differentiated (DIF) and immune-reactive (IMR) subtypes, whereas C1 had higher fractions of mesenchymal (MES) and proliferative (PRO) subtypes ( Figure 4F and Supplementary Figure 7C ). C1 samples exhibited significantly higher deserts and excluded immune phenotypes, whereas C2 samples showed a higher infiltration phenotype ( Figure 4G ) [ 75 ]. These findings validate our initial observations and provide new insights into the biological and clinical significance of the C1/C2 classification. The distinct molecular and cellular characteristics associated with each type of cancer suggest potential therapeutic implications, particularly regarding immunotherapy responses and targeted intervention strategies. 3.4. GPER1 signaling activity correlates with MES and DIF subtypes Understanding the molecular mechanisms that determine ovarian cancer subtypes is crucial for developing effective targeted therapies. Previous studies have demonstrated the importance of estrogen signaling in ovarian cancer progression; [ 76 ] however, the specific pathways and receptors that mediate subtype-specific effects remain poorly understood. We investigated the factors determining the four consensus subtypes [ 33 ], focusing on estrogen responsiveness, which is crucial for the C1/C2 classification. Analysis of estrogen receptor (ESR) activity, which mediates estrogen responsiveness, did not refine subtype classification based on the C1/C2 type nor did it show differences in molecular subtypes ( Supplementary Figure 8 ). We then analyzed the activity of GPER, which is highly expressed in the C1 cluster. GPER-high and low subtypes were defined based on GPER1 signaling enrichment scores (ES) above or below the threshold of -0.09890369, respectively. C1-GPER-high samples were enriched in the mesenchymal (MES) subtype, C1-GPER-low samples were enriched in proliferative (PRO), C2-GPER-high samples were enriched in immunoreactive (IMR), and C2-GPER-low samples were enriched in differentiated (DIF) subtype (Supplementary Figures 9A and 4B). Similar patterns were observed in the GSE51088 and GSE26193 ovarian cancer cohorts ( Supplementary Figure 10 ). GPER activity significantly correlated with molecular subtypes: low GPER activity associated with higher DIF scores, whereas high GPER activity associated with higher MES scores ( P = 1.5 × 10⁻⁹ and P = 4.9 × 10⁻¹⁰, respectively; Supplementary Figure 9C ). Functional enrichment analysis revealed that samples with high GPER activity had higher scores for RNA processing and stem cell-related genes (e.g., FOXM1 and Oct4), whereas samples with low GPER activity had higher scores for extracellular matrix organization and immune system processes ( Supplementary Figure 9D ). Experimental validation using four ovarian cancer cell lines revealed no significant difference in GPER1 expression. However, PA-1 cells exhibited markedly high GPER1 signaling activity ( Supplementary Figure 9E ). PA-1 cells with high GPER1 activity had a lower DIF subtype score but a significantly higher MES score than other cell lines ( Supplementary Figure 9F ). The inhibition of GPER1 by G15 in PA-1 cells decreased GPER1 activity, increased DIF score, and decreased MES score ( Supplementary Figure 9G ). These findings suggest that GPER1 activation may drive the MES subtype [ 74 ]. These results reveal a previously unrecognized role of GPER1 signaling in determining the molecular subtypes of ovarian cancer, particularly in driving the mesenchymal phenotype. The identification of GPER1 as a key regulator of subtype switching has important therapeutic implications, suggesting that targeting GPER1 signaling may provide a novel strategy for modulating the tumor phenotype. 3.5. ScRNA-seq defines the estradiol-responsive cell fraction Although estrogen has been implicated in cancer immune evasion [ 77 ], the specific cellular mediators of estrogen-induced immunosuppression remain poorly defined. Previous studies have largely focused on bulk tissue responses or isolated cell populations, making it difficult to identify the key cellular fractions that directly mediate the immunosuppressive effects of estrogen. To address these knowledge gaps and identify specific estradiol-responsive cell populations, we performed scRNA-seq on four dissociated samples with or without estradiol ( Figure 5A and Supplementary Table 7). Among the 18,399 filtered cells, we identified six major cell types (Supplementary Figures 11A-B, and Supplementary Table 8). Each cluster expressed significantly higher levels of cell-type specific markers ( Supplementary Figure 12 ). Download figure Open in new tab Figure 5. Estradiol-responsive cell fractions contribute to the immunosuppressive TME A. A schematic representation of the scRNA-seq-based detection of estradiol-responsible cell fractions B. UMAP plot of malignant cells and fibroblast subtypes (left). C1-C2 scores based on percentile rank difference between C1 and C2 upregulated DEG ESs in subtypes (center). Estrogen response as the fold difference between estradiol-treated and untreated samples (right). C, D. Antigen presentation (C) and oxidative phosphorylation (D) ESs across malignant cell subtypes. Box plots show the quartiles and 1.5 IQR whiskers. The median is marked by a central line. Significance was determined using the Mann-Whitney U test. E, F. CAF.S1 (E) and TGF-β-CAF (F) ESs across fibroblast subtypes. Box plots show the quartiles and 1.5 IQR whiskers. The median is marked by a central line. P values from the Mann-Whitney U test. G, H. MAL.PDCD5 (G) and FB.TNFSF10 (H) cell fractions across TCGA tumor immune phenotypes were calculated using CIBERSORTx. Box plots show the quartiles and 1.5 IQR whiskers. The median is marked by a central line. Significance was assessed using the Mann-Whitney U test. I. Heatmap of the differential interaction strength between malignant cells and fibroblasts in estradiol-treated vs. untreated groups. The top and right bar plots show the total incoming and outgoing signals, respectively. Red indicates increased signaling and blue indicates decreased signaling in the treated group. J. Kaplan-Meier analysis of overall survival in FB.TNFSF10-high and -low TCGA patient groups. Subgroups were determined using CIBERSORTx, with mean FB.TNFSF10 expression as a threshold. P-values from the log-rank test. K. Box plots show FB.TNFSF10 cell fractions in the indicated patient-derived CAFs under 100 nM estradiol-treated and untreated conditions for 14 days. Box plots show the quartiles and 1.5 IQR whiskers. The median is marked by a central line. L. Heatmaps displaying ssGSEA scores for FB.TNFSF10 marker sets (All, Gene, and Network markers) in the indicated CAFs under estradiol-treated and untreated conditions. Using PRIMUS [ 42 ], we identified three malignant (MAL), five fibroblast (FB), and two macrophage subtypes ( Figure 5B , Supplementary Figures 11C-E and 13; Supplementary Tables 9 and 10) [ 42 ]. All FB populations exhibited CAF-positive patterns, with FB.KRT19 showing a relatively strong antigen-presenting CAF signature ( Supplementary Figure 14 ) [ 78 , 79 ]. The MAL.PDCD5, FB.COL4A2, FB.TNFSF10, and FB.KRT19 subtypes had higher C1 than C2 scores, and MAL.PDCD5, FB.MMP3, FB.TNFSF10, and FB.KRT19 clusters showed quantitative increases upon estradiol treatment ( Figure 5B and Supplementary Figure 15 ). MAL.PDCD5 had a significantly lower enrichment for antigen presentation and oxidative phosphorylation functions, suggesting immunosuppressive properties ( Figures 5C and 5D ) [ 80 ]. Of estradiol-responsive CAFs, FB.TNFSF10 had higher scores for immunosuppressive CAF types CAF.S1 and TGF-β CAF, whereas FB.KRT19 showed lower scores ( Figures 5E and 5F ) [ 47 , 81 ] In TCGA samples, the desert immune phenotype [ 46 ] had a significantly larger MAL.PDCD5 fraction, whereas the excluded type had a larger FB.TNFSF10 fraction ( Figures 5G and 5H ). Cell-cell interaction analysis using CellChat[ 48 ] revealed that MAL.PDCD5 and FB.TNFSF10 showed the largest increase in interaction strength after estradiol treatment ( Figure 5I and Supplementary Figure 16 ). The FB.TNFSF10-high group showed a significantly shorter overall survival ( P = 0.016; Figure 5J ). Given the significant association between the FB.TNFSF10 fraction and patient outcomes, coupled with its distinctive tumor immune escape phenotype, we performed in-depth characterization of this cellular subset. We sought to validate whether the FB.TNFSF10 signature identified by scRNA sequencing can be translated for bulk transcriptome analysis. CIBERSORTx [ 32 ] analysis confirmed substantial CAF fractions (> 0.6) across the four established CAF cell lines, with consistent proportions maintained regardless of estradiol treatment ( Supplementary Figure 17 ). Three out of the four CAF cells exhibited increased FB.TNFSF10 fractions following estradiol treatment, with O010 cells showing the most pronounced response (0.12 to 0.20, representing a ∼1.7-fold increase; Figure 5K ). Single-sample Gene Set Enrichment Analysis (ssGSEA) of the FB.TNFSF10 gene set revealed a remarkable upregulation in response to estradiol treatment, which was consistently observed across both gene and network markers ( Figure 5L ). These findings provide the first detailed cellular map of estradiol responsiveness in ovarian cancer, revealing specific cell populations that mediate hormone-dependent changes in the tumor microenvironment. 3.6. Estradiol-responsive cell fractions correlate with pancancer immunosuppression Although estrogen signaling has traditionally been studied in hormone-dependent cancers, emerging evidence suggests a broader role in modulating immune responses across various cancer types [ 17 , 77 ]. However, the extent to which estrogen-responsive cell populations influence the immune landscape beyond hormone-dependent cancers remains unclear. To elucidate the relationship between estradiol and the TME, we conducted a cell deconvolution analysis of 11,069 TCGA transcriptome samples. The normalized proportions of FB.TNFSF10 and MAL.PDCD5 cells varied across the 33 TCGA cancer types ( Supplementary Figure 18 ). We found significant negative correlations between MAL.PDCD5 and various immune cell populations, including dendritic cells, CD8+ T cells, and B cells, across the pancancer cohort ( Figure 6A , left). FB.TNFSF10 exhibited positive correlations with TGFβ-CAF and CAF.S1 fractions, consistent with patterns observed in TCGA ovarian cancer cohort ( Figure 6A , right) [ 82 ]. Download figure Open in new tab Figure 6. Pan-cancer immune phenotypes associated with estradiol-responsive fractions A. A heat map of correlations between immune phenotypes and FB.TNFSF10 and MAL.PDCD5 expression across 33 TCGA cancer types. B, C. Spearman correlation analysis of PDL1 expression with FB.TNFSF10 (B) and MAL.PDCD5 (C) cell abundance in TCGA pan-cancer dataset. D. Distribution of fibroblast (FB) subtypes across the TCGA samples. FB subtype classification was based on the predominant fibroblast population identified by CIBERSORTx analysis for each patient. E, F. Enrichment scores (ESs) for TGF-β CAF (E) and CAF.S1 (F) signatures stratified by FB subtype. G, H. X-cell-derived regulatory T cell (Treg) scores (G) and stromal scores (H) across patient FB subtypes. Box plots show the median (center line), first to third quartiles (box), and 1.5× interquartile range (whiskers). Statistical significance was assessed using the Mann-Whitney U test. I. Predicted fibroblast and T/NK cell subtypes in xenium ovarian cancer samples (top). The Xenium data were divided into 56 tiles (7×8 grids, labeled S1-1 through S7-8). Pearson correlation coefficients (PCC) between normalized cell fractions of fibroblast subtypes and T/NK cells across tiles are shown (bottom). Asterisks indicate a significant correlation ( P < 0.05). J. Distribution of four fibroblast subtypes and T/NK cells in three representative tiles (S2-7, S5-2, and S7-2), where FB.TNFSF10 cells numbered < 100. Gray dots represent fibroblast subtypes and red dots indicate T/NK cells. Analysis of PD-L1 expression [ 83 ], revealed a significant inverse correlation with MAL.PDCD5 fraction (Spearman R = -0.331, P = 2.98 × 10⁻²⁸¹, Figure 6b ). Conversely, FB.TNFSF10 demonstrated a positive correlation with PD-L1 expression (Spearman R = 0.194, P = 4.33 × 10⁻⁹⁴, Figure 6C ). These findings prompted us to focus on subsequent ICI treatment analyses of FB.TNFSF10. Next, we classified patients into five categories based on their predominant fibroblast (FB) subtype, with the distribution of these FB types varying across the 33 TCGA cancer types ( Figure 6D ). Consistent with our previous observations, FB.TNFSF10-type patients exhibited the highest ES for TGFβ-CAF and CAF.S1 markers among the five FB subtypes (median ES = 0.257 and 0.312, respectively; Figs. 6E and 6F ) [ 47 , 81 ]. Furthermore, this subtype ranked second and third in xCell scores for Tregs and stromal content, respectively (median ES = 0.0149 and 0.110, respectively; Figures 6G and 6H ). To investigate the spatial relationships between FB subtypes, we analyzed 5k Xenium data from an ovarian cancer sample ( Figure 6I ). We quantified the cell-type fractions across a 7×8 grid of tiles ( Supplementary Figure 19A ) and discovered that FB.TNFSF10 cells showed mutual exclusivity with T and NK cells ( Figure 6J ). Across all 56 tiles, the normalized fraction of FB.TNFSF10 cells displayed a significant negative correlation with T/NK cell fractions (PCC r= -0.317, P = 0.0185; Figure 6I and Supplementary Figure 19B ). Consistent with previous results, FB.TNFSF10 cells exhibited elevated TGFβ CAF signature scores in the Xenium data (median ES = 0.520; Supplementary Figures 19C-D). These findings provide compelling evidence that estrogen-mediated signaling shapes the immunosuppressive landscape across diverse cancer types through conserved cellular mechanisms. Furthermore, high-resolution spatial transcriptomic analysis revealed consistent patterns of immune cell exclusion and stromal modulation associated with FB.TNFSF10 cell population. This spatial validation provided crucial contextual evidence for the functional role of FB.TNFSF10 cells within the complex tissue architecture, highlighting their potential utility as histopathological markers. 3.7. Estradiol-responsive CAF fraction associates with ICI therapy resistance Although ICIs have shown remarkable efficacy against various cancers, understanding their mechanisms of resistance remains a critical challenge. Our identification of estradiol-responsive cellular populations with immunosuppressive properties raises the possibility that these cells may influence immunotherapy outcomes. The strong association between FB.TNFSF10 and immunosuppressive features suggests its potential role in treatment resistance. In addition, as previously discussed, the MAL.PDCD5 fraction showed an inverse correlation with PD-L1 expression, which is a target of ICI therapy ( Figure 6B ). We focused our subsequent analysis on the FB.TNFSF10 fraction and analyzed the distribution of FB.TNFSF10 between the responder and non-responder groups across the 18 ICI-treated cohorts. FB.TNFSF10-type patients exhibited significantly lower response rates compared to other FB-type patients (Fisher’s exact test P = 3.96 × 10⁻³ in the total cohort; Figure 7A ). This pattern was remarkably consistent across 15 of the 18 cohorts, with FB.TNFSF10 subtypes consistently underrepresented in responders compared with non-responders ( Figure 7B ). Download figure Open in new tab Figure 7. Effect of estradiol-related TME subtypes on ICI treatment response. A. Differences in FB subtype ratios between ICI responders and non-responders across the 18 cohorts. Asterisks indicate p values from Fisher’s exact test comparing the corresponding cell types between responders and non-responders. B. Comparison across 15 cohorts showing lower proportions of FB.TNFSF10-types in responders than in non-responders. C. Kaplan-Meier analysis of progression-free survival in three patient groups from the Liu et al. cohort: FB.TNFSF10-high, FB.TNFSF10-medium, and FB.TNFSF10-low. The first and third quartiles of the FB.TNFSF10 proportions were used as thresholds. The p-value was determined using the log-rank test. D. Schematic representation of estradiol-modified TME influencing ICI treatment response. To validate these findings in terms of clinical outcomes, we analyzed progression-free survival in a cohort study by Liu et al. [ 56 ]. Patients with low FB.TNFSF10 expression (fourth quartile) demonstrated significantly longer progression-free survival than those with medium (second and third quartiles) or high (first quartile) expression levels (log-rank test, p = 0.034; Figure 7C ). These results establish that FB.TNFSF10 is a key mediator of immunotherapy resistance. In tumors exposed to estrogen, including estradiol, increased distribution of estradiol-responsive cell fractions contributes to an immunosuppressive microenvironment through two mechanisms: MAL.PDCD5-mediated inhibition of immune cell infiltration through reduced antigen presentation and oxidative phosphorylation, and FB.TNFSF10-driven immunosuppression through enhanced CAF.S1 and TGF-β CAF populations, and increased Treg infiltration ( Figure 7D ). The consistency of these findings across multiple cohorts suggests a conserved mechanism of resistance, with implications for patient stratification and therapeutic intervention. These results provide a rationale for developing combination strategies that target estrogen-responsive stromal populations to enhance immunotherapy efficacy. 4. Discussion Our biofabricated tumor model system addresses critical unmet needs in cancer research by revealing estradiol-responsive cellular fractions that drive immunosuppression across diverse cancer types. More importantly, this work establishes a transformative biofabrication platform that could revolutionize how we model complex tissue environments and develop therapeutic strategies. The biofabrication approach employed in this study offers several key advantages over traditional culture systems. First, our patient-derived 3D culture platform maintains the cellular heterogeneity and architectural complexity characteristic of native tumors while enabling systematic manipulation of growth factor environments. This approach allowed us to identify specific estradiol-responsive populations that would be difficult to characterize in conventional 2D culture systems [ 3 ]. Second, our system provides a more comprehensive model of the TME, enabling investigation of tumor-stromal interactions that are critical for understanding immunotherapy resistance [ 16 ]. The discovery that FB.TNFSF10 correlates with immunotherapy resistance across multiple cancer cohorts represents a significant advancement in biomarker development. This finding is particularly noteworthy because it transcends traditional boundaries of hormone-dependent cancers, suggesting a fundamental biological principle that could revolutionize patient stratification for immunotherapy [ 84 ]. The consistency of this correlation across various malignancies indicates a conserved mechanism of estradiol-mediated immune modulation that can be effectively modeled using our biofabrication platform. Our spatial transcriptomic analysis of engineered tissue constructs revealed compelling evidence of immune cell exclusion mediated by FB.TNFSF10 cells. This spatial validation provides crucial contextual evidence for the functional role of these cells within the complex tissue architecture, highlighting the importance of 3D tissue models in understanding cellular interactions. The mutual exclusivity between FB.TNFSF10 cells and T/NK cells observed in our biofabricated models suggests active immune exclusion mechanisms that would be difficult to identify in traditional culture systems. Our systematic growth factor profiling approach represents a paradigm shift in biofabrication methodology that extends far beyond cancer research. The platform’s ability to systematically optimize culture conditions for patient-specific requirements makes it applicable to engineering diverse tissue constructs, from organ-on-chip models to regenerative medicine applications. The 128-combination screening approach can be adapted for any tissue type, providing a standardized framework for biofabrication across multiple research domains. This systematic approach enables the development of reproducible protocols that can be implemented across different laboratories and manufacturing facilities, addressing a critical need in the biofabrication field for consistent methods that can be translated from research to clinical applications. The integration of 3D culture systems with comprehensive molecular profiling offers unprecedented opportunities for pharmaceutical applications and high-throughput drug screening. The principles underlying our growth factor-responsive culture systems can be adapted for engineering other complex tissue constructs that require precise control of cellular heterogeneity and signaling environments. Potential applications include modeling neurodegenerative diseases, cardiovascular disorders, and inflammatory conditions. The platform’s ability to maintain patient-specific characteristics while enabling systematic experimental manipulation provides a powerful tool for understanding tissue development and disease progression across multiple biological systems, ultimately advancing personalized medicine approaches beyond oncology. 5. Conclusion The potential for translating our biofabrication platform to clinical applications is significant. The ability to rapidly generate patient-specific tumor models for drug testing could revolutionize personalized medicine approaches in oncology. Furthermore, the identification of biomarkers using our 3D culture systems could lead to more accurate predictive tests for immunotherapy responses, ultimately improving patient outcomes and reducing healthcare costs. Data and code availability Raw FASTQ files of the whole-exome, bulk RNA, and scRNA sequencing data were deposited in the European Genome-Phenome Archive (EGAS50000000422) and are available from the corresponding authors upon reasonable request. TCGA gene expression data (STAR-Counts) were downloaded using the R package TCGAbiolinks, and clinical data were obtained from cBioPortal (https://www.cbioportal.org/study/summary?id=ov_tcga_pub/). The public ovarian cancer bulk RNA-seq data reanalyzed in this study are available in the Gene Expression Omnibus under accession codes GSE51088 and GSE26193. The batch-corrected transcriptomes of 11,096 samples across 33 TCGA tumor types are available in the PancanAtlas ( https://gdc.cancer.gov/about-data/publications/pancanatlas ). All relevant packages and software information are provided in the Star Methods and Supplementary Table 12. No custom codes were used in this study. Declaration of interests The authors declare that they have no conflicts of interest. Supplementary figures and legends Download figure Open in new tab Supplementary Figure 1. Correlations between responsiveness to the seven GFs. A. Normalized cell growth of samples from clusters C1 and C2 in the presence of NAC, FGF10, EGF, or Noggin. B. Pearson correlation coefficients between the indicated GFs. Download figure Open in new tab Supplementary Figure 2. Top 50 differentially expressed genes (DEGs) used in C1/C2 clustering. The heatmap shows the z-score-normalized expression of the top 25 upregulated and downregulated DEGs. Download figure Open in new tab Supplementary Figure 3. GF responsiveness of four HGSOC consensus subtypes A–D. Growth scores of DIF (A), PRO (B), IMR (C), and MES (D) subtype samples incubated in the presence or absence of the indicated GFs. Box plots span the first to third quartiles and whiskers represent the 1.5 interquartile range. The central line indicates the median value. P -values were calculated using the Mann–Whitney U test. Download figure Open in new tab Supplementary Figure 4. Mutation signatures associated with the GF responses. A. Oncoplot of 20 frequently mutated genes and BRCA1/2 in the 16 EOC samples. Mutation type, C1/C2 cluster, pathological cancer type, and consensus subtype of each sample are shown. B. Volcano plot showing pairs of mutations and GFs that exhibited significant differences in growth scores between mutant and wild-type cells under the corresponding GF conditions. The x-axis represents the difference in cell growth scores, and the y-axis shows the significance of cell growth score differences, as assessed using the Mann–Whitney U test. Dot size indicates mutation frequency across samples. Download figure Open in new tab Supplementary Figure 5. GF responsiveness according to the mutation status. Growth scores of samples harboring mutations (red line) in the indicated genes or of WT cells (gray line) cultured in the presence or absence of the indicated GFs. Box plots span the first to third quartiles and whiskers represent the 1.5 interquartile range. The central line indicates the median value. P -values were calculated using the Mann–Whitney U test. Download figure Open in new tab Supplementary Figure 6. Mutational L-landscape of C1/C2-type samples in the TCGA cohort. Alterations in the top 20 most frequently mutated genes in TCGA-OV cohort. The mutation types and C1/C2 statuses are shown. Download figure Open in new tab Supplementary Figure 7. Clinical characteristics of C1/C2-type samples in TCGA cohort. A. Seven clinical characteristics of C1/C2-type samples obtained from cBioPortal: age at diagnosis, tumor mutation burden, ESR1 expression (log2[TPM+1]), aneuploidy score, winter hypoxia score, MSI score (MANTIS), and mutation count. B. Kaplan–Meier analysis of overall survival (OS) and progression-free interval (PFI) in C1- and C2-type samples. Statistical significance was calculated using log-rank tests. C. Classification probability of C1- and C2-type samples according to HGSOC consensus subtypes. Classification probability was calculated using the R package ConsensusOV. Box plots span the first to third quartiles and whiskers represent the 1.5 interquartile range. The central line indicates the median value. P -values were calculated using the Mann–Whitney U test. Download figure Open in new tab Supplementary Figure 8. ESR signatures show no association with consensus molecular subtypes in the TCGA-OV cohort. a. Sankey plot illustrating the relationships between C1/2 subtypes, ESR signatures, and HGSOC consensus subtypes. b. Distribution of consensus subtypes across the four groups: C1 and ESR-high, C1 and ESR-low, C2 and ESR-high, and C2 and ESR-low. The threshold for ESR-high/low subtypes was defined as the mean value of the ESR-mediated signalling (R-HSA-8939211) enrichment score (ES). c. Classification probability of ESR-high (red line) and ESR-low (blue line) subtypes to HGSOC consensus subtypes. Box plots span the first to third quartiles and whiskers represent the 1.5 interquartile range. The central line indicates the median value. Statistical significance was calculated using the Mann–Whitney U test. d. Functional enrichment analysis results of DEGs between ESR-high (red) and ESR-low (sky blue) subtypes. Red bars represent enriched functions of the upregulated DEGs in ESR-high subtypes, whereas sky-blue bars represent enriched functions in ESR-low subtypes. Curated gene sets from the Gene Ontology, Reactome, and Transcription Factor databases were used. Download figure Open in new tab Supplementary Figure 9. GPER signatures coordinate MES subtypes A. Sankey plot showing the relationships between C1/2-type samples based on GPER1 signaling (R-HSA-9634597) and HGSOC consensus subtypes. B. Stacked bar plot showing the distribution of consensus subtypes across four groups: C1&GPER-high, C1&GPER-low, C2&GPER-high, and C2&GPER-low. The threshold for GPER high/low subtypes was defined as the mean GPER1 signaling (R-HSA-9634597) ES value. C. Classification probability of GPER-high (red) and GPER-low (sky blue) subtypes for HGSOC consensus subtypes. Box plots span from the first to third quartiles, and the whiskers represent the 1.5 interquartile range. The central line indicates the median value. Statistical significance was calculated using the Mann–Whitney U test. D. Functional enrichment analysis results of DEGs between GPER-high (red) and GPER-low (sky blue) subtypes. Curated gene sets from the Gene Ontology, Reactome, and TRANSFAC databases were used. E. Bar plots displaying GPER1 expression levels [log2(TPM+1)] (top) and ESs of the GPER1 signaling pathway (R-HSA-9634597, bottom) in the indicated ovarian cancer cell lines. F. Bar plots show subtype scores for DIF (top) and MES (bottom) across the indicated ovarian cancer cell lines. Subtype scores represent the classification probabilities calculated using the R package ConsensusOV. G. ES of the GPER1 signaling pathway in PA-1 cells treated or not (control) with G15 (GPER1 inhibitor) (50 μM, 12 h, top). Stacked bar plot showing relative subtype scores between the DIF and MES subtypes in the control and G15-treated cells (bottom). Download figure Open in new tab Supplementary Figure 10. Association of GPER signatures with mesenchymal subtype-related genes in two ovarian cancer cohorts (GSE51088 and GSE26193). A, C Sankey plots illustrating the relationships among C1/2 type samples, GPER signatures, and HGSOC consensus subtypes in the GSE51088 (A) and GSE26193 (C) cohorts. B, D Distribution of consensus subtypes across the four groups in GSE51088 (B) and GSE26193 (D) cohorts: C1 and GPER-high, C1 and GPER-low, C2 and GPER-high, and C2 and GPER-low. The threshold for the GPER-high/low subtypes was defined as the mean value of GPER1 signalling (R-HSA-9634597) ES. Download figure Open in new tab Supplementary Figure 11. General single-cell characteristics of four ovarian cancer samples. A. Distribution of cell populations across four ovarian cancer samples in the presence or absence of estradiol based on single-cell RNA sequencing. B. Uniform Manifold Approximation and Projection plot showing cell populations across four ovarian cancer samples. C–E. Subtypes of malignant cells (C), fibroblasts (D), and macrophages (E) were identified using PRIMUS analysis. The subtyping data are summarized in UMAP plots and color-coded according to the cell subtypes (left). The top five markers for each subtype are shown as dotted plots (right panels). Color bars indicate the average expression of genes in the subtypes, and the dot size represents the percentage of cells expressing these genes. Markers were selected based on the fold change between cells. Download figure Open in new tab Supplementary Figure 12. Expression pattern of major cell markers that define the distinct cell fractions. a–f. Violin plots showing representative marker expression of the indicated cell types: fibroblasts (A), malignant cells (B), macrophages (C), dendritic cells (D), B cells (E), and T cells (F). Download figure Open in new tab Supplementary Figure 13. Representative gene sets of major cell clusters in single-cell analysis. A–C. Heatmaps (left) and dot blots (right) show the normalized ESs of the indicated cell subtypes for malignant cells (A), fibroblasts (B), and macrophages (C). Dot color indicates the adjusted P -value of enrichment, and dot size represents the overlap between representative gene sets and curated gene sets from Gene Ontology. Download figure Open in new tab Supplementary Figure 14. Analysis of CAF subtype signatures in single-cell fibroblast fractions. A. ESs of five fibroblast subtypes for CAF markers and subtype markers. B, C. Gene expression levels of CAF-positive and -negative markers. D–F. Gene expression levels in three CAF subtypes: myofibroblastic CAFs (myCAF), inflammatory CAFs (iCAF), and antigen-presenting CAFs (apCAF). Download figure Open in new tab Supplementary Figure 15. Relative fractions of cell subtypes in TCGA-OV samples. The indicated cell fractions were determined by CIBERSORTx cell deconvolution of C1-(red) and C2-(green)-type samples. Download figure Open in new tab Supplementary Figure 16. Ligand–receptor interaction between malignant cells and fibroblasts. Differential communication probability of ligand–receptor pairs that increased (left) or decreased (right) in estradiol-treated (sky blue) versus non-treated (pink) conditions between MAL.PDCD5 and TNFSF10 (A) and between MAL.TIMP1 and FB.TNFSF10 (B). Download figure Open in new tab Supplementary Figure 17. Analysis of CAF fractions in established CAF lines Bar plots show CAF fractions calculated using CIBERSORTx in established CAFs under estradiol-treated and estradiol-untreated conditions. Download figure Open in new tab Supplementary Figure 18. The cell deconvolution results of FB.TNFSF10 and MAL.PDCD5 cells across 33 TCGA cancer types A, B. The normalized proportions of the two subtypes, FB.TNFSF10 (A) and MAL.PDCD5 (B), were measured across 33 TCGA cancer types using CIBERSORTx. Download figure Open in new tab Supplementary Figure 19. Spatial associations between FB fractions and immune infiltration. A. The Xenium dataset was divided into a 7×8 grid (56 tiles, labeled S1-1 through S7-8), and cell type proportions were quantified for each tile. Stacked bar plots show the relative abundance of different FB subtypes and T/NK cells. B. Correlation analysis between normalized FB subtype proportions and T/NK cell fractions across the tiles. Pearson’s correlation coefficients (PCC) were calculated for each fibroblast subtype. Black dots represent individual tiles (n = 56) and the red line indicates the linear regression fit with a 95% confidence interval (shaded area). C , TGF-β CAF enrichment scores in a Xenium data (left). Box plots show the quartiles and 1.5 IQR whiskers. The median is marked by the central line. Significance was assessed using the Mann-Whitney U test (right). Acknowledgments This work was supported by the Creative-Pioneering Researchers Programme in SNU; the Cooperative Research Program of Basic Medical Science and Clinical Science from SNUCM (800-20240351), the National Research Foundation (NRF) of Korea funded by the Korean government (MIST) (NRF2020M3A9D803800912, NRF2022R1A5A102641311, RS-2023-00272547, RS-2023-00222910, and RS-2023-00208519); the Ministry of Science and ICT (2024-22030007-30); the Global Physician-Scientist Training Program from the Korea Health Industry Development Institute (KHIDI, RS-2024-00440302); and the Commercialization Promotion Agency for R&D Outcomes (COMPA). Footnotes All authors declare no conflict of interest. REFERENCES [1]. ↵ Moroni L , Burdick J A , Highley C , Lee S J , Morimoto Y , Takeuchi S and Yoo J J 2018 Biofabrication strategies for 3D in vitro models and regenerative medicine Nat Rev Mater 3 21 – 37 OpenUrl CrossRef PubMed [2]. ↵ Baptista L S , Mironov V , Koudan E , Amorim E A , Pampolha T P , Kasyanov V , Kovalev A , Senatov F and Granjeiro J M 2024 Bioprinting Using Organ Building Blocks: Spheroids, Organoids, and Assembloids Tissue Eng Part A 30 377 – 86 OpenUrl PubMed [3]. ↵ Jensen C and Teng Y 2020 Is It Time to Start Transitioning From 2D to 3D Cell Culture? Front Mol Biosci 7 33 OpenUrl PubMed [4]. ↵ Pampaloni F , Reynaud E G and Stelzer E H 2007 The third dimension bridges the gap between cell culture and live tissue Nat Rev Mol Cell Biol 8 839 – 45 OpenUrl CrossRef PubMed Web of Science [5]. ↵ Clevers H 2016 Modeling Development and Disease with Organoids Cell 165 1586 – 97 OpenUrl CrossRef PubMed [6]. Tuveson D and Clevers H 2019 Cancer modeling meets human organoid technology Science 364 952 – 5 OpenUrl Abstract / FREE Full Text [7]. Murphy S V and Atala A 2014 3D bioprinting of tissues and organs Nat Biotechnol 32 773 – 85 OpenUrl CrossRef PubMed [8]. ↵ Whitesides G M 2006 The origins and the future of microfluidics Nature 442 368 – 73 OpenUrl CrossRef PubMed Web of Science [9]. ↵ Gjorevski N , Sachs N , Manfrin A , Giger S , Bragina M E , Ordonez-Moran P , Clevers H and Lutolf M P 2016 Designer matrices for intestinal stem cell and organoid culture Nature 539 560 – 4 OpenUrl CrossRef PubMed [10]. ↵ Drost J and Clevers H 2018 Organoids in cancer research Nat Rev Cancer 18 407 – 18 OpenUrl CrossRef PubMed [11]. ↵ McGranahan N and Swanton C 2017 Clonal Heterogeneity and Tumor Evolution: Past , Present, and the Future Cell 168 613 – 28 OpenUrl PubMed [12]. ↵ Burrell R A , McGranahan N , Bartek J and Swanton C 2013 The causes and consequences of genetic heterogeneity in cancer evolution Nature 501 338 – 45 OpenUrl CrossRef PubMed Web of Science [13]. ↵ Quail D F and Joyce J A 2013 Microenvironmental regulation of tumor progression and metastasis Nat Med 19 1423 – 37 OpenUrl CrossRef PubMed [14]. ↵ Hanahan D and Coussens L M 2012 Accessories to the crime: functions of cells recruited to the tumor microenvironment Cancer Cell 21 309 – 22 OpenUrl CrossRef PubMed Web of Science [15]. ↵ Li Y , Liu Q , Jing X , Wang Y , Jia X , Yang X and Chen K 2025 Cancer-Associated Fibroblasts: Heterogeneity, Cancer Pathogenesis, and Therapeutic Targets MedComm (2020) 6 e70292 OpenUrl [16]. ↵ Sahai E , Astsaturov I , Cukierman E , DeNardo D G , Egeblad M , Evans R M , Fearon D , Greten F R , Hingorani S R , Hunter T , Hynes R O , Jain R K , Janowitz T , Jorgensen C , Kimmelman A C , Kolonin M G , Maki R G , Powers R S , Pure E , Ramirez D C , Scherz-Shouval R , Sherman M H , Stewart S , Tlsty T D , Tuveson D A , Watt F M , Weaver V , Weeraratna A T and Werb Z 2020 A framework for advancing our understanding of cancer-associated fibroblasts Nat Rev Cancer 20 174 – 86 OpenUrl CrossRef PubMed [17]. ↵ Rothenberger N J , Somasundaram A and Stabile L P 2018 The Role of the Estrogen Pathway in the Tumor Microenvironment Int J Mol Sci 19 [18]. ↵ Kopper O , de Witte C J , Lohmussaar K , Valle-Inclan J E , Hami N , Kester L , Balgobind A V , Korving J , Proost N , Begthel H , van Wijk L M , Revilla S A , Theeuwsen R , van de Ven M , van Roosmalen M J , Ponsioen B , Ho V W H , Neel B G , Bosse T , Gaarenstroom K N , Vrieling H , Vreeswijk M P G , van Diest P J , Witteveen P O , Jonges T , Bos J L , van Oudenaarden A , Zweemer R P , Snippert H J G , Kloosterman W P and Clevers H 2019 An organoid platform for ovarian cancer captures intra- and interpatient heterogeneity Nat Med 25 838 – 49 OpenUrl CrossRef PubMed [19]. ↵ Sulaiman R , De P , Aske J C , Lin X , Dale A , Koirala N , Gaster K , Espaillat L R , Starks D and Dey N 2023 Patient-Derived Primary Cancer-Associated Fibroblasts Mediate Resistance to Anti-Angiogenic Drug in Ovarian Cancers Biomedicines 11 [20]. ↵ Garcia M , Juhos S , Larsson M , Olason P I , Martin M , Eisfeldt J , DiLorenzo S , Sandgren J , Diaz De Stahl T , Ewels P , Wirta V , Nister M , Kaller M and Nystedt B 2020 Sarek: A portable workflow for whole-genome sequencing analysis of germline and somatic variants F1000Res 9 63 OpenUrl [21]. ↵ Li H and Durbin R 2010 Fast and accurate long-read alignment with Burrows-Wheeler transform Bioinformatics 26 589 – 95 OpenUrl CrossRef PubMed Web of Science [22]. ↵ Chen Z , Yuan Y , Chen X , Chen J , Lin S , Li X and Du H 2020 Systematic comparison of somatic variant calling performance among different sequencing depth and mutation frequency Sci Rep 10 3501 OpenUrl CrossRef PubMed [23]. ↵ Kim S , Scheffler K , Halpern A L , Bekritsky M A , Noh E , Kallberg M , Chen X , Kim Y , Beyter D , Krusche P and Saunders C T 2018 Strelka2: fast and accurate calling of germline and somatic variants Nat Methods 15 591 – 4 OpenUrl CrossRef PubMed [24]. ↵ McLaren W , Gil L , Hunt S E , Riat H S , Ritchie G R , Thormann A , Flicek P and Cunningham F 2016 The Ensembl Variant Effect Predictor Genome Biol 17 122 OpenUrl CrossRef PubMed [25]. ↵ Deng Z L , Munch P C , Mreches R and McHardy A C 2022 Rapid and accurate identification of ribosomal RNA sequences via deep learning Nucleic Acids Res 50 e60 OpenUrl CrossRef PubMed [26]. ↵ Dobin A , Davis C A , Schlesinger F , Drenkow J , Zaleski C , Jha S , Batut P , Chaisson M and Gingeras T R 2013 STAR: ultrafast universal RNA-seq aligner Bioinformatics 29 15 – 21 OpenUrl CrossRef PubMed Web of Science [27]. ↵ Li B and Dewey C N 2011 RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome BMC Bioinformatics 12 323 OpenUrl CrossRef PubMed [28]. ↵ Zhang Y , Parmigiani G and Johnson W E 2020 ComBat-seq: batch effect adjustment for RNA-seq count data NAR Genom Bioinform 2 lqaa078 OpenUrl [29]. ↵ Kolberg L , Raudvere U , Kuzmin I , Vilo J and Peterson H 2020 gprofiler2 -- an R package for gene list functional enrichment analysis and namespace conversion toolset g:Profiler F1000Res 9 [30]. ↵ Malta T M , Sokolov A , Gentles A J , Burzykowski T , Poisson L , Weinstein J N , Kaminska B , Huelsken J , Omberg L , Gevaert O , Colaprico A , Czerwinska P , Mazurek S , Mishra L , Heyn H , Krasnitz A , Godwin A K , Lazar A J , Cancer Genome Atlas Research N , Stuart J M , Hoadley K A , Laird P W , Noushmehr H and Wiznerowicz M 2018 Machine Learning Identifies Stemness Features Associated with Oncogenic Dedifferentiation Cell 173 338 – 54 e15 OpenUrl CrossRef PubMed [31]. ↵ Aran D , Hu Z and Butte A J 2017 xCell: digitally portraying the tissue cellular heterogeneity landscape Genome Biol 18 220 OpenUrl CrossRef PubMed [32]. ↵ Newman A M , Steen C B , Liu C L , Gentles A J , Chaudhuri A A , Scherer F , Khodadoust M S , Esfahani M S , Luca B A , Steiner D , Diehn M and Alizadeh A A 2019 Determining cell type abundance and expression from bulk tissues with digital cytometry Nat Biotechnol 37 773 – 82 OpenUrl CrossRef PubMed [33]. ↵ Chen G M , Kannan L , Geistlinger L , Kofia V , Safikhani Z , Gendoo D M A , Parmigiani G , Birrer M , Haibe-Kains B and Waldron L 2018 Consensus on Molecular Subtypes of High-Grade Serous Ovarian Carcinoma Clin Cancer Res 24 5037 – 47 OpenUrl Abstract / FREE Full Text [34]. ↵ Desbois M and Wang Y 2021 Cancer-associated fibroblasts: Key players in shaping the tumor immune microenvironment Immunol Rev 302 241 – 58 OpenUrl PubMed [35]. ↵ Liberzon A , Birger C , Thorvaldsdottir H , Ghandi M , Mesirov J P and Tamayo P 2015 The Molecular Signatures Database (MSigDB) hallmark gene set collection Cell Syst 1 417 – 25 OpenUrl PubMed [36]. ↵ Heaton H , Talman A M , Knights A , Imaz M , Gaffney D J , Durbin R , Hemberg M and Lawniczak M K N 2020 Souporcell: robust clustering of single-cell RNA-seq data by genotype without reference genotypes Nat Methods 17 615 – 20 OpenUrl CrossRef PubMed [37]. ↵ Huang Y , McCarthy D J and Stegle O 2019 Vireo: Bayesian demultiplexing of pooled single-cell RNA-seq data without genotype reference Genome Biol 20 273 OpenUrl CrossRef PubMed [38]. ↵ Kang H M , Subramaniam M , Targ S , Nguyen M , Maliskova L , McCarthy E , Wan E , Wong S , Byrnes L , Lanata C M , Gate R E , Mostafavi S , Marson A , Zaitlen N , Criswell L A and Ye C J 2018 Multiplexed droplet single-cell RNA-sequencing using natural genetic variation Nat Biotechnol 36 89 – 94 OpenUrl CrossRef PubMed [39]. ↵ Hao Y , Stuart T , Kowalski M H , Choudhary S , Hoffman P , Hartman A , Srivastava A , Molla G , Madad S , Fernandez-Granda C and Satija R 2024 Dictionary learning for integrative, multimodal and scalable single-cell analysis Nat Biotechnol 42 293 – 304 OpenUrl CrossRef PubMed [40]. ↵ Hafemeister C and Satija R 2019 Normalization and variance stabilization of single-cell RNA-seq data using regularized negative binomial regression Genome Biol 20 296 OpenUrl CrossRef PubMed [41]. ↵ Korsunsky I , Millard N , Fan J , Slowikowski K , Zhang F , Wei K , Baglaenko Y , Brenner M , Loh P R and Raychaudhuri S 2019 Fast, sensitive and accurate integration of single-cell data with Harmony Nat Methods 16 1289 – 96 OpenUrl CrossRef PubMed [42]. ↵ Staples J , Qiao D , Cho M H , Silverman E K , University of Washington Center for Mendelian G , Nickerson D A and Below J E 2014 PRIMUS: rapid reconstruction of pedigrees from genome-wide estimates of identity by descent Am J Hum Genet 95 553 – 64 OpenUrl CrossRef PubMed [43]. ↵ Zhang K , Erkan E P , Jamalzadeh S , Dai J , Andersson N , Kaipio K , Lamminen T , Mansuri N , Huhtinen K , Carpen O , Hietanen S , Oikkonen J , Hynninen J , Virtanen A , Hakkinen A , Hautaniemi S and Vaharautio A 2022 Longitudinal single-cell RNA-seq analysis reveals stress-promoted chemoresistance in metastatic ovarian cancer Sci Adv 8 eabm1831 OpenUrl CrossRef PubMed [44]. ↵ Su C , Xu Z , Shan X , Cai B , Zhao H and Zhang J 2022 Cell-type-specific co-expression inference from single cell RNA-sequencing data bioRxiv [45]. ↵ Langfelder P and Horvath S 2008 WGCNA: an R package for weighted correlation network analysis BMC Bioinformatics 9 559 OpenUrl CrossRef PubMed [46]. ↵ Hornburg M , Desbois M , Lu S , Guan Y , Lo A A , Kaufman S , Elrod A , Lotstein A , DesRochers T M , Munoz-Rodriguez J L , Wang X , Giltnane J , Mayba O , Turley S J , Bourgon R , Daemen A and Wang Y 2021 Single-cell dissection of cellular components and interactions shaping the tumor immune phenotypes in ovarian cancer Cancer Cell 39 928 – 44 e6 OpenUrl CrossRef PubMed [47]. ↵ Kieffer Y , Hocine H R , Gentric G , Pelon F , Bernard C , Bourachot B , Lameiras S , Albergante L , Bonneau C , Guyard A , Tarte K , Zinovyev A , Baulande S , Zalcman G , Vincent-Salomon A and Mechta-Grigoriou F 2020 Single-Cell Analysis Reveals Fibroblast Clusters Linked to Immunotherapy Resistance in Cancer Cancer Discov 10 1330 – 51 OpenUrl Abstract / FREE Full Text [48]. ↵ Jin S , Guerrero-Juarez C F , Zhang L , Chang I , Ramos R , Kuan C H , Myung P , Plikus M V and Nie Q 2021 Inference and analysis of cell-cell communication using CellChat Nat Commun 12 1088 OpenUrl CrossRef PubMed [49]. ↵ Liu J , Lichtenberg T , Hoadley K A , Poisson L M , Lazar A J , Cherniack A D , Kovatich A J , Benz C C , Levine D A , Lee A V , Omberg L , Wolf D M , Shriver C D , Thorsson V , Cancer Genome Atlas Research N and Hu H 2018 An Integrated TCGA Pan-Cancer Clinical Data Resource to Drive High-Quality Survival Outcome Analytics Cell 173 400 – 16 e11 OpenUrl CrossRef PubMed [50]. ↵ Gao F , Wang W , Tan M , Zhu L , Zhang Y , Fessler E , Vermeulen L and Wang X 2019 DeepCC: a novel deep learning-based framework for cancer molecular subtype classification Oncogenesis 8 44 OpenUrl CrossRef PubMed [51]. ↵ Kim S T , Cristescu R , Bass A J , Kim K M , Odegaard J I , Kim K , Liu X Q , Sher X , Jung H , Lee M , Lee S , Park S H , Park J O , Park Y S , Lim H Y , Lee H , Choi M , Talasaz A , Kang P S , Cheng J , Loboda A , Lee J and Kang W K 2018 Comprehensive molecular characterization of clinical responses to PD-1 inhibition in metastatic gastric cancer Nat Med 24 1449 – 58 OpenUrl CrossRef PubMed [52]. ↵ Huang A C , Orlowski R J , Xu X , Mick R , George S M , Yan P K , Manne S , Kraya A A , Wubbenhorst B , Dorfman L , D’Andrea K , Wenz B M , Liu S , Chilukuri L , Kozlov A , Carberry M , Giles L , Kier M W , Quagliarello F , McGettigan S , Kreider K , Annamalai L , Zhao Q , Mogg R , Xu W , Blumenschein W M , Yearley J H , Linette G P , Amaravadi R K , Schuchter L M , Herati R S , Bengsch B , Nathanson K L , Farwell M D , Karakousis G C , Wherry E J and Mitchell T C 2019 A single dose of neoadjuvant PD-1 blockade predicts clinical outcomes in resectable melanoma Nat Med 25 454 – 61 OpenUrl CrossRef PubMed [53]. ↵ Van Allen E M , Miao D , Schilling B , Shukla S A , Blank C , Zimmer L , Sucker A , Hillen U , Foppen M H G , Goldinger S M , Utikal J , Hassel J C , Weide B , Kaehler K C , Loquai C , Mohr P , Gutzmer R , Dummer R , Gabriel S , Wu C J , Schadendorf D and Garraway L A 2015 Genomic correlates of response to CTLA-4 blockade in metastatic melanoma Science 350 207 – 11 OpenUrl Abstract / FREE Full Text [54]. ↵ Riaz N , Havel J J , Makarov V , Desrichard A , Urba W J , Sims J S , Hodi F S , Martin-Algarra S , Mandal R , Sharfman W H , Bhatia S , Hwu W J , Gajewski T F , Slingluff C L , Jr. , Chowell D , Kendall S M , Chang H , Shah R , Kuo F , Morris L G T , Sidhom J W , Schneck J P , Horak C E , Weinhold N and Chan T A 2017 Tumor and Microenvironment Evolution during Immunotherapy with Nivolumab Cell 171 934 – 49 e16 OpenUrl CrossRef PubMed [55]. ↵ Chen P L , Roh W , Reuben A , Cooper Z A , Spencer C N , Prieto P A , Miller J P , Bassett R L , Gopalakrishnan V , Wani K , De Macedo M P , Austin-Breneman J L , Jiang H , Chang Q , Reddy S M , Chen W S , Tetzlaff M T , Broaddus R J , Davies M A , Gershenwald J E , Haydu L , Lazar A J , Patel S P , Hwu P , Hwu W J , Diab A , Glitza I C , Woodman S E , Vence L M , Wistuba , II , Amaria R N , Kwong L N , Prieto V , Davis R E , Ma W , Overwijk W W , Sharpe A H , Hu J , Futreal P A , Blando J , Sharma P , Allison J P , Chin L and Wargo J A 2016 Analysis of Immune Signatures in Longitudinal Tumor Samples Yields Insight into Biomarkers of Response and Mechanisms of Resistance to Immune Checkpoint Blockade Cancer Discov 6 827 – 37 OpenUrl Abstract / FREE Full Text [56]. ↵ Liu D , Schilling B , Liu D , Sucker A , Livingstone E , Jerby-Arnon L , Zimmer L , Gutzmer R , Satzger I , Loquai C , Grabbe S , Vokes N , Margolis C A , Conway J , He M X , Elmarakeby H , Dietlein F , Miao D , Tracy A , Gogas H , Goldinger S M , Utikal J , Blank C U , Rauschenberg R , von Bubnoff D , Krackhardt A , Weide B , Haferkamp S , Kiecker F , Izar B , Garraway L , Regev A , Flaherty K , Paschen A , Van Allen E M and Schadendorf D 2019 Integrative molecular and clinical modeling of clinical outcomes to PD1 blockade in patients with metastatic melanoma Nat Med 25 1916 – 27 OpenUrl CrossRef PubMed [57]. ↵ Hugo W , Zaretsky J M , Sun L , Song C , Moreno B H , Hu-Lieskovan S , Berent-Maoz B , Pang J , Chmielowski B , Cherry G , Seja E , Lomeli S , Kong X , Kelley M C , Sosman J A , Johnson D B , Ribas A and Lo R S 2016 Genomic and Transcriptomic Features of Response to Anti-PD-1 Therapy in Metastatic Melanoma Cell 165 35 – 44 OpenUrl CrossRef PubMed [58]. ↵ Gide T N , Quek C , Menzies A M , Tasker A T , Shang P , Holst J , Madore J , Lim S Y , Velickovic R , Wongchenko M , Yan Y , Lo S , Carlino M S , Guminski A , Saw R P M , Pang A , McGuire H M , Palendira U , Thompson J F , Rizos H , Silva I P D , Batten M , Scolyer R A , Long G V and Wilmott J S 2019 Distinct Immune Cell Populations Define Response to Anti-PD-1 Monotherapy and Anti-PD-1/Anti-CTLA-4 Combined Therapy Cancer Cell 35 238 – 55 e6 OpenUrl CrossRef PubMed [59]. ↵ Prat A , Navarro A , Pare L , Reguart N , Galvan P , Pascual T , Martinez A , Nuciforo P , Comerma L , Alos L , Pardo N , Cedres S , Fan C , Parker J S , Gaba L , Victoria I , Vinolas N , Vivancos A , Arance A and Felip E 2017 Immune-Related Gene Expression Profiling After PD-1 Blockade in Non-Small Cell Lung Carcinoma, Head and Neck Squamous Cell Carcinoma, and Melanoma Cancer Res 77 3540 – 50 OpenUrl Abstract / FREE Full Text [60]. ↵ Lee J S , Nair N U , Dinstag G , Chapman L , Chung Y , Wang K , Sinha S , Cha H , Kim D , Schperberg A V , Srinivasan A , Lazar V , Rubin E , Hwang S , Berger R , Beker T , Ronai Z , Hannenhalli S , Gilbert M R , Kurzrock R , Lee S H , Aldape K and Ruppin E 2021 Synthetic lethality-mediated precision oncology via the tumor transcriptome Cell 184 2487 – 502 e13 OpenUrl CrossRef PubMed [61]. ↵ Cho J W , Hong M H , Ha S J , Kim Y J , Cho B C , Lee I and Kim H R 2020 Genome-wide identification of differentially methylated promoters and enhancers associated with response to anti-PD-1 therapy in non-small cell lung cancer Exp Mol Med 52 1550 – 63 OpenUrl CrossRef PubMed [62]. ↵ Jung H , Kim H S , Kim J Y , Sun J M , Ahn J S , Ahn M J , Park K , Esteller M , Lee S H and Choi J K 2019 DNA methylation loss promotes immune evasion of tumours with high mutation and copy number load Nat Commun 10 4278 OpenUrl CrossRef PubMed [63]. ↵ Hwang S , Kwon A Y , Jeong J Y , Kim S , Kang H , Park J , Kim J H , Han O J , Lim S M and An H J 2020 Immune gene signatures for predicting durable clinical benefit of anti-PD-1 immunotherapy in patients with non-small cell lung cancer Sci Rep 10 643 OpenUrl CrossRef PubMed [64]. ↵ Miao D , Margolis C A , Gao W , Voss M H , Li W , Martini D J , Norton C , Bosse D , Wankowicz S M , Cullen D , Horak C , Wind-Rotolo M , Tracy A , Giannakis M , Hodi F S , Drake C G , Ball M W , Allaf M E , Snyder A , Hellmann M D , Ho T , Motzer R J , Signoretti S , Kaelin W G , Jr . ., Choueiri T K and Van Allen E M 2018 Genomic correlates of response to immune checkpoint therapies in clear cell renal cell carcinoma Science 359 801 – 6 OpenUrl Abstract / FREE Full Text [65]. ↵ Phillips D , Matusiak M , Gutierrez B R , Bhate S S , Barlow G L , Jiang S , Demeter J , Smythe K S , Pierce R H , Fling S P , Ramchurren N , Cheever M A , Goltsev Y , West R B , Khodadoust M S , Kim Y H , Schurch C M and Nolan G P 2021 Immune cell topography predicts response to PD-1 blockade in cutaneous T cell lymphoma Nat Commun 12 6726 OpenUrl CrossRef PubMed [66]. ↵ Mariathasan S , Turley S J , Nickles D , Castiglioni A , Yuen K , Wang Y , Kadel E E , III , Koeppen H , Astarita J L , Cubas R , Jhunjhunwala S , Banchereau R , Yang Y , Guan Y , Chalouni C , Ziai J , Senbabaoglu Y , Santoro S , Sheinson D , Hung J , Giltnane J M , Pierce A A , Mesh K , Lianoglou S , Riegler J , Carano R A D , Eriksson P , Hoglund M , Somarriba L , Halligan D L , van der Heijden M S , Loriot Y , Rosenberg J E , Fong L , Mellman I , Chen D S , Green M , Derleth C , Fine G D , Hegde P S , Bourgon R and Powles T 2018 TGFbeta attenuates tumour response to PD-L1 blockade by contributing to exclusion of T cells Nature 554 544 – 8 OpenUrl CrossRef PubMed [67]. ↵ Robinson M D and Oshlack A 2010 A scaling normalization method for differential expression analysis of RNA-seq data Genome Biol 11 R25 OpenUrl CrossRef PubMed [68]. ↵ Robinson M D , McCarthy D J and Smyth G K 2010 edgeR: a Bioconductor package for differential expression analysis of digital gene expression data Bioinformatics 26 139 – 40 OpenUrl CrossRef PubMed Web of Science [69]. ↵ Zhang Z , Murtagh F , Van Poucke S , Lin S and Lan P 2017 Hierarchical cluster analysis in clinical research with heterogeneous study population: highlighting its visualization with R Ann Transl Med 5 75 OpenUrl PubMed [70]. ↵ Hernandez-Vega A M , Del Moral-Morales A , Zamora-Sanchez C J , Pina-Medina A G , Gonzalez-Arenas A and Camacho-Arroyo I 2020 Estradiol Induces Epithelial to Mesenchymal Transition of Human Glioblastoma Cells Cells 9 [71]. ↵ Tang Y , Xu Q , Hu L , Yan X , Feng X , Yokota A , Wang W , Zhan D , Krishnamurthy D , Ochayon D E , Wen L , Huo L , Zeng H , Luo Y , Huang L F , Wunderlich M , Zhang J , Vivier E , Zhou J , Waggoner S N and Huang G 2021 Tumor Microenvironment-Derived R-spondins Enhance Antitumor Immunity to Suppress Tumor Growth and Sensitize for Immune Checkpoint Blockade Therapy Cancer Discov 11 3142 – 57 OpenUrl Abstract / FREE Full Text [72]. ↵ Mungenast F and Thalhammer T 2014 Estrogen biosynthesis and action in ovarian cancer Front Endocrinol (Lausanne) 5 192 OpenUrl PubMed [73]. O’Donnell A J , Macleod K G , Burns D J , Smyth J F and Langdon S P 2005 Estrogen receptor-alpha mediates gene expression changes and growth response in ovarian cancer cells exposed to estrogen Endocr Relat Cancer 12 851 – 66 OpenUrl Abstract / FREE Full Text [74]. ↵ Yan Y , Liu H , Wen H , Jiang X , Cao X , Zhang G and Liu G 2013 The novel estrogen receptor GPER regulates the migration and invasion of ovarian cancer cells Mol Cell Biochem 378 1 – 7 OpenUrl CrossRef PubMed [75]. ↵ Desbois M , Udyavar A R , Ryner L , Kozlowski C , Guan Y , Durrbaum M , Lu S , Fortin J P , Koeppen H , Ziai J , Chang C W , Keerthivasan S , Plante M , Bourgon R , Bais C , Hegde P , Daemen A , Turley S and Wang Y 2020 Integrated digital pathology and transcriptome analysis identifies molecular mediators of T-cell exclusion in ovarian cancer Nat Commun 11 5583 OpenUrl CrossRef PubMed [76]. ↵ Langdon S P , Herrington C S , Hollis R L and Gourley C 2020 Estrogen Signaling and Its Potential as a Target for Therapy in Ovarian Cancer Cancers (Basel) 12 [77]. ↵ Svoronos N , Perales-Puchalt A , Allegrezza M J , Rutkowski M R , Payne K K , Tesone A J , Nguyen J M , Curiel T J , Cadungog M G , Singhal S , Eruslanov E B , Zhang P , Tchou J , Zhang R and Conejo-Garcia J R 2017 Tumor Cell-Independent Estrogen Signaling Drives Disease Progression through Mobilization of Myeloid-Derived Suppressor Cells Cancer Discov 7 72 – 85 OpenUrl Abstract / FREE Full Text [78]. ↵ Ohlund D , Handly-Santana A , Biffi G , Elyada E , Almeida A S , Ponz-Sarvise M , Corbo V , Oni T E , Hearn S A , Lee E J , Chio , II , Hwang C I , Tiriac H , Baker L A , Engle D D , Feig C , Kultti A , Egeblad M , Fearon D T , Crawford J M , Clevers H , Park Y and Tuveson D A 2017 Distinct populations of inflammatory fibroblasts and myofibroblasts in pancreatic cancer J Exp Med 214 579 – 96 OpenUrl Abstract / FREE Full Text [79]. ↵ Elyada E , Bolisetty M , Laise P , Flynn W F , Courtois E T , Burkhart R A , Teinor J A , Belleau P , Biffi G , Lucito M S , Sivajothi S , Armstrong T D , Engle D D , Yu K H , Hao Y , Wolfgang C L , Park Y , Preall J , Jaffee E M , Califano A , Robson P and Tuveson D A 2019 Cross-Species Single-Cell Analysis of Pancreatic Ductal Adenocarcinoma Reveals Antigen-Presenting Cancer-Associated Fibroblasts Cancer Discov 9 1102 – 23 OpenUrl Abstract / FREE Full Text [80]. ↵ Vardhana S A , Hwee M A , Berisa M , Wells D K , Yost K E , King B , Smith M , Herrera P S , Chang H Y , Satpathy A T , van den Brink M R M , Cross J R and Thompson C B 2020 Impaired mitochondrial oxidative phosphorylation limits the self-renewal of T cells exposed to persistent antigen Nat Immunol 21 1022 – 33 OpenUrl CrossRef PubMed [81]. ↵ Ghahremanifard P , Chanda A , Bonni S and Bose P 2020 TGF-beta Mediated Immune Evasion in Cancer-Spotlight on Cancer-Associated Fibroblasts Cancers (Basel) 12 [82]. ↵ Togashi Y , Shitara K and Nishikawa H 2019 Regulatory T cells in cancer immunosuppression - implications for anticancer therapy Nat Rev Clin Oncol 16 356 – 71 OpenUrl CrossRef PubMed [83]. ↵ Jiang X , Wang J , Deng X , Xiong F , Ge J , Xiang B , Wu X , Ma J , Zhou M , Li X , Li Y , Li G , Xiong W , Guo C and Zeng Z 2019 Role of the tumor microenvironment in PD-L1/PD-1-mediated tumor immune escape Mol Cancer 18 10 OpenUrl CrossRef PubMed [84]. ↵ Sharma P , Hu-Lieskovan S , Wargo J A and Ribas A 2017 Primary, Adaptive, and Acquired Resistance to Cancer Immunotherapy Cell 168 707 – 23 OpenUrl CrossRef PubMed View the discussion thread. Back to top Previous Next Posted September 09, 2025. Download PDF Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following Systematic Growth Factor Profiling Platform for 3D Tumor Models Reveals Estradiol-Responsive Cellular Mechanisms of Immunotherapy Resistance Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share Systematic Growth Factor Profiling Platform for 3D Tumor Models Reveals Estradiol-Responsive Cellular Mechanisms of Immunotherapy Resistance Kwanghwan Lee , Minsung Kim , Si On Lim , Dong-Ju Shin , Yun Shin , Jung-Joo Choi , Maria Lee , Hyun Ju Kang , Jeong-Won Lee , Jin-Ku Lee bioRxiv 2025.09.04.674374; doi: https://doi.org/10.1101/2025.09.04.674374 Share This Article: Copy Citation Tools Systematic Growth Factor Profiling Platform for 3D Tumor Models Reveals Estradiol-Responsive Cellular Mechanisms of Immunotherapy Resistance Kwanghwan Lee , Minsung Kim , Si On Lim , Dong-Ju Shin , Yun Shin , Jung-Joo Choi , Maria Lee , Hyun Ju Kang , Jeong-Won Lee , Jin-Ku Lee bioRxiv 2025.09.04.674374; doi: https://doi.org/10.1101/2025.09.04.674374 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Cancer Biology Subject Areas All Articles Animal Behavior and Cognition (7640) Biochemistry (17706) Bioengineering (13902) Bioinformatics (41978) Biophysics (21465) Cancer Biology (18611) Cell Biology (25528) Clinical Trials (138) Developmental Biology (13387) Ecology (19920) Epidemiology (2067) Evolutionary Biology (24332) Genetics (15615) Genomics (22519) Immunology (17747) Microbiology (40424) Molecular Biology (17194) Neuroscience (88662) Paleontology (667) Pathology (2839) Pharmacology and Toxicology (4827) Physiology (7650) Plant Biology (15160) Scientific Communication and Education (2046) Synthetic Biology (4302) Systems Biology (9826) Zoology (2271)