We got RFX1 overexpression plasmids (GP-PL-RFX1, GenePharma, Shanghai, China) (Supplementary Table 1). We used Lipofectamine 2000 (cat. no. 11668019, Invitrogen; Thermo Fisher Scientific, Inc., Carlsbad, CA, USA) for transfection following the manual. To perform gene silencing, we put small interfering RNA (siRNA) (GP-siR-RFX1, GenePharma, Shanghai, China) for RFX1 into HP75 and AtT20 cells (authenticated by STR profiling and tested negative for mycoplasma, sourced from the Cell Bank of the Chinese Academy of Sciences) with Lipofectamine 2000. Also, we used RT-qPCR and Western blot to check if the overexpression and knockdown worked.

We seeded cells (1 $\times$ 10⁴/well) into 96-well plates and let them grow for 24 hours. Then, we added CCK-8 (cat. no. CK04, Dojindo Molecular Technologies, Inc., Kumamoto, Japan) reagent. After 2 hours at 37 °C, we used a microplate reader (CCP-96H, Servicebio, Wuhan, Hubei, China) to measure absorbance at 450 nm to see cell viability.

We seeded 1000 cells in each well of 6-well plates and kept them for 14 days. Then, we fixed the colonies with 4% paraformaldehyde and stained them with crystal violet (Sigma-Aldrich; Merck KGaA, Darmstadt, Germany). We used ImageJ (v1.53, National Institutes of Health, Bethesda, MD, USA) to count the results.

We used Transwell inserts coated with 100 µL of Matrigel Matrix (cat. no. 354234, Corning, Inc., Corning, NY, USA). We put them at 37 °C for 30 minutes to make the gel solid. Then, we put the inserts into a 24-well plate to make the upper chamber. We seeded 5 $\times$ 10⁴ cells in serum-free medium in the top and 600 µL of complete medium in the bottom. After 48 hours, we fixed the cells on the lower membrane with 4% paraformaldehyde. Finally, we stained them with crystal violet (Sigma-Aldrich; Merck KGaA, Darmstadt, Germany) and counted them under a microscope.

We used Trizol (cat. no. 15596026, Invitrogen; Thermo Fisher Scientific, Inc., Carlsbad, CA, USA) to get total RNA. Then, 1 µg of RNA was turned into complementary DNA using a reverse transcription kit (RR036A, Takara Bio, Inc., Kusatsu, Shiga, Japan). We did qPCR with PowerUp SYBR Green Master Mix (A25742, Applied Biosystems (Thermo Fisher Scientific, Inc., Austin, TX, USA)) on a Quantitative real-time PCR system (QuantStudio 1, Thermo Fisher Scientific, Inc., Waltham, MA, USA). The settings were: 95 °C for 5 minutes, then 40 cycles of 95 °C for 10 seconds and 60 °C for 30 seconds. We calculated gene expression using the 2^{-(Δ⁢Δ⁢CT)} method and used GAPDH as a control. The specific primer sequences used in this study are listed in Supplementary Table 2.

We used RIPA lysis buffer (P0013B, Beyotime Biotechnology, Shanghai, China) with inhibitors to lyse the cells. We used SDS-PAGE to separate the proteins and moved them to PVDF membrane (ISEQ00010, Merck Millipore, Billerica, MA, USA) membranes using a semi-dry system. After using blocking buffer, we put the membranes with primary antibodies (Anti-RFX1 primary antibody (1:1000 dilution, 26859-1-AP, Proteintech, Wuhan, Hubei, China); Anti-GAPDH primary antibody (1:50,000 dilution, 81640-5-RR, Proteintech); Anti-Cyclin D1 primary antibody (1:2000 dilution, 26939-1-AP, Proteintech); Anti-CDK2 primary antibody (1:1000 dilution, 10122-1-AP, Proteintech)) at 4 °C overnight. Then, we used HRP-conjugated goat anti-rabbit IgG secondary antibody (1:50,000 dilution, ZB-2301, Beijing Zhongshan Golden Bridge Biotechnology Co., Ltd., Beijing, China) for 1 hour at room temperature. We saw the protein bands using NcmECL Ultra (NCM Biotech, Suzhou, China) and a ChemiDoc XRS+ imaging system (ChemiDoc XRS+, Bio-Rad Laboratories, Inc., Hercules, CA, USA). We used ImageJ (v1.53) to measure the data.

We used the AddModuleScore function in Seurat to measure cuproptosis activity. We looked for differential gene expression using the Wilcoxon rank-sum test (p 0.05). For genes linked to RFX1 in normal cells, we did GO and KEGG enrichment with clusterProfiler.

