Here, ssGSEA is used to calculate the enrichment score of each sample and obtain the correlation between PCSK9 expression and pathway score. R software GSVA package was used to analyze, choosing parameter as method = ‘ssgsea’. The correlation between genes and pathway scores was analyzed by Spearman correlation. ClusterProfiler package (version: 3.18.0) in R software was employed to analyze the Gene Ontology (GO) function of potential targets and enrich the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway, p <0.05 or FDR <0.05 is considered to be a meaningful pathway (enrichment score with −log10 (P) of more than 1.3).

To analyze the stemness features, we performed the spearman correlation between mRNAsi of various tumors and PCSK9 expression. This method refers to the OCLR algorithm constructed by Malta, which contains 11774 different gene profiles (22).

The Drug Sensitivity Genomics Project (GDSC) and The Cancer Therapeutics Response Portal (CTRP)are two databases, which combine drug sensitivity and genome data sets to promote the new therapeutic biomarkers for cancer therapy. Through these two databases, we investigated the role of PCSK9 expression in cancer therapeutic response. Pearson correlation analysis was performed to obtain the correlation between PCSK9 mRNA expression and drug IC50. P-value was adjusted by FDR.

A total of 25 neuroblastoma (NB) patients from Children’s Hospital of Soochow University (Suzhou, China) were consecutively enrolled in a study between January 2016 and December 2019 and followed up until December 2021. NB specimens and adjacent normal tissues were collected at the time, frozen in liquid nitrogen, and stored at −80^°C until use. No patients received chemotherapy or radiotherapy or any treatment for the tumor before surgery or tissue biopsy. The sample collection and related experiments met the ethical requirements of the Children’s Hospital of Suzhou University, Suzhou, China. The Clinical Research Ethics Committee approved this study at the Children’s Hospital of Suzhou University, Suzhou, China, and all patients or their parents signed informed consent forms.

Proteins were isolated from tissues using RIPA lysis buffer (Biotime, Shanghai, China), separated using SDS-PAGE, and transferred onto a nitrocellulose membrane (Bio-Rad Laboratories, CA, USA). Then, we blocked proteins with 5% skim milk for 30 min and incubated the membranes with diluted primary antibodies. Primary antibodies for GAPDH (ab8245, 1:10000 dilution; Abam, Cambridge, UK), PCSK9 (ab181142, 1:1000 dilution; Abam, Cambridge, UK), CD11b (ab133357, 1:1000 dilution; Abam, Cambridge, UK), CD45 (ab40763, 1:5000 dilution; Abam, Cambridge, UK), CD68 (ab213363, 1:1000 dilution; Abam, Cambridge, UK), and BSA-1 (ab219724, Abam, Cambridge, UK) were purchased from Abcam (Cambridge, UK). After incubating with horseradish peroxidase-conjugated secondary antibodies, the immune complexes were detected with an ECL detection kit (Millipore, Billerica, MA, USA) and quantified using a Gel-Pro Analyzer (Media Cybernetics Corporation, USA). Total RNA was isolated from patients with NB by using Trizol and then converted into cDNA by reverse transcription (GoScriptTM Reverse Transcription system, USA). RT-qPCR was used to assay the mRNA expression of PCSK9, immune-related genes, and pathway-related genes. PCSK9 primers were purchased from Sino Biological Inc. The primers used for qRT-PCR were as follows: PCSK9, 5′-GCT GAGCTGCTCCAGTTTCT-3′ (forward) and 5′-AAT GGCGTAGACACCCTCAC-3′ (reverse); and GAPDH, 5′-AAGGTGAAGGTCGGAGTCAAC-3′ (forward) and 5′-GGGGTCATTGATGGCAACAATA-3′ (reverse).

