ZBED6–IGF2–PIK3C3 autophagy axis drives ccRCC progression: A multi-omics integration study

TL;DR

This study identifies a new regulatory pathway involving ZBED6, IGF2, and PIK3C3 that drives kidney cancer progression and offers a six-gene model for predicting patient outcomes.

Contribution

The discovery of the ZBED6–IGF2–PIK3C3 regulatory axis and its role in autophagy-driven ccRCC progression.

Findings

ZBED6 activates PIK3C3 via repression of IGF2 to promote pro-tumorigenic autophagy.

A six-gene prognostic signature stratifies ccRCC patients with distinct survival outcomes.

ZBED6 promotes ccRCC cell proliferation, migration, and invasion in functional assays.

Abstract

Autophagy supports clear cell renal cell carcinoma (ccRCC) progression, yet its upstream regulatory mechanisms remain to be fully defined. Integrating bulk, single-cell, and spatial transcriptomics, we identify a regulatory axis wherein the transcription factor ZBED6 activates the expression of the autophagy-initiating kinase PIK3C3 via the repression of IGF2, thereby driving pro-tumorigenic autophagy. Spatial analysis confirms the co-localization of ZBED6 and PIK3C3 in tumor tissues. Using genes associated with this axis, we develop a six-gene prognostic signature that stratifies patients with distinct survival outcomes and differential responses to immunotherapy and targeted therapy. Functional assays show that ZBED6 promotes ccRCC cell proliferation, migration, and invasion. This work elucidates a pathway governing autophagy in ccRCC and provides a framework for prognostic assessment and precision therapy.

•Identification of the ZBED6-PIK3C3 regulatory axis in ccRCC•Discovery of the IGF2-mediated ZBED6-PIK3C3 regulatory axis•Analysis based on multi-omics data•Construction of a 6-gene prognostic model based on PIK3C3

Identification of the ZBED6-PIK3C3 regulatory axis in ccRCC

Discovery of the IGF2-mediated ZBED6-PIK3C3 regulatory axis

Analysis based on multi-omics data

Construction of a 6-gene prognostic model based on PIK3C3

Renal cancer; Omics; Transcriptomics; Bioinformatics; Machine learning

Figures

The technical route and logical sequence of the research

Identify the potential regulatory relationship between ZBED6 and PIK3C3(A) Venn diagram shows the intersection of differentially expressed genes (DEG) analysis and PIK3C3 highly correlated gene collections.(B) Correlation analysis between ZBED6 and PIK3C3.(C) Survival curves of ZBED6 and PIK3C3. The correlation between ZBED6 and PIK3C3 and age (D), survival status (E), gender (F), T stage (G), M stage (H), and stage (I) was analyzed, respectively.(J) Heatmap shows the correlation between high and low expression of ZBED6 and PIK3C3 and clinical features.(K) Boxplot shows the expression of ZBED6

Preliminary verification of the regulatory relationship between ZBED6 and PIK3C3(A) Scatterplot shows the correlation between ZBED6 and PIK3C3 in ArrayExpress (left) and ICGC-RECU (right).(B) Bar graph illustrates the expression levels of ZBED6 and PIK3C3 across different cell lines.(C) Western blot (WB) analysis validates differential protein expression of ZBED6 and PIK3C3 in renal tubular epithelial cell lines and clear cell renal cell carcinoma (ccRCC) cell lines.(D) WB analysis demonstrates changes in ZBED6 and PIK3C3 expression in the 769-P cell line following ZBED6 knockdown and overexpr

Joint verification using single-cell and spatial transcriptomics(A) The t-SNE plot shows the 10 cell types identified by the marker genes.(B) Expression distribution of ZBED6 in t-SNE plot. (C) Expression distribution of PIK3C3 in t-SNE plot.(D) Distribution of PIK3C3-related pathway scores in t-SNE plot (whole cell).(E) Boxplot shows the distribution of PIK3C3-related pathway scores across 10 cell types.(F) Distribution of PIK3C3-related pathway scores (tumor cells) under the t-SNE plot.

Identify the set of prognostic-related genes closely associated with PIK3C3(A) Hierarchical clustering treemap. The bottom heatmap represents the PIK3C3-related pathway score for each sample.(B) Cluster dendrogram of WGCNA analysis.(C) Module-trait heatmap shows that the blue module is closely related to PIK3C3-related pathway score.(D) Scatterplot of gene significance (GS) and module membership (MM) in the blue module. (E) Venn diagram shows the intersection genes of the renal cancer differential gene and the blue module gene set.(F) Univariate Cox regression analysis of intersecting genes an

Construction and validation of machine learning prognosis models(A) A total of 101 prediction models was constructed through a 10-fold cross-validation framework, and the C-index of each model on all validation datasets was further calculated. Kaplan-Meier curves for TCGA cohort (B), ArrayExpress cohort (C), TCGA DSS cohort (D), and RECA-cohort (E). ROC curves for ArrayExpress cohort (F), TCGA DSS cohort (G), and RECA cohort (H).

Validation of the clinical relevance of core genes(A) Heatmap shows the expression of six model genes between high and low risk groups and the clinical data of patients with kidney cancer. Boxplots (B and C) show the correlation of different T-stages with risk score.(D) Comparison of the distribution of clinical data between the two groups of high and low risk.(E) Comparison of the distribution of T stages between the two groups at high and low risk.(F) Diagnostic ROC curves of prognostic models in predicting metastasis of kidney cancer.(G) Kaplan-Meier curve in survival differences between th

The transformation of clinical prognosis prediction modelsUnivariate (A) and multivariate (B) Cox regression analyses for prognostic models.(C) Nomogram, calibration curve, and DCA curve.(D) C-index of nomogram and other clinical data for 1–10 years survival prediction.

Enrichment analysis of PIK3C3-related prognostic genes(A) Analysis of the characteristic gene set in the low-risk GSVA group.(B) Analysis of the characteristic gene set in the GSVA high-risk group.(C) Correlation of risk score with 50 signaling pathways.(D) Risk score and the relationship between the signaling pathways are plotted, interacting. K-M curve analysis of protein secretion (E) and P53 pathway (F). ∗ indicates *p* < 0.05, ∗∗ indicates *p* < 0.01, ∗∗∗ indicates *p* < 0.001, and ∗∗∗∗ indicates *p* < 0.0001.

Immune infiltration analysis of PIK3C3-related prognostic genes(A) Immune scores, matrix scores, and estimated scores of the high-risk group and the low-risk group.(B) Immune cell infiltration in the high- and low-risk groups was scored.(C) Correlation between 6 model genes and 22 immune cells.(D) Correlation analysis of risk score with 22 immune cells. ∗ indicates *p* < 0.05, ∗∗ indicates *p* < 0.01, ∗∗∗ indicates *p* < 0.001, and ∗∗∗∗ indicates *p* < 0.0001.

Immune response analysis of prognosis-related genes related to PIK3C3(A) Distribution of CR/PR and PD/SD in patients treated with atezolizumab between high- and low-risk groups.(B) Correlation between risk score and CR/PR and PD/SD.(C) Correlation between risk score and CR, PR, PD, and SD. CR complete response: all target lesions disappeared, no new lesions appeared, and tumor markers were normal for at least 4 weeks. PR: partial response was achieved by ≥ 30% reduction in the sum of the maximum diameters of the target lesions for at least 4 weeks. SD: stable disease, the sum of the maximum di

Predictive analysis of potential targeted drugs for PIK3C3-related prognostic genes(A) Drug susceptibility prediction of 8 targeted drugs.(B) The top five drugs calculated by the XSum analysis all had the potential to reverse the phenotype of high-risk populations.(C and D) Results from Spearman’s correlation analysis and differential drug response analysis for 4 CTRP-derived compounds and 8 PRISM-derived compounds.

Experimental verification of ZBED6-IGF2-PIK3C3 in regulating the phenotype of renal cancer(A) The CCK-8 assay shows the proliferation viability of renal cancer cells after ZBED6 knockout or overexpression compared with the NC group.(B and C) The wound healing assay and Transwell migration assay showed the migration ability of renal cancer cells treated with NC or ZBED6 overexpression. The image is magnified 40X.(D) The Transwell invasion assay showed the invasion ability of renal cancer cells treated with NC or ZBED6 overexpression. The image is magnified 40X.(E) The heatmap of qPCR experiment

Prediction of drugs with high binding potential for ZBED6 through molecular docking(A) Drugs with potential binding to ZBED6 and their *p*-value ranking.(B and C) The molecular docking diagram shows the best predicted binding conformation of BML- 277/Tubastatin A and ZBED6 protein.

Full text

Introduction

Renal cell carcinoma (RCC) remains a lethal urological malignancy characterized by high metastatic potential and resistance to conventional therapies.1 Clear cell renal cell carcinoma (ccRCC) is the predominant pathological subtype of RCC, accounting for more than 80% of all cases.2 Despite advances in targeted and immune therapies, persistent challenges in overcoming treatment resistance underscore the need for a deeper understanding of the molecular mechanisms driving ccRCC progression.3 Autophagy, a conserved lysosomal degradation pathway, plays a critical and context-dependent role in cancer biology.4 In ccRCC, accumulating evidence indicates that autophagy functions as a pro-survival mechanism, enabling tumor cells to endure metabolic stress, chemotherapy, and targeted agents, thereby contributing to disease aggressiveness and therapeutic resistance. However, the precise regulatory networks governing autophagy initiation in ccRCC, particularly the upstream signaling axes, remain incompletely elucidated.5

