Raw RNA sequencing data from 33 different cancer tissues and normal tissues were obtained from the TCGA and GTEx databases for pan-cancer analysis. The TCGA-KIRC dataset’s corresponding clinical data was also downloaded. The raw data was acquired in HTSeq-FPKM format, and converted to TPM format, followed by log2 transformation.

In this study, 16 pairs of KIRC tumor tissues and paired normal tissues (derived from patients who underwent surgery in Shandong Provincial Hospital from 2019 to 2021) were collected for quantitative polymerase chain reaction (qPCR) and immunohistochemistry (IHC) to verify the expression level of C1QTNF1. All the patients were aware of the intent of the study and provided written informed consent. This study followed the Declaration of Helsinki and was approved by the Hospital Medical Ethics Committee of Shandong Provincial Hospital Affiliated to Shandong First Medical University.

KIRC patients were separated into C1QTNF1 high expression and low expression groups based on the median C1QTNF1 expression in the TCGA-KIRC dataset. Differential analysis between the high and low expression groups was conducted using the R package “DESeq2”. A threshold of p values <0.05 and |logFC| > 1.5 was set. The Spearman correlation coefficient method was employed to analyze the correlation between C1QTNF1 and the top 10 differentially expressed genes (DEGs). The results were presented using volcano plots and heat maps created using the R packages “ggplot2” and “pheatmap”.

