Knockdown TNF family prognosis index crucial gene PDE4B promoted PANoptosis of ovarian carcinoma cell:Based in vitro and in vivo experiments.

Qianqian Yu: Writing – original draft, Software, Methodology, Investigation, Formal analysis, Data curation. Yunxiao Wang: Writing – original draft, Visualization, Formal analysis, Data curation. Ting Fu: Writing – original draft, Methodology, Formal analysis, Data curation. Dongyu Han: Writing – review & editing, Formal analysis, Data curation. Linlin Wang: Writing – review & editing, Validation, Project administration. Lin Zhao: Writing – review & editing, Validation, Supervision, Resources, Project administration, Methodology, Funding acquisition, Conceptualization. Yongle Xu: Writing – review & editing, Supervision, Resources, Methodology, Formal analysis, Data curation.

This study received no external funding.

We used R software to calculate and analyze the enrichment scores of TNF family genes, with each tumor sample corresponding to a TNF_Sig. We determined β=7 as the optimal soft threshold by assessing Scale Independence and Mean Connectivity and constructed a co-expression network ( Fig. 1 A). Genes were grouped into modules, and a Cluster Dendrogram was generated, where the MEgrey module represented unclassified genes ( Fig. 1 B). Subsequently, we generated a heatmap to display the correlation between different gene modules and Age, Stage, Alive, Dead, OS.time, and TNF_Sig. Most modules exhibited a correlation with TNF_Sig, with the MEyellow module showing the strongest correlation (R = 0.79, p 0.2 and GS > 0.2 criteria ( Fig. 1 D). Enrichment analysis of GO_BP was performed on genes from modules excluding the MEgrey module. The top five enriched pathways for each module were displayed. Among the gene co-expression modules identified, the MEpink module exhibited significantly higher enrichment in extracellular matrix organization, compared to external encapsulating structure organization, and extracellular structure organization pathways. Additionally, the cytoplasmic translation pathway in the MEgreenyellow module exhibited significantly higher expression compared to other pathways ( Fig. 1 E). Differential gene expression analysis was performed comparing cancer tissues to normal tissues. Then, volcano plot revealed SNX20, IL21R, SLAMF8, CCR5, PTPRC, LCP2, CSF2RB, P2RY10, PTPN22, and RHOH as upregulated genes in tumor tissues, while AC098476.1 , TRIM43, CCDC181, COL2A1, PRRT1B, GDF10, AC08484866.1, FGF17, TRIM9, and AC107057.1 were downregulated genes ( Fig. 2 A). Subsequently, we intersected highly expressed genes in tumor tissues with hub genes identified through WGCNA, resulting in 91 differentially expressed hub TNF family genes ( Fig. 2 B). We conducted enrichment analyses separately for the 91 intersected genes in CC, BP, and MF, and displayed the top 10 pathways for each analysis. The establishment of organelle localization pathway was significantly enriched in BP analysis, while the collagen-containing extracellular matrix pathway was significantly enriched in CC analysis, and GTPase regulator activity and nucleoside-triphosphatase regulator activity pathways were significantly enriched in MF analysis ( Fig. 2 C). Furthermore, univariate Cox analysis were performed on the intersected genes. The forest plot revealed GFPT2, PDE4B, PODNL1, COL8A1, VCAN, TNC, PTGIS, CSF1R, TGFBI, and COL5A2 as risk factors, while C2, CD3D, and SLAMF7 were identified as protective factors ( Fig. 2 D). Fig. 1 Discovery of TNF family related genes by WGCNA. (a) The distribution and trends of scale free topology model fit and mean connectivity along with soft threshold. (b) The clustering of gens among different modules. The gray modules represent unclassified genes. (c) Average correlation between multiple modules and TNF family related genes. The color of the cell indicates the strength of the correlation, and the number in parentheses indicates the p-value for the correlation test. (d) Correlation between module membership and gene significance in the yellow module. Dots in color were regarded as the hub genes of the corresponding module. (e) Top five enriched GO-BP terms of module genes of each module. Fig 1 Fig. 2 Identification of prognostic TNF family related genes. (a) Volcano plot illustrating the upregulated genes (colored in red) and downregulated genes (colored in blue) in OV cancer tissues. (b) 91 common TNF family related genes were obtained from WGCNA and high and low score groups. (c) The top ten enriched GO-BP/CC/MF terms by the 91 common TNF family related genes. (d) Forest plot showing the prognostic TNF family related genes identified by the univariate cox regression analysis. Red indicates risk factor whereas blue indicates protective factors. Fig 2 Discovery of TNF family related genes by WGCNA. (a) The distribution and trends of scale free topology model fit and mean connectivity along with soft threshold. (b) The clustering of gens among different modules. The gray modules represent unclassified genes. (c) Average correlation between multiple modules and TNF family related genes. The color of the cell indicates the strength of the correlation, and the number in parentheses indicates the p-value for the correlation test. (d) Correlation between module membership and gene significance in the yellow module. Dots in color were regarded as the hub genes of the corresponding module. (e) Top five enriched GO-BP terms of module genes of each module. Identification of prognostic TNF family related genes. (a) Volcano plot illustrating the upregulated genes (colored in red) and downregulated genes (colored in blue) in OV cancer tissues. (b) 91 common TNF family related genes were obtained from WGCNA and high and low score groups. (c) The top ten enriched GO-BP/CC/MF terms by the 91 common TNF family related genes. (d) Forest plot showing the prognostic TNF family related genes identified by the univariate cox regression analysis. Red indicates risk factor whereas blue indicates protective factors. A prognostic model generated by machine learning was devised and rigorously validated, leveraging the TCGA-OV dataset as the foundational training cohort and subsequently employing two GEO datasets as independent validation cohorts to ensure robustness and reproducibility. LASSO analysis was employed to determine the optimal parameter, λ=0.040, as depicted in Fig. 3 A. Furthermore, a multivariable Cox regression analysis was conducted to obtain coefficients for 14 model genes, which were then illustrated using a forest plot ( Fig. 3 B). The plot of cumulative risk factors unveiled the variation in risk scores among distinct cohorts, portraying their distribution patterns ( Fig. 3 C). Subsequently, according to the value of median scores, we stratified all cohorts into low-risk and high-risk groups. Fig. 3 Models Construction of prognostic TNF family related genes. (a) The selection of genes based on the optimal parameter λ that was obtained in the LASSO regression analysis. (b) Lollipop chart of the coefficients of prognostic model genes determined by the multiCox regression analysis. (c) Distribution of patients based on the median risk score in the TCGA-OV training set. Fig 3 Models Construction of prognostic TNF family related genes. (a) The selection of genes based on the optimal parameter λ that was obtained in the LASSO regression analysis. (b) Lollipop chart of the coefficients of prognostic model genes determined by the multiCox regression analysis. (c) Distribution of patients based on the median risk score in the TCGA-OV training set. The KM survival analysis indicated markedly better survival status in low-risk groups of both datasets compared to the high-risk groups, with progressively widening gaps in survival rates over time. Additionally, the analysis results of ROC indicated that the AUCs for all three datasets were greater than 0.6 at 1, 3, and 5 years, suggesting satisfactory discriminative performance of the constructed model at these time points ( Fig. 4 A-C). Fig. 4 Performance evaluation of model prognosis. (a) Survival differences between two groups in the TCGA-OV (up panel). Time-dependent ROC analysis of the model in the TCGA-OV (bottom panel). (b) Survival differences between two groups in the GSE102073 (up panel). Time-dependent ROC analysis of the model in the GSE102073 (bottom panel). (c) Survival differences between two groups in the GSE102073 (up panel). Time-dependent ROC analysis of the model in the GSE102073 (bottom panel). Fig 4 Performance evaluation of model prognosis. (a) Survival differences between two groups in the TCGA-OV (up panel). Time-dependent ROC analysis of the model in the TCGA-OV (bottom panel). (b) Survival differences between two groups in the GSE102073 (up panel). Time-dependent ROC analysis of the model in the GSE102073 (bottom panel). (c) Survival differences between two groups in the GSE102073 (up panel). Time-dependent ROC analysis of the model in the GSE102073 (bottom panel). The violin plot revealed that the risk score for individuals aged 65 and above was slightly higher than that for those below 65, although not statistically significant ( Fig. 5 A). Additionally, Stage II exhibited the lowest risk score, while Stage I had the highest, with marginal differences observed among Stage I, Stage III, and Stage IV ( Fig. 5 B). Furthermore, individuals of American Indian or Alaska Native descent exhibited the highest risk score, whereas Asians had the lowest ( Fig. 5 C). The application of Pearson correlation analysis unveiled a salient positive association linking the model's genetic signature to the risk score, particularly among individuals categorized as high-risk ( Fig. 5 D). In the TCGA-OV dataset, elevated TMB levels were observed in high-risk group. Moreover, results exhibited that the expression level of immune checkpoint genes and TNF family genes was relatively higher in high-risk group ( Fig. 6 A). However, immune cell infiltration analysis using EPIC, MCP-Counter, quanTIseq, TIMER, and xCell algorithms revealed no significant differences in the abundance status of immune cell between different risk groups ( Fig. 6 B). Fig. 5 Clinical analysis of the model. (a) The distribution of risk scores between the two age populations. (b) The distribution of risk scores across stages. (c) The distribution of risk scores across race. (d) Correlation between the risk score and the model genes. Fig 5 Fig. 6 Analysis of immune cell infiltration related to the model. (a) Differences in clinical data, TNF family genes, and check points between high and low-risk groups. TMB is presented in the form of a bar chart and density plot respectively (the heatmap shows the results with statistical differences). (b) The differences in the abundance of immune cell infiltration algorithms in IOBR R package including EPIC, MCP-Counter, quanTIseq, TIMER, and xCell between the high-risk and low-risk groups. Fig 6 Clinical analysis of the model. (a) The distribution of risk scores between the two age populations. (b) The distribution of risk scores across stages. (c) The distribution of risk scores across race. (d) Correlation between the risk score and the model genes. Analysis of immune cell infiltration related to the model. (a) Differences in clinical data, TNF family genes, and check points between high and low-risk groups. TMB is presented in the form of a bar chart and density plot respectively (the heatmap shows the results with statistical differences). (b) The differences in the abundance of immune cell infiltration algorithms in IOBR R package including EPIC, MCP-Counter, quanTIseq, TIMER, and xCell between the high-risk and low-risk groups. A comprehensive examination of single-cell sequencing datasets was undertaken, utilizing the UMAP method for dimensionality reduction on integrated single-cell sequencing data from three ovarian cancer datasets, yielding 16 distinct cell clusters ( Fig. 7 A-B). Utilizing TISCH-specific markers for each cell subtype, we annotated nine major cell subtypes: CD4Tconv, CD8Tex, DC, Endothelial, Fibroblasts, Malignant, Mono/Macro, Myofibroblasts, and Plasma ( Fig. 7 C). A heatmap illustrated that each cell subtype exhibited high expression of certain marker genes, with CLDN4, EPCAM, WFDC2, CD24, and SLPI significantly overexpressed in the Malignant subtype ( Fig. 7 D). Using the AddModuleScore, we calculated the TNF-related Signature for each cell and visualized it on the UMAP plot ( Fig. 7 E). Employing a resolution of 0.2, we further performed dimensionality reduction clustering on malignant cells, identifying six subgroups related to malignancy ( Fig. 8 A). The UMAP plot revealed relatively high expression levels of VCAN and C2 in malignant cells, while SLAMF7 and CSF1R exhibited low expression ( Fig. 8 B). We constructed developmental trajectories for each subgroup of malignant cells ( Fig. 8 C-D), demonstrating an increase in expression with pseudo-time for C2, whereas VCAN exhibited a decrease in distribution with pseudo-time ( Fig. 8 E). Fig. 7 The highly activated TNF-related signature in scRNA-seq datasets of OV. (a, b) UMAP visualization of single cells from public three OV scRNA-seq cohorts. A total of 16 cell clusters were identified. (c) 9 major cell types were manually annotated. (d) Heatmap illustrating the expression values of cell type-specific markers. (e) The signature genes expression at single cell level determined by AddModuleScore () function in Seurat. Fig 7 Fig. 8 The dynamic expression profiles of the signature genes revealed by the pseudotime analysis. (a) A total of six malignant subpopulations were identified under the single-cell resolution of 0.2. (b) UMAP visualization of the expression levels of signature genes at the single-cell level. (c, d) Pseudotime analysis uncovers the potential developmental pathways of the six subpopulations. (e) The dynamic expression profiles of the signature genes along with the trajectory. Fig 8 The highly activated TNF-related signature in scRNA-seq datasets of OV. (a, b) UMAP visualization of single cells from public three OV scRNA-seq cohorts. A total of 16 cell clusters were identified. (c) 9 major cell types were manually annotated. (d) Heatmap illustrating the expression values of cell type-specific markers. (e) The signature genes expression at single cell level determined by AddModuleScore () function in Seurat. The dynamic expression profiles of the signature genes revealed by the pseudotime analysis. (a) A total of six malignant subpopulations were identified under the single-cell resolution of 0.2. (b) UMAP visualization of the expression levels of signature genes at the single-cell level. (c, d) Pseudotime analysis uncovers the potential developmental pathways of the six subpopulations. (e) The dynamic expression profiles of the signature genes along with the trajectory. We identified differential genes between two risk groups. Results from Gene Set Enrichment Analysis (GSEA) indicated upregulation of KEGG pathways including Focal Adhesion, ECM Receptor Interaction, Hedgehog Signaling Pathway, Axon Guidance, and Pathways in Cancer in the low-risk group, while in the high-risk group, Hallmark pathways such as Apical Junction, Epithelial-Mesenchymal Transition, Myogenesis, UV Response Down, Early Estrogen Response, Angiogenesis, and Hypoxia were upregulated ( Fig. 9 A). Cancer signatures exhibited a notable enrichment within the high-risk group ( Fig. 9 B). Mantel plots revealed a high degree of mutual correlation between risk score, anti-cancer immunity cycle, and therapeutic-related pathways ( Fig. 9 C-D). Drug sensitivity analysis revealed that Acetalax demonstrated reduced IC50 values in the high-risk group, suggesting superior therapeutic efficacy in the former. Conversely, other drugs displayed enhanced efficacy in the low-risk group ( Fig. 10 ). Fig. 9 Biological functions. (a) Significantly enriched KEGG pathways in the high-risk and low-risk groups. The ridge plots in the left to zero indicates the upregulated pathways in low-risk and vice versa. (b) Significantly enriched cancer hallmarks in the high-risk and low-risk groups. (c, d) The correlation between the risk score and the anti-cancer immunity cycles (c) and the therapeutic-related pathways (d). Fig 9 Fig. 10 Therapeutic sensitivity between two risk groups. Fig 10 Biological functions. (a) Significantly enriched KEGG pathways in the high-risk and low-risk groups. The ridge plots in the left to zero indicates the upregulated pathways in low-risk and vice versa. (b) Significantly enriched cancer hallmarks in the high-risk and low-risk groups. (c, d) The correlation between the risk score and the anti-cancer immunity cycles (c) and the therapeutic-related pathways (d). Therapeutic sensitivity between two risk groups. In this study, we analyzed protein samples extracted from 20 cases of ovarian cancer and adjacent tissue samples. The results presented in Fig. 11 demonstrate high expression of PDE4B in ovarian cancer tissues. To explore the function of PDE4B in ovarian cancer, we examined its expression in ovarian epithelial cells and ovarian carcinoma cells using RT-PCR. Results indicated notably higher PDE4B expression in ovarian carcinoma cells compared to ovarian epithelial cells ( Fig. 11 B). We further tested PDE4B protein levels via Western blot, confirming the RT-PCR results ( Fig. 11 C). The above results show that the PDE4B gene high expression in ovarian cancer. Fig. 11 Expression status of PDE4B. (a) Expression of PDE4B in matched ovarian cancer and adjacent tissues. (b-c) mRNA and protein expression levels of PDE4B in cancer cells and normal cells. (d-e) Verification of PDE4B knockdown efficiency in Caov-3 and HEY-A8 cells. Fig 11 Expression status of PDE4B. (a) Expression of PDE4B in matched ovarian cancer and adjacent tissues. (b-c) mRNA and protein expression levels of PDE4B in cancer cells and normal cells. (d-e) Verification of PDE4B knockdown efficiency in Caov-3 and HEY-A8 cells. In order to find out the impact of PDE4B inhibition on aggressiveness of carcinoma cells, we utilized siRNA to downregulate PDE4B expression in HEY-A8 and Caov-3 cells. The efficiency of two siRNA sequences was confirmed by Western blot ( Fig. 11 D, E). The EdU assay showed a significant reduction in EdU-positive cells in the siPDE4B group compared to the control ( Fig. 12 A). CCK-8 assays revealed a significant decrease in cell proliferation within 4 days of PDE4B knockdown ( Fig. 12 B, C). Wound healing assays demonstrated a marked reduction in cell migration with PDE4B silencing ( Fig. 12 D, E), and Transwell migration assays indicated significant suppression of migration in both Caov-3 and HEY-A8 cells with PDE4B knockdown ( Fig. 12 F). Thus, our findings suggest that silencing PDE4B attenuates proliferation and migration in ovarian carcinoma cells. Fig. 12 Expression and Function of PDE4B in Ovarian Cancer Cells. (a-c) Analysis of cancer cell proliferation following PDE4B knockdown using EdU and CCK-8 assays. (d-e) Wound healing assay results and statistical analysis for Caov-3 and HEY-A8 cells at 48 hours. (f) Transwell migration efficiency and statistical analysis at 24 hours. (g-i) Expression of specific molecules involved in pyroptosis, apoptosis, and necroptosis during PANoptosis following PDE4B knockdown. Fig 12 Expression and Function of PDE4B in Ovarian Cancer Cells. (a-c) Analysis of cancer cell proliferation following PDE4B knockdown using EdU and CCK-8 assays. (d-e) Wound healing assay results and statistical analysis for Caov-3 and HEY-A8 cells at 48 hours. (f) Transwell migration efficiency and statistical analysis at 24 hours. (g-i) Expression of specific molecules involved in pyroptosis, apoptosis, and necroptosis during PANoptosis following PDE4B knockdown. PANoptosis, a distinctive form of inflammatory cell death, involves interactions between pyroptosis, apoptosis, and necroptosis. There is growing interest in its process and function. In this study, we knocked down PDE4B, an oncogene actively expressed in cancer, to explore the function in promoting tumor cell death. Compared to the NC group, the PDE4B knockdown group exhibited higher expression level of GSDMD, activated caspase 1 (CASP1) (p20), and GSDME N-terminal fragments ( Fig. 12 G). The PDE4B knockdown group also exhibited higher expression of CASP3, CASP7, and CASP8, associated with apoptosis pathways ( Fig. 12 H). Additionally, phosphorylation of MLKL, involved in necroptosis, was significantly enhanced after PDE4B knockdown ( Fig. 12 I). In summary, our findings indicate that the knockdown of PDE4B facilitates ovarian cancer cell death via PANoptosis, suggesting a promising therapeutic strategy for ovarian cancer treatment.

The "TCGAbiolinks" R package was employed to retrieve bulk transcriptomic data of ovarian cancer TCGA-OV and corresponding information of clinical characteristics from The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/ ). Utilizing the extensive resources of the Genotype-Tissue Expression (GTEx) database( www.org/home/index.html ), we acquired bulk transcriptomic data pertaining to normal ovarian tissues. Additionally, Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/ ) was used to procure two comprehensive bulk transcriptomic datasets, GSE26712 and GSE102073 , specifically focused on ovarian cancer. Single-cell sequencing datasets of ovarian cancer ( GSE115007 , GSE118828 , and GSE130000 ) were acquired from Tumor Immune Single-cell Hub (TISCH, https://comp-genomics.org ). All publicly available databases utilized in this study permit unrestricted access and usage, without requiring additional ethical approval. In adherence to pertinent regulations, our methodologies for data acquisition and subsequent analysis were rigorously executed. The list of TNF family genes, comprising 29 TNFRSF genes and 18 TNFSF, was obtained from a literature [ 11 ]. We computed the enrichment scores of TNF family genes using the "GSVA" package, with each tumor sample corresponding to a TNF Signature Score (TNF_Sig). WGCNA was implemented using the "WGCNA" R package. Genes with low expression levels or negligible differences across all samples were excluded. Subsequently, correlation matrices and adjacency matrices were constructed. We jointly selected an appropriate soft threshold of 7 based on Scale independence and Mean connectivity, constructed a co-expression network, divided genes into modules, and generated a Cluster Dendrogram. A heatmap was generated to demonstrate the associations linking distinct gene modules with clinical parameters, including Stage, Age, Alive, Dead, OS.time, and TNF_Sig. Finally, the MEyellow module was identified, which exhibited the highest correlation with TNF_Sig. Hub genes within the MEyellow module were selected based on GS (gene significance) > 0.2 and MM (module membership) > 0.2 criteria. Apart from the unclassified genes in the MEgrey module, genes in other modules underwent enrichment analysis of Biological Process in Gene Ontology (GO_BP). The top five enriched pathways in each module were presented. We conducted a differential gene expression analysis comparing ovarian cancer tissues with normal tissues utilizing R package "limma", with criteria of |logFC| > 1 & adjusted P-value(adj.p.val.) < 0.01 for selecting differentially expressed genes, visualized through volcano plots. Subsequently, the intersection of genes highly expressed in tumor tissues and hub genes identified by WGCNA was obtained, yielding differentially expressed hub TNF family genes. The intersected genes underwent GO- MF, CC, and BP enrichment analyses, with the top 10 pathways in each analysis presented. Additionally, we utilized the R software package "ezcox" to conduct univariate Cox analysis on the intersected genes and visualized results using forest plots, further refining the selection of prognosis-related genes. We selected TCGA-OV as the role of training set to construct the model, with both GEO datasets chosen for model validation. The Least Absolute Shrinkage and Selection Operator (LASSO) was utilized to obtain parameter λ that is optimal for further filtering the prognostic genes obtained. We conducted an analysis of multivariate Cox regression to determine the coefficient of each gene in the model and presented them in a lollipop plot to build the prognostic signature. The risk score was derived through the aggregation of products obtained by multiplying individual gene expression levels with their corresponding coefficients. Risk score = ∑ i = 1 n [ E x p g e n e i * β i ] Here, the E x p g e n e i designates the quantitative value of model gene expression, whereas β i signifies associated coefficient, specific to respective model gene. Furthermore, we divided the two groups(low risk and high risk) according to the values of median scores in each dataset. Survival analysis was independently conducted on two sets to examine survival disparities between the two risk cohorts. Additionally, an evaluation of the model's predictive capability was undertaken via ROC curve analysis, spanning across the time points of 1, 3, and 5 years. Typically, an AUC value greater than 0.6 is considered indicative of a model with acceptable discriminatory ability, while values closer to 1.0 suggest excellent performance. To ensure transparency and reproducibility of our work, we specify that an AUC threshold of 0.7 or higher is used to indicate a model with acceptable discriminatory ability in our study. We conducted separate analyses on the correlations between clinical information (age, stage, race) and risk score and visualization. In order to explore the relationship between the risk score and the model genes, a Pearson correlation analysis was employed. In addition, for the TCGA-OV dataset, we conducted the analysis of status in TMB, immune checkpoint genes, and TNF family genes between different risk groups. Finally, we utilized built-in TME (Tumor Microenvironment) analysis algorithms (EPIC, MCP-Counter, quanTIseq, TIMER, and xCell) from the "IOBR" package ( https://github.com/IOBR/IOBR ) to analyze the infiltration abundance status of immune cell between different risk groups. To ensure the veracity and robustness of downstream analyses, we implemented the Seurat framework for single-cell sequencing data, integrating rigorous quality control(QC) measures and data purification. These QC thresholds were: nFeature_RNA < 9000 and percent.mt < 25. We utilized the harmony method for batch integration of the data across multiple samples. Through the application of the Uniform Manifold Approximation and Projection (UMAP) technique, a dimensionality reduction process was executed on the consolidated data stemming from three ovarian cancer-focused single-cell sequencing datasets, ultimately yielding sixteen distinct cell clusters. Annotation of nine major cell subtypes was performed based on specific markers provided by the TISCH database, followed by visualization and heatmap demonstration of the expression differences of marker genes within each cell subtype. The AddModuleScore function within the Seurat package were utilized to compute the TNF-related signature for each cell and visualized it on the UMAP plot. Malignant cells were further extracted for dimensionality reduction clustering, identifying six malignant-related cell subtypes at a resolution of 0.2. Following this, we depicted the expression patterns of chosen model genes across various cell subtypes utilizing UMAP. Furthermore, employing the monocle2 package, we constructed developmental trajectories for each malignant subtype and visualized the expression changes of selected model genes along these trajectories over pseudo-time. We utilized the "limma" package to identify differential genes between two risk groups (based on the criteria |logFC| > 1 & adj.P.val. < 0.01). The fold change is calculated by dividing the expression level of a gene in one condition by its expression level in another condition. Taking the logarithm of this ratio, usually to the base 2, allows for easier interpretation and comparison of changes across multiple genes, as it converts large ratios into more manageable numbers. Subsequently, we performed Gene Set Enrichment Analysis (GSEA) using Kyoto Encyclopedia of Genes and Genomes(KEGG) and list of cancer hallmark gene obtained from MSigDB ( https://www.gsea-msigdb.org/gsea/msigdb/ ), enriching and visualizing significantly enriched pathways in two risk groups. We generated Mantel plots to illustrate correlation between the values of risk scores and both anti-cancer immunity cycle and therapeutic-related pathways. The relevant gene sets were sourced from literature [ 12 ]. Additionally, leveraging the "oncopredict" tool, we conducted drug sensitivity analysis using six selected drugs and visualized the drug sensitivity status between different risk groups. All the ovarian cancer tissues and adjacent normal tissues were obtained from Suzhou Hospital, Affiliated Hospital of Meddical School, Nanjing University. Patients included in the study underwent surgery as their primary treatment, while those who received radiotherapy, chemotherapy, or immunosuppression were excluded. During surgical procedures, resected tissue specimens were acquired from a cohort of 20 patients, all of whom provided their informed consent prior to the collection. Our study design was approved by hospital's ethics committee. Caov-3, IOSE-80, and HEY-A8 cells were obtained from the American Type Culture Collection (USA), and SK-OV-3 cell line were sourced from Wuhan Procell Life Science & Technology Co., Ltd(China). IOSE-80, HEY-A8, and Caov-3 cells were cultured in Dulbecco's Modified Eagle Medium (DMEM; Gibco, USA) added with fetal bovine serum(10 %), whereas SK-OV-3 cells were cultured in DMEM with high glucose (Gibco, USA). Four types of cells were cultured in a 37°C incubator with CO 2 (5 %) for the experiments. Small interfering RNA (siRNA) targeting PDE4B was used, with sequences purchased from GenePharma (China), and transfection procedure was performed using LipoFiter 3.0 (HANBIO, China). The specific siRNA sequences were: si-NC: UUCUCCGAACGUGUCACGUTT; siPDE4B-1: TTGGAATTGTATCGGCAATC; siPDE4B-2: TCCTAAAGACATTCAGAAT. In this study, the Trizol reagent (Takara, Japan) was utilized to extract total RNA of cells, and cDNA was reverse transcripted and synthesized by the cDNA synthesis reagent kit (Takara, Japan). RT-PCR was conducted using a mixture of cDNA and RT-PCR SYBR Green (Takara, Japan). The primers for cDNA amplification were: GAPDH-Forward: 5′-GTCAAGGCTGAGAACGGGAA-3′; GAPDH-Reverse: 5′-AAATGAGCCCCAGCCTTCTC-3′; PDE4B-Forward: 5′-AACGCTGGAGGAATTAGACTGG-3′; PDE4B-Reverse: 5′-GCTCCGGTTCAGCATTCT-3′. Cell lysates were generated by incubating cells in radioimmunoprecipitation assay (RIPA) buffer for 30 min. We quantified protein concentrations using BCA kit (Vazyme Biotech, China). For sample preparation, 5 × loading buffer (Beyotime, China) was mixed with the protein at a 1:4 ratio and heated at 95°C for 10 min. The protein was loaded onto 10 % SDS-PAGE gels and transferred to PVDF membranes. Following transfer, we soaked membranes for 2 hours in skim milk(5 %) which was dissolved in Tris-buffered saline with 0.1 % Tween(TBST). Then, all the membranes were exposed to primary antibodies and left for 10 hours at 4°C for incubation. The primary antibodies used in our study were: GAPDH (1:4000; Proteintech, USA), PDE4B (1:4000; Abcam, USA), caspase-1 (1:1000; CST, USA), cleaved caspase-1 (1:1000; CST, USA), caspase-3 (1:1000; Abcam, USA), caspase-7 (1:1000; CST, USA), caspase-8 (1:1000; CST, USA), cleaved caspase-3 (1:1000; Abcam, USA), cleaved caspase-7 (1:1000; CST, USA), cleaved caspase-8 (1:1000; CST, USA), GSDMD (1:4000; Proteintech, USA), GSDME (1:5000; Proteintech, USA), pMLKL (1:4000; Proteintech, USA), and MLKL (1:10,000; Proteintech, USA). Subsequent to the primary antibody incubation, the membranes underwent a thorough washing process with TBST solution, followed by a two-hour incubation period with the designated secondary antibodies. Blots were visualized using the SuperFemto ECL Chemiluminescence kit (Vazyme Biotech Co., Ltd, China) and the Amersham Imager 600 (GE Healthcare). We employed BeyoClick EdU-488 cell proliferation kit (Beyotime, China) to evaluate the impact of PDE4B on cellular proliferation. We seeded cells into 12-well plates at a density of 20,000 cells/well and cultured overnight at 37 °C. According to the protocol of manufacturer, we treated cells with medium containing EdU(10 µM) and incubated for 2 hours. Subsequently, they were fixed with 4 % paraformaldehyde (PFA) at 23 °C for 15 min. Then, cells underwent a 30-minute treatment with the Click reaction solution, conducted within a dark environment to ensure optimal conditions for the reaction. We captured photos using a fluorescence microscope (Nikon Corporation, Japan). Two types of cancer cells (5 × 10 3 /well) were cultured in a plate with 96 wells overnight. Then, all the medium was then removed, and 90 µL of medium mixed with CCK-8 reagent(15 µL) was transfered to each group. We used spectrophotometer to measure the cell viability after 2 hours of incubation. All the cancer cells were transfered in plates with 6 wells and they reached nearly 100 % confluence. Scratching the layer of cells with 200 µL pipette tip to generate wounds. The cells which floated in the liquid were removed by PBS washing, then cells which still attach to plate were cultured in serum-free medium for 48 hours. Data analysis was conducted utilizing Image J software. We conducted the evaluation procedure on the ability of cancer cells migration using Transwell chambers equipped with membranes (8.0-μm pore, Procell Life Science & Technology, China). The plate with 24 wells was filled with complete medium(500 μL), while cells at a density of 1 × 10 4 cells/chamber were seeded into chambers with serum-free medium. One day after the incubation, successfully migrated cells were underwent the procedure of fixing with 4 % PFA and dyed with crystal violet(0.1 %). We counted the stained cells under a microscope. In this study, our statistical analyses were comprehensively executed utilizing the R programming environment (version 4.1.3). Analysis related to the Single-cell sequencing data was conducted by R software package "Seurat" and "SCP pipeline." The "IOBR" package was employed for immune infiltration analysis. We conducted drug sensitivity analysis using "oncoPredict". We performed enrichment analysis using the "clusterProfiler" package. Univariate Cox analysis was conducted using "ezcox." Differential expression analysis was accomplished by "limma" package. A statistical significance threshold of p<0.05 was established, marking the cut-off for inferring meaningful differences in the data (** p-value < 0.01 denotes strong statistical significance; *** p-value < 0.001 signifies highly significant statistical difference; **** p-value < 0.0001 demonstrates extremely significant statistical significance.).

This study utilized multiple public databases to analyze bulk transcriptomic and single-cell data of ovarian cancer and normal tissues. Through calculating GSVA scores of TNF family genes, screening hub genes via WGCNA, and analyzing differentially expressed genes, genes with significant regulatory roles in ovarian cancer were identified. A 14-gene prognostic model constructed showed favorable performance in both training and validation sets, providing a reliable basis for clinical assessment. Additionally, single-cell analysis revealed cellular subpopulations of ovarian cancer and their expression differences throughout pseudotime. Exploration of biological functions and drug sensitivity analysis provided insights for the formulation of new therapeutic approaches. In summary, our study contributes insights for improving prognosis, guiding treatment, and promoting drug development for ovarian cancer patients. We further emphasize the role of PDE4B, which has been shown to play a significant role in the progression and treatment response of ovarian cancer. By elucidating the mechanisms through which PDE4B influences the disease, our findings underscore its potential as a therapeutic target and a prognostic biomarker. Future research should focus on the development of PDE4B inhibitors and their clinical application to enhance patient outcomes.

Among gynecological malignancies, ovarian cancer ranks third in terms of mortality, trailing behind cervical and endometrial cancer, and is among the most prevalent of such cancers. The prognosis of this disease is characterized by its poor outlook, marked by a high recurrence rate post-standard first-line treatment, and an approximate 5-year survival rate of 45 %. Most patients also develop chemoresistance, and advanced ovarian cancer exhibits high invasiveness. Hence, delving into the underlying mechanisms of ovarian cancer-associated genes, particularly those that contribute to its pronounced metastatic potential and recurrent nature, alongside the identification of pivotal biomarkers and crucial target genes within the biological processes, holds paramount significance for advancing the diagnosis, therapeutic strategies, and prognostic outcomes of ovarian cancer. We calculated the enrichment scores of TNF family genes and obtained 91 differentially expressed hub TNF family genes through WGCNA and differential expression analysis. Then, we conducted analysis of enrichment for BP, CC, MF. We selected TCGA-OV dataset as the training set to participant in the formulation of signarture, and two GEO datasets as the validation set. LASSO and multivariable Cox regression analysis were utilized to perform establishment procedure of the prognostic model, from which 14 genes were screened and identified, including GFPT2, PDE4B, PODNL1, TGFBI, CSF1R, PTGIS, SFRP2, COL5A2, TRAC, SLAMF7, VCAN, GBP1P1, C2, and TRBV28. Studying the hexosamine biosynthetic pathway, particularly focusing on GFPT2, a pivotal enzyme within it, offers insights into the biology process of colorectal carcinoma cells. This investigation entails augmentation of the p65 glycosylation and process of the NF-κB pathway activation [ 13 ]. Downregulation of GFPT2 expression can reduce the nuclear localization of β-catenin as well as the expression levels of Slug and Zeb1 in ovarian cancer cells [ 14 ]. Research findings highlight GFPT2 as a high-risk gene in ovarian cancer, correlating with an unfavorable prognosis [ 15 ]. PDE4B, as a constituent of the phosphodiesterase family, participates in the degradation of cyclic nucleotides like cGMP and cAMP. In addition to its role in regulating intracellular signaling pathways, PDE4B has been implicated in various cancer types. Specifically, PDE4B can influence cancer cell proliferation, apoptosis, and migration. For instance, in breast cancer cells, PDE4B has been shown to promote cell survival and growth by modulating the levels of cAMP, which in turn affects the activity of protein kinase A (PKA) and downstream signaling pathways. Furthermore, PDE4B expression has been associated with resistance to chemotherapeutic agents, suggesting that it may play a role in cancer cell resistance mechanisms. Understanding the specific contributions of PDE4B to cancer progression could lead to the development of targeted therapies aimed at modulating its activity, consequently attenuating the transduction for signal mediated by these pivotal second messengers [ 16 ]. Research indicates the upregulation of PDE4B in malignant tumors, particularly in the context of bladder cancer, presenting it as a promising therapeutic target. Its expression correlates with aggressive clinicopathological features and an unfavorable prognosis [ 17 ]. PODNL1, an extracellular protein, is expressed in various tissues including the tibial nerve, coronary arteries, and bone marrow mesenchymal stem cells, contributing to the formation of the extracellular matrix. As a member of the small leucine-rich proteoglycan (SLRP) family, which encompasses 17 genes, it falls under Class V SLRP [ 18 ]. Previous studies suggest that PODNL1 is a prognostic biomarker for ovarian cancer and gliomas [ 19 ]. TGFBI gene encodes an RGD-containing protein that can bind to collagen types I, II, and IV and is a critical hub for HRG. Confirmed as a HIF-2α responsive gene, it facilitates drug resistance of cisplatin in ovarian cancer through PI3K/Akt pathway activation [ 20 ]. TGFBI expression is closely associated with ovarian cancer macrophages, contributing to the immune-suppressive microenvironment of ovarian cancer [ 21 ]. Belonging to the III receptor tyrosine kinase family, CSF1R undergoes homodimerization upon binding with either CSF1 or its newly identified ligand, IL-34, initiating downstream receptor signaling activation [ 22 ]. In pancreatic cancer models, inhibition of CSF1/CSF1R can reprogram tumor-infiltrating macrophages, enhancing the efficacy of T cell checkpoint immunotherapy [ 23 ]. Cytochrome P450 superfamily member PTGIS participates in two pivotal PGI signaling pathways: one involving direct activation through cell surface engagement with PGI receptor PTGIR, and the other mediated by nuclear activation of peroxisome proliferator-activated receptor β/δ (PPARβ/δ) [ 24 ]. PTGIS, identified as an anti-metastatic agent, suppresses ovarian cancer cell invasion in vitro by downregulating MMP2/MMP9. Its expression correlates with diminished survival rates among breast and ovarian cancer patients [ 25 ]. SFRP2, belonging to the secreted frizzled-related protein (SFRP) family, functions as a canonical regulatory factor modulating the intricate dynamics of the WNT signaling cascade [ 26 ]. SFRP2 serves to impede tumor metastasis and holds promise as a candidate biomarker for the early detection of ovarian cancer [ 27 ]. COL5A2 belongs to the collagen type V family and is located at 2q32.2. This gene encodes a low-abundance fibrillar collagen α chain [ 28 ]. Studies have shown strong correlation and precise predictability of COL5A2 in patients with renal metastases from gastric cancer [ 29 ]. TRAC expresses the α and β chains of T cell receptors [ 30 ]. TRAC is associated with the prognosis of endometrial cancer [ 31 ]. SLAMF7, as a SLAM family receptors member, performs a crucial function in immune activation and suppression [ 32 ]. SLAMF7 exhibits dual roles in promoting and inhibiting tumors in human cancers. Evidence implicates SLAMF7 in tumor metastasis and invasion in multiple myeloma, clear cell renal cell carcinoma, and lymphoma, supporting its tumor-promoting effects. However, in lymph node metastatic breast cancer, high expression of SLAMF7 serves as a favorable prognostic indicator [ 33 ]. VCAN, a substantial proteoglycan constituent of the extracellular matrix, has been extensively investigated for its overexpression in tumor tissues. Its upregulation has been implicated in facilitating various tumor-promoting processes, including migration, invasion, proliferation, adhesion, and angiogenesis [ 34 ]. In advanced serous ovarian cancer, elevated VCAN expression within the stroma of tumor is related with poorer prognosis, while in vitro investigations indicate a promotional role of VCAN in ovarian cancer cell invasion [ 35 ]. GBP1P1 is a pseudogene. Reports indicate overexpression of the pseudogene GBP1P1 in endometriosis, while underexpression is found in HCC [ 36 ]. Compared to normal tissues, GBP1P1 expression is significantly higher in tumors [ 37 ]. Component C2, a serum glycoprotein, operates within the classical pathway of the complement system. C2 is part of both the CP (classical pathway) and LP (lectin pathway), both of which are associated with diseases driven by autoantibodies or ischemia-reperfusion [ 38 ]. TRBV28 is expected to participate in cell surface receptor signaling pathways and is anticipated to be part of the T cell receptor complex. The V region within the variable domain of the T cell receptor (TR) β chain, crucial for antigen recognition, is referenced [ 39 ]. We divided the risk groups into high and low based on the median score of each dataset. The results of KM survival analysis indicate poorer prognosis in the high-risk group. The results of ROC demonstrate good diagnostic efficacy of our model at 1, 3, and 5 years. We analyzed the correlation between clinical information and risk score, finding that age lacks statistical significance, while among stage categories, Stage I exhibits the highest risk score. In terms of race, individuals of American Indian or Alaska Native descent have the highest risk score, while Asians have the lowest. Upon conducting a Pearson correlation analysis, we observe a pronounced positive linkage between the assigned risk score and the expression patterns of the model genes, underscoring their concerted role in predicting risk. The high-risk group demonstrates elevated expression levels of TMB, TNF family genes, and immune checkpoint genes. Evaluation of immune cell infiltration utilizing five algorithms reveals no significant variance in the abundance status of immune cell between different risk groups. Using UMAP, we reduced the dimensionality of integrated single-cell sequencing data from three ovarian cancer datasets to identify 16 cell clusters. Nine major cell subtypes were annotated based on TISCH-specific markers, and differential expression of marker genes for each cell subtype was analyzed. We calculated TNF-related signatures for each cell and visualized them on the UMAP plot. We further clustered malignant cells to identify six malignant-related cell subtypes, with higher expression levels of model genes VCAN and C2 in malignant cells, and lower expression of SLAMF7 and CSF1R. We constructed developmental trajectories for each malignant subtype, showing an increase in C2 expression with pseudotime and a decrease in VCAN expression. Differentially expressed genes were identified among the two risk groups, followed by GSEA enrichment analysis conducted on pathways enriched within each group. Cancer markers were significantly enriched in the high-risk group. A Mantel plot displays a high degree of mutual correlation between risk score and anti-cancer immunity cycle as well as therapeutic-related pathways. Drug sensitivity analysis suggests most drugs have better efficacy in the low-risk group, but Acetalax shows better therapeutic effects in the high-risk group. Finally, we confirmed the expression of PDE4B and its role in ovarian cancer through in vitro wet experiments. PDE4B exhibited a higher expression level in ovarian cancer, promoting cancer cell proliferation, invasion and migration, and inhibiting PANoptosis. The study's classification of risk groups into high and low based on the median score may not fully capture the complexity and variability of patient risk. Different datasets may have different risk score distributions, and the median may not always be the most appropriate threshold for dividing groups. Although the study conducted a comprehensive analysis of various factors, including clinical information, race, gene expression patterns, and immune cell infiltration, the sample size and diversity of the datasets used may limit the generalizability of the results to a broader population. Additionally, the study's findings on drug sensitivity and the role of PDE4B in ovarian cancer require further validation in larger and more diverse patient cohorts.

Ovarian cancer ranks seventh among malignant tumors and eighth as a cause of cancer-related deaths in women globally. As a prevalent gynecologic malignancy, ovarian cancer exhibits the third highest fatality rate, trailing behind cervical and uterine cancers [ 1 ]. The scientific literature predominantly acknowledges the existence of three fundamental categories of ovarian cancer: epithelial carcinoma, germ cell tumor, and sex cord-stromal tumor, with epithelial carcinoma being the most prevalent, while the latter two collectively represent only about 5 % of all ovarian cancers. Four distinct histological subtypes constitute the primary classification of epithelial ovarian cancer: mucinous, serous, endometrioid, and clear cell carcinoma. Within the serous category, tumors are subdivided into high-grade serous carcinoma (HGSC) or low-grade serous carcinoma (LGSC) [ 2 ]. Globally, ovarian cancer affects 239,000 patients annually, resulting in 152,000 deaths [ 3 ]. The Centers for Disease Control and Prevention (CDC) reports a prevalence peak among white females in terms of the incidence rate of the condition, with 11.3 cases per 100,000 individuals. The next highest incidence rates by race after white women are observed among Hispanic women [ 2 ]. Ovarian cancer is rare among young women but shows a sharp increase in incidence after the age of 50. Factors influencing ovarian cancer development include genetic, environmental, and lifestyle factors [ 4 ]. Additionally, factors such as pregnancy, breastfeeding, and oral contraceptive use play roles in reducing the hazard of this disease. The heightened prevalence of oral contraceptives usage in recent times has been accompanied by a discernible decline in ovarian cancer incidence [ 5 ]. The current standard of care for first-line treatment involves administering a platinum-taxane chemotherapy regimen following debulking surgery. After first-line treatment, cancer recurrence occurs in 60 %-70 % of optimally debulked patients (1 cm residual disease), with a 5-year survival rate of approximately 45 %. While maintenance regimens utilizing bevacizumab or PARP inhibitors have demonstrated an extension in progression-free survival (PFS), they fail to elicit a commensurate prolongation in overall survival (OS), underscoring the imperative for the development of more efficacious maintenance therapeutic strategies [ 6 ]. Predictions suggest a substantial increase in the mortality rate of this cancer by 2040 [ 7 ]. Most patients experience recurrence and develop chemotherapy resistance, and advanced ovarian cancer exhibits aggressive behavior and rapid dissemination. Hence, it is essential to investigate the underlying mechanisms that orchestrate the initiation and progression of ovarian cancer, particularly those associated with its elevated rates of metastasis and recurrence. This exploration aims to pinpoint vital biomarkers and elucidate pivotal target genes crucial for diagnosing, treating, and prognosticating ovarian cancer. Numerous studies have investigated the mechanisms of action of ovarian cancer-related biomarkers and constructed prognostic models [ 8 , 9 ]. Previous research has indicated the potential application of circulating free DNA methylation patterns in ovarian cancer detection and prognostic assessment [ 10 ]. However, the extent to which TNF family genes, as immune-focused therapeutic targets, contribute to therapeutic efficacy and prognostic insights in ovarian cancer remains an area of unexplored scientific territory. We analyzed bulk transcriptomic data of ovarian cancer and normal ovarian tissues as well as single-cell sequencing data achieved from publicly available databases. We calculated the GSVA scores of TNF family genes, performed WGCNA, identified hub genes, and conducted GO enrichment analysis. Additionally, we performed differential gene expression analysis and intersected the obtained tumor tissue highly expressed genes and hub genes for GO enrichment analysis in BP, CC, and MF. We used LASSO and multivariate Cox regression analysis to construct a prognostic model, comprising 14 genes (GFPT2, PDE4B, PODNL1, TGFBI, CSF1R, PTGIS, SFRP2, COL5A2, TRAC, SLAMF7, VCAN, GBP1P1, C2, TRBV28). Both the training and validation sets exhibited good discriminatory performance and prognostic assessment capability. The model also underwent clinical information analysis and evaluation of immune cell infiltration. In single-cell sequencing analysis, we dimensionally reduced and annotated cell subtypes, analyzed the expression of relevant genes in each subtype, and explored expression differences over pseudo-time. Biological functional analysis and drug sensitivity analysis were also performed. Our study provides new perspectives and a robust theoretical framework for the treatment, prognosis, and drug development of ovarian cancer patients.

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

The datasets analyzed for this study can be found in the GEO website ( https://www.ncbi.nlm.nih.gov/geo/ ), TCGA website( https://portal.gdc.cancer.gov/ ) and MSigDB database ( https://www.gsea-msigdb.org/gsea/index.jsp ).