The initiation of autophagy is tightly regulated by the class III phosphatidylinositol 3-kinase (PIK3C3) complex, with PIK3C3 serving as its indispensable catalytic core.6 PIK3C3 generates phosphatidylinositol 3-phosphate, a lipid second messenger essential for autophagosome nucleation.7 Although PIK3C3 is recognized as a major upstream regulator of autophagy induction, the transcriptional and epigenetic mechanisms governing its expression and activity in ccRCC remain poorly defined. Identifying key upstream regulators of PIK3C3 may reveal novel therapeutic vulnerabilities for disrupting pro-tumorigenic autophagy in ccRCC and other cancers.8 Leveraging multi-omics integration, we hypothesized that combining bulk transcriptomic, single-cell, and spatial transcriptomic data could uncover a novel regulatory axis of autophagy initiation in ccRCC. To this end, we systematically analyzed TCGA-KIRC: https://portal.gdc.cancer.gov RNA-seq datasets to identify upstream regulatory factors whose expression correlated with autophagy activity and ccRCC progression. To further characterize cell type-specific expression patterns within the tumor microenvironment (TME) and map the spatial co-localization of candidate regulators and autophagy markers, we re-analyzed single-cell RNA sequencing (scRNA-seq) and spatial transcriptomic data from GEO. This integrated computational approach identified ZBED6 (zinc finger BED-type containing 6), a transcription factor with emerging yet incompletely characterized roles in cancer, as a top upstream regulator strongly associated with PIK3C3 expression. ZBED6, known for its DNA-binding and transcriptional regulatory functions, has been implicated in cell growth.9 Previous studies have shown that it influences cancer cell growth and invasion in brain glioma and colon cancer, but its association with autophagy regulation in ccRCC has not been previously reported.10^,^11 Meanwhile, our results clearly demonstrated that IGF2, a downstream target of ZBED6, is significantly inhibited by ZBED6 in renal cell carcinoma. Moreover, IGF2 exerts a notable, indirect inhibitory effect on the PI3K pathway, a key regulator of autophagy, suggesting that IGF2 may play a crucial role in the ZBED6–PIK3C3 regulatory network.12

In this study, we provide compelling evidence establishing a novel ZBED6–IGF2–PIK3C3 regulatory axis that critically governs autophagy initiation in ccRCC. We demonstrate that ZBED6 positively regulates PIK3C3 expression and promotes the initiation of autophagy in tumor cells. Importantly, our results show that autophagy is indispensable for ccRCC cell survival, and the ZBED6–PIK3C3 axis serves as a major conduit for maintaining this pro-tumorigenic process (Figure 1). These findings identify ZBED6 as a previously unrecognized key regulator of autophagy initiation through its indirect control of PIK3C3. The discovery of the ZBED6–PIK3C3 axis not only enhances our understanding of autophagy dysregulation in ccRCC but also underscores its potential as a therapeutic target for disrupting the pro-survival autophagy network in this refractory malignancy.Figure 1. The technical route and logical sequence of the research

Results

Screening and correlation analysis of class III phosphatidylinositol 3-kinase upstream genes

First, four candidate genes were identified based on differentially expressed gene (DEG) analysis with an FDR <0.05 and |log_2_FC| > 1 in the TCGA-KIRC cohort, combined with the PIK3C3 high-correlation gene set, which was also found in the TCGA-KIRC cohort. After taking the intersection of the two screening gene sets, we found 4 common genes: AC108449.2, AC018752.1, NBEAL1, and ZBED6. ZBED6, considered a repressor clearly associated with tumor disease progression, was identified as the target gene (Figure 2A). The scatterplot demonstrated a positive correlation between ZBED6 and PIK3C3 (Figure 2B) in TCGA-KIRC. Survival analysis indicated that elevated expression levels of ZBED6 and PIK3C3 were associated with a favorable prognosis, suggesting their potential roles as protective genes in patients with kidney cancer (Figure 2C). To further explore the relationship between ZBED6 and PIK3C3 expression and clinical features, boxplots were generated to depict their associations with age, survival status, sex, T stage, M stage, and stage I vs. stage IV (Figures 2D–2I), further confirming that ZBED6 and PIK3C3 act as protective factors in kidney cancer. A statistically significant difference was observed between the early disease state (T1 and Stage I) and the advanced disease state. A heatmap was generated to visualize the correlations between ZBED6 and PIK3C3 expression levels and various clinical parameters (Figure 2J). In addition, boxplots comparing tumor and normal tissues revealed that ZBED6 waw expressed at lower levels in tumor tissues than in normal tissues, but PIK3C3 was not (Figure 2K).Figure 2. Identify the potential regulatory relationship between ZBED6 and PIK3C3(A) Venn diagram shows the intersection of differentially expressed genes (DEG) analysis and PIK3C3 highly correlated gene collections.(B) Correlation analysis between ZBED6 and PIK3C3.(C) Survival curves of ZBED6 and PIK3C3. The correlation between ZBED6 and PIK3C3 and age (D), survival status (E), gender (F), T stage (G), M stage (H), and stage (I) was analyzed, respectively.(J) Heatmap shows the correlation between high and low expression of ZBED6 and PIK3C3 and clinical features.(K) Boxplot shows the expression of ZBED6 and PIK3C3 in tumor tissues and normal tissues (normal = 72, tumor = 539). ∗ indicates p < 0.05, ∗∗ indicates p < 0.01, ∗∗∗ indicates p < 0.001, and ∗∗∗∗ indicates p < 0.0001.

In vitro validation of ZBED6 and class III phosphatidylinositol 3-kinase expression levels and their correlation with the autophagy pathway