The relationship between C1QTNF1 expression in KIRC samples and the abundance of 24 immune cells in tumors was evaluated using the ssGSEA algorithm in the “GSVA” R package, including neutrophils, cytotoxic cells, dendritic cells (DCs), CD8⁺ T-cells, plasmacytoid DC (pDC), natural killer (NK) cells, mast cells, T gamma delta (Tgd), type 17 Th (Th17) cells, immature DCs (iDCs), eosinophils, NK CD56dim cells, regulatory T-cells (TReg), T effector memory (Tem), T-cells, T central memory (Tcm), B cells, type 1 Th (Th1) cells, macrophages, NK CD56bright cells, activated DC (aDC), T follicular helper (TFH), T helper cells, and type 2 Th (Th2) cells. Correlation analysis between C1QTNF1 and immune checkpoints was conducted via the “ggplot2” R package and Spearman correlation coefficient. When p < 0.05, the correlation was significant. TISIDB (http://cis.hku.hk/TISIDB/) is a portal containing a variety of tumor immunology data resources, which can be used to explore the relationship between genes and immunological features (Ru et al., 2019). We used TISIDB to determine the association of C1QTNF1 with TIL infiltration as well as immune checkpoints.

The GEPIA2 database (http://gepia2.cancer-pku.cn/#index) is a free and open online database that can be used by users to determine the expression level and prognostic impact of target molecules in different tumors (Tang et al., 2019). The Kaplan-Meier Plotter online site assesses the relationship between genes and survival outcomes in a variety of cancers, and users can look for the prognostic value of specific genes analyzed (Lánczky and Győrffy, 2021). The GEPIA2 database was used to analyze the prognostic significance of C1QTNF1 expression in pan-cancer. For the TCGA-KIRC dataset, prognostic analysis was performed using R packages “survival” and “survminer”. KIRC patients were divided into a C1QTNF1 high expression group and a low expression group according to the median value of C1QTNF1 expression. The prognostic value of the C1QTNF1 expression level in KIRC was evaluated by the Kaplan-Meier curve and log-rank test. The Kaplan-Meier Plotter was used to verify the relationship between C1QTNF1 expression level and the survival outcome of KIRC patients. Based on Cox regression analysis, the independent risk factors affecting the prognosis of KIRC were analyzed, and then the factors obtained by multivariate Cox regression were included in the prognostic nomogram to evaluate the prognosis of KIRC patients. The discrimination of the nomogram was then quantified using the concordance index (C-index), and the performance of the nomogram was evaluated using a calibration plot. The R package “rms” was used to draw the nomogram and calibration plots. “pROC” and “timeROC” were used to draw the diagnostic ROC curve and the time-dependent survival ROC curve. The prognosis analysis results of KIRC subgroups were visualized by the R package “ggplot2.”

GeneMENIA (http://www.genemania.org) is a flexible, user-friendly online tool that allows users to upload their own data sets to predict the priority and function of genes (Warde-Farley et al., 2010). Genes co-expressed with C1QTNF1 were predicted by the GeneMANIA online tool. The Comparative Toxicogenomic Database (CTD, http://ctdbase.org/) is a friendly, innovative database containing toxicological information on chemicals, genes, phenotypes, diseases, and exposures that provides potential molecular mediators to develop testable hypotheses (Davis et al., 2023). CTD was used to predict the chemicals interacting with C1QTNF1. ENCORI (http://starbase.sysu.edu.cn/) is a comprehensive online database that allows users to explore the miRNA-mRNA and miRNA-lncRNA interaction networks and predict the function of miRNAs and other ncRNAs (Li et al., 2014). The authors predicted the upstream potential miRNAs and lncRNAs that interacted with C1QTNF and hsa-miR-27b-3p through the RNA Interactomes Encyclopedia (ENCORI) database. The interaction network was visualized by Cytoscape software.

Total RNA was extracted from cells and tissues with the SteadyPure Universal RNA Extraction Kit II (Accurate Biotechnology, Hunan, China). Reverse transcription of RNA into cDNA was performed with Evo M-MLV RT Premix (Accurate Biotechnology, Hunan, China) according to the manufacturer’s instructions. qPCR assays were performed using the SYBR^® Green Premix Pro Taq HS qPCR Kit (Accurate Biotechnology, Hunan, China) and amplified in a LightCycler 480II (Roche). The primers used in this study were as follows: β-actin-F: AAG​TGT​GAC​GTG​GAC​ATC​CGC, β-actin-R: CCG​GAC​TCG​TCA​TAC​TCC​TGC​T, C1QTNF1-F: AAG​TTC​TAC​TGC​TAC​GTG​CCC, C1QTNF1-R: TGT​GCA​GGT​AGG​TCT​CCT​TCT.

Paraffin blocks containing tissue were cut into 4 μm sections and placed on slides. After deparaffinization and rehydration, they were incubated in 3% hydrogen peroxide to block endogenous peroxidase. Slides were next incubated with the C1QTNF1 (Proteintech, 12209-1-AP, Wuhan, China) antibody at a 1:50 dilution at 4°C overnight. Finally, the slides were washed with secondary antibody at room temperature for 30 min and stained with DAB and hematoxylin.

786-O cells were purchased from the Cell Bank of the Chinese Academy of Sciences. 786-O cells were placed in RPMI-1640 medium (Gibco; Thermo Fisher Scientific, Inc., Waltham, MA, United States) supplemented with 10% fetal bovine serum (FBS) and 1% penicillin/streptomycin. The above cells were cultured in a 5% CO₂ cell incubator at 37°C. siRNA targeting C1QTNF1 was synthesized by Sangon Bio (Shanghai, China) and transfected using INTERFERin^® (Polyplus, Illkirch, France) according to the manufacturer’s instructions.

The treated cells were seeded in 96-well plates with 2000 cells per well, and 10 μL CCK8 reagent was added to the corresponding wells every 24 h, respectively. The cells were incubated in a 5% CO₂ incubator at 37°C for 1 h. The optical density (OD) of each well was measured using a microplate reader at a wavelength of 450 nm. The appropriate cells were seeded in 6-well plates, and when the cells grew to the appropriate density, the monolayer cells were scratched with the suction tip of 200 μL pipetter, and the cells were photographed by microscope 24 h after injury.

All statistical analyses were performed using R (V 3.6.3). Wilcoxon rank-sum test was used to investigate the levels of C1QTNF1 in both tumor and noncancer samples. The chi-square test, or Fisher exact test, was used to compare and analyze the statistical significance between the two groups of categorical variables. Univariate and multivariate analyses were performed using the Cox proportional hazards model to evaluate the independent prognostic factors of KIRC patients. Only p values of less than 0.05 (two-sided) were considered to indicate statistical significance. For data on C1qTNF1-related functions, statistical analysis was performed using GraphPad Prism 8.0.

The authors evaluated the expression levels of C1QTNF1 mRNA in 33 different cancer tissues and normal tissues using TCGA and GTEx data. Among them, C1QTNF1 expression levels were increased in several tumor tissues, including Lymphoid Neoplasm Diffuse Large B-cell Lymphoma (DLBC), Glioblastoma multiforme (GBM), Head and Neck squamous cell carcinoma (HNSC), and KIRC (Figure 1A). Meanwhile, we investigated the expression levels of C1QTNF1 in paired tumor tissues and normal tissues in the TCGA dataset and found that the expression levels of C1QTNF1 were increased in Colon adenocarcinoma (COAD), HNSC, and KIRC (Figure 1B). Kaplan-Meier survival curves showed that high expression of C1QTNF1 was significantly associated with a worse prognosis in Bladder Urothelial Carcinoma (BLCA), Brain Lower Grade Glioma (LGG), KIRC, and Uveal Melanoma (UVM) (Figure 1C).

By comparing the expression of C1QTNF1 in unpaired and paired KIRC tumor tissues and normal tissue samples, we found that the expression level of C1QTNF1 was significantly increased (Figures 2A, B). The ROC curve showed that C1QTNF1 could effectively diagnose KIRC (Figure 2C). Survival analysis suggested that KIRC patients with high C1QTNF1 expression had a poor prognosis (Figures 2D–F). In addition, the prognostic analysis of C1QTNF1 using the Kaplan-Meier Plotter online website also obtained the same results (Supplementary Figure S1).

As shown in Figures 3A–I, C1QTNF1 expression was significantly correlated with TNM stage, Pathologic stage, age, sex, Overall Survival (OS) events, Disease Specific Survival (DSS) events, and Progress Free Interval (PFI) events (Table 1). Logistic regression analysis also indicated that NCAPG2 expression was significantly correlated with poor clinicopathological prognosis, including T stage (T3, T4 vs. T1, T2), N stage (N1 vs. N0), M stage (M1 vs. M0), Pathological stage (III and IV vs. I and II), and age (>60 vs.≤60) (Table 2).

We analyzed the gene expression profile between the C1QTNF1 high expression group and the low expression group, and 357 differentially expressed genes (DEGs) were obtained. These included 170 downregulated genes and 187 upregulated genes (Figure 4A). The heat map shows the expression of the top five genes in the high versus low expression groups (Figure 4B). In addition, we performed Gene Ontology (GO) functional enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis of differentially expressed genes between high and low-expression groups (Tables 3, 4). The main biological processes included cornification, acute-phase response, keratinization, isoprenoid metabolic process, terpenoid metabolic process, extracellular structure organization, monovalent inorganic cation homeostasis, anion transmembrane transport, hormone metabolic process, and diterpenoid metabolic process (Figure 4C). The most enriched cellular components were collagen-containing extracellular matrix, apical part of cell, Golgi lumen, high-density lipoprotein particle, endoplasmic reticulum lumen, basolateral plasma membrane, apical plasma membrane, vacuolar proton-transporting V-type ATPase complex, intermediate filament, and plasma lipoprotein particle (Figure 4D). The most enriched molecular functions were serine-type endopeptidase activity, receptor ligand activity, anion transmembrane transporter activity, heparin binding, serine-type peptidase activity, serine hydrolase activity, extracellular matrix structural constituent, inorganic anion transmembrane transporter activity, structural constituent of cytoskeleton, and glycosaminoglycan binding (Figure 4E). The KEGG pathway was enriched in Retinol metabolism, Collecting duct acid secretion, Steroid hormone biosynthesis, Complement, and coagulation cascades, Drug metabolism - cytochrome P450, Metabolism of xenobiotics by cytochrome P450, Synaptic vesicle cycle, Chemical carcinogenesis, Rheumatoid arthritis, and Linoleic acid metabolism (Figure 4F).

We explored the relationship between C1QTNF1 and immune cell infiltration based on the ssGSEA method. C1QTNF1 and NK cells, pDC, Tem, Th2 cells, Th1 cells, Macrophages, DC, TReg, Mast cells, B-cells, Tgd, NK CD56dim cells, TFH, Cytotoxic cells, aDC, iDC, T-cells, CD8 T-cells, and NK CD56 bright cells were positively correlated, while Th17 cells were negatively correlated (Figure 5A). In addition, we analyzed the infiltration levels of several immune cells with the highest correlation between the C1QTNF1 high expression group and the low expression group and found that the infiltration levels of NK cells, pDC, Tem, and Th17cells were statistically significant between the C1QTNF1 high expression group and the low expression group (Figures 5B–E). Finally, we also found that two immune checkpoints, PDCD1 and CTLA4, were positively correlated with the expression level of C1QTNF1 (Figures 5F, G). In addition, we found a strong association between C1QTNF1 and TIL infiltration through the TISIDB database, and it was also positively correlated with two immune checkpoints, PDCD1 and CTLA4 (Supplementary Figure S2). This suggests that C1QTNF1 may be closely related to the immune infiltration of tumors.

Analysis of the prognostic impact of C1QTNF1 on different clinical subgroups showed that C1QTNF1 was a risk factor for stages T1 and T2, T3 and T4, N0, M0, and Stage I and Stage II (HR > 1, p < 0.05) (Figure 6A). Time-dependent ROC curve analysis showed that the AUC values for predicting the 1-, 3-, and 5-year survival rates of KIRC patients based on C1QTNF1 expression levels were all greater than 0.5 (Figure 6B). Univariate and multivariate Cox analyses were performed by combining multiple clinical factors with C1QTNF1 to finally determine the independent factors predicting the prognosis (Table 5). The predictors obtained above were used to construct a prediction model, and a nomogram was drawn to assess the risk probability (Figure 6C). In addition, the calibration plots showed that the model had good predictive value for the survival time of patients at 1, 3, and 5 years (Figure 6D).

To verify the expression level of C1QTNF1 in KIRC, the authors examined the expression of C1QTNF1 in clinical samples using qRT-PCR. The results confirmed that C1QTNF1 was elevated in KIRC tissues (Figure 7A). Immunohistochemical staining results also showed that the expression of C1QTNF1 protein was increased in KIRC tissues (Supplementary Figure S3). Subsequently, we knocked down C1QTNF1 in KIRC cell lines, and qRT-PCR showed that the expression of C1QTNF1 was significantly reduced after knockdown (Figure 7B). CCK-8, wound healing assay, and Transwell assay confirmed that knockdown of C1QTNF1 inhibited the proliferation, invasion, and migration of KIRC cells (Figures 7C–E).

Using GeneMANIA, we generated a gene-gene network of 20 potential C1QTNF1 interactions, They are AVPR2, SGTA, HEYL, C1QTNF6, MAFF, SOD3, BDKRB2, PRPF31, GADD45G, MFGE8, CST3, IRF4, ABLIM1, GSN, FSTL3, SF3A1, CMKLR1, ITGA7, SYNPO, and PLA2G 5 (Supplementary Figure S4A). Chemical agents associated with the C1QTNF1 gene were predicted from the CDC database (Supplementary Figure S4B) (Supplementary Table S1). Increasing evidence has confirmed that ncRNAs can regulate gene expression in various ways in tumors (Qi et al., 2015). To determine whether C1QTNF1 is regulated by certain ncRNAs, we predicted the upstream miRNAs that might bind to C1QTNF1 and finally found 28 miRNAs (Supplementary Figure S4C) (Supplementary Table S2). Due to the mechanism by which upstream miRNAs negatively regulate C1QTNF1 expression at the post-transcriptional level, there should be a negative relationship between C1QTNF1 and upstream miRNAs. Therefore, the correlation between C1QTNF1 and 28 miRNAs was examined in the ENCORI database. The results suggested a negative correlation between hsa-miR-27b-3p and C1QTNF1 (Supplementary Figure S4D). We then determined the expression level of hsa-miR-27b-3p in the TCGA-KIRC database. The results showed that the expression level of hsa-miR-27b-3p in KIRC was lower than that in adjacent normal tissue controls (Supplementary Figure S4E). We also explored the correlation between the hsa-miR-27b-3p expression level and the prognosis of KIRC patients. The results showed that high expression of hsa-miR-27b-3p was significantly associated with a good prognosis of KIRC (Supplementary Figure S4F). Combined with correlation analysis, expression analysis, and survival analysis, we suggest that hsa-miR-27b-3p may be the most likely upstream regulatory miRNA of C1QTNF1 in KIRC.

Next, ENCORI was used to predict the lncRNAs upstream of hsa-miR-27b-3p. Finally, 143 possible lncRNAs were obtained (Supplementary Table S3). The competing endogenous RNA (ceRNA) hypothesis suggests that lncRNA competitively binds tumor suppressor miRNAs to reduce the inhibitory effect of miRNAs on target mRNAs. It can be seen that there is a negative correlation between lncRNA and target miRNA in the ceRNA network, while there is a positive correlation between lncRNA and target mRNA. Therefore, we evaluated the expression correlation of hsa-miR-27b-3p/C1QTNF1 and 143 lncRNAs in the ENCORI database. The results suggested that CYTOR, AC040970.1, AC016717.2, AC010980.2, and LINC02381 were negatively correlated with the expression of hsa-miR-27b-3p, but positively correlated with the expression of C1QTNF1 (Supplementary Figures S5A, B). Subsequently, we performed expression analysis and prognostic analysis of the above 5 lncRNAs in KIRC (Supplementary Figures S5C, D). CYTOR and AC040970.1 may be the two most potential lncRNAs upstream of the hsa-miR-27b-3p/C1QTNF1 axis in KIRC.