Transcriptome analysis indenties KRT23 as a immunotherapeutic target in cervical cancer

Background: Cervical cancer (CC) is one of the most common malignancies in females worldwide. Traditional treatments have been used widely, but the prognosis remains poor. Therefore, new strategies are needed to improve outcomes. Immunotherapy has been used to treat various types of solid tumors. Subtypes of the tumor microenvironment (TME) are associated with the response to immunotherapy, so understanding complexity of TME is pivotal for immunotherapy. Methods: In this study, we used two methods, “ssGSEA” and “xCell”, to estimate the immune prole in CC. R packages “Corrplot” was used to analysis the correlation of immune cells. “ConsensusClusterPlus” was used to cluster CC based on inltration of immune cells. “Limma” was used to identify differentially expressed genes (DEGs), and “clusterprole” was used to perform enrichment analysis. “survival”and “glmnet”was used to construct perdition model. RT-PCR was used to detect Keratin, type I cytoskeletal 23 (KRT23) expression (cid:0) Cytometric bead arrays and ELISA were used to detect CCL5expression. Transwell assay was used to detect migration of CD8+T cells. Results: We subdivided CC into “hot” and “cold” tumors, in which hot tumors had inltration of more immune cells and longer survival. Enrichment analyses of DEGs revealed that the number of activated immune signaling pathways was higher in hot tumors. KRT23 showed high expression in cold tumors and its expression was negatively correlated with inltration of immune cells. In vitro experiments, knockdown of KRT23 expression promoted secretion of CCL5, and promoted recruitment of CD8 + T cells. We also constructed a model based on DEGs that had high ecacy for predicting the survival of CC and patients receiving immunotherapy. Conclusion: Our study provides deep insights in inltration of immune cells into CC. KRT23 may act as an immunotherapeutic target. Our model can predict the prognosis of CC patients and may guide immunotherapy.

In recent years, increasing numbers of studies have focused on the crucial role of immunotherapy in CC.
Notably, several studies have reported that cisplatin-based chemotherapy can increase PD-L1 expression in CC (15)(16)(17). Combinations of immune-checkpoint inhibitors with chemotherapy, radiotherapy, or other novel approaches may improve the results of CC treatment. Even though immunotherapy has achieved remarkable e cacy, accumulated data in recent years have demonstrated that many patients experience minimal or no clinical bene t if provided with identical treatment. This phenomenon could be attributed to the complexity and uniqueness of the tumor microenvironment(TME).
The TME is a complex, plastic, and dynamic system "sculpted" by tumor cells and surrounding cells (18,19).Cells from the innate immune system and adaptive immune system are important components of tumor stroma, can be reprogrammed according to the TME, and may be involved (positively or negatively) in the survival and progression of tumor cells (20,21). For example, tumor-associated macrophages (TAMs) are usually the largest population of myeloid cells in ltrating in most solid tumors (22). TAMs display a high degree of functional plasticity if exposed to various microenvironmental conditions, and can be classi ed as "M1-like" (pro-in ammatory and usually anti-tumor) or "M2-like" (anti-in ammatory and pro-tumor) (23,24). Accumulating evidence suggests the critical roles of the TME on promoting tumor progression. However, how the TME affects the e cacy of immunotherapy is not known. Immunotherapy harnesses or restores the immune system to kill tumor cells, but this requires the in ltration of immune cells to the tumor site. Studies have demonstrated that the type of the TME is associated with the clinical e cacy of immunotherapy. A "hot" tumor that has su cient tumor-in ltrating lymphocytes and antigenpresenting cells can respond robustly to immunotherapy. A"cold" tumor lacking immune cells, in general, cannot elicit an effective response to immunotherapy (25).Therefore, understanding and distinguishing the unique classes of the TME is useful for predicting and guiding immunotherapy.
We undertook a comprehensive analysis to explore the in ltration of immune cells in CC using two methods and constructed a prediction model. We observed that CC patients with in ltration of a higher number of immune cells survived longer. To uncover the underlying mechanisms of in ltration of immune cells, we subdivided tumors into hot and cold types and ascertained the differentially expressed genes (DEGs) between them. In addition, our model performed well in predicting the overall survival (OS) of patients with CC and the OS of patients receiving immunotherapy.