Analysis of external ArrayExpress: https://www.ebi.ac.uk/biostudies/arrayexpress and ICGC: https://dcc.icgc.org databases confirmed a significant positive correlation between ZBED6 and PIK3C3 expression (Figure 3A). Evaluation of ZBED6 and PIK3C3 expression patterns across various renal carcinoma cell lines using the CCLE: https://sites.broadinstitute.org/ccle database demonstrated consistent co-expression (Figure 3B). Western blot analysis revealed markedly elevated PIK3C3 expression in 769-P, Caki-1, and ACHN cell lines compared with the normal renal proximal tubular epithelial cell line HK-2. However, there was no such significant high expression of ZBED6 (p > 0.05, ns) (Figure 3C and Table 1).Figure 3. Preliminary verification of the regulatory relationship between ZBED6 and PIK3C3(A) Scatterplot shows the correlation between ZBED6 and PIK3C3 in ArrayExpress (left) and ICGC-RECU (right).(B) Bar graph illustrates the expression levels of ZBED6 and PIK3C3 across different cell lines.(C) Western blot (WB) analysis validates differential protein expression of ZBED6 and PIK3C3 in renal tubular epithelial cell lines and clear cell renal cell carcinoma (ccRCC) cell lines.(D) WB analysis demonstrates changes in ZBED6 and PIK3C3 expression in the 769-P cell line following ZBED6 knockdown and overexpression.(E) Subcellular localization of ZBED6 in cells based on the Human Protein Atlas (HPA) database.(F) Immunohistochemical (IHC) results from HPA show PIK3C3 expression in ccRCC tumor tissues versus normal renal tissues.(G) Altered expression of autophagy markers after ZBED6 overexpression. ∗ indicates p < 0.05, ∗∗ indicates p < 0.01, ∗∗∗ indicates p < 0.001, and ∗∗∗∗ indicates p < 0.0001.Table 1. Experimental orginal dataCCK80hNC10.1908NC20.1938NC30.1759OE10.2104OE20.2071OE30.2171SI10.1923SI20.2043SI30.199024hNC10.3625NC20.3667NC30.3452OE10.3741OE20.3885OE30.3714SI10.2926SI20.3155SI30.283348hNC10.5399NC20.5642NC30.5055OE10.5529OE20.5633OE30.5310SI10.4251SI20.4159SI30.417272hNC10.7084NC20.7128NC30.7476OE10.7591OE20.7760OE30.7734SI10.5473SI20.6017SI30.5679qPCR-1ZBED6-OE-8GENEZBED6VPS34ASAH1ATP1A1DSTNEIF1BPGK1SCP22-ΔΔT'10.9416960170.7973766884.2378523775.07473797578.4299803231.6324486526.784627614.5736099482-ΔΔT'21.0163049320.7755723813.714925644.29701059283.4785324633.4360688928.31183464.8009946672-ΔΔT'31.0448771530.8254961173.3480784524.66971121595.2292808933.9028190226.599612684.702191625ZBED6-OE-8GENEZBED6-OEVPS34-OEASAH1-OEATP1A1-OEDSTN-OEEIF1B-OEPGK1-OESCP2-OE2-ΔΔT'126.415875730.9159452911.418754033.514533809181.43806356.623669225.693532222.7959398682-ΔΔT'231.196955381.15936379111.183759473.211691527185.250462456.2325409427.537645992.5906845042-ΔΔT'323.479557811.02337389212.669900883.588381635201.318135864.1480423629.310249292.367448977Transwell-migrationCell numbernc1.tif796nc2.tif828nc3.tif781oe1.tif1116oe2.tif1085oe3.tif993Transwell-invasionCell numbernc1.tif713nc2.tif726nc3.tif871oe1.tif1167oe2.tif994oe3.tif985WB data-1ZBED6-ALL CELL LINESpro area valuegap area valuerationormalizeEXP1-HK-2166871312290.1871894111.040047638EXP1-769-P254881675140.2287663241.271054136EXP1-Caki-1222241325970.2054524781.141519509EXP1-ACHN251081007920.262531631.458658374EXP2-HK-2261281391390.2048483711.138163016EXP2-769-P346291619440.2109813381.172238543EXP2-Caki-1265371311810.1800607961.000440167EXP2-ACHN27202837050.2672863591.485076239EXP3-HK-2265491665350.147906940.821789346EXP3-769-P294621517920.1394281280.774680013EXP3-Caki-1189531026790.1124920620.625019881EXP3-ACHN19786576340.1717788220.95442449VPS34-ALL CELL LINESEXP1-HK-2418961312290.31925871.773841027EXP1-769-P1070481675140.6390391253.550580822EXP1-Caki-1710731325970.5360076022.978124872EXP1-ACHN672941007920.6676521953.709558596EXP2-HK-2557531391390.4007000192.22633912EXP2-769-P1274211619440.786821374.371677347EXP2-Caki-1671161311810.5116289712.842674167EXP2-ACHN57749837050.6899109973.833231267EXP3-HK-2511671665350.3072447231.707089877EXP3-769-P1200821517920.791095714.395426112EXP3-Caki-1576611026790.5615656563.120128597EXP3-ACHN33435576340.5801263143.22325392ZBED6-OENC176041130100.0672860810.609415501OE1389461411990.2758234832.498155688NC254042293260.0235647070.213427465OE2580262282590.2542112252.302411712NC3151892039770.0744642780.674429014OE3706522218590.3184545142.884268399VPS34-OENC1512182233580.2293090022.076870256OE11180752254220.523795374.74405721NC2237892053030.1158726371.049467885OE2765312274940.3364088723.046882477NC3270991783200.1519683711.376389879OE3905512251800.4021271873.642098594ZBED6-SINC1633681221190.5189036920.747969188SI1388491025450.378848310.546087582NC2999761240970.8056278561.161265226SI234868916580.3804141480.548344646NC3989621006610.9831215661.417111981SI325864884390.2924501630.421549729VPS34-SINC1739801583440.467210630.673456676SI1407751630160.2501288220.360546002NC2956721288390.7425701841.070371296SI2372431350930.2756841580.397382518NC31160521799070.6450666180.929825634SI3382641284970.2977812710.429234208ULK1NC1103133292590.0313218470.792547746NC2135173169580.0426460291.079087519NC3146943295100.0445934871.128364735OE1384203380370.1136561972.875882854OE2315053038590.1036829582.623526469OE3274482448260.1121122762.836816476NC162102329190.0266616290.614365278NC2122232444570.0500006141.1521667NC3119882239540.053528851.233468022OE1173432724420.0636575861.466865003OE2211912428070.0872750792.011084099OE3207472044810.101461752.337988283NC125362466320.0102825260.710344125NC241132853960.0144115540.995588313NC344312365450.0187321651.294067562OE195441901830.0501832453.466791453OE296142442640.0393590542.719027676OE366562017690.0329882192.278913523p62NC1585421058060.5532956541.016401335NC2797701646640.4844410440.889915762NC31185391991030.5953652131.093682903OE11174201697170.6918576221.270938973OE21181621628610.7255389571.332811416OE3736711071790.6873641291.26268445NC1710341567330.4532166170.714049037NC21185721606450.7380995361.16288601NC31459382047320.7128245711.123064953OE11084571574700.6887470631.085130507OE21016761703050.5970229880.940617961OE3672641612170.4172264710.657346066NC1518421433540.3616362291.070262674NC2515431434040.3594251211.063718896NC3471651611800.2926231540.86601843OE1502271506690.3333598820.986578803OE2581281698180.342295871.013024866OE3748861681070.445466281.318357769LC3NC123257978400.2377044150.614426141NC2575681282690.4488068041.160090495NC3773561631620.474105491.225483364OE12223681718401.2940409683.344879383OE22472771541181.604465414.147274623OE3185997894932.0783413235.372164569NC1649691064440.6103584980.91022081NC2835391205040.6932466971.033831055NC3773951093030.7080775461.055948134OE11835631534051.1965907241.78446238OE22217571185671.8703096142.789171834OE31493891381681.0812127271.612400461NC1875741184970.7390398070.868203549NC21179811337360.8821932761.03637629NC31273331365570.9324531151.09542016OE11938531358281.4271946871.676628893OE22210421566351.4111916241.657828936OE31855281248781.4856740181.745328794Wound healing assay area0h24hrateNC1214930020875920.028710743NC2239250122429830.062494436NC3231247521418890.073767716ZBED6-OE1228805816023770.299678155ZBED6-OE2235574719689420.164196325ZBED6-OE3233476819814690.151320816Wound healing assay length0h24hrateNC17336720.083219645NC28377860.0609319NC38187460.08801956ZBED6-OE1815662.66666670.186912065ZBED6-OE2814.66666676820.162847791ZBED6-OE3769.66666675550.278908618qPCR-2IGF2NCZBED6-SIZBED6-OE2-ΔΔT'10.8577723234.815713070.0194736442-ΔΔT'20.69817793432.823715890.0247173292-ΔΔT'31.66978999140.522932250.022902782WB data-2IGF2pro area valuegap area valuerationormalizeNC-1107398860311.2483639621.112140504ZBED6-SI-1164771735952.2388885111.994577441ZBED6-OE-162932955130.6588841310.586985648NC-288976844051.054155560.939124431ZBED6-SI-2136403725611.8798390321.674708011ZBED6-OE-261261916410.6684889950.595542414NC-391501859211.0649433780.948735066ZBED6-SI-3121970766081.592131371.418395466ZBED6-OE-350323972070.5176890550.461198003PIK3C3NC-184808627421.351694241.064815115ZBED6-OE-1100029569041.7578553351.38477392ZBED6-OE+pro-IGF2-170697611211.1566728290.911184405NC-287438726051.2042972250.948701156ZBED6-OE-2106669640301.6659222241.312352389ZBED6-OE+pro-IGF2-275316794780.9476333070.746510741NC-381209648501.2522590590.986483729ZBED6-OE-394516543071.740401791.371024657ZBED6-OE+pro-IGF2-375363656511.1479337710.904300096IGF2NC-1302831670110.1813233860.946801483ZBED6-OE-1180221789980.1006826890.525726558ZBED6-OE+pro-IGF2-11205811778570.6779660063.540079606NC-2418501835580.2279933321.19049412ZBED6-OE-2190511943070.0980458760.511958125ZBED6-OE+pro-IGF2-2976881950860.5007432622.614690107NC-3264751602430.1652178250.862704396ZBED6-OE-3165221632930.1011800870.528323781ZBED6-OE+pro-IGF2-3763981538910.4964422872.592232058ULK1NC-198986708341.3974362590.965946224ZBED6-OE-1148192650282.2788952451.575234822ZBED6-OE+pro-IGF2-176825757961.0135759140.700611438NC-2100882663601.5202230261.050819801ZBED6-OE-2138666692142.0034386111.38483165ZBED6-OE+pro-IGF2-280115722881.1082752320.766070201NC-394156661931.4224464820.983233975ZBED6-OE-3113090626161.8060879011.248417433ZBED6-OE+pro-IGF2-376211738971.0313138560.712872389P62NC-115578828230.188087850.916913288ZBED6-OE-111738855980.1371293720.668494767ZBED6-OE+pro-IGF2-128897899380.3212991171.566307606NC-215662832000.1882451920.917680319ZBED6-OE-210928810400.1348469890.657368332ZBED6-OE+pro-IGF2-223773856950.2774140851.352371568NC-317976751940.2390616271.165406392ZBED6-OE-39860705630.1397332880.68118865ZBED6-OE+pro-IGF2-323307781780.298127351.453347088LC3NC-1907961066460.8513774541.112458055ZBED6-OE-11356381045511.2973381411.695175574ZBED6-OE+pro-IGF2-1563621151650.4894021620.639480614NC-2661441133960.5833009980.762174157ZBED6-OE-2125566948211.324242521.730330361ZBED6-OE+pro-IGF2-2525011141890.4597728330.600765252NC-3918911066940.8612574281.125367788ZBED6-OE-3122970970371.2672485751.655858841ZBED6-OE+pro-IGF2-339676938250.4228723690.552549014

Following ZBED6 knockdown (Table S1) and overexpression in 769-P cells, Western blot analysis showed coordinated changes in both ZBED6 and PIK3C3 protein levels, confirming a strong positive correlation between them (Figure 3D and Table 1). Subcellular localization analysis using the Human Protein Atlas (HPA) database indicated predominant nucleoplasmic localization of ZBED6 (Figure 3E). HPA data further revealed significant upregulation of PIK3C3 protein in ccRCC tissues compared with normal kidney tissues (Figure 3F).

Moreover, Western blot analysis of autophagy markers in ZBED6-overexpressing (OE) versus negative control groups showed significant increases in ULK1 and LC3 expression, accompanied by decreased p62 levels (Figure 3G and Table 1), indicating the activation of autophagic flux following ZBED6 overexpression.

Expression characteristics of class III phosphatidylinositol 3-kinase in the single-cell landscape and spatial transcriptome analysis