We did all stats in R software (v4.3.1, R Core Team, Vienna, Austria). We showed data as mean $\pm$ standard deviation. We used the Bartlett test for variance homogeneity and Spearman correlation to see how variables relate. A two-tailed p 0.05 was the sign of statistical significance.

Analysis of 19 CRGs in invasive versus non-invasive PitNET samples revealed nine significantly dysregulated genes: ATP7B (log2FC = 0.225, p 2.403693 $\times$ 10^-7), LIAS (log2FC = 0.097, p 2.680955 $\times$ 10^-4), PDHA1 (log2FC = 0.119, p 0.003273), MTF1 (log2FC = 0.063, p 0.008561), and DLST (log2FC = 0.195, p 0.007547) were upregulated in invasive tumors, whereas ATP7A (log2FC = –0.096, p 0.005471), FDX1 (log2FC = –0.177, p 3.266440 $\times$ 10^-4), DLAT (log2FC = –0.073, p 0.025148), and GLS(log2FC = –0.439, p 1.065999 $\times$ 10^-5) were downregulated (Supplementary Fig. 1A,B and Supplementary Fig. 2A). Correlation analysis identified coordinated shifts, such as the synergistic effect between PDHA1 and DBT, suggesting cuproptosis-related metabolic rewiring (Supplementary Fig. 2B,C). This dysregulation coincided with significant immune microenvironment remodeling. Cell-type Identification By Estimating Relative Subsets Of RNA Transcripts (CIBERSORT) analysis showed higher CD8+ T cell infiltration and reduced levels of M2 macrophages and naïve CD4+ T cells in invasive PitNETs (Supplementary Fig. 2D,E). Notably, FDX1 expression was positively correlated with M0 macrophage levels, while DLST negatively associated with plasma cell counts, indicating that CRG status may modulate the tumor immune landscape (Supplementary Fig. 2F).

Consensus clustering was performed based on the expression profiles of the nine differentially expressed cuproptosis-related genes. The optimal number of clusters was determined to be k = 2, as indicated by the highest clustering stability (Supplementary Fig. 1E–H). The cumulative distribution function (CDF) curves fluctuated within a minimal consensus index range of 0.03 to 0.55, demonstrating robust clustering stability Additionally, variations in the area under the curve (AUC) between k-1 and k further supported the selection of k = 2, as observed in the CDF curve analysis across values ranging from 2 to 6 (Supplementary Fig. 1I). Consequently, invasive PitNET samples were categorized into two distinct subtypes, designated as Cluster 1 and Cluster 2. Principal component analysis (PCA) further validated the substantial divergence between these two clusters, highlighting distinct transcriptional patterns associated with the differentially expressed cuproptosis-related genes (Supplementary Fig. 1J). Notably, Cluster 1 exhibited higher proportions of CD8⁺ T cells and resting mast cells, whereas Cluster 2 was enriched with follicular helper T cells and activated mast cells (Supplementary Fig. 2H). Enrichment analysis further revealed that starch and sucrose metabolism as well as mismatch repair pathways were prominently enriched in Cluster 2 (Supplementary Fig. 2I). These distinct signatures suggest that CRG-based subtyping effectively captures the metabolic and immunogenic heterogeneity of invasive PitNETs.

Six-week-old BALB/c nude mice were obtained from GemPharmatech Co., Ltd. and maintained in a specific pathogen-free (SPF) environment. All mice were allowed a 1-week acclimation period after being transferred to the experimental facility before being used in the study. For the subcutaneous xenograft assay, the mice were randomly stratified into distinct groups, and 2 × 10⁷ target cells were subcutaneously injected into the right axillary region of each mouse. The tumor size and body weight of the mice were measured every 2 to 3 days using a digital caliper. At the conclusion of the experiment, the tumors were excised and weighed.

Weighted Gene Co-expression Network Analysis (WGCNA) was performed to construct co-expression networks for both invasive and non-invasive PitNET samples. Initially, variance computations were conducted across all genes within the training dataset, followed by further analysis of the top 25% of genes exhibiting the highest variance. Co-expressed gene modules were identified by setting the soft-threshold power to 6, achieving a scale-free R² value of 0.9 (Supplementary Fig. 3A). The dynamic tree-cutting algorithm partitioned the genes into six distinct co-expression modules, each assigned a unique color, with the topological overlap matrix heatmap further illustrating module connectivity (Supplementary Fig. 3B–D).

To evaluate the relationship between gene modules and clinical traits, correlation analysis was performed across the six identified modules. The turquoise module exhibited the strongest association with invasive PitNET samples (Supplementary Fig. 3E), suggesting its potential involvement in invasive PitNETs pathogenesis. Moreover, correlation analysis revealed a significant positive association between turquoise module genes and module-related genes, reinforcing their functional relevance in invasive PitNET (Supplementary Fig. 3F).

