A Pyroptosis-Related LncRNA signature in Prognosis and Immune Microenvironment infiltration Prediction of Ovarian Cancer

Simple Summary To evaluate the role of pyroptosis-related lncRNAs (PRLs) in clinical outcome and immue microenviroment infiltration in ovarian cancer, 32 prognostic PRLs were screened and a 7-PRLs related signature with prediction accuracy of overall survival (OS) and immune microenvironment infiltration for ovarian cancer were developed based on expression profiles in TCGA dataset. Patients with high-risk score showed worse outcome, lower infiltration abundance and immunotherapy response. Additionally, lncRNA TYMSOS was proved to serve as a novel tumor biomarker. It might promote cell proliferation, invasion and migration via inhibiting the process of pyroptosis by regulating the expression of GPX4. Our findings provided a Abstract LncRNA and pyroptosis play important roles in cancer development and tumor immune microenviroment. However, pyroptosis-related lncRNAs (PRLs) in ovarian cancer have not been identified and its impact on prognosis and immune response are not fully understood. In TCGA-RNA-seq cohort (n=377), 32 prognostic PRLs were screened by using pearson correlation analysis and univariate cox analysis. Subsequently, a 7-PRLs related signature with high prediction accuracy of overall survival (OS) for ovarian cancer were developed. Multivariate Cox regression analysis suggested that it might serve as an independent prognostic indicator and the risk score showed significantly associated with prognosis. Except to clinical outcome, the signature was also associated with tumor immune microenvironment infiltration. Patients with high risk score exhibited lower infiltration abundance of MHC class Ⅰ cells, Type Ⅰ IFN response and immunotherapy response. Furthermore, it was revealed that TYMSOS was highly expressed and its high expression was associated with worse OS in ovarian cancer. TYMSOS deletion in ovarian cancer cell lines inhibited the cell proliferation, invasion and migration. TYMSOS might promote tumor progression via inhibiting the process of pyroptosis by regulating the expression of GPX4. In conclusion, the prognostic PRLs-related signature constructed in our study can be served as efficient biomarkers for prognostic prediction and immune microenvironment infiltration in ovarian cancer.