Using marker genes specific to distinct cell types, cells were categorized into ten major clusters: B cells, endothelial cells, epithelial cells, macrophages, monocytes, NK cells, RCC cells, smooth muscle cells, T cells, and unclassified cells. The t-SNE dimensionality reduction plot depicts each cell type in distinct color clusters (Figure 4A). Figures 4B and 4C illustrate the spatial distribution of ZBED6 and PIK3C3 expression within the single-cell landscape. To assess PIK3C3 signaling pathway activity across various cell types, the AddModuleScore function in the Seurat package was used to calculate the expression scores of the PIK3C3 high-correlation gene set in all cells, integrating these data into a unified single-cell transcriptomic profile (Figure 4D). Among the ten identified cell types, B cells and T cells exhibited higher PIK3C3 signaling pathway activity, whereas epithelial cells and smooth muscle cells showed lower activity, as quantified by the violin plot (Figure 4E). Further analysis restricted to tumor cells revealed that PIK3C3 high-correlation gene set, was generally elevated in RCC tumor cells (Figure 4F).Figure 4. Joint verification using single-cell and spatial transcriptomics(A) The t-SNE plot shows the 10 cell types identified by the marker genes.(B) Expression distribution of ZBED6 in t-SNE plot. (C) Expression distribution of PIK3C3 in t-SNE plot.(D) Distribution of PIK3C3-related pathway scores in t-SNE plot (whole cell).(E) Boxplot shows the distribution of PIK3C3-related pathway scores across 10 cell types.(F) Distribution of PIK3C3-related pathway scores (tumor cells) under the t-SNE plot.

The expression domain of PIK3C3 was spatially colocalized with that of ZBED6. Unlike conventional analyses that rely primarily on quantitative data, co-localization analysis emphasizes the visualization and interpretation of spatial expression patterns. Using public 10x Genomics spatial transcriptomic data of RCC tumor tissues, we identified spatial co-localization between the expression regions of PIK3C3 and ZBED6 (Figures 4G and 4H). As illustrated, the expression domains of PIK3C3 and ZBED6 exhibited a high degree of overlap, further supporting their potential regulatory relationship.

Identification of class III phosphatidylinositol 3-kinase-related genes by weighted gene co-expression network analysis combined with single-cell RNA sequencing

The ssGSEA algorithm was used to calculate enrichment scores for each sample-gene set pair, generating PIK3C3 signaling pathway activity scores for each single-cell sample as phenotypic data for weighted gene co-expression network analysis (WGCNA). Figure 5A presents the dendrogram of DEGs identified from scRNA-seq data used for co-expression network construction. The upper section of the treemap shows hierarchical clustering results, whereas the lower heatmap depicts PIK3C3 pathway activity scores for each sample, calculated using the ssGSEA algorithm. WGCNA was performed with an optimal soft-thresholding power of 7 (R2 = 0.874), a minimum module size of 60 genes, and a module dissimilarity threshold of 0.25, resulting in nine distinct gene modules (Figure 5B). Among these, the blue module showed the strongest association with the PIK3C3 pathway activity score (correlation = 0.66, p = 4 × 10^-^66; Figure 5C). Furthermore, the scatterplot of gene significance versus module membership for the blue module demonstrated a strong correlation (correlation = 0.76, p = 3.3 × 10^−24^; Figure 5D), indicating that genes within this module were highly correlated with PIK3C3 pathway activity. Subsequently, the genes within the blue module were intersected with differentially expressed renal cancer genes, yielding 23 overlapping genes (Figure 5E). Log rank survival analysis of these 23 genes identified 18 with significant prognostic value (p < 0.05), as shown in Figure 5F. Finally, copy number variation analysis of 14 interacting genes revealed that EIF1B exhibited the highest amplification frequency, with copy number gains exceeding 13% (Figure 5G).Figure 5. Identify the set of prognostic-related genes closely associated with PIK3C3(A) Hierarchical clustering treemap. The bottom heatmap represents the PIK3C3-related pathway score for each sample.(B) Cluster dendrogram of WGCNA analysis.(C) Module-trait heatmap shows that the blue module is closely related to PIK3C3-related pathway score.(D) Scatterplot of gene significance (GS) and module membership (MM) in the blue module. (E) Venn diagram shows the intersection genes of the renal cancer differential gene and the blue module gene set.(F) Univariate Cox regression analysis of intersecting genes and the correlation between genes. The circle size represents the p-value of the univariate Cox regression analysis. (G) Copy number variation (CNV) frequencies of intersecting genes.

Development of a kidney cancer prognostic model using 101 machine learning algorithms

To construct an accurate prognostic model for the PIK3C3 signaling pathway, machine learning algorithms were applied to analyze 18 key genes identified through the log rank survival test. Ten algorithms were employed within a 10-fold cross-validation framework in the training cohort to generate 101 predictive models. The C-index was calculated for each model across all cohorts, and the mean C-index values were summarized on the right side of the plot (Figure 6A). Among these models, the Lasso-SuperPC algorithm achieved the highest average C-index of 0.651 and was therefore selected for final model construction. The resulting prognostic model comprised six genes: ASAH1, ATP1A1, DSTN, EIF1B, PGK1, and SCP2. Risk scores for each patient were calculated by applying the regression coefficients derived from the Cox proportional hazards model to the expression levels of the six genes. Based on the median risk score, samples from all cohorts were subsequently categorized into high-risk and low-risk groups.Figure 6. Construction and validation of machine learning prognosis models(A) A total of 101 prediction models was constructed through a 10-fold cross-validation framework, and the C-index of each model on all validation datasets was further calculated. Kaplan-Meier curves for TCGA cohort (B), ArrayExpress cohort (C), TCGA DSS cohort (D), and RECA-cohort (E). ROC curves for ArrayExpress cohort (F), TCGA DSS cohort (G), and RECA cohort (H).

Kaplan-Meier survival curves were generated for four cohorts to evaluate the prognostic performance of the constructed model. Across the TCGA test set, ArrayExpress, and RECA cohorts, samples classified as low-risk by the model exhibited significantly longer overall survival, whereas the high-risk group showed markedly poorer prognosis with less consistent statistical significance (Figures 6B–6E). The TCGA cohort demonstrated the most pronounced prognostic separation, likely because it served as the model’s training dataset. To further validate the predictive accuracy of the model, ROC curves were plotted for 1-, 2-, and 3-year survival in the ArrayExpress, TCGA test, and RECA cohorts. The area under the curve (AUC) values for the TCGA test cohort were 0.69, 0.65, and 0.65 at 1, 2, and 3 years, respectively. In the ArrayExpress cohort, the AUC values were 0.87, 0.88, and 0.88, while the RECA cohort yielded corresponding AUCs of 0.53, 0.50, and 0.52 (Figures 6F–6H). These findings preliminarily suggest that the constructed prognostic model demonstrates moderate yet consistent predictive capability across independent datasets.

Performance evaluation of prognostic models

A heatmap integrating clinical data visualized the expression profiles of the six model genes across risk groups in kidney cancer, revealing lower expression levels in high-risk patients, which were associated with poorer outcomes (Figure 7A). Boxplots showed significantly higher risk scores in advanced T stages (T3–4 vs. T1–2, p < 0.05; Figures 7B and 7C). Fan plot analysis (Figure 7D) and T-stage distribution (Figure 7E) further confirmed more advanced disease status and a higher proportion of T2–4 stages among high-risk patients. The prognostic model demonstrated modest diagnostic performance for metastasis prediction (AUC = 0.621; Figure 7F). Consistent survival differences between high- and low-risk groups were observed in Kaplan-Meier analysis (Figure 7G), collectively indicating that the activation of the PIK3C3 signaling pathway may serve as a poor prognostic indicator in kidney cancer.Figure 7. Validation of the clinical relevance of core genes(A) Heatmap shows the expression of six model genes between high and low risk groups and the clinical data of patients with kidney cancer. Boxplots (B and C) show the correlation of different T-stages with risk score.(D) Comparison of the distribution of clinical data between the two groups of high and low risk.(E) Comparison of the distribution of T stages between the two groups at high and low risk.(F) Diagnostic ROC curves of prognostic models in predicting metastasis of kidney cancer.(G) Kaplan-Meier curve in survival differences between the high- and low-risk groups. ∗ indicates p < 0.05, ∗∗ indicates p < 0.01, ∗∗∗ indicates p < 0.001, and ∗∗∗∗ indicates p < 0.0001.

Univariate and multivariate Cox regression analyses were performed on kidney cancer samples to determine whether the prognostic model served as an independent prognostic factor (Figures 8A and 8B). The univariate analysis identified T stage, M stage, overall stage, and risk score as significant prognostic indicators (p < 0.001). Furthermore, multivariate analysis confirmed that the risk score was an independent prognostic factor (p < 0.001), underscoring its robust predictive value in patients with kidney cancer. To improve clinical applicability, a nomogram was developed to estimate 1-, 3-, and 5-year survival probabilities in patients with kidney cancer (Figure 8C). Calibration curves demonstrated strong concordance between the nomogram-predicted and observed outcomes, while decision curve analysis showed a higher net clinical benefit compared with other individual variables (Figure 8C). Additionally, the C-index verified the nomogram’s stable and superior predictive accuracy, outperforming other clinical parameters in forecasting 1–10-year survival (Figures 8D and 8E).Figure 8. The transformation of clinical prognosis prediction modelsUnivariate (A) and multivariate (B) Cox regression analyses for prognostic models.(C) Nomogram, calibration curve, and DCA curve.(D) C-index of nomogram and other clinical data for 1–10 years survival prediction.

Functional enrichment analysis of prognostic gene sets