Additionally, WGCNA was applied to identify key gene modules specifically associated with cuproptosis-related clusters. Following the determination of optimal soft-thresholding parameters (soft threshold = 4, R² = 0.9), a scale-free network was successfully constructed (Supplementary Fig. 3G). The analysis identified six major co-expression modules, with their Topological Overlap Matrix (TOM)-based heatmap visualization further delineating inter-module relationships (Supplementary Fig. 3H–J). Correlation analysis of module-clinical trait relationships indicated that invasive PitNETs-related clusters were most strongly associated with the green module (Supplementary Fig. 3K). Furthermore, gene-module correlation analysis demonstrated that genes within the green module exhibited strong connectivity and functional relevance within the selected module (Supplementary Fig. 3L).

We found 15 hub genes by looking at the overlap between genes from invasive PitNETs and cuproptosis clusters (Fig. 2A). To make a predictive framework, we built four machine learning models: Random Forest (RF), Support Vector Machine (SVM), Generalized Linear Model (GLM), and Extreme Gradient Boosting (XGB). We used the DALEX R package (v2.5.3, ModelOriented, Warsaw, Poland) to check model performance and investigated the residual distribution in the validation set. Among the four, RF and SVM had the lowest residuals. This means they had the best accuracy (Fig. 2B).

Fig. 2.

Optimization and validation of machine-learning-based predictive models and nomogram construction in invasive pituitary neuroendocrine tumors (PitNETs). (A) Identification of overlapping genes between disease-WGCNA and cluster-WGCNA. (B) Box plot comparing the residuals of the four candidate models. (C) Key performance characteristics of the evaluated machine learning models. (D) Reverse cumulative distribution curves for the four models. (E) ROC curves assessing the performance of Random Forest (RF), Support Vector Machine (SVM), Extreme Gradient Boosting (XGB), and Generalized Linear Model (GLM) models. (F) Nomogram for predicting the risk of invasive PitNET clusters based on the Random Forest (RF) model. (G) Calibration curve assessing the predictive accuracy of the nomogram. (H) Decision curve analysis evaluating the clinical utility of the nomogram. (I) ROC curve of the five-gene RF model, demonstrating its predictive performance.

Next, we ranked the top 15 genes in each model using the root mean square error (Fig. 2C). We used ROC analysis to check how well the models classified the data. The RF model had the highest AUC (0.902). This shows it has a very good predictive ability (Fig. 2D,E). By looking at both residuals and AUC, the RF model was the best at finding invasive PitNET patients. So, we picked the top five genes from the RF model—WLS, WFDC2, RFX1, TAGLN3, and TGFBR3L—for more study.

We made a nomogram to help estimate risk for individual invasive PitNET patients (Fig. 2F). We checked the accuracy of the nomogram using calibration curves and decision curve analysis (DCA). The calibration curves showed that the predicted and actual results were very similar (Fig. 2G). Also, the DCA showed that the nomogram gives a net clinical benefit (Fig. 2H). In the test data, the ROC curve for the five genes (WLS, WFDC2, RFX1, TAGLN3, and TGFBR3L) showed a good AUC of 0.767 (Fig. 2I).

We investigated the expression of the five candidate genes in invasive and non-invasive PitNET samples. RFX1 was the only gene that was significantly different between the two groups. It was much lower in the invasive PitNETs (Fig. 3A).

Fig. 3.

RFX1 overexpression inhibits proliferation and invasion of PitNET cells. (A) Relative mRNA expression levels of five key genes in invasive and non-invasive PitNET samples. (B) Validation of RFX1 overexpression and knockdown in constructed cell lines. (C) CCK-8 assay assessing cell viability following RFX1 Overexpression or knockdown. (D) Colony formation assay evaluating the proliferative capacity of cells with RFX1 overexpression or knockdown. (E) Transwell invasion assay measuring the invasive Capacity of cells with RFX1 overexpression or knockdown. (F) Representative images of xenograft tumors from nude mice. (G) Tumor weight measured at the end of the experiment. (H) Tumor volume measured over time. Scale bar, 100 µm. *p 0.05, **p 0.01, ***p 0.001, ****p 0.0001; ns, not significant (p $\ge$ 0.05).