Introduction
Ovarian cancer (OC) is the third most malignant cancer of all female reproductive systems and patients with OC suffer from an extremely high recurrence and mortality rate [1]. In 2020, there were 21,750 new cases and 13,940 deaths of ovarian cancer in the United States, and it was the fifth most prevalent cause of death among all women who died of cancer [2]. Due to the absence of obvious early physical signs, the median age of OC patients at diagnosis is 63 years, and over 70% of EOC cases are reported to be newly diagnosed at terminal stages with a five-year survival rate of approximately 48% [3]. The standard of therapy for OC is tumor cytoreductive surgery in conjunction with platinum-based chemotherapy. However, the vast majority of patients recur within two years. It indicates that recurrence and drug resistance are the major challenges that need to be addressed [4]. Given the limitations of current OC treatment, new targets for therapy are desired to enhance the clinical outcome of OC. In this light, there is an overwhelming urgency for robust new prognostic models to render targeted therapies more plausible.
At present, multiple evidence have demonstrated that pyroptosis played importanted roles in cancer [5]. Pyroptosis, which is morphologically characterized by swollen cell lysis, rupture of cell membranes and the release of cell contents, is a programmed cell death referring to the Gasdermin family-induced and caused by the inflammasome, and ultimately activating a cascade of enlarged inflammatory responses [6][7][8]. As is know to us, there are two main forms of pyroptosis: i) the caspase-1 dependent classical pathway; ii) the caspase4/5/11 reliant non-classical pathway. In the non-classical inflammasome pathway, bacterial lipopolysaccharide (LPS) identifies and initiates caspase-4/5/11 cleavage of gasdermin D (GSDMD) and induces pyroptosis. In a typical inflammasome pathway, the inflammasome recruits and combines with ASC (CARD-containing apoptosis-associated spot-like protein), resulting in ASC gathering, which in turn recruits procaspase-1 and activates caspase-1. caspase-1 is involved in the cleavage and maturation of pro-IL-18/1b and the cleavage of GSDMD. The C-terminal fragment of GSDMD (GSDMD-CT) remains in the cytoplasm. Meanwhile, the released N-terminal fragment of GSDMD (GSDMD-NT) is released, leading to building a porosity in the plasma membrane. Subsequently, the excretion of IL-18/1b causes the influx of water, cell swelling and permeability lysis. [9][10][11]. Recent findings demonstrated that the pyroptosis-related gasdermins playes an emerging role and might serve as potential new therapeutic targets in various diesases, including autoimmune and inflammatory diseases, infectious diseases, deafness, and cancer [12][13][14][15][16]. In cancer, it is suggested that the process of pyroptosis inhibits tumorigenesis and progression, as well as serves as a pro-inflammatory signal, and establishes a microenvironment suitable for tumor cell growth [17][18][19][20]. Yet, the exact function of the pyroptosis was poorly investigated in OC.
Long non-coding RNAs (lncRNAs) are identified as a group of RNA molecules over 200bp in length that are not translated into proteins. They are reported to participate in various biological processes, such as epigenetic modifications, inheritance stamping, chromatin organization, and protein amendment [21]. Several studies have indicated that lncRNAs might be engaged in ovarian cancer and pyroptosis [22][23][24][25]. While, previous studies have primarily focused on the utility of protein-coding genes of pyroptosis, lncRNAs associated with pyroptosis have barely been reported in OC [26][27][28]. Therefore, the identification of pyroptosis-related lncRNAs (PRLs) is essential for deciphering the underlying motifs of pyroptosis in OC and investigating new therapeutic targets.
In our study, PRLs were firstly screened by using the pearson correlation analysis. Furthermore, a prognostic signature was constructed based on the PRLs. The prognostic signature significantly predicted the clinical outcome of OC patients in high-risk group and low-risk group with a high diagnostic accuracy. In addition, it was correlated with the immune microenvironment infiltration and immunotherapy response. Finally, we validated the effects of TYMSOS on cell proliferation, invasion and migration in OC cell lines.

Clinical characteristics of the study patients
The mRNA and LncRNA expression profiles of 377 OC subjects were screened in the TCGA-RNA-seq dataset. All patients were randomly separated into training cohort (n=189) and test cohort (n=188) cohorts and there was no difference between the two cohorts in age, Stage, Grade, tumor residual, lymphatic invasion, subdivision and chemotherapy (Table1). The flow chart of this study was presented in Figure 1A.

Identification of prognostic PRLs in OC patients
Firstly, we compared the expression of the 33 PRGs between OC tissuses and normal ovarian tissues using the TCGA datasets and GTEx datasets. The results suggested that 31 among 33 PRGs were either upregulated or downregulated in OC, compared to normal ovarian tissues ( Figure 1B). Then, we identified 14826 lncRNAs in the TCGA-RNA-seq dataset, based on the lncRNA annotation file from the GENCODE website.
To extract the potential PRLs, the Pearson correlation analysis was performed to assess the correlation between PRGs and lncRNAs. Subsequently, we obtained 2792 PRLs in the TCGA-RNA-seq dataset (Supplement Table 2). To assess the prognostic value of PRLs, we initially used univariate Cox regression analysis and identified 32 LncRNAs with a significant value(p<0.05). The forest plot demonstrated that 32 RPLs were notably associated with the prognosis of OC patients ( Figure 1C). Besides, the Cytoscape software (3.8.2) was applied to briefly display the interactive relationships between prognosis PRLs and PRGs ( Figure 1D).