To further elucidate the molecular mechanisms linking the prognostic model to renal cancer outcomes, functional enrichment analysis was performed. Gene set variation analysis using the HALLMARK gene sets in GSEA: https://www.gsea-msigdb.org/gsea/msigdbdatabse revealed that the low-risk group was enriched in pathways associated with the negative regulation of oxidative phosphorylation, heme metabolism, pancreatic β-cell function, androgen response, adipogenesis, bile acid metabolism, fatty acid metabolism, and protein secretion (Figure 9A). In contrast, the high-risk group exhibited enrichment in pathways promoting allograft rejection, epithelial-mesenchymal transition, inflammatory response, interferon-γ response, and TNF-α signaling via NF-κB (Figure 9B). Correlation analysis showed that the risk score was positively associated with the p53, Wnt/β-catenin, IL6/JAK/STAT3, allograft rejection, and DNA repair pathways, while being negatively correlated with pathways related to protein secretion, androgen response, heme metabolism, fatty acid metabolism, bile acid metabolism, and adipogenesis (Figures 9C and 9D). These findings suggest that the prognostic gene set is functionally linked to cancer-associated metabolic and cell death processes. To validate these associations, Kaplan-Meier survival analyses were performed for the two pathways exhibiting the strongest correlations with the risk score: protein secretion (the strongest negative correlation) and the p53 pathway (the strongest positive correlation). Patients with higher p53 pathway activity, which was positively correlated with risk score, exhibited poorer survival, although the association did not reach statistical significance (p = 0.09). Conversely, higher activity of the protein secretion pathway, which was negatively correlated with risk score, was associated with improved survival outcomes (Figures 9E and 9F).Figure 9. Enrichment analysis of PIK3C3-related prognostic genes(A) Analysis of the characteristic gene set in the low-risk GSVA group.(B) Analysis of the characteristic gene set in the GSVA high-risk group.(C) Correlation of risk score with 50 signaling pathways.(D) Risk score and the relationship between the signaling pathways are plotted, interacting. K-M curve analysis of protein secretion (E) and P53 pathway (F). ∗ indicates p < 0.05, ∗∗ indicates p < 0.01, ∗∗∗ indicates p < 0.001, and ∗∗∗∗ indicates p < 0.0001.

Analysis of immune microenvironment differences between risk subgroups

The composition of immune and stromal cells within tumor tissues plays a critical role in determining patient prognosis. These non-tumor cellular components constitute the major elements of the TME and are essential for accurate tumor characterization and prognostic assessment. Using the ESTIMATE algorithm, we calculated immune, stromal, and total ESTIMATE scores across risk subgroups, revealing significantly higher scores in the high-risk group (Figure 10A). Immune cell infiltration analysis indicated that the high-risk subgroup exhibited elevated proportions of activated B cells, activated CD4^+^ T cells, activated CD8^+^ T cells, and myeloid-derived suppressor cells (Figure 10B). The correlations between the six model genes and 22 immune cell types are presented in Figure 10C. Furthermore, correlation analysis between the risk score and immune cell subsets demonstrated positive associations with regulatory T cells (Tregs), follicular helper T cells, and CD8^+^ T cells, and negative associations with M2 macrophages, resting mast cells, and eosinophils (Figure 10D). Although CD8^+^ and CD4^+^ T cells were elevated in the high-risk group, this increase may not represent effective antitumor immunity, as these T cells are likely functionally exhausted. In addition, a trend toward higher Treg cell levels was observed in high-risk patients, although this was not statistically significant, further suggesting the presence of an immunosuppressive tumor microenvironment.Figure 10. Immune infiltration analysis of PIK3C3-related prognostic genes(A) Immune scores, matrix scores, and estimated scores of the high-risk group and the low-risk group.(B) Immune cell infiltration in the high- and low-risk groups was scored.(C) Correlation between 6 model genes and 22 immune cells.(D) Correlation analysis of risk score with 22 immune cells. ∗ indicates p < 0.05, ∗∗ indicates p < 0.01, ∗∗∗ indicates p < 0.001, and ∗∗∗∗ indicates p < 0.0001.

Immunotherapy and targeted drug susceptibility analysis in different risk subgroups

The marked tumor heterogeneity of renal cancer renders immune cell infiltration assessment alone insufficient for accurately characterizing antitumor immune activation and dysfunction. Tumor immunotherapy aims to restore antitumor immunity by reactivating and sustaining the tumor-immune cycle, thereby restraining tumor growth. Comprehensive evaluation of the stages of antitumor immune response provides a deeper understanding of renal cancer immune dynamics and enhances the development of effective therapeutic strategies.

To improve the prognostic model’s predictive capacity for immunotherapy response, we incorporated the atezolizumab-treated IMvigor210 cohort to predict immune checkpoint inhibitor responsiveness. Analysis revealed a higher proportion of progressive disease and stable disease outcomes in high-risk patients (Figures 11A–11C), indicating poorer responses to immunotherapy. This observation was supported by violin plots stratified according to risk score. SubMap analysis further demonstrated differential associations with checkpoint therapies: high-risk patients exhibited stronger similarity to anti-PD-1 response signatures, whereas low-risk patients were more closely aligned with anti-CTLA-4 efficacy profiles (Figure 11D). Immune target gene profiling showed increased CD274 expression in high-risk patients, whereas TIGIT, CTLA4, HAVCR2, and LAG3 were more highly expressed in the low-risk group (Figure 11E).Figure 11. Immune response analysis of prognosis-related genes related to PIK3C3(A) Distribution of CR/PR and PD/SD in patients treated with atezolizumab between high- and low-risk groups.(B) Correlation between risk score and CR/PR and PD/SD.(C) Correlation between risk score and CR, PR, PD, and SD. CR complete response: all target lesions disappeared, no new lesions appeared, and tumor markers were normal for at least 4 weeks. PR: partial response was achieved by ≥ 30% reduction in the sum of the maximum diameters of the target lesions for at least 4 weeks. SD: stable disease, the sum of the maximum diameters of the target lesions decreased without reaching PR, or increased without reaching PD. PD: progressive disease, the sum of the maximum diameters of the target lesions increased by at least 20% or new lesions appeared.(D) Submap algorithm to compare correlations between patients in different risk subgroups and those receiving immune checkpoint therapy.(E) Differences in the expression of 6 immune target genes in the high and low risk groups. ∗ indicates p < 0.05, ∗∗ indicates p < 0.01, ∗∗∗ indicates p < 0.001, and ∗∗∗∗ indicates p < 0.0001.

For therapeutic prediction, we assessed susceptibility to eight first-line clinical drugs. High-risk patients displayed higher sensitivity to axitinib, sunitinib, metformin, and temsirolimus, whereas low-risk patients were more responsive to sorafenib, imatinib, pazopanib, and tipifarnib (Figure 12A). XSum analysis of the top 300 differentially expressed genes identified five phenotype-reversing compounds based on Connectivity Map scoring: PHA.00816795, MS.275, STOCK1N.35696, STOCK1N.35874, and NU.1025 (Figure 12B). Complementary analyses of CTRP: https://portals.broadinstitute.org/ctrpderived (Figure 12C) and PRISM: https://depmap.org/portal/prismderived (Figure 12D) compounds revealed differential drug sensitivities, where lower boxplot values indicated increased sensitivity. Collectively, these results establish an integrated framework linking molecular risk stratification with precision therapeutic selection in renal cancer.Figure 12. Predictive analysis of potential targeted drugs for PIK3C3-related prognostic genes(A) Drug susceptibility prediction of 8 targeted drugs.(B) The top five drugs calculated by the XSum analysis all had the potential to reverse the phenotype of high-risk populations.(C and D) Results from Spearman’s correlation analysis and differential drug response analysis for 4 CTRP-derived compounds and 8 PRISM-derived compounds.

Overexpression of ZBED6 promotes invasion, metastasis, and proliferation of clear cell renal cell carcinoma cells

A series of in vitro functional assays was conducted to systematically evaluate the biological role of ZBED6 in ccRCC. The CCK-8 proliferation assay showed that ZBED6 overexpression significantly enhanced the proliferative capacity of 769-P cells compared with the negative control, whereas ZBED6 silencing markedly suppressed cell proliferation (Figure 13A and Table 1). Wound healing and Transwell migration assays revealed that ZBED6 overexpression substantially increased the migratory potential of 769-P cells (Figure 13B, 13C and Table 1). Similarly, the Transwell invasion assay demonstrated that ZBED6 overexpression significantly promoted the invasive capacity of 769-P cells through Matrigel (Figure 13D and Table 1). Collectively, these findings suggest that ZBED6 functions as a potential oncogenic driver and therapeutic target in ccRCC, with its expression level bearing prognostic significance.Figure 13. Experimental verification of ZBED6-IGF2-PIK3C3 in regulating the phenotype of renal cancer(A) The CCK-8 assay shows the proliferation viability of renal cancer cells after ZBED6 knockout or overexpression compared with the NC group.(B and C) The wound healing assay and Transwell migration assay showed the migration ability of renal cancer cells treated with NC or ZBED6 overexpression. The image is magnified 40X.(D) The Transwell invasion assay showed the invasion ability of renal cancer cells treated with NC or ZBED6 overexpression. The image is magnified 40X.(E) The heatmap of qPCR experiments shows the changes in the content of highly correlated genes in renal cancer cells after ZBED6 overexpression.(F) The heatmap shows the correlation between ZBED6, PIK3C3, and highly correlated genes.(G and H) qPCR and WB experiments verified the correlation between ZBED6 and IGF2. Overexpression of ZBED6 inhibited the expression of IGF2, while silencing ZBED6 promoted the expression of IGF2.(I) The WB results demonstrate the expression changes of IGF2 and PIK3C3 in the NC group, the ZBED6-OE group, and the ZBED6-OE + pro-IGF2 group.(J) The WB experiment shows the differences in the protein expression of autophagy markers ULK1, p62, and LC3. ∗ indicates p < 0.05, ∗∗ indicates p < 0.01, ∗∗∗ indicates p < 0.001, and ∗∗∗∗ indicates p < 0.0001.