Gene expression data from the TCGA and GTEx databases were analyzed by using t-tests. The Kruskal-Wallis test was used to evaluate the difference among various tissues, and the Wilcoxon test was used to determine the gene expression differences between normal and tumor tissues. In PrognoScan, univariate Cox regression analysis was used to analyze the survival time of patients with the HR and p value. In GEPIA and Kaplan-Meier Plotter, log rank test was used to compare survival rate of patients stratified according to the different expression levels of PCSK9. Other online analysis websites of GEPIA2, cBioportal and GSCALite were also used. Spearman correlation analysis was calculated between the expression of PCSK9 and the level of infiltrating immune cells, the level of immunosuppressant or immunostimulant factors and infiltration scores of six immune infiltrations. Correlations were considered statistically significant when p < 0.05 for all statistical analyses. The experimental data were analyzed using SAS 9.3 statistical software. Statistical analysis was performed using the t test. Differences were considered statistically significant at p < 0.05. OS was defined as the time interval between the date of surgery and the date of progression or death. Survival analyses were conducted by Kaplan–Meier curves (p values from log-rank test) by R version 3.5.3 and the survival package. HRs were calculated using the R package.

Abnormal PCSK9 expression has been reported in various cancer types. Previous studies on PCSK9 expression in cancer used several research methods, such as DNA microarrays, but were limited to relatively small sample sizes and limited numbers of cancer types. This study has provided a more comprehensive analysis of PCSK9 expression in cancer. Since PCSK9 has a potential role as an important new target for cancer diagnosis and prognosis, we analyzed the PCSK9 mRNA levels across different cancers from TIMER, GEPIA, and UALCAN databases. Data from these databases indicated that PCSK9 mRNA expression had inter-tumor heterogeneity, with some tumors having very high levels of PCSK9 (BRCA, CESC, CHOL, COAD, ESCA, HNSC, LIHC, READ, SKCM, STAD, THCA, and UCEC). In contrast, others were characterized by low levels of PCSK9 expression (LUAD, PRAD, KIRP, KIRC, PCPG, LUSC, and GBM) ( Figure 1A ). Different stages and subtypes of a tumor may exhibit differential expression of PCSK9. Thus, we further assessed PCSK9 expression in different clinical stages and subtypes from GEPIA and UALCAN. Briefly, differential expression of PCSK9 was obtained from GEPIA to correlate with clinical subtypes of tumors, including BRCA, COAD, ESCA, HNSC, STAD, LUSC, OV, and UCEC ( Figure 1B ). As shown in Figure 1C , PCSK9 expression was elevated in some tumors, including BLCA, BRCA, CESC, HNSC, READ, STAD, THCA, COAD, ESCA, UCEC, and LIHC. Whereas lower PCSK9 expression in later stages was observed only in KIRC, KIRP, and LUAD.

It has been widely acknowledged that genomic mutations are closely associated with tumorigenesis. To determine the genomic mutations of PCSK9 in tumors, we reviewed the genetic alterations of the PCSK9 gene in cancer patients using the cBioPortal database. Notably, patients with ovarian epithelial tumors had the highest frequency of PCSK9 genetic alterations (5%), including amplification of copy numbers and deep deletions ( Figure 2A ). In addition, several cancer types (e.g., non-small cell lung cancer, cervical squamous cell carcinoma, ESCA, LIHC, SARC, BLCA, and BRCA) had PCSK9 mutations (amplifications or deep deletions). Tumors with dominant PCSK9 mutations included cervical adenocarcinoma, esophageal squamous cell carcinoma, esophagogastric adenocarcinoma, renal non-clear cell carcinoma, COAD, HNSC, KIRC, PAAD. Deep deletions were more common in various neuroepithelial tumors, PCPG, LGG, and GBM. These results revealed the highly heterogeneous inheritance and expression changes of PCSK9 in different types of cancer.

We next differentiated the distribution of mutations and collected the types, sites, and case number of genetic alterations of PCSK. As shown in Figures 2B, C , there were 93 missense and 6 truncation mutations in PCSK9. The site S91T/E92Afs*78/S91L, in the inhibitor_I9 domain, was confirmed to have the highest abundance of mutations. We future searched the cBioPortal website for PCSK9-related tumor gene mutations and identified six types of tumors with a total of eleven protein changes: Melanoma (D321N, E84K, G176E, E405K, S249G), Uterine Endometrioid Carcinoma (V79M), Colorectal Adenocarcinoma (R272L), Diffuse Large B-Cell Lymphoma (Q454H), Serous Ovarian Cancer (R93H, S5G), and Esophageal Adenocarcinoma (C526F). Subsequently, we conducted a search on the PolyPhen-2 website (http://genetics.bwh.harvard.edu/pph2/) and found that in five tumor types (Melanoma, Uterine Endometrioid Carcinoma, Diffuse Large B-Cell Lymphoma, Serous Ovarian Cancer, and Esophageal Adenocarcinoma), the protein changes were predicted to be either “probably damaging” or “possibly damaging.” However, the mutation in Colorectal Adenocarcinoma exhibited benign characteristics. Moreover, when we further correlated the prognosis of these six tumor types, we observed that three of them (Melanoma, Serous Ovarian Cancer, and Uterine Endometrioid Carcinoma) displayed worse prognoses, suggesting a potential association with inherent mutations ( Supplementary Table 1 ).

DNA methylation is an epigenetic modification that can alter gene expression. The alteration of DNA methylation may be an essential factor in tumorigenesis. Through the UALCAN database, we explored PCSK9 methylation between tumors and normal tissues. The results showed that the methylation of PCSK9 was upregulated in various tumors, including BRCA, CHOL, HNSC, KIRC, KIRP, LGG, LUAD, LUSC, PRAD, SARC, and THCA ( Figure 2D ). Accordingly, these results indicated that PCSK9 may mediate tumorigenesis by regulating DNA damage or methylation status in human cancers. Furthermore, similar results were concluded from GSCALite data, which revealed a relatively high mutation frequency in diverse types of cancer. We also investigated the frequency of CNV changes of PCSK9. The results showed that tumors, including TGCT, BRCA, COAD, SKCM, PCPG, PRAD, LIHC, BLCA and KIRC, were significantly correlated with the CNV changes ( Figure 2E ).

The PrognoScan database was used to explore the relationship between PCSK9 expression and prognosis in each cancer, results are summarized in Figure 3 . Notably, PCSK9 expression was significantly associated with five cancer types: brain, colorectal, lung, head and neck, and ovarian cancer ( Figure 3A ). Among them, PCSK9 showed a detrimental relationship in three cancer types, including colorectal cancer [OS: total = 177, HR = 1.46, Cox p = 0.048], lung cancer [OS: total = 56, HR = 1.43 Cox p = 0.028], lung cancer [RFS: total = 56, HR = 1.30, Cox p = 0.011] and head and neck cancer [RFS: total = 28, HR = 3.26, Cox p = 0.048]. However, PCSK9 exhibited a protective effect in two other cancer types, brain cancer [OS: total = 67, HR = 0.01, Cox p = 0.030] and ovarian cancer [PFS (progression free survival): total = 110, HR = 0.89, Cox p = 0.044] ( Figure 3A ).

To evaluate the relationship between PCSK9 protein expression levels and the prognosis of tumor patients, the Kaplan-Meier Plotter database was used to study the PCSK9 protein expression level in 21 cancer tissues ( Figure 3B ). In some cancer tissues, such as BLCA (OS, HR = 1.68, p = 0.0005), THYM (OS, HR = 4.26, p = 0.02), KIRC (OS, HR = 2.56, p = 0.002), LIHC (OS, HR = 1.64, p = 0.005), SARC (OS, HR = 1.9, p = 0.008), LUAD (OS, HR = 1.43, p = 0.015), OV (OS, HR = 1.38, p = 0.014), PAAD (OS, HR = 2.00, p = 0.003), KIRP (RFS, HR = 3.44, p = 0.001), LIHC (RFS, HR = 1.46, p = 0.025), TGCT (RFS, HR = 3.44, p = 0.032), and PAAD (RFS, HR = 3.21, p = 0.006), high expression of PCSK9 correlated with poor OS and RFS. Other cancers exhibited a protective role for PCSK9 with low expression, such as BRCA (OS, HR = 0.58, p = 0.001; RFS, HR = 0.64, p = 0.047), UCS (OS, HR = 0.61, p = 0.022; RFS, HR = 0.56, p = 0.044), ESCA (RFS, HR = 0.26, p = 0.013), and HNSC (RFS, HR = 0.26, p = 0.007).

Similar work was also performed using the GEPIA database. High expression of PCSK9 was associated with a poor OS prognosis in BLCA (p = 0.002), UVM (p = 0.012), LUAD (p = 0.028), SKCM (p = 0.0002), KIRC (p = 0.014), KIRP (p = 0.001), LIHC (p = 0.023), and BRCA (p = 0.009) ( Figure 3C ). The data showed that high expression of PCSK9 was associated with an adverse DFS prognosis in PAAD (p = 0.031), KIRC (p = 0.014), KICH (p = 0.046), BLCA (p = 0.009), and LUAD (p = 0.026) ( Figure 3C ).

To verify the correlation between PCSK9 expression and multiple clinicopathological characteristics, we used LIHC as an example. The results indicated that PCSK9 was associated with a detrimental prognosis with five patient or tumor characteristics: female (OS: p = 0.01; PFS: p = 0.02), Asian race [OS: p = 0.01; PFS: p = 0.01], no hepatitis virus infection [OS: p = 0.03; PFS: p = 0.02], pathology stage 2 [OS: p = 0.01] and stage 3 [OS: p = 0.04], grade 2 [OS: p = 0.01], and AJCC stage 2 (OS: p = 0.01) and stage 3 (PFS: p = 0.03) ( Figure 4A ). A nomogram prediction model was generated by integrating the above clinicopathological parameters and PCSK9 expression levels. The calibration curve is further evaluated to show that the predictions made by nomogram is in considerable agreement with the actual survival ( Figures 4B–D ).

Cancer progression involves the gradual loss of differentiated phenotypes and stemness features. Our correlation analysis revealed that the expression of PCSK9 was positively correlated with the mRNAsi in COAD, LUAD and LIHC. As for tumors such as SKCM, BLCA and LGG, there is no correlation between the characteristics of stemness and the expression of PCSK9 ( Figures 5A, B ).

Taking BLCA as an example, we found that abnormal expression of PCSK9 was associated with tumor invasion characteristics, cell hypoxia manifestations, and EMT-related markers ( Figure 6A ). In BRCA, PCSK9 expression was involved in tumor inflammation signature, reactive oxygen species and transforming growth factor beta (TGF-β) ( Figure 6B ). The KEGG pathway analysis and GO enrichment analysis were performed to demonstrate the primary biological pathways and potential targets of major potential PCSK9 mRNA. The results showed that PCSK9 positively regulated cellular adhesion, cholesterol and fatty acid metabolism, PI3K-Akt signaling pathway and immune-related functions in BRCA, LIHC and LUAD ( Figures 7A–C ).

Tumor-infiltrating lymphocytes are independent predictors of sentinel node status and cancer prognosis. The correlation between the expression levels of PCSK9 protein and immune cells in pan-cancer tissues were analyzed using the TIMER database. PCSK9 expression levels significantly correlated with immune cells (CD8+ T cells in 14 types of cancer, dendritic cells (DCs) in 11 types of cancer, and macrophages in 13 types of cancer). We further studied whether PCSK9 expression correlated with the infiltration of different immune cell subtypes using the xCell online tool. The results showed that PCSK9 expression was significantly correlated with subtypes of infiltrating immune cells in various tumors, including HNSC, TGCT, ESCA, COAD, STAD, LUSC, and LIHC ( Figure 8A ). CD8+ T cells, DCs, and M2 macrophages were the immune cell subtypes most positively associated with PCSK9 expression in these different cancers.

Immune checkpoint inhibitors (ICIs) are novel tumor immunotherapy agents that play an essential role in tumor immunotherapy. Subsequently, we analyzed the correlation between PCSK9 expression and immune checkpoint gene expression levels using the R software package in the SangerBox database. The correlations between 46 immune checkpoint genes and PCSK9 protein expression levels were calculated, and a significant relationship was found in many cancer types with many of the 46 genes, such as THCA (40 of 46), BRCA (37 of 46), LUAD (37 of 46), BLCA (35 of 46), and TGCT (33 of 46) ( Figure 8B ). Furthermore, the results showed that the expression of PCSK9 in BRCA, BLCA, PRAD, THCA, LIHC, PCPG, and UVM was negatively correlated with immune checkpoint genes. Furthermore, a positive correlation was found in CESC, CHOL, LUAD, STAD, and TGCT. Notably, these immune checkpoint genes contained a broad spectrum of immune regulators, including signaling chemokines, immune stimulators, immune inhibitors, and major histocompatibility complex molecules, such as those related to regulatory T (Treg) cells (chemokine receptor 8, CCR8; forkhead box protein p3, FOXP3; signal transducer and activator of transcription 5B, STAT5B), B cells (CD19), macrophages (CD68), interleukin-10 (IL10), Th17 (signal transducer and activator of transcription 3 (STAT3), neutrophils (integrin subunit alpha m, ITGAM), natural killer (NK) cells killer cell immunoglobulin-like receptor 2DL (KIR2DL), DCs (major histocompatibility complex class II DR beta 1 (HLA-DPB1); histocompatibility complex class II DR alpha (-DRA; integrin subunit alpha X(ITGAX), Th1 cells (interferon gamma (IFNy); tumor necrosis factor (TNF), Th2 cells (signal transducer and activator of transcription 6 (STAT6; signal transducer and activator of transcription 5A (STAT5A), and exhausted T cells (programmed cell death protein 1 (PDCD1); CTAL4; hepatitis A virus cellular receptor 2 (HAVCR2). Additionally, our results showed a significant correlation between PCSK9 expression with PDCD1 and CTLA4, which referred to T cell exhaustion and led to the loss of T cell function in patients with common chronic infections and cancer.

Immune infiltration plays a vital role in the TME. To evaluate the association between PCSK9 expression and immune infiltration, we further investigated the relationships between the PCSK9 expression and immune cells in three different cancers (BLCA, BRCA, and LIHC) from the TIMMER database ( Figure 8C ). The immune cells included B cells, CD8+ T cells, CD4+ T cells, macrophages, neutrophils, and DCs. In BLCA, the expression level of PCSK9 significant and positively correlated with CD4+ T cells (R = 0.126, p = 8.55E−05), CD8+ T cells (R = 0.11, p = 6.00E−04), macrophages (R = 0.09, p = 4.88E−03), neutrophils (R = 0.158, p = 1.01E−06), and DCs (R = 0.162, p = 5.23E−07). In BLCA, the correlation between the expression level of PCSK9 with immune cells was significant, including CD8+ T cells (R = 0.246, p = 1.82E−06), CD4+ T cells (R = 0.086, p = 9.93E−02), macrophages (R = 0.167, p = 1.35E−03), neutrophils (R = 0.231, p = 8.38E−06), and DCs (R = 0.356, p = 2.50E−12). In LIHC, we also found a significant correlation of PCSK9 with B cells (R = 0.136, p = 1.17E−02), CD4+ T cells (R = 0.132, p = 1.40E−02), macrophages (R = 0.143, p = 7.98E−03), neutrophils (R = 0.098, p = 7.03E−02), and DCs (R = 0.131, p = 1.60E−02). There was no correlation between PCSK9 and B cells in BRCA or BLCA (p > 0.05), or LIHC with CD8+ T cell (p > 0.05). In addition, PCSK9 expression strongly correlated with tumor-associated macrophages and DCs in these three cancers. These findings strongly suggested that PCSK9 affected the immune microenvironment by interacting with immune cell infiltration in various cancers.

Considering the heightened sensitivity of gamma delta T cells, particularly gamma 9 delta 2 T cells, towards alterations in the cholesterol pathway, it would be significant to include these cells as a specific control in the immune infiltration analysis. Thus, we further evaluated the immune cell infiltration scores associated with tumor samples obtained from the TCGA database. By utilizing the “corr.test” function in R software, we ultimately identified a significant correlation between PCSK9 expression and gamma delta T cell immune infiltration scores. Specifically, we observed a significant association between PCSK9 expression and gamma delta T cell infiltration in five types of tumors (TGCT, COAD, PRAD, LUSC and LIHC). These results indicated incorporating gamma delta T cells into the analysis would provide additional insights into the potential impact of PCSK9 and its role in modulating immune responses within the tumor microenvironment ( Supplementary Figure 1 ).

To further unravel the potential predictive value of PCSK9 gene alterations for ICI treatment, we then investigated the relationship between PCSK9 alterations and six common immune infiltrates (B cells, CD4+ T cells, CD8+ T cells, macrophages, neutrophils, and DCs) across multiple cancer types. We demonstrated that PCSK9 is frequently altered across different cancer types. Furthermore, we investigated the relationship between the PCSK9 CNV and immune infiltration in various cancers ( Figure 8D i–v). For six types of immune cells, we verified a correlation of the PCSK9 CNV with immune infiltration. CD8+ and CD4+ T cells were associated with PCSK9 deletions in LIHC, BLCA, BRCA, and THCA. Deletion of PCSK9 was significantly related to infiltration of CD4+ T cells (p = 0.001), neutrophils (p = 0.001), and DCs (p = 0.05) in BRCA. These results suggested a possible mechanism by which immune cells may be affected by PCSK9 in the TME, which may help direct future immunotherapy treatments.

GDSC showed that high expression of PCSK9 could make cancer more sensitive to IPA-3 (target PAK1), (5Z)-7-Oxozeaenol (target TAK1), Nutlin-3a (target MDM2), Navitoclax (target BCL2) and resistant to Docetaxel (target microtubule stabilizer), Epothilone B (target microtubule stabilizer), OSU-03012 (target PDK1) ( Figure 9A ). While in CTRP database, IC50 of QW-BI-011, CCT036477, CIL70, PRIMA-1, PRIMA-1-Met, teniposide, ML210, BRD-K92856060, BRD-K26531177 and avrainvillamide showed top 10 significant positive correlations with the expression level of PCSK9 ( Figure 9B ). Collectively, these results may provide new ideas for developing potential drugs relating to the expression of PCSK9 for treating cancers.

To confirm PCSK9 expression in NB, we measured mRNA and protein levels of PCSK9 in paired NB and adjacent non-tumor tissues using qRT-PCR (n = 18; Figure 10A ) and western blotting (n = 18; Figure 10A ). Our findings revealed that PCSK9 expression in tumor tissues was significantly higher at both the transcriptional and protein levels than in adjacent normal tissues (p < 0.001).

Immune cells play an essential role in tumor immune tolerance. Immunotherapy (such as disialoganglioside GD2) has been incorporated into first-line treatment regimens of relapsed NB to significantly improve patient outcomes. Remarkably, immune cells, such as NK cells and macrophages, are thought to be the main effectors of anti-GD2 antibody potency in NB tumors. Thus, to investigate the role of immune cells in the NB TME, we used western blotting to measure the levels of CD11b, CD45, and CD68 expression in NB tissues ( Figure 10B ). The results showed that CD11b (PDC#2, PDC#4, and PDC#5) was positive for immune infiltration of NB, which represented the potential involvement of NK cells in immune infiltration ( Figure 10B ).

Furthermore, we investigated the PCSK9 prognostic value in NB patients. Upregulated PCSK9 was associated with a lower PFS in the NB cohort, which was consistent with our database survival analysis (HR = 1.51, Figure 10C ). Our findings showed that PCSK9 was negatively associated with NK cell infiltration in NB, implying a potential role of PCSK9 modulation in the NB microenvironment.