Construction and validation of the prognostic PRLs signature
To establish an optimal prognostic signature for predicting clinical outcome in OC patients, the LASSO-Cox regression analysis was performed to screen out the most robust model from the candidated lncRNAs (Figure 2A-B). In summary, we constructed a prgnostic PRLs signature and the risk score for each patient was assigned based on the coefficients of each lncRNA in the TCGA cohort ( Figure 2C). The following is the were divided into high-risk group and low-risk group depending on the best cut-off of risk score. Survival curves were obtained using the Kaplan-Meier analysis and the log-rank test was conducted to compare OS between the two groups. Kaplan-Meier survival curves exhibited that the OS of patients with lower risk was substantially longer than the high-risk group both in the training cohort and validation cohorts (test and sum cohort) ( Figure 2D

Stratification analysis of the prognostic PRLs signature
To better evaluate the prediction ability of the prognostic PRLs signature, the stratification analysis was performed. Compared to patients with high risk, patients with low risk had better OS in patients aged ≤55 and aged >55 subgroup ( Figure Figure 4E, p=0.078). However, high risk patients in both subgroups had the tendency for worse OS, in contrst to low risk patients. Altogether, these results suggest that the prognostic PRLs retained its ability to predict OS in varous subgroups and it could be served as a potential predictor for OC patients.

Modeling the prognostic nomogram
In addition, whether the prognostic PRLs signature was an independent prognostic facter was assessed by using univariate and multivariate Cox regression analysis. As shown in the forest plot, the red color presented the risk factors (HR>1) and the blue one indicated the protective factors (HR<1). The results depicted that risk score (HR:  Figure 5B). Thus, Age, FIGO stage and risk score were applied in the construction of a nomogram model ( Figure 5C).
Calibration plots indicated that the actual vs predicted rates of 3-and 5-year OS showed perfect concordance ( Figure 5D). The diagram verified that the nomogram has a reliable and robust ability to predictive the prognosis for OC patients.

Identify the biological function of prognostic PRLs signature in OC
Gene Set Enrichment Analysis (GSEA) was employed to find the key pathways and biological functions that different in the high-and low-risk groups. Firstly, DEGs were screened and |logFC|>2, p<0.05 was identified as statistical significance ( Figure 6A).
GSEA results suggested that the DEGs were mainly enriched in the inflammatory response pathway, p53 pathway, TGF-β signaling and TNF α signaling via NF-kB and so on( Figure 6B). Then, KEGG analysis and GO analysis were conducted and the outcomes displayed that the DEGs were mainly enriched in cell adhesion molecules, MAPK signaling pathway, NF-kB signaling pathway, PI3K-AKT signaling pathway, Wnt signaling pathway, Primary immunodeficiency, and plentiful immune-related biological process ( Figure 6C-D). Thus, these results suggested the signature may be involved in the above pathways and influence the survival of OC.

Association between PRLs signature and immune infiltration
Considering the aforementioned enrichment analysis and the stronger association between pyroptosis and immune status, we explored the relationships between risk score and immune microenvironment. To explain the immune cell and stromal cell infiltration situation, wecalculated the ESTIMATEscore, ImmuneScore, PurityScore and StromalScore. The correlation analysis implied that the risk score was positively relevant with the ESTIMATEscore (Figure 7A, r=0.14, p=0.006), and StromalScore (Figrue 7C, r=0.20, p=9.64e-05) but negatively correlated with the PurityScore ( Figure   7D, r=-0.14, p=0.005). whereas, there was no correlation between risk score and ImmuneScore ( Figure 7B, r=0.215, p=0.06). After that, the distribution proportion of different subpopulations of adaptive immunity cells and innate immunity cells in highrisk groups and low-risk groups were analyzed by using the ssGSEA method. The distribution proportion of MHC class Ⅰ cells and Type Ⅰ IFN response were significantly lower in the high-risk group than the low-risk group ( Figure 7E). In addition, the potential response of individual patients to immunotherapy was appraised by using the TIDE algorithm. The results showed that patients with low risk scores were more sensitive to immunotherapy than patients with high risk scores( Figure 8F), which might be associated with higher expression of PD-L1, CTAL4 and LAG3 in low risk group ( Figure 7G-I). Taken together, we speculated that patients with high risk score might recruit multiple immune cells and stromal cells and escape immune surveillance.

Comparison of chemotherapy response between patients with different risk scores
Thus far, the predominant therapeutics for OC patients were still chemotherapy and different molecular subtypes of OC should direct clinical decisions. Therefore, we contrasted the sensitivity to a variety of anticancer drugs between the high-and lowrisk groups to imply potential treatment modalities (supplement Table 3). The results displayed that the IC50s of veliparib and metformin, two food and drug administration (FDA)-approved drugs, were higher in patients with higher risk. which means the shrinkage risk occurs with the growing sensitivity to veliparib and metformin ( Figure   8). Unfortunately, there was no significantly difference in IC50 of cisplatin and paclitaxel between two groups (supplement Figure 1). All these evidence indicated that veliparib and metformin might be a good choice for OC patients with low risk score

Discussion
Ovarian cancer is one of the most common malignancies in the world and associated with a high mortality rate. Consequently, there is of great importance to identify reliable and effective biomarkers for the OC prognosis. In prior research, the lncRNA signatures for prognostic prediction have been valicated in many categories of cancers [29], and there is even a database Lnc2Cancer 3.0, which includes comprehensive data on experimentally supported long non-coding RNAs (lncRNAs) and circular RNAs (circRNAs) associated with human cancers [30]. In the same way, based on the differentially expressed lncRNAs and illness etiology, several lncRNA-associated signatures have correspondingly been industrialized to forecast the outcome of OC patients [31,32]. Nevertheless, the mutual interactions of PRLs in prognosis of OC patients remain unclear. In this study, we firstly reported a prognostic PRLs signature, providing a promising strategy for prognosis and immune response in OC patient, which have important clinical implications for guiding individual treatment and improving the effectiveness of the immune response.
Taking advantage of the TAGA and GTEx databases, we compared the expression of 33 PRGs in OC samples and normal ovary tissues by using the R package "limma".
Surprisingly, we discovered that all PRGs except two gene (CASP1 and CASP4) were differently expressed in tumor samples. Afterwards, we performed a correlation examination of the 31 PRGs with lncRNAs. Subsequently, 32 lncRNAs assocaited with OS were picked out and defined as prognostic PRLs. Utimately, several prognostic PRLs were identified to play important roles in cancer. For instance, Overexpression of TOPORS-AS1 was supposed to inhibit ovarian cancer cell proliferation and restrain cell aggressive behavior in vivo and in vitro [33]. LINC01281 was confirmed as immune-associate lncRNAs for Predicting Prognosis in Cervical Cancer [34]. In conclusion, it is still unknown whether these lncRNAs are involved in ovarian oncogenesis and development. Subsequently, we performed LASSO-Cox regression to construct a prognostic PRLs signature. The KM plotter analysis and AUC curves suggested that the prognostic PRLs signature might effectively predict the clinical outcome of OC patients. In further, stratification analysis indicated that the prognostic PRLs signature still retained its prognostic ability to predict OS for patients without considering other clinical features. Multivariate cox regression analysis showed that risk score was an independent risk factor for prognosis of OC patients. For better clinical applicability of the signature, we created a nomogram.
GSEA analysis was then carried out between the two groups. The outcomes demonstrated that the model had a robust association with cell proliferation and immunity. What's more, GSDME expression was reported to enhance both the number and functions of tumor-infiltrating natural killer and CD8+ T lymphocytes to increase the phagocytosis of tumor cells [35]. Natural killer cells and cytotoxic T lymphocytes kill gasdermin B (GSDMB)-positive cells through pyroptosis [36]. In addition, lncRNAs was also validated to participate in immune pathways through pyroptosis [37]. Therfore, we then analyzed the different distributions of immune function and immune cells between high-and low-risk group. The results indicated that the risk score was associated with immune microoenvironment and immune response. Various immune cells and immune reponse, especially MHC class Ⅰ cells and Type Ⅰ IFN response, were differently distributed in the high-risk group and low-risk group. All the evidence reveals that patients with high risk score might recruit multiple immune cells and stromal cells and escape immune surveillance. Also, we conducted a comparison of the sensitivity of 138 common anti-cancer drugs between the high-risk and low-risk groups.
It suggested that patients with high risk acquired a large number of drug-resistance in metformin and veliparib, which acted as the first-line Chemotherapy and Maintenance Therapy in OC [38].

Data acquisition and preprocessing
We obtained all datasets shown in this study including TCGA and GTEx datasets were publicly accessible and downloaded from the University of California, Santa Cruz TCGA data were randomly split into two equal cohorts: training cohort and test cohort, which was also applied to validate as a sum cohort. There was no different in pathological features and treatment. The PRGs (pyroptosis-related genes) were obtained from the previously published literature and MSigDB dataset (www.gseamsigdb.org/gsea/msigdb/) [22,39].

Identification of PRLs
The lncRNA annotation file was acquired from the GENCODE website for the annotation of lncRNAs. Consequently, 14826 lncRNAs were obtained from the TCGA-RNA-Seq cohort [40]. Pearson correlation analysis was used to screen PRLs.
Those lncRNAs with r>0.4 and p<0.001 were considered as the PRLs [41]. To determine the prognostic value of PRLs, we further conducted univariate Cox regression analysis by using the "survival" package, and the hazard ratios (HR) with 95% confidence intervals (CIs) were examined. p<0.05 indicated that pyroptosis-related lncRNAs were significantly associated with overall survival (OS) and considered as prognostic PRLs.
PRLs with HR>1 were considered to be risk factors, whereas HR<1 were considered to be protective factors.

Construction of prognostic PRLs signature
A risk signature was constructed by performing the LASSO-Cox regression on the prognostic-related lncRNAs using the "glmnet" package [42]. Through 1000 crossvalidation, a panel of genes and their LASSO coefficients were obtained. The risk score for the signature was calculated using the following formula: (n, is the number of the gene; β, LASSO coefficient; X, the expression of each prognosticPRLs in each sample). Based on the best cut-off risk score determined by the "surv_cutpoint" function of the "survminer" R package, patients were divided into high-risk and low-risk groups. Kaplan-Meier method with the long-rank test was performed to reveal the difference of OS between the high-risk and low-risk group by using the "survival" package. Besides, the time-dependent ROC curve and area under the curve (AUC) were applied to evaluate the prediction accuracy of the signature. All the time-dependent ROC curves were calculated by the "SurvivalROC" package and drawn by the "ggplot2" package.

Nomogram construction based on clinical features and risk score
Univariate and multivariate COX regression were performed to select the prognostic risk factors. The nomogram model was constructed using the "RMS" package to predict the 3, 5-year survival probability. The calibration curves were used to assess the concordance of the observed and predicted rates of 3, 5-year overall survival [43].

Estimation of tumor-infiltration, immunotherapy and chemotherapy response
Firstly, All microenvironment scores, including Estimate score, Immune score, Purity score, and stromal score were calculated by using the ESTIMATE algorithm (https://bioinformatics.mdanderson.org/public-software/estimate/) [44]. The infiltrating immune cells scores and the activity of immune-related pathways were calculated by performing the ssGSEA analysis with "gsva" package [45]. Tumor Immune Disfunction and Exclusion (TIDE) algorithm (http://tide.dfci.harvard.edu/), which is commonly utilized to accurately predict the outcome of patients treated with immune checkpoint blockade (ICB), were employed to evaluate the immunotherapy response [46]. The chemotherapy response of each patient was evaluated by using the Genomics of Drug Sensitivity in Cancer database (GDSC, https://www.cancerrxgene.org). The halfmaximal inhibitory concentration (IC50) of all drugs commonly used to treat tumors was calculated and represented the drug response. The R package 'pRRopheticRredic' was used with tenfold cross-validation and other parameters by default [47].

Function Enrichment Analysis
Gene set enrichment analysis (GSEA) was performed to identify the potential molecular mechanisms or potential functional pathways associated with prognostic PRLs signature. The TCGA samples were divided into a high-risk group and a low-risk group Gene Ontology (GO) enrichment analyses were performed to identify enriched pathways between the high-risk group and the low-risk group by using the "clusterProfiler" R package. |NES| >1 and false discovery rate (FDR) <0.05 were considered statistically significant

Cell Culture
The human ovarian cancer cell lines A2780 and SKOV3 were both cultured in RMPI-1640 medium, supplemented with 10% fetal bovine serum (FBS) and 1% penicillin/streptomycin. They were cultured in a sterile incubator maintained at 37°C with 5% CO2. The cells in logarithmic growth phase were collected for subsequent experiments.

RNA Extraction and Quantitative PCR
Total RNA was extracted using TRIZOL reagent according to the manufacturer's protocol (Invitrogen, 15596-026). Reverse transcription of cDNA was performed using the PrimeScript RT kit (Takara, RR047A, Japan). qPCR assays were performed using LightCycler480 detector (Roche, USA). The relative mRNA levels were calculated using the comparative Ct method with GAPDH as the reference gene. All the primers are listed in supplement Table 1.

Colony Forming Assay
The A2780 and SKOV3 cells were transfected with or without lncRNA-targeted siRNAs for about 48h. Then, 200 cells were plated in the six-well culture plates and cultured for about 2 weeks. The cellular colonies were counted by staining with crystal violet.

Cell Counting Kit-8 assay
The A2780 and SKOV3 cells were measured for cell survival after transfection with or without lncRNA-targeted siRNAs for 24h, 48h and 72h, respectively. Following the protocols of the manufacturer, cell viability was detected by using Cell Counting Kit-8 (Proteintech: PF00004).

Transwell migration and invasion
A total of 2*10 5 ovarian cells transfected as above were seeded in the upper chamber of Transwell plates with serum-free medium. The cells were then incubated for 48 hours.
To perform invasion experiments, the upper chambers were covered with a mixture of RPIM-1640 and Matrigel. Finally, the cells at the top of the lumen were removed with a cotton swab, whereas the cells across the membrane were stained with 0.5% crystal violet, observed and counted under 100x magnification.

Statistical analysis
Significant quantitative differences between groups were analyzed using the two-tailed Students' t-test, whereas differences among groups were analyzed using the one-way ANOVA. Kaplan-Meier curves and log-rank test were used to calculate the overall survival rate. All statistical analyses were performed using R sofware (version 4.0.2). * means p<0.05, ** means p<0.01, ***means p<0.001. p<0.05 was considered statistically significant.

Conclusions
In conclusion, our study constructed for the first time a novel prognostic signature based on lncRNAs associated with pyroptosis. It effectively predicted the clinical out of OC patients. Except to prognosis, the tumor immue microenvironment and immune response were significantly different in high-risk and low-risk group.

Data availability statement
Publicly available datasets were analyzed in this study. The datasets analyzed for this study were obtained from TCGA, GEO and GTEx databases.

Conflicts of Interest
The authors declare that there are no conflicts of interest.

Supplement materials
Supplement Table 1. A list of Primer sequence used in the study. Table 2. A list of all the PRLs screened by Pearson correlation analysis; Supplement Table 3. The association between risk score and drug sensitivity of all the drugs; Supplement Figure 1. The difference IC50 of cisplatin and paclitaxel between two groups.