qPCR analysis revealed significant expression changes in six genes highly correlated with PIK3C3 following ZBED6 overexpression in renal cancer cells. The mRNA levels of PIK3C3, ASAH1, ATP1A1, DSTN, and EIF1B were markedly upregulated, while SCP2 was downregulated, and PGK1 expression remained largely unchanged (Figure 13E and Table 1). Correlation analysis using TCGA datasets confirmed a significant positive correlation between these six genes and the expression levels of both ZBED6 and PIK3C3 (Figure 13F).

To elucidate the intermediate mechanism through which ZBED6 regulates PIK3C3, we first verified the expression relationship between ZBED6 and IGF2 using qPCR and Western blot (WB) analyses. The results showed that IGF2 expression was markedly decreased in the ZBED6-OE group and increased in the ZBED6-silencing group, indicating a negative regulatory relationship between ZBED6 and IGF2 (Figures 13G and 13H). To further investigate this interaction, rescue experiments were performed. In the ZBED6-OE group, exogenous IGF2 protein (pro-IGF2) was supplemented. WB results demonstrated that PIK3C3 protein levels were significantly increased, whereas IGF2 protein levels were decreased in the ZBED6-OE group. However, following exogenous IGF2 supplementation, IGF2 protein expression increased, and PIK3C3 protein levels decreased (Figure 13I).

Evaluation of autophagy-related markers revealed that ZBED6 overexpression significantly upregulated ULK1 and LC3 while downregulating p62, indicating the activation of autophagic flux. Conversely, after exogenous IGF2 supplementation, ULK1 and LC3 expression levels decreased, and p62 expression increased, reversing the autophagy-activating effect. Collectively, these findings demonstrate that ZBED6 promotes PIK3C3 upregulation by suppressing IGF2 expression, thereby activating the autophagy pathway. IGF2 thus serves as a key mediator in the ZBED6–PIK3C3–autophagy regulatory axis (Figure 13J).

Identification of compounds with potential binding affinity to ZBED6

Based on Enrichr’s online analysis results, we identified the top ten compounds most likely to bind ZBED6 from the Proteomics Drug Atlas 2023 dataset. Among these, BML-277 and Tubastatin A exhibited the highest predicted binding affinities and were most likely to interact with and inhibit ZBED6 protein activity (Figure 14A). To further explore their potential interactions, molecular docking simulations between ZBED6 and these two compounds were performed using AutoDock software. The results revealed the two most probable binding conformations (Figure 14B, 14C and Table S2).Figure 14. Prediction of drugs with high binding potential for ZBED6 through molecular docking(A) Drugs with potential binding to ZBED6 and their p-value ranking.(B and C) The molecular docking diagram shows the best predicted binding conformation of BML- 277/Tubastatin A and ZBED6 protein.

Discussion

By systematically integrating multi-omics data with experimental validation, this study identifies the transcription factor ZBED6 as a key upstream regulator of the autophagy core kinase PIK3C3 in ccRCC for the first time. Furthermore, a six-gene prognostic signature with strong predictive performance was constructed based on this regulatory axis. These findings not only deepen our understanding of the mechanisms underlying autophagy initiation in ccRCC but also provide valuable insights for prognostic assessment and the development of targeted therapeutic strategies for this malignancy.

The most significant finding of this study is the identification of a functional regulatory axis between ZBED6 and PIK3C3 in ccRCC. As a class III phosphatidylinositol 3-kinase, PIK3C3 serves as the core catalytic component required for autophagosome formation. Although its downstream effects have been extensively characterized, the upstream transcriptional regulation of PIK3C3 in ccRCC has remained poorly understood. For the first time, our study identifies ZBED6 as a regulator of this pivotal autophagic node. This discovery is particularly important because ZBED6 is a transcription factor derived from a domesticated DNA transposon and is known to exert essential biological functions in placental mammals.13 Notably, the mechanism by which ZBED6 regulates PIK3C3 may operate independently of its classical regulatory pathway.11 Under physiological conditions such as muscle development, ZBED6 suppresses IGF2 transcription primarily by binding to a specific motif within the third intron.11^,^13 In accordance with this, our results demonstrated that the ZBED6-mediated regulation of PIK3C3 was accompanied by a consistent inverse change in IGF2 expression in ccRCC models. This finding suggests that ZBED6 may enhance PIK3C3 expression by inhibiting the suppressive activity of IGF2 on PIK3C3 transcription in ccRCC. Such tissue-specific regulatory divergence underscores the complexity and plasticity of intracellular signaling networks in cancer cells. Furthermore, spatial transcriptomics analysis (p < 0.001) revealed the co-localization of ZBED6 and PIK3C3 expression, supporting the potential of a functional interaction between these proteins within the tumor microenvironment—a promising direction for future research.14 However, the observed spatial co-expression pattern indicated a significant negative correlation between their expression levels across two ccRCC samples. Taken together, these results suggest that ZBED6 may influence PIK3C3 regulation through the intermediate action of IGF2. The precise mechanism of this IGF2*–*PIK3C3 regulatory axis warrants further elucidation using experiments such as chromatin immunoprecipitation sequencing (ChIP-seq).

Building on the ZBED6/PIK3C3 regulatory network, we developed a prognostic signature comprising six genes—ASAH1, ATP1A1, DSTN, EIF1B, PGK1, and SCP2—using WGCNA and multiple machine learning algorithms.15^,^16 The strong predictive capability of this signature derives from its solid biological foundation, as these genes are closely associated with the intersection between autophagy and tumor progression. For instance, PGK1 (phosphoglycerate kinase 1), a key glycolytic enzyme, is known to be modulated by autophagy and, in turn, influence autophagic flux17; DSTN participates in actin cytoskeleton remodeling, a process that is critically linked to both autophagosome formation and cell migration.18 Collectively, these genes delineate an aggressive tumor phenotype characterized by metabolic reprogramming and cytoskeletal plasticity. From a translational perspective, this prognostic marker effectively stratified high-risk patients with shorter overall survival (HR = 2.98, p < 0.0001) and predicted patient responses to immune checkpoint inhibitors, particularly anti-PD-1 therapy.19 Patients in the high-risk group exhibited marked resistance to immunotherapy, reflected by a 38% reduction in CR/PR rate (p = 0.007), yet appeared more sensitive to tyrosine kinase inhibitors such as axitinib20^,^21 and sunitinib.22 These findings highlight the clinical utility of this model as a tool for precise risk stratification and individualized treatment selection in ccRCC. Nevertheless, this model demonstrated suboptimal performance (AUC ≈0.5) in the RECA validation cohort. This limitation may be attributable to multiple factors, including cohort heterogeneity, inconsistencies in sample collection and processing, and batch effects arising from distinct sequencing platforms. These observations underscore the need to validate genomic signatures in larger, prospective, multicenter cohorts and to apply advanced batch-correction methods such as Harmony23 or ComBat24 to enhance model robustness and generalizability prior to clinical implementation.

Our study uncovered an intriguing phenomenon that merits further investigation: higher ZBED6 and PIK3C3 mRNA expression levels appeared to correlate with better patient prognosis in portions of the public cohorts, which seems to contradict the oncogenic roles confirmed by our in vitro functional experiments. Several potential explanations may account for this discrepancy. First, there is a confounding effect of cellular composition. RNA extracted from bulk tumor tissue represents an aggregate of multiple cellular populations, including malignant cells, immune cells, and stromal cells. The elevated expression of ZBED6 and PIK3C3 may therefore reflect active physiological autophagy occurring in tumor-infiltrating immune cells or residual normal renal tubular epithelial cells.25^,^26 Such autophagic activity could contribute to antitumor immunity or tissue repair, thereby appearing statistically as a protective factor.27 Second, a tumor stage-dependent effect may exist. Autophagy plays a well-documented “double-edged sword” role in cancer biology. In the early stages of tumorigenesis, autophagy may suppress malignant transformation by clearing damaged organelles and maintaining genomic stability. However, in advanced disease, it can facilitate tumor survival and metastasis by enabling cells to withstand metabolic stress, hypoxia, and nutrient deprivation.28 Our functional assays were conducted primarily in established ccRCC cell lines, which more closely reflect the behavior of late-stage tumors. Third, post-transcriptional and post-translational regulation may contribute to the apparent inconsistency. We observed a decoupling phenomenon between mRNA and protein expression levels of ZBED6 and PIK3C3 (i.e., low RNA abundance accompanied by high protein expression) in ccRCC cells.29^,^30 Furthermore, ZBED6 can partially regulate Igf2 expression by modulating miR483 expression. ZBED6 reduces the interaction between Igf2 and miR483 by inhibiting the Igf2 binding site; meanwhile, miR483 (particularly miR4835p) promotes Igf2 transcription and downstream growth-promoting signals such as PI3KAkt, which play crucial roles in muscle development and tumor progression. In transgenic mice, when the ZBED6-Igf2 axis is disrupted, miR483 expression is significantly upregulated, and CRISPR/Cas9-mediated knockout of miR483 results in downregulated Igf2 expression and reduced cell proliferation rates.31 This observation strongly implies the presence of post-transcriptional regulatory mechanisms, such as microRNA-mediated translational repression or alterations in protein degradation pathways.32 Therefore, mRNA expression levels alone may not reliably represent the functional activity of this signaling pathway. Such context-dependent dual functionality is a well-recognized challenge in cancer biology and underscores the critical need for the precise design of targeted therapeutic strategies.33^,^34 To explore potential therapeutic agents acting on the ZBED6–PIK3C3 regulatory axis, we performed molecular docking analyses to evaluate the binding affinity between compounds identified through Enrichr’s online prediction platform and PIK3C3. The results revealed a high probability of drug-PIK3C3 interaction. Nonetheless, this remains a preliminary, in silico investigation without validation in cellular or in vivo models. These findings are intended to provide a foundation for future experimental studies and drug development efforts.

Taken together, this study proposes and validates a novel ZBED6-driven PIK3C3 signaling pathway that regulates pro-survival autophagy and promotes tumor progression in ccRCC. The resulting six-gene prognostic signature provides a valuable tool for patient stratification and individualized treatment decision-making. Moving forward, several critical issues warrant urgent investigation. (1) Verification of the direct regulatory mechanism: It is imperative to determine whether ZBED6 directly binds to the promoter or enhancer regions of the PIK3C3 gene using ChIP-qPCR and ChIP-seq assays, to confirm its role as a direct transcriptional regulator of PIK3C3. (2) Exploration of the RNA-protein uncoupling mechanism: The molecular basis underlying the observed discordance between ZBED6/PIK3C3 mRNA and protein levels in ccRCC cells requires elucidation. This may involve regulatory mechanisms such as microRNA interference, RNA-binding proteins, or alterations in the ubiquitin-proteasome pathway. (3) In vivo functional validation: Conditional knockout or transgenic mouse models with ZBED6 overexpression should be employed to establish the causal role of the ZBED6–PIK3C3 axis in ccRCC initiation, progression, and therapy resistance. (4) The development of therapeutic targeting strategies: Future studies should assess whether the pharmacological inhibition of this axis—such as through the design of small-molecule inhibitors that specifically disrupt ZBED6–DNA binding—can effectively suppress ccRCC progression, either alone or in combination with immunotherapy or targeted therapy. Collectively, the continued exploration of this regulatory axis may not only elucidate a previously unrecognized oncogenic pathway but also provide a promising avenue for overcoming therapeutic resistance in ccRCC.

By integrating multi-omics data, this study is the first to demonstrate that the transcription factor ZBED6 activates PIK3C3 expression through the regulation of IGF2, thereby initiating pro-survival autophagic flux and promoting the malignant progression of ccRCC. Spatial transcriptomic analysis confirmed the significant co-localization of these molecules within the tumor microenvironment. Based on genes associated with this pathway, we constructed and validated a six-gene prognostic signature capable of effectively distinguishing high-risk patients characterized by shorter overall survival, reduced responsiveness to anti-PD-1 immunotherapy, and increased sensitivity to targeted agents. Functional assays further demonstrated that ZBED6 interference markedly suppressed the invasive and migratory capabilities of ccRCC cells. The principal limitation of this study lies in the absence of definitive experimental verification of the indirect transcriptional cascade through which ZBED6 regulates PIK3C3 via IGF2. Key regulatory interactions have yet to be confirmed by ChIP assays, which represent a critical direction for future research aimed at elucidating the precise transcriptional mechanisms underlying this pathway.

Limitations of the study

Although this study revealed the role of the ZBED6-IGF2-PIK3C3 axis in ccRCC through multi-omics integration and experimental validation, several limitations remain. First, the key regulatory mechanism—whether ZBED6 directly binds to the promoter/enhancer regions of PIK3C3 to regulate its transcription—has not been conclusively confirmed through experiments such as chromatin immunoprecipitation (ChIP). Second, the predictive performance of the six-gene prognostic model constructed in this study (AUC ≈0.5) in external validation cohorts (e.g., RECA) was suboptimal, suggesting that the model may be influenced by cohort heterogeneity, batch effects, or technical variations, and its generalizability requires further validation in larger, prospective multicenter cohorts. Additionally, we observed that the mRNA expression levels of ZBED6/PIK3C3 were associated with favorable patient outcomes, which appears to contradict the pro-oncogenic effects revealed by in vitro functional experiments. This discrepancy may stem from the heterogeneous cell populations in the tumor microenvironment, the “double-edged sword” effect of autophagy at different stages of cancer, or post-transcriptional regulatory mechanisms, but the underlying causes have not been fully elucidated in this study. Finally, the potential therapeutic compounds screened for this axis were only analyzed through computational docking, lacking pharmacodynamic validation at cellular or in vivo levels. These limitations highlight key areas for future research breakthroughs.

Resource availability

Lead contact

Guangzhen Wu: Department of urology, First affiliated hospital of Dalian Medical University. No.222 Zhongshan Road, Xigang District, Dalian City, Liaoning Province, China. Email: wuguang0613@hotmail.com.

Materials availability

The study did not generate new unique reagents, or there are restrictions on availability. All supplementary material can be found via DOI: https://doi.org/10.57760/sciencedb.27556.

Data and code availability

Data availability

All rawdata are publicly available in scienceDB and can be accessed directly via DOI: https://doi.org/10.57760/sciencedb.27556.

Code availability

All raw-code are publicly available in scienceDB and can be accessed directly via DOI: https://doi.org/10.57760/sciencedb.27556.

All other items availability

Any additional information required to reanalyze the data reported in this article is available from the lead contact upon request.

STAR★Methods

Key resources table

REAGENT or RESOURCESOURCEIDENTIFIERAntibodiesAnti-ZBED6 AntibodyImmunowayCat Num: YT6573Anti-Vps34 AntibodyBioworldCat Num: AP1010Anti-LC3 AntibodyBioworldCat Num: AP1013Anti-p62 AntibodyBioworldCat Num: BS65004Anti-ULK1 AntibodyBioworldCat Num: BS91404Anti-IGF2 AntibodyProteintechCat Num: 84643-1-RR; RRID:AB_3672102Anti-GAPDH AntibodyImmunowayCat Num: PT0582RAnti-β-actin AntibodyProteintechCat Num: 81115-1-RR; RRID:AB_2923704HRP-conjugated goat anti-rabbit IgG secondary antibodyABclonalCat Num: AS014PrimerZBED6-PrimerGENERAL BIOLforward 5'-GAAGGGTTTGCGAATTAAGGGG-3', reverse 5'-GGGTCATTGGAAGCTAACAAAGC-3'Vps34-PrimerGENERAL BIOLforward 5'-CCTGGAAGACCCAATGTTGAAG-3', reverse 5'-CGGGACCATACACATCCCAT-3'ASAH1-PrimerGENERAL BIOLforward 5'-AGATGTCATGTGGATAGGGTTCC-3', reverse 5'-GGGGCCAATATCTTGGTCTTG-3'ATP1A1-PrimerGENERAL BIOLforward 5'-CTGTGGATTGGAGCGATTCTT-3', reverse 5'-TTACAACGGCTGATAGCACCA-3'DSTN-PrimerGENERAL BIOLforward 5'-ATTTTGTGGGAATGCTTCCTGA-3', reverse 5'-GCATCCTTGGAGCTTGCATAG-3'EIF1B-PrimerGENERAL BIOLforward 5'-CGGAACGGCAGAAAGACACT-3', reverse 5'-ACCTCTCCGTATTCAGGATGTT-3'PGK1-PrimerGENERAL BIOLforward 5'-TGGACGTTAAAGGAAGCGG-3', reverse 5'-GCTCATAAGGACTACCGACTTGG-3'SCP2-PrimerGENERAL BIOLforward 5'-CCACTGATAAGCATGTTGACCT-3', reverse 5'-TGGAACTGGGAATACGGGTTTA-3'Cell culturefetal bovine serum (FBS)Wuhan PricellaCat# 164210-50streptomycin-penicillinWuhan PricellaCat# PB180120McCoy's 5A MediumWuhan PricellaCat# PM150710RPMI-1640 MediumWuhan PricellaCat# PM150110Minimum Essential Medium (MEM)Wuhan PricellaCat# PM150410Dulbecco's Modified Eagle's Medium (DMEM)Wuhan PricellaCat# PM150210Main chemicals, peptides, and recombinant proteinsCCK-8APExBIOCat# K1018ZBED6 oeRNA plasmidApplied Biological Materials3×Flag tag (pLenti-GIII-CMV-C-term-3xFlag-CBH-GFP-2A-Puro, CL097, 2937bp)Lipo2K Transfection ReagentAPExBIOZBED6 siRNASuzhou Haixing Biosciences#1, forward 5'-GCUUCAUUGUUUCUGACAATT-3', reverse 5'-UUGUCAGAAACAAUGAAGCTT-3', #2, forward 5'-GGCUAAAGACUGUUUGAUATT-3', reverse 5'-UAUCAAACAGUCUUUAGCCTT-3', #3, forward 5'-CAACGACAUCUGCAAGCAATT-3', reverse 5'-UUGCUUGCAGAUGUCGUUGTT -3'.2×RealStar Universal SYBR qPCR Mix KitGenStarCat# A308-05TRIGene Plus Total RNA Extraction Reagent and Auxiliary KitGenStarCat# P108-05StarScript Pro All-in-one RT Mix with gDNA RemoverGenStarCat# A240-10Ultra Sensitive ECL Chemiluminescence KitSevenbioCat# SW134-01MatrigelKeygentecCat# KGL5102-1crystal violetBeyotimeCat# C0121recombinant human pro IGF-II proteinProteintechCat# HZ-1161Software and algorithmsRR Foundationhttps://www.r-project.org/GraphPad PrismGraphPad SoftwareVersion 9.0ImageJNIHVersion 1.54gDeposited dataTCGANCIhttps://www.cancer.gov/tcgaCCLEBroadinstitutehttps://sites.broadinstitute.org/ccleArrayExpressEMBL-EBIhttps://www.ebi.ac.uk/biostudies/arrayexpressICGCICGC-ARGOhttps://www.icgc-argo.org/GEONCBIhttps://www.ncbi.nlm.nih.gov/geo/Code and datascienceDBhttps://doi.org/10.57760/sciencedb.27556