Materials And Methods
Ethics statement CC specimens were obtained after surgical treatment in the First A liated Hospital of Zhengzhou University. Specimens were frozen in the biobank of Zhengzhou University. All participants provided written informed consent for their specimens to be used in this study. The study protocol was approved by the ethics committee of First A liated Hospital of Zhengzhou University.

Cell culture
A human cervical cell line (HeLa) was purchased from the Institute of Biochemistry and Cell Biology of the Chinese Academy of Sciences (Shanghai, China). Cells were sustained in RPMI1640 medium with 5% fetal bovine serum and an atmosphere of 5% CO 2 in a humidi ed incubator at 37°C.

Acquisition and normalization of data
Level-2 mRNA-sequencing data (fragment per kilobase of transcript per million mapped reads) of CC were downloaded from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/) and transformed to transcripts per millionfor further analyses. The clinical data of CC were downloaded from Xena within University of California Santa Cruz (http://xena.ucsc.edu/). The GSE78220 dataset was downloaded from the Gene Expression Omnibus database (www.ncbi.nlm.nih.gov/geo/). A dataset of patients with metastatic urothelial cancer treated with anti-PDL1 agents was downloaded on online website supplied in the article(http://resea rch-pub.gene.com/IMvigor210CoreBiologies/)

Estimation of the immune pro le
The immune pro le (i.e., the number and type of immune cells)was estimated by the packages "ssGSEA" and "xCell" within R (R Institute for Statistical Computing, Vienna, Austria). For xCell, we selected samples with p <0.05 and we included onlyimmune cells for further analyses. The Immune Score, Stromal Score and tumor purity were calculated by the R package "ESTIMATE".

Identi cation and functional annotation of DEGs
Tumor samples were divided into two subtypes: "cold" and "hot". DEGs were calculated by the R package "Limma" and visualized by volcano plots constructed by the R package "ggplot2". DEGs with log fold change>1 and p<0.05 were selected for annotation using the Kyoto Encyclopedia of Genes and Genomes (KEGG; https://www.genome.jp/) and Gene Ontology (GO; http://geneontology.org/)databases by the R package "clusterpro le". A protein-protein interaction (PPI) network was constructed usingSearch Tool for the Retrieval of Interacting Genes/Proteins (STRING;www.string-db.org/) and visualized by Cytoscape3.6.1 (https://cytoscape.org/).

Correlation and survival analyses
The R package "corrplot" was used to analyze the correlation of immune cells. The correlation of keratin, type I cytoskeletal 23 (KRT-23) and C-X-C motif chemokine ligand 9 (CXCL9), and CXCL10 and C-C motif chemokine ligand 5 (CCL5), in TCGA dataset was analyzed through cbioportal(www.cbioportal.org/).Correlation in tumor tissues from patients was conducted by Prism7 (GraphPad, San Diego, CA, USA). For survival analyses, samples were divided into four clusters of hot and cold tumors. The R package "survival" was used to assess the survival difference using the logrank test.

Construction and validation of prediction model
The DEGs between hot and cold tumor were utilized to perform unicox analysis using R package "survival" to select survival-related DEGs with p<0.05 and followed with LASSO regression to optimize the gene sets using R package "glmnet". The risk formula was calculated by multicox in R package "survival".
The R package" survivalROC" was used to analyze the sensitivity and speci city of the model.

Small interfering RNA(siRNA) transfection
Knockdown of KRT23 expression was achieved using the jetPRIME ® Transfection Reagent kit(Polyplustransfection, Illkirch-Graffenstaden, France).HeLa cells (1×10 5 )were seeded in six-well plates with RPMI1640 medium. Before transfection, siRNA of KRT23 was diluted to 20 µM according to manufacturerinstructions.Then, 200 µL of transfection buffer and 4 µL of jetPRIME reagents were mixed and incubated for 10 sat room temperature.After that, 50nM of siRNA was added and incubationallowed for 15 min at room temperature. siRNA e cacy was analyzed by RT-qPCR after 48 h. The sequence of siRNA synthesized by Gene Pharma(Shanghai, China)is listed in Table2.

Transwell™ assay
Migration of CD8 + T cells was analyzed through the Transwell assay. CD8 + Tcells (2×10 4 ) isolated by microbeads from healthy donors were activated with CD3/CD28 beads and seeded in the upper chamber of the Transwell apparatus with serum-free medium (Millipore, Billerica, MA, USA).HeLa cells (2×10 4 ) were seeded in the lower chamber with RPMI1640 medium. The number of CD8 + T cells was calculated using ow cytometry.

Enzyme-linked immunosorbent assay (ELISA)
Tumor cells were transfected with siRNA for 48 h.Then, supernatants were collected andcentrifugation(1500rpm,5 minute) to remove debris. The CCL5 concentration was measured by the LEGEND MAX™ Human CCL5 (RANTES) ELISA kit according to manufacturer(Biolegend, San Diego, CA, USA) instructions. Brie y, standard dilutions and samples were prepared, followed by addition of 50 μL of Assay Buffer B to each well. Then, 50 μL of the standard or samplewas added to the appropriate well. This action was followed by incubation of the plate at room temperature for 2 hwith agitation at 200 rpm.
Then, 100 μL of Human CCL5 Detection Antibody solution was added, followed by 100 μL of Avidin-HRP A solution, to each well. Results were read at an optical density of 450 nm.

Detection of multiple chemokines
We used the LEGENDplex ™ kit(Biolegend) to detect the chemokines secreted by tumor cells.This kit can detect 13 types of chemokines. First, 25 µL of assay buffer was added to the standard or sample in each tube. Second, we added 25 µLof mix beads (A and B) and permitted incubation at room temperature for 2 hwith agitation at 500rpm. Third, we added 25 µLof antibodies to each tube and allowed incubation at room temperature for 1 h with agitationat 500rpm. Next, we added 25 µL ofSA-PE to each tube and washed with washing buffer. The uorescence intensity was detected by a ow cytometer and analyzed by LEGENDplexv8.0.

Statistical analyses
Statistical analyses were undertaken using by Prism 7(GraphPad) and R 3.6.3.Two-tailed unpaired t-tests and the Wilcoxon test were used to compare the difference in data between two groups. Spearman's rank correlation coe cientwas used to evaluate correlation.p < 0.05 was considered signi cant.

In ltration pattern of immune cells in adjacent tissue and tumor tissue
We carried out a multistep analysis to explore the in ltration of immune cells into CC (Fig. 1). First, we estimated the number of immune cells in each sample by ssGSEA and xCell, which have various algorithms and focus on different immune cells between adjacent tissue and tumor tissue. The numberof each cell type in ltrated into the TME was different, which revealed the complexity of the TME.ssGSEA and xCellshowed consistent results. In general, the number of the most adoptive immune cells in tumor tissue was higher than that in adjacent tissue, such as activated CD4 + T cells, effector memory CD4 + T cells, type-17 T-helper (Th17) cells,and Th2 cells, which indicated an activated immune response in tumor tissue. The number of CD8 + T cells tendedto behigher in tumor tissues, but not signi cantly. Cells from the innate immunesystem showed highin ltration in normal tissues (Fig2A,B). Tumor tissues had a lowerImmune Score, but there was no signi cant difference in the Stromal Score (Fig.2C,D). We also compared the difference in immune cells in patients who received radiotherapy.Pro-B cells and Th1cells tended to accumulate in tumor tissue after radiotherapy(sFig1A,B). These results revealeddistinct patterns of in ltration between adoptive immune cells and innate immune cells.

Characterization of immune clusters in CC tissue
Activation of an e cacious anti-tumor immune response requires the synergistic action of multiple cells.
To explore the relationships between different cell types, we undertook a correlation analysis of in ltrating cells in tumor tissues. Most in ltrating cells showed a high correlation with each other, especially activated CD8 + ,CD4 + T, dendritic, and B cells. We observed a higher correlation of immunosuppressive cells and immune cells, such as regulatory T cells, myeloid-derived suppressor cells, and M2 macrophages, which suggested immune suppression induced by tumor cells after activation of the immune system. Cells of the innate immune system, such as monocytes, neutrophils, and natural killer cells, showed a weak association with other cells, which demonstrated a unique anti-tumor immune process (Fig.3A,B).
Next, we performed consensus clustering of all samples based on the proportions of immune cells to identify the subtypes of in ltrating immune cells. The consensus matrix heatmap showed four clearly identi ed groups estimated by two methods, respectively (Fig.3C,D). We observed a gradual increase in in ltration of immune cells in tumor tissue from groups 1 to 4, in which group1 and group 2 lacked in ltration of immune-related cells, group 3 had modest in ltration, and group 4 showed abundant in ltration of immune cells (sFig2A,B).In accordance with these results, group 4 had the highest Immune Score (sFig.2C,D). To further characterize the clusters of CC cells, we took the intersection of each group of the two methods, denoted as clusters 1 to 4(sFig.3A).In accordance with the results stated above, cluster4 had a high Immune Score (Fig4A).Next, we analyzed expression of genes involved in the immune response, immune tolerance, and antigen presentation in the four clusters. Expression of immune checkpoint-related genes (CD276, CD274, CD40,CTLA4, HAVCR2,LAG3, PDCD1), antigen presentationrelated genes (B2M, HLA-B, HAL-C, HLA-DQA1, TAP1, TAP2, HLA-DQA2), cytokine-related genes (GZMB, GZMH, IFNG, PRF1, TNF) and chemokine-related genes(CCL5, CXCL10, CXCL13, CXCL9) increased gradually from cluster 1 to cluster 4 (Fig.4B, sFig.3C, D). Survival analyses revealed that cluster 4 had the longest survival relative to that in clusters 1,2, and 3 in terms of OS and progression-free interval ( Fig.4C,D).

Survival status and signaling alterations between hot tumors and cold tumors
To further explore the mechanisms of immune-cell in ltration, we rede ned cluster1, cluster2, and cluster3 as cold tumors and cluster4 as hot tumors based on immune cells and survival status. Hot tumors had longer OS and progression-free interval compared with those for cold tumors (Fig. 5A,B).Next, we analyzed the difference between the two groups at the transcriptional level. Hot tumors and cold tumors showed different transcription patterns according to volcano plots (Fig. 5C). Finally, 657 messenger (m)RNAs with upregulated expression in hot tumors and 55 mRNAs with upregulated expression in cold tumors were identi ed. To further explore DEGs function, enrichment analyses using GO and KEGG databases were done. The GO database revealed that DEGs in hot tumors were enriched primarily in "T cell activation", "regulation of lymphocyte activation", "leukocyte cell-cell adhesion", "regulation of T cell activation" and "leukocyte proliferation"; none of these GO terms were enriched in cold tumors (Fig 5D).Analyses of enrichment of DEGs using the KEGG database revealed that DEGs in cold tumors were signi cantly enriched in "apical part of cell", "actin-based cell projection" and "apical plasma membrane" (Fig. 5E). In hot tumors, DEGs were enriched mainly in "cytokine-cytokine receptor interaction", "chemokine signaling pathway" and "cell adhesion molecules", which indicated an activated immune response in hot tumors (Fig.5F). These results suggested that immunity was activated comprehensively in hot tumors, particularly the T cell-mediated immune response. PPI networks also revealed that the "hub genes" in DEGs of hot tumors were mainly immune-related chemokines and cytokines(sFig. 4A,B).

Inhibition of KRT23expression promotes in ltration of CD8 + T cells
The results stated above revealed a correlation between in ltration of many immune cells and longer survival. Therefore, promoting in ltration of immune cells in cold tumors may enhance antitumor immunity and prolong survival. KRT23 showed the highest expression in cold tumors as compared with hot tumors. Enrichment analyses using the KEGG and GO databases revealed that KRT23-related genes were negatively correlated with the immune response (Fig. 6A,B). Knockdown of KRT23 expression in HeLa cells inhibited cell proliferation (Fig. 6C, D), which suggested an important role of KRT23. To explore how KRT23 affected in ltration of immune cells, we measured expression of KRT23 and CD8 + T cellrelated chemokines in tumor samples(CCL5, CXCL9,CXCL10). KRT23 expression was negatively correlated with these chemokines, and this result was con rmed using TCGA database (Fig. 6E,sFig. 5).Furthermore, knockdown of KRT23 expression promoted CCL5 expression at the mRNA level. Cytometric bead arrays and ELISA con rmed that inhibition of KRT23 expression increased CCL5 secretion (Fig. 6G, H)and recruitment of CD8 + T cells (Fig.6I).

Construction and validation of a prediction model based on DEGs
Next, we used DEGs to construct a prediction model. We undertook a univariate Cox analysis followed by a lasso regression analysis (sFig.6A). To optimize the model, we carried out a multivariate Cox analysis and, nally, identi ed 11 genes to construct our model(sFig.6B). Heatmaps revealed expression of these genes in high and low risk groups in a training cohort and testing cohort. Survival analyses showed that patients with a high risk of CC had shorter survival in the training cohort and testing cohort (Fig. 7A, B). To explore the accuracy of our model, we analyzed receiver operating characteristic (ROC)curves in the training cohort and testing cohort at 1,3, and 5 years. Our model had a higher area under the ROC curve (AUC) in the training cohort and testing cohort (Fig. 7C, D).Our model was constructed by DEGs in hot tumors and cold tumors, so we hypothesized that this model may predict the response to immunotherapy. Hence, we used two external cohorts of patient sreceiving immunotherapy. Patients with a higher risk of CC had shorter survival in the two cohorts (Fig. 7E,F).These results suggested that our model could not only predict the survival of patients with CC, it could also be used to predict the response to immunotherapy.

Discussion
Most cases of CC are caused by HPV infection (26,27). Chemotherapy for CC is limited. The optimal regimen against recurrent CC or mCC is a combination of cisplatin, paclitaxel, and bevacizumab. This regimen is associated with an overall response rate of 48% and median survival of 17 months (28).The side-effects caused by radiotherapy restrict its clinical application in CC (29).Therefore, new and e cacious strategies for CC are needed urgently.
Immunotherapy has shown sustainable clinical response and is rst-line treatmentfor various types of tumors (30)."Cancer immunotherapy" is a general term describing harnessing of the immune system of a patient to elicit antitumor effects (31). Antibodies against PD-1 and PD-L1 are used commonly for cancer immunotherapy. They work by releasing the "inhibitory brakes" of T cells, resulting in robust activation of the antitumor immune response (32).
The major risk factor for CC is HPV infection (33).High-risk HPV types, such as 16 and 18, are more likely to persist and integrate into the host genome to enable excess expression of the oncoproteins E6 and E7. These oncoproteins interfere with the immune response (34).Therefore, therapies targeting the HPV have been attempted, but the effects have been suboptimal (12).However, the retained viral antigens in CC make immunotherapy an attractive option because they could be recognized as foreign. Indeed, several clinical trials have used antibodies against PD-1 or PD-L1 (31,35). E cacious immunotherapy is reliant on the in ltration of lymphocytes and antigen-presenting cells.In general, the TME can be divided into two broad categories: "T cell-in amed" and "non-T cell-in amed" (36). Several methods have been used to estimate the immune pro le in the TME: ssGSEA, CIBERSORT, TIMER, MCP-counter, and xCell (37)(38)(39)(40)(41). ssGSEA and MCP-Counter use speci c cell-maker genes and score the immune pro le through expression of these genes. CIBERSORT focuses on the ratios of each cell type using Nu-support vector regression. xCell is an integration of these methods and expands the cells that can be evaluated to 64 types. In order to more accurately re ect the level of immune cells in the TME of CC, we used two methods to score the immune cells. Results from the comparison between tumor and adjacent tissues as well as correlation analysis of estimated immune cells showed constant ndings of the two methods, suggesting that these two methods can be used to estimate immune levels. Based on the consensus clustering, the CCcan be divided into to 4 clusters, and clusters with higher immune in ltration had a better survival. In line with previous studies, Wang J et.al, also reported that CD4 + T cells are independent prognostic factor of CC estimated by CIBERSORT [42] . Meanwhile, in ltration of immune cells also correlated with the response of chemotherapy [43] . We further re-divided the 4 clusters into 2 subtypes and donated as "hot" and "cold" tumor based on the immune levels, which "hot" tumor is T cell-in amed and "cold" tumor is Non-T cellin amed. Pathway enrichment analysis con rmed that hot tumor is a state of active immune response.
Cold tumors are characterized by in ltration of few immune cells. They are the most challenging to eradicate and are invariably associated with a poor prognosis (44). Several strategies have been used to covert cold tumors to hot tumors:radiotherapy, chemotherapy, targeted therapy, and adoptive-cell therapy (45)(46)(47)(48)(49). We compared the difference between hot tumors and cold tumors, and identi ed KRT23 as most upregulated gene in cold tumors. Keratin is the main component of epithelial cells, and malignant tumor cells originate from these epithelial cells. KRT23 is a newly identi ed gene in the KRT family (50,51).Studies have reported that KRT23 overexpression promotes the migration of ovarian cancer cells via epithelial-mesenchymal transition (52). KRT23 promotes proliferation of colorectal tumor cells through increasing expression of telomerase reverse transcriptase (52). Although the oncogenic role of KRT23 has been explored, how KRT23 affects the immune response is not known. We observed that KRT23 expression was negatively correlated with the immune response. Knockdown of KRT23 expression in tumor cells results in increased secretion of CCL5(a potent chemokine that recruits CD8 + T cells)and inhibits the proliferation of tumor cells. Our results indicate that inhibition of KRT23 expression not only inhibits tumor growth but also enhances the anti-tumor response. Hence, a potential combination strategy of targeting KRT23 and anti-PD1 and PD-L1 could be a rational approach against CC.

Conclusion
In summary, we undertook a comprehensive analysis of the in ltration of immune cells in CC. We identi ed hot tumors and cold tumors of CC, and observed the former to have a favorable outcome. We demonstrated KRT23 to be a negative regulator of the immune response, and that knockdown of KRT23 expression promotedCCL5 secretion. In addition, we constructed a prediction model based on DEGs between two types of CC. This model performed well in predicting the survival of CC patients receiving immunotherapy. We provided new insights into the in ltration of immune cells into CC and highlighted KRT23 as a potential target to enhance immunotherapy against CC.

Declarations
Xia Li designed and performed the experiments, wrote the manuscript, analyzed the data. Yanmei Cheng prepared gures 1-6 and Yanyan Jia performed s gures 1-6. Huirong Shi supported, designed the study, analyzed the data and revised the manuscript. All authors read and approved the nal manuscript.

Funding
This study had no grants to support.

Data availability statement
The data sets used and/or analyzed during the current study are available from the corresponding author on reasonable request.The data that support the ndings of this study are openly available on the online website UCSCXena (https://xenabrowser.net/). GSE77280 are available from on the GEO database (https://www.ncbi.nlm.nih.gov/geo/). and patients with metastatic urothelial cancer treated with anti-PDL1 agents were downloaded as online website supplied in the article(http://resea rchpub.gene.com/IMvigor210CoreB iologies/) Ethics approval and consent to participate All participants provided written informed consent for their specimens to be used in this study. The study protocol was approved by the ethics committee of First A liated Hospital of Zhengzhou University.

Consent for publication
Not applicable.

Competing interests
The authors declare that they have no competing interests  Figure 1 Multiple step analysis of this study