To further investigate the functional significance of RFX1, we did overexpression and knockdown tests (Fig. 3B). The tests showed that making more RFX1 stopped PitNET cells from growing and invading. But, taking away RFX1 made the tumor cells grow and spread more (Fig. 3C–E. Supplementary Fig. 4A,B). Mice were anesthetized by intraperitoneal injection of 50 mg/kg pentobarbital sodium (cat. no. P3761, Sigma‑Aldrich, St. Louis, MO, USA) before tumor implantation. Also, the in vivo xenograft assay showed that RFX1 overexpression slowed down tumor growth in mice (Fig. 3F–H).

Copper is an important micronutrient for cell growth and metabolism. To see how copper affects PitNET cells, we treated them with copper ions and carriers. The results showed that tumor cells are sensitive to copper. This effect depended on the dose of Elesclomol (ELES), a copper ionophore, from 0 to 200 nM. Also, using 5 µM tetrathiomolybdate (TTM), which is a copper chelating agent, stopped the effects of ELES (Fig. 4A) [17, 18]. These findings show a link between cuproptosis and PitNET biology.

Fig. 4.

RFX1 overexpression enhances cuproptosis in PitNET cells. (A) The inhibitory effect of ELES-Cu on PitNET viability was dose-dependent, and this effect was reversible by tetrathiomolybdate (TTM) treatment. (B) PitNET cell viability decreased in a dose-dependent manner after ELES-Cu treatment. Restoring RFX1 expression promoted cuproptosis, whereas knocking down RFX1 inhibited cuproptosis. (C) Western blot analysis of CYCLIND1 and CDK2 in PitNET cells treated with ELES-Cu for 24 h, with or without RFX1 overexpression or knockdown. (D) Comparison of PitNET cell viability with or without ELES-Cu treatment for 24 h, and with or without RFX1 overexpression or knockdown. (E) Western blot analysis of CYCLIND1 and CDK2 in PitNET cells with or without ELES-Cu treatment, and with or without RFX1 overexpression or knockdown. ***p 0.001.

We also investigated how RFX1 regulates cuproptosis. Our results showed that cells with more RFX1 were more sensitive to Elesclomol-Cu and had lower cell viability. But, cells with low RFX1 survived better when exposed to copper (Fig. 4C). Western blot tests for CYCLIND1 and CDK2 also showed changes in these cell cycle proteins (Fig. 4B, Supplementary Material-original images of western blot).

When we used the same dose of Elesclomol-Cu, the cells with more RFX1 had more cell death. In contrast, cells with RFX1 knockdown had less cell death than the control group (Fig. 4D). Western blot data for CYCLIND1 and CDK2 supported these results again (Fig. 4E).

We combined 20,591 cells from normal (n = 7) and PitNET (n = 47) groups. We found 11 clusters, including tumor cells and seven normal subtypes (Fig. 5A). The Uniform Manifold Approximation and Projection (UMAP) plot showed clear groups for tumor and normal cells (Fig. 5A, right). Tumor clusters mostly came from PitNET patients, while normal clusters came from controls (Fig. 5B). Marker analysis confirmed the cell types, such as POU1F1 for tumor cells and POMC for normal cells (Fig. 5C).

Fig. 5.

Integrated single-cell profiling links RFX1 loss to cuproptosis dysregulation in PitNET. (A) Uniform Manifold Approximation and Projection (UMAP) plots: Left, cell type annotation (11 clusters); Right, tumor/normal origin. (B) Stacked bar plot showing sample source proportions across clusters. (C) Dot plot displaying top 3 marker genes per cluster. Tumor cells unified as “Tumor cell”. (D) Violin plots of cuproptosis scores across clusters. (E) Left: Violin plot comparing RFX1 expression (normal vs. tumor); Right: UMAP feature plot of RFX1. (F) Enriched pathways for top 3 markers per cluster (GO/KEGG). (G) Enrichment results for RFX1-correlated genes in normal cells, highlighting mitochondrial metabolism pathways.

We calculated cuproptosis scores using 19 genes. The scores were lower in tumor clusters than in normal cells (Fig. 5C,D). RFX1 is a transcription factor linked to mitochondrial metabolism. It was almost absent in tumor cells but high in normal cells (Fig. 5E). Feature plots showed RFX1 was mainly in normal clusters (Fig. 5E, right). Multiplex immunofluorescence on clinical samples showed that RFX1 protein levels were low when Topoisomerase II alpha (TOP2A) and marker of proliferation Ki-67 (Ki-67) were high (Supplementary Fig. 4C). In normal cells, we found 98 genes linked to RFX1. These genes are involved in mitophagy, the TCA cycle, and oxidative phosphorylation (Fig. 5F). Pathway enrichment showed that RFX1-linked genes in normal cells relate to mitochondrial autophagy and ATP-Binding Cassette (ABC) transporters (Fig. 5G). These processes are part of copper balance and cell death.