Experimental model and study participant details

Public data acquisition and new data & original code

This study did not generate new unique reagents. Data and code availability: All public datasets analyzed in this study were obtained from their original publications, including CCLE,35 TCGA,36 ArrayExpress,37 ICGC,38 and GEO39^,^40 (accession numbers: [GSE156632](GSE156632) and [GSE210041](GSE210041)). All the codes used for the bioinformatics analyses, along with the initial data and newly generated data, have been packaged and uploaded to scienceDB. DOI: https://doi.org/10.57760/sciencedb.27556. All the original images of the cell validation experiments, including western blot, cell scratch, and transwell experiments, can be found in the "Original image data" and "Original image data-revision" sections of the supplemental information.

Cell lines

The human renal cancer cell lines 769-P, Caki-1, and ACHN, as well as the normal renal proximal tubular epithelial cell line HK-2, were obtained from Wuhan Pricella Biotechnology Co., Ltd. All cell lines were authenticated using short tandem repeat profiling and confirmed to be free of mycoplasma contamination. Cells with low passage numbers (passage numbers <3) in the logarithmic growth phase were used for all experiments. Cells were cultured under standard conditions (37°C, 5% CO_2_) in the following media, each supplemented with 10% fetal bovine serum and 1% penicillin–streptomycin (all reagents from Wuhan Pricella Biotechnology Co., Ltd.): Caki-1 in McCoy’s 5A Medium, 769-P in RPMI-1640 Medium, ACHN in Dulbecco’s Modified Eagle Medium, and HK-2 in Minimum Essential Medium.

Public sequencing data subjects

The use of de-identified patient data obtained from public databases such as TCGA, ICGC, and ArrayExpress complied with the respective data-use policies and ethical guidelines of those databases.

Method details

Bioinformatic analysis of public datasets

Gene expression matrices and clinicopathological data were obtained from the CCLE, TCGA, ArrayExpress, and ICGC databases. The ggplot241 R package was used to generate scatterplots showing correlations between two genes, while the pheatmap R package42 employing Spearman’s correlation analysis was used to visualize correlations among multiple genes. A gene co-expression network was constructed using the WGCNA R package43 to identify phenotype-associated gene modules. Utilizing 11 machine learning algorithms and 10-fold cross-validation, 101 prognostic models for kidney renal clear cell carcinoma were developed, with model performance evaluated using Harrell’s concordance index, leading to the identification of an optimal six-gene signature.44 Survival analysis was performed using the survival and survminer R packages,45^,^46 receiver operating characteristic (ROC) analysis using the pROC R package,47 and a nomogram was constructed using the rms R package.48 Gene set enrichment analysis was employed to explore expression differences in specific gene sets between high- and low-risk groups.49 The tumor immune microenvironment was characterized using the estimate R package50 and single-sample gene set enrichment analysis (ssGSEA).51 The pRRophetic algorithm, integrated with GDSC and TCGA gene expression profiles, was used to predict half-maximal inhibitory concentration values.52^,^53 The Seurat R package was used to process scRNA-seq data, including quality control, normalization, dimensionality reduction, clustering, and cell-type annotation using the SingleR R package.54 Spatial transcriptomic data from 10x Genomics were analyzed using R packages including Seurat,55 which involved quality control, normalization, clustering, differential expression analysis, identification of spatially variable genes, and cell-type annotation.

Cell culture and transfection

Cell culture was performed as described in the Experimental Model section. The ZBED6 lentiviral vector containing a 3×FLAG tag was transfected into 769-P cells using the Lipo2K Transfection Reagent, and cells were collected 48 h after transfection. ZBED6-targeting small interfering RNA was transfected into 769-P cells using the GP-transfect-Mate transfection reagent, and subsequent assays were performed 48 h later.

Quantitative PCR (qPCR)

Total RNA was extracted using the TRIGene Plus Total RNA Extraction Reagent, and complementary DNA was synthesized using the StarScript Pro All-in-One RT Mix with gDNA Remover. qPCR was performed using the 2× RealStar Universal SYBR qPCR Mix Kit on a quantitative PCR instrument. Each sample was analyzed in three technical replicates. Relative gene expression levels were calculated using the 2^–ΔΔCt^ method and normalized to GAPDH or β-actin as internal controls. Primer sequences are provided in Section 2.14.

Western blotting

Total protein was extracted from tissues or cells. Equal amounts of protein were separated by SDS–PAGE and transferred onto PVDF membranes. Membranes were blocked with 5% skim milk for 2 h at 37°C and then incubated overnight at 4°C with primary antibodies against ZBED6, PIK3C3, LC3, p62, ULK1, GAPDH, and β-actin. After washing, membranes were incubated with an HRP-conjugated goat anti-rabbit IgG secondary antibody for 1.5 h at room temperature. Protein bands were visualized using an Ultra Sensitive ECL Chemiluminescence Kit and quantified with ImageJ software.

Cell proliferation assay (CCK-8)

Cells were seeded in 96-well plates at a density of 2 × 10^3^ cells per well. At the indicated time points (0, 24, 48, and 72 h), 10 μL of CCK-8 reagent mixed with 100 μL of serum-free medium was added to each well. After incubation for 2 h at 37°C, absorbance was measured at 450 nm using a microplate reader.

Wound healing assay

Cells were seeded in 6-well plates and cultured until reaching full confluence. A sterile 200 μL pipette tip was used to create a linear wound. Detached cells were gently removed by washing with PBS, and the medium was replaced with low-serum medium. Wound images were captured at regular intervals under a light microscope, and the migration area was quantified using ImageJ software.

Transwell migration and invasion assays

For the migration assay, 2 × 10^4^ cells suspended in serum-free medium were seeded into the upper chamber of a Transwell insert. The lower chamber contained medium supplemented with 10% fetal bovine serum as a chemoattractant. After 48 h of incubation, cells that migrated to the lower surface were fixed, stained with crystal violet, and counted. For the invasion assay, the upper chamber membrane was pre-coated with Matrigel for 3 h before seeding 5 × 10^4^ cells, and subsequent steps were performed as in the migration assay. Cell counting was conducted manually in three random fields at 100× magnification for each replicate.

Molecular docking analysis

Ligands for molecular docking were obtained based on Enrichr’s online analysis results, with PIK3C3, the downstream molecule of the ZBED6–PIK3C3 regulatory axis, serving as the receptor. Rigid docking was performed using the Easy2MD platform56 to validate receptor–ligand binding.

Quantification and statistical analysis

All bioinformatic analyses were performed using R software. Data from cell-based experiments were analyzed using GraphPad Prism. For two-group comparisons of continuous variables with large data from the TCGA-KIRC cohort and the associated validation set cohort, Wilcoxon-test was used. Statistical analysis of data from experiments was performed using ANOVA. Comparisons among multiple groups were performed using one-way analysis of variance. Correlation analyses were performed using Pearson’s or Spearman’s correlation coefficients, as appropriate. Survival differences were assessed using the log rank test. The predictive performance of prognostic models was evaluated using Harrell’s concordance index (C-index), while the performance of diagnostic or predictive models was assessed using the area under the ROC curve. Statistical significance was defined as a false discovery rate (FDR) or p-value less than 0.05. Exact n values are provided in the respective figure legends. All quantitative data are presented as the mean ± standard deviation. In all statistical test results, ∗ indicates p < 0.05, ∗∗ indicates p < 0.01, ∗∗∗ indicates p < 0.001, and ∗∗∗∗ indicates p < 0.0001.

Acknowledgments

We thank GEO database and the TCGA database for open access to data, and we thank Director Qin of the Department of Pathogenic Microbiology, 10.13039/501100012465Dalian Medical University, for providing key experimental instruments. This research was supported by the Scientific Research Fund of 10.13039/501100007620Liaoning Provincial Education Department (No. JYTMS20230577), the Joint Program of Science and Technology Program of Liaoning Province, the Zhejiang Tsinghua Yangtze River Delta Research Institute Key Laboratory of Multidimensional Omics and Molecular Enzymology 2025 Open to the Outside World Fund projects (No. 2025YB007), the Dalian Science and Technology Talent Innovation Support Plan (No. 2024RQ034), and the Dalian Life and Health Guidance Program Project (No. 2024ZDJH01PT068).

Author contributions

G.W. conceived the experimental design, and X.Q. and G.L. organized the data and wrote the article. X.C. and X.R. provided the experimental platform and financial support. J.S., Y.H., and Y.L. contributed to the writing of the article. Y.L. assisted and guided the experimental operation part. Q.W. and F.C. provided subsequent proofreading of the article.

Declaration of interests

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Declaration of generative AI and AI-assisted technologies in the writing process

During the preparation of this work, the authors used DEEPSEEK to improve the readability and language quality of the article. Following its use, the authors thoroughly reviewed and edited the content and take full responsibility for the integrity and accuracy of the published article.