The gene expression data of LEPROT and clinical follow-up information in pan-cancer atlas studies were extracted from The Cancer Genome Atlas (TCGA) database (Tomczak et al., 2015) using the cBio cancer genomics portal (cBioPortal) cBioPortal server (https://www.cbioportal.org/) for further analysis (Cerami et al., 2012; Gao et al., 2013). Apart from the expression profiles of TCGA, the Genotype-Tissue Expression (GTEx) database, which builds on molecular profiling of nondiseased tissue sites; the NCBI’s Gene Expression Omnibus (GEO), which stores curated gene expression data sets and might serve as a complement of TCGA data sets; and the European Genome-phenome Archive (EGA) are also included in the current study. The abbreviations of the types of cancers examined in our analyses are as follows: ACC: adrenocortical carcinoma; BLCA: bladder urothelial carcinoma; BRCA: breast invasive carcinoma; CESC: cervical squamous cell carcinoma; CHOL: cholangiocarcinoma; COAD: colon adenocarcinoma; CRC: colorectal carcinoma; DLBC: lymphoid neoplasm diffuse large B cell lymphoma; ESCA: esophageal carcinoma; GBM: glioblastoma multiforme; LGG: brain lower grade glioma; HNSC: head and neck squamous cell carcinoma; KICH: kidney chromophobe; KIRC: kidney renal clear cell carcinoma; KIRP: kidney renal papillary cell carcinoma; LAML: acute myeloid leukemia; LIHC: liver hepatocellular carcinoma; LUAD: lung adenocarcinoma; LUSC: lung squamous cell carcinoma; MESO: mesothelioma; OV: ovarian serous cystadenocarcinoma; PAAD: pancreatic adenocarcinoma; PCPG: pheochromocytoma and paraganglioma; PRAD: prostate adenocarcinoma; READ: rectum adenocarcinoma; SARC: sarcoma; SKCM: skin cutaneous melanoma; STAD: stomach adenocarcinoma; TGCT: testicular germ cell tumors; THCA: thyroid carcinoma; THYM: thymoma; UCEC: uterine corpus endometrial carcinoma; UCS: uterine carcinosarcoma; and UVM: uveal melanoma; STES: stomach and esophageal carcinoma cohort; KIPAN, Pan-kidney cohort (KICH + KIRC + KIRP). Plots in the study were conducted in Graphpad Prism 9.0, R version 4.0.5, or downloaded from online servers and then compiled and edited by Affinity Designer 1.8.5. All data or analyses are in accordance with the latest version in July 2021.

The Human Protein Atlas (HPA, https://www.proteinatlas.org/) is a comprehensive resource for the exploration of the human proteome, and it contains information of mRNA expression and protein level in cells, tissues, and organs (Uhlén et al., 2015). The LEPROT expression level in normal tissues was achieved from HPA based on the GTEx database. An overview plot to compare LEPROT mRNA expressions in tumor tissues and corresponding normal tissues across TCGA cancer types was achieved by the “Cancer Exploration” > “Gene_DE” module of TIMER2.0 (http://timer.comp-genomics.org/) (Li et al., 2020). TIMER is currently using the expression data of tumor tissues and the corresponding normal tissues in TCGA of the latest version. Differences in LEPROT mRNA expression between the tumor and the corresponding normal tissues were also shown by boxplots generated by the UALCAN (http://ualcan.path.uab.edu/analysis.html), a web portal for analyses, including expression distributions according to data from the TCGA and Clinical Proteomic Tumor Analysis Consortium (CPTAC) data set (Chandrashekar et al., 2017; Chen et al., 2019a). Expression of LEPROT (presented as log2 [counts per million (CPM)]) among molecular subtypes across pan-cancer were checked by the “Subtype” > “Molecular subtypes” module of TISIDB portal (http://cis.hku.hk/TISIDB/index.php) (Ru et al., 2019) with all 17 available cancer types analyzed, including ACC, BRCA, COAD, ESCA, GBM, HNSC, KIRP, LGG, LIHC, LUSC, OV, PCPG, PRAD, READ, SKCM, STAD, and UCEC, and those significant ones presented (p-value < .05). Mutation features of LEPROT in the TCGA pan-cancer samples were analyzed using the cBioPortal. “TCGA Pan Cancer Atlas Studies” was selected before “LEPROT” was queried in “Query by gene” to get the results of the gene alteration patterns across pan-cancer types in the “Cancer Types Summary” module.

LinkedOmics (http://www.linkedomics.org/) is a web portal for comprehensive multigroup data analysis for multiple cancer types (Vasaikar et al., 2018). The relationship of LEPROT DNA methylation (methylation preprocessor, presented as beta values (value-0.5) with its mRNA expression (presented as RNA-Seq by Expectation-Maximization (RSEM), Log2 (value + 1)) across cancer types was analyzed using the LinkedOmics with Methylation450 and high-output RNA sequencing (HiSeq) data. Pearson correlation was conducted for a correlation coefficient R and a p-value for each cancer type.

Partial correlation analysis between LEPROT and expression or mutation status of certain genes was conducted using the “Exploration” > “Gene_Corr”/“Gene_Mutation” module of TIMER2.0. The heatmap of the correlations shows the coefficient of Spearman correlation with red indicating a statistically significant positive association and blue indicating a statistically significant negative association (significant when p < .05); gray denotes a nonsignificant result. The scatterplots were based on log2 [TPM]. The cBioportal was used to identify the top 500 LEPROT-correlated genes and to generate scatterplots of the correlations of LEPROT and genes of interest across cancer types, via selecting each “PanCancer Atlas” study before “LEPROT” was queried in “Query by gene” to obtain the top correlated genes with LEPROT in the “coexpression” module. The intersection of the top 500 correlated genes among different cancer types was obtained and visualized using the Upset tool by R package “UpSetR” (Lex et al., 2014).

The stromal and immune scores were obtained by the Estimate algorithm (https://bioinformatics.mdanderson.org/estimate/) (Yoshihara et al., 2013). A higher stromal or immune score indicates a larger number of stromal or immune components; the ESTIMATE score is the sum of the stromal and immune scores. The relationship between stromal/immune/ESTIMATE scores and LEPROT expression is indicated by the Pearson correlation coefficient and visualized by R package “ggplot2” (Wickham, 2016). The relationship between LEPROT and infiltrating CD4⁺, CD8⁺, regulatory T cells, and CAFs are demonstrated by the “Immune” > “Gene” module of TIMER2.0 with results conducted by different algorithms, including TIMER (Li et al., 2016), xCell (Aran et al., 2017), MCP-counter (Becht et al., 2016), CIBERSORT (Newman et al., 2015), EPIC (Racle et al., 2017), and quanTIseq (Finotello et al., 2019). Overall correlations of LEPROT expression and gene signatures (including signatures of specific immune cells or CAFs) using expression data from all cancer samples were generated by the “Correlation Analysis” module of the web server GEPIA2 (http://gepia2.cancer-pku.cn/) (Tang et al., 2019).

Gene Set Enrichment Analysis (GSEA) software (Subramanian et al., 2005) was established to determine the statistically significant gene sets in the comparison between different subgroups. GSEA (Version 4.10) was obtained from the Broad Institute (http://software.broadinstitute.org/gsea/index.jsp). Expression profiles of cancer studies were pre-ranked according to the Pearson correlation coefficient with LEPROT and input with the annotation file “H: hallmark gene sets” and “C5: gene ontology sets”. The cutoff values were predefined as FDR <0.25. GeneMANIA (http://www.genemania.org) (Warde-Farley et al., 2010) was used to construct protein–protein interaction networks of certain proteins and related proteins and enrich the co-regulated genes into functions. The proteins were connected based on “physical interactions,” “co-expression,” “predicted co-localization,” “genetic interactions,” and “pathway shared” protein domains.”

Survival analysis was conducted using R package “ggplot2” (Wickham, 2016) to get forest plots showing the hazard ratio and p-value for overall survival (OS) or disease-specific survival (DSS) in patients across cancer types in TCGA, based on LEPROT high versus low expression, in the Sangerbox tools, a free online platform for data analysis (http://www.sangerbox.com/tool). Kaplan–Meier analysis of OS or recurrence-free survival (RFS) was obtained using the Kaplan–Meier plotter (KM-plotter) (www.kmplot.com) (Nagy et al., 2021), using “mRNA” > “Start KM Plotter for pan-cancer” module with the setup of “split patients by auto select best cutoff” to divide LEPROT mRNA into high and low groups.

Gene expression profile and corresponding clinical responses to immune checkpoint inhibitors (CPIs) in GSE67501 and GSE79691 data sets were downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/gds). The associations of LEPROT expression and patient responses to CPIs were investigated by comparing LEPROT-low (lower than the median) and LEPROT-high (higher than the median) groups using the Mann–Whitney U test. A simple logistic regression was also conducted to generate a receiver operating curve (ROC). Graphpad Prism 9.0 was used to conduct statistics.

As LEPROT was expressed widely in human tissues (Uhlén et al., 2015), we first explored the level of LEPROT in different tissues, especially in normal and tumor tissues across organs. By exploring the HPA and GTEx databases, LEPROT was found to be highly expressed in female genital organs, including the endometrium, cervix, and placenta and also in the kidney (Supplementary Figure S1A). The mRNA level of LEPROT was various in different tissues with relatively low tissue and low cell-type specificity (Supplementary Figure S1A).

Moreover, a comparison between tumor tissues and adjacent normal tissues showed the LEPROT mRNA levels were significantly lower in the tumor tissues of BLCA, BRCA, COAD, KICH, KIRC, KIRP, LUAD, LUSC, PRAD, READ, THCA, and UCEC (p < .001) and higher in CHOL (p < .001), GBM (p < .01), and HNSC (p < .05) (Figure 1A) compared with the corresponding normal tissues. Thirteen cancers (ACC, DLBC, LAML, LGG, MESO, OV, PCPG, SARC, SKCM, TGCT, THYM, UCS, UVM, UCEC, and PCPG) were removed from the comparison analysis due to a lack of corresponding normal tissue or a normal tissue sample of less than three cases. Similar results were also observed in the UALCAN (Supplementary Figure S1B). Further evaluation of LEPROT mRNA levels was performed in different molecular subtypes across pan-cancer types. LEPROT was found to be distinctly expressed in different molecular subtypes of LGG, BRCA, OV, HNSC, STAD, LUSC, PRAD, and PCPG (p < .05) (Figure 1B). Specifically, LEPROT was expressed in a significantly lower level in the proliferative subtype of OV (p < .05) and in a significantly higher level in the mesenchymal subtype of HNSC (p < .05) (Figure 1B), suggesting that LEPROT might be related to more aggressive subtypes.

To investigate whether the alterations of expression or functional roles of LEPROT resulted from its mutations, the mutation feature of LEPROT in TCGA pan-cancer atlas was examined and demonstrated in Supplementary Figure S2A. The pan-cancer mutation spectra of LEPROT indicated a low mutant frequency of LEPROT across all cancer types with the maximal alteration frequency of LEPROT in SARC (only 3.14%). There was no significant difference between the expression of LEPROT in wild-type and mutated cases except for the limited LEPROT-mutated cases (n = 2) in COAD that had a lower level of LEPROT compared with LEPROT wild-type cases (p < .05) (Supplementary Figure S2B). These results indicate LEPROT was aberrantly expressed in cancer tissues, but it did not result from its mutations in most cancer types.

To interpret the cause of low LEPROT expression across cancer types, DNA methylation of LEPROT, the epigenetic modifications that altered gene expression, were investigated. First, correlations of LEPROT expression and methylation level were examined by LinkedOmics. It revealed LEPROT methylation reduced its expression in most cancer types except for ACC, BLCA, HNSC, KICH, MESO, SKCM, and UCEC (Supplementary Figure S3). The negative correlations between the expression of LEPROT and its methylation were most significant in OV (Pearson R = -0.745, p = .02) and CHOL (Pearson R = -0.523, p = .001). In addition, as key factors in regulating DNA methylation levels (Cui and Xu, 2018), the correlation of the DNA methyltransferase (DNMT) family (DNMT1, DNMT2, DNMT3A, DNMT3B) and LEPROT expression were also investigated. Unexpectedly, the expression of LEPROT was positively correlated with that of DNMT1 and DNMT3A but not DNMT3B in most cancers (Figure 1C) and was consistently positive-correlated with DNMT2 in all cancers. Therefore, we further examined the correlation between DNMTs and methylation level of LEPROT, and the results show no significant correlations, represented by DLBC and OV (Supplementary Figure S3B, C). These results suggest that the aberrant expression of LEPROT might be regulated by DNA methylations, but the methylation level of LEPROT was regulated independent of DNMTs.

Tumorigenesis is often driven by mutation or compromising of tumor-suppressor genes, which restricted the cell cycle and repaired DNA (Lee and Muller, 2010; van der Weyden et al., 2012). Although, generally, mutations of tumor suppressor genes did not result in the alteration of LEPROT expression (Supplementary Figure S4A), the expression of tumor-suppressor genes is shown to be highly related to LEPROT expression. Tumor-suppressor genes BAX, CDKN2A, and MSH5 were negatively correlated with LEPROT in multiple cancer types, and most tumor-suppressor genes were suggested to be positively correlated with LEPROT in a gene-consistent way across cancer types (Figure 1C). Notably, as tumor-suppressor genes, mismatch repair (MMR) genes constitute the MMR system that recognizes and repairs DNA mistakes (Marti et al., 2002; Putnam, 2016), six of which (MLH1, MLH3, MSH2, MSH3, MSH6, and PMS2) was overall positively correlated with LEPROT in mRNA expression. These results suggest that the aberrant expression LEPROT may be accompanied by abnormal tumor suppressor gene functions, especially compromised DNA mismatch repair.

Oncogenes are genes that are abnormally expressed and have the potential to drive cancer. To understand the role of the aberrantly expressed LEPROT in cancer development, we explored the correlation between the mRNA expression level of LEPROT and commonly reported oncogenes across various cancer types by TIMER2. As shown in Figure 1C, LEPROT expression was positively correlated with most oncogenes in a broad spectrum of cancer types, suggesting a role of LEPROT in promoting cancer progression. However, mutations of oncogenes did not result in the alteration of LEPROT expression (Supplementary Figure S4B), which indicated that mutations of those oncogenes did not drive expression alterations of LEPROT.

As a member of the leptin signaling pathway, LEPROT was not only involved in the regulation of intracellular signaling pathways but, more importantly, was involved in the interaction of tumor cells with the extracellular TME, which played a regulatory role in tumor progressing through reversible changes in signal transduction or gene expression programs. Thus, we next estimated the association of LEPROT and TME components. Stromal score, which was calculated based on gene signature, was a direct quantitative indicator of stromal cell infiltration in the microenvironment. Thus, the relationship between LEPROT expression and stromal scores was first estimated. Pearson correlation analysis showed LEPROT expression was positively correlated with the stromal score in 23/33 cancers (69.7%; they were BLCA, BRCA, CESC, COAD, DLBC, ESCA, GBM, HNSC, KICH, KIRC, LGG, LIHC, LUAD, LUSC, OV, PAAD, PCPG, PRAD, READ, SARC, STAD, TGCT, and UCEC (p < .05)) without any significant negative correlation among TCGA cancer types (Figure 2A).

To better understand which stromal components were positively associated with LEPROT, we checked the relationship between LEPROT mRNA level and immune cell infiltration (TII). As the most important compositions of the stromal, the ratio of TII reflected by immune score was calculated. It was positively correlated with LEPROT expression in COAD, DLBC, GBM, HNSC, KICH, LGG, LIHC, LUAD, LUSC, OV, PAAD, PCPG, PRAD, READ, STAD, and UVM and negatively associated with LEPROT only in TGCT and KIRP (Figure 2A). These results aroused our interest because the composition of TII cells could serve as biomarkers for predicting patient responses to treatment and survival in terms of chemotherapy and immunotherapy (Zhang et al., 2019). As the immune score reflecting TII as a whole, the composition of TII cells was further examined, especially the distinct T cells in cancers. TIMER2, which offers many existing algorithms (Li et al., 2020; Aran et al., 2017; Newman et al., 2015; Racle et al., 2017; Li et al., 2017) for estimating tumor-infiltrating immune cell populations, showed a remarkable correlation between LEPROT expression and CD4⁺, CD8⁺, and regulatory T cells (Figure 2B). Generally, the correlation was presented in a type-dependent manner and algorithm-dependent way. But the positive relationship of LEPROT and memory CD4⁺ T cells resting, and the negative relationship of LEPROT and Th1 seemed to be strong and consistent across most cancer types. Moreover, CAFs, a prominent component of the microenvironment in most types of solid tumors, which are also considered to be involved in TII (Barrett and Puré, 2020), were also shown to consistently positively correlate with LEPROT expression regardless of the cancer types (Figure 2B). The LEPROT mRNA level was overall positively correlated with the signatures of effect T cells, resident memory T cells, Th1 cells, and CAFs across all cancer types using whole TCGA cancer data sets (Figure 2D).

Considering that LEPROT was highly associated with TII, we next investigated its association with expressions of immune checkpoint molecules (ICMs). ICMs and the composition of the TME are shown to influence each other intensively and both influence the effect of immune therapy (Jiang et al., 2019; Wang et al., 2019a; Wang et al., 2019b; Zhang et al., 2020b; Lee et al., 2019). ICMs, including programmed death 1 (PD-1), programmed cell death ligand 1 (PD-L1), and cytotoxic T lymphocyte-associated antigen 4 (CTLA-4), which played key roles in immune checkpoint inhibition, were included in the correlation analysis. As shown in Figure 2C, PD-L1 was positively correlated with LEPROT in all cancer studies, and PD-1 and CTLA-4 were limited related to LEPROT (Figure 2C). Apart from ICMs, cumulative evidence suggests profound and complex roles of toll-like receptor (TLR) signals in TME (Angrini et al., 2020). And TLR1-10 showed an overall positive correlation with LEPROT expression (Figures 2C,E). All results above suggest that LEPROT was involved in the interactions of tumor cells with TME, especially with TII and CAFs.

To gain mechanistic insight into LEPROT, GSEA was performed on five TCGA studies (BLCA, BRCA, COAD, KIRC, and STAD). Pre-ranked lists of genes according to the Pearson correlation coefficient with LEPROT were adjusted for the investigation of enriched genes involving hallmarks, gene ontology (GO), and pathways. Not surprisingly, inflammatory/immune response markers and extracellular component-related pathways were enriched in these cancer types. Items significantly enriched in hallmarks were inflammatory or immune response, including IL6_JAK_STAT3 signaling, inflammatory response, complement, TNF-alpha signaling via NFKB, and IL2_STAT5 (Table 1). Items significantly enriched in GO analysis were extracellular matrix structural constituent, collagen binding of molecular functions (MF), cell matrix adhesion, cell substrate adhesion, cell substrate junction organization, and regulation of cell substrate junction organization in biological process (BP) (Table 2). These enrichment items once again illustrated the interactions of LEPROT with the TME, especially with the extracellular matrix and immune response.

To further identify the key genes that interacted with LEPROT in tumorigenesis and tumor progression, LEPROT co-regulated genes were examined using the cBioportal. Intersection analyses were performed among 28 non-germ cell solid cancers with the top 500 LEPROT-related genes (Supplementary Table S1, Figure 3A), from which IL6ST turned out to be the only intersected gene. Pearson correlation analysis was performed on the level of LEPROT and the top 10 intersected genes. The genes JAK1, JAK2, and STAT3 were also included in the analysis as they were the downstream molecules of IL6ST (Sanz-Moreno et al., 2011; Looyenga et al., 2012; Pilati and Zucman-Rossi, 2015) (Figure 3B). Results show that the level of LEPROT was consistently positively correlated with that of JAK1, JAK2, and STAT3 (Figure 3B), which suggested that JAK1/2/STAT3 might be the downstream LEPROT and mediate cancer cell proliferation and TME alterations. Moreover, to better understand the way they interacted with each other, GeneMANIA was used to construct protein–protein interaction networks of proteins encoded by LEPROT, genes highly co-regulated with IL6ST and their correlated proteins (Figure 3C). Functional enrichment analysis for these proteins suggest they were enriched in functions involving “functions cytokine receptor binding,” “cellular response to interleukin-6,” “tyrosine phosphorylation of STAT protein,” “receptor signaling pathway via JAK-STAT,” and “receptor signaling pathway via STAT.” These correlation analyses suggest that the interactions of LEPROT with the TME might be mediated by inflammatory signaling pathway IL6ST/JAK1/2/STAT3 (Figure 3C).

Prognostic markers were biomarkers that were heterogeneously expressed in cancer samples with distinct patient outcomes. Identifying prognostic markers would help to group patients and guide precise drug discovery. Common prognostic markers were tumor suppressor, proto-oncogene, and TME components. Considering LEPROT was highly correlated with these common prognostic markers, we were curious if LEPROT could also provide prognostic value in patients. Thus, Cox analyses were performed for OS and DSS of LEPROT across pan-cancer types. The forest plots indicate that higher LEPROT was linked to worse OS and DSS in PAAD, CESC, LGG (p < .0001), and KICH (p < .05), worse OS in STAD (p < .05), worse DSS in BLCA (p < .05), and longer OS and DSS in SKCM (p < .05) and KIRC (p < .0001) (Figures 4A,B). Kaplan–Meier analysis showed an association of high LEPROT expression with better prognosis (both OS and PFS) in HNSC (p < .01), PAAD (p < .05), and STAD (p < .01), longer OS in KIRC (median OS HR = 0.54, p = 6e-5) and THCA (p < .05), and longer DFS in UCEC (p < .01) (Figures 4C,D). These indicate LEPROT might serve as potential prognostic markers for several cancer types.

Considering that the LEPROT expression correlated with TII and immune checkpoints, we investigated whether it could imply patient responses to immune checkpoint inhibitors (CPIs). In KIRC (GSE67501), patients with higher LEPROT (above the median expression) had a significant (p = .022) lower response rate (0%, 0 in 5) to the PD-1 inhibitor nivolumab than those with low LEPROT (67%, 4 in 6) (Figure 4E, left panel). Although in SKCM (GSE79691), patients with higher LEPROT had a response rate of 80% (4 in 5) compared with a response rate of 40% (2 in 5) for patients with low LEPROT (Figure 4F, left panel), there was no significance between the two groups. The predictive role of LEPROT expression in CPI therapy was shown by ROC curves, and the area under the curve (AUC) was 0.875 in GSE67501 (Figure 4E, right panel) and 0.708 in GSE79691 (Figure 4F, right panel). These indicate that the LEPROT may influence CPI therapy response in different directions and has the potential to predict patient responses to CPI therapy.