Transcriptome RNA-sequencing and clinical data were from the TCGA-UVM dataset (https://portal.gdc.cancer.gov/) (24). Additionally, gene expression and clinical information files of GSE22138 were downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) (5).

Based on the RNA-seq data of 80 samples from the TCGA-UVM dataset, the “ESTIMATE” R package v1.0.13 was used to calculate ESTIMATE, immune, and stromal scores in the TME for all samples (25). The samples were divided into two groups based on the medians of these three scores. Thereafter, Kaplan–Meier analysis and the log-rank test (significance level: p < 0.05) were used to compare overall survival (OS) in two groups of patients using the R packages “survminer” v0.4.9 and “survival” v3.2-13.

Afterward, the “limma” package v3.48.3 in R was used to assess differential gene expression in the high and low immune/stromal score groups. The criteria for selecting differentially expressed genes (DEGs) included |log2(fold change)| >1 and a false discovery rate (FDR)-adjusted p <0.05. Results from these analyses were plotted as volcano plots using the “ggplot2” package v3.3.5 in R. Co-upregulated and co-downregulated DEGs in the high and low immune/stromal score groups were regarded as TME- and prognosis-related genes. These genes were subjected to the Gene Ontology (GO) functional enrichment and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis, with the results visualized using the R packages “clusterProfiler” v4.0.5 and “enrichplot” v1.12.2 (p <0.05, q <0.05).

Univariate Cox regression analyses (significance level: p <0.05) were performed to assess the relationships between the levels of expression of these potential candidate genes and OS in patients with UM, followed by least absolute shrinkage and selection operator (LASSO) regression analysis using the “glmnet” package v4.1-2 in R to avoid overfitting. This was followed by multivariate Cox regression analysis (significance level: p <0.05) to determine whether there was any independent prognostic effect of these candidate genes on OS.

Patients in the TCGA and GEO cohorts were divided into those with high and low levels of expression of PRRX1 based on the median levels of PRRX1 in each dataset. OS in the TCGA cohort and metastasis-free survival (MFS) in the GEO cohort were analyzed by the Kaplan–Meier method and compared in each pair of groups by the log-rank test (significance level: p <0.05), using the “survival” package v3.2-13 in R. The predictive power of the model was evaluated using receiver operating characteristic (ROC) and Cox regression (significance level: p <0.05) analyses, as calculated using the R packages “survival” v3.2-13, “survminer” v0.4.9 and “timeROC” v0.4. The prognostic impact of PRRX1 expression in several types of cancer was investigated by pan-cancer analysis. RNA-seq data and clinical information come from patients with 33 types of tumors were downloaded from the TCGA database (https://portal.gdc.cancer.gov/). Univariate Cox regression (significance level: p <0.05) analysis was used to assess the relationship between OS and PRRX1 expression, with the results displayed as forest plots using the “forestplot” package v2.0.1 in R. In addition, because cytogenetic studies have found that the loss of chromosome 3 (monosomy 3) and a gain of chromosome 8q are associated with UM metastasis (26) and the status of chromosomes 3 and 8q has been determined in every patient in the TCGA-UVM project ( Table S1 ) (24), Spearman rank correlation analysis (significance level: p <0.05) was performed to assess the correlation between the number of copies of these chromosomes and PRRX1 expression.

To determine the biological states and processes that correlate with PRRX1 expression in UM, RNA-seq data from the TCGA-UVM dataset were analyzed by gene set enrichment analysis (GSEA) v4.1.0 software using the Hallmark gene sets (http://www.broadinstitute.org/gsea), with significance set at FDR <0.25 and p <0.05.

The level of expression of immune checkpoint genes may correlate with responses to treatment with immune checkpoint inhibitor therapies (27). The immune checkpoint genes expressed in UM include CTLA4, PD1 (PDCD1), PDL1 (CD274), TIGIT, and LAG3 (2, 11, 24). Therefore, the correlation between PRRX1 expression and the expression of these critical immune checkpoint molecules was subsequently explored using the “ggplot2” v3.3.5 and “reshape2” v1.4.4 packages in R. In addition, the correlation between PRRX1 expression and tumor mutation burden (TMB) was analyzed by Spearman’s correlation analysis (significance level: p <0.05). To quantify the composition of TME in UM, the abundance of TME-related cell populations was estimated from gene expression data using the “MCPcounter” package v1.2.0 in R.

The human UM cell line MuM-2B purchased from iCell Bioscience Inc. (iCELL-h148;Shanghai, China) was chosen for its highly aggressive phenotype (28, 29). Cells were cultured in RPMI-1640 medium with 10% FBS and 1% penicillin–streptomycin (Gibco, USA) at 37°C with 5% CO₂ condition. Experiments were performed with cells in the logarithmic growth period.

PRRX1 small interfering (si) RNA (si-PRRX1) and negative control (si-NC) oligonucleotide were purchased from RiboBio (Ribobio, Guangzhou, China) with the sequences GGAATAGGACAACCTTCAA (Si1), ACACTATCCTGATGCTTTT (Si2), and GTTCCGCAGGAATGAGAGA (Si3). Transient transfection with siRNA at a final concentration of 100 nM was performed with Lipofectamine 2000 (Invitrogen, USA) according to manufacturer’s protocol.

Total protein from cells was harvested using radio immunoprecipitation assay buffer (RIPA; Beyotime, China), and was quantified by BCA Assay (Beyotime, China). In brief, protein samples were separated by SDS–polyacrylamide gel electrophoresis and transferred onto polyvinylidene fluoride (PVDF) membranes, which was then incubated with primary antibodies and HRP-conjugated secondary antibodies ( Table S2 ). Band intensity was quantified by Quantity One software (Bio-Rad) and defined as the ratio of target protein to β-actin. The assay was repeated three times.

The migration of UM cells was evaluated with a wound healing scratch assay. Cells were seeded in 6-well plates (5 × 10⁵/well) until reached full confluence. Cells were serum-starved overnight and an artificial scratch wound was created at the center of the well and photographed, the scratch areas were photographed at 0, 24 and 48 h. Migration index was calculated as the follows: [(original scratch width-scratch width at 24/48 h)/(original scratch width)] × 100%. The assay was repeated three times.

Transwell invasion assays were performed to assess cell invasion. 24-well Millicell hanging cell culture inserts (8.0 µm, Millipore, USA) and Matrigel (BD Biosciences, USA) were used as per manufacturer’s protocol. Approximately 2 × 10⁵ cells in serum-free medium (100 ul, 2 × 10⁶ cells/ml) were added per well to the upper chamber, 500 ul complete medium with 10% FBS in the lower chamber served as a chemoattractant. After 48 h of incubation, cells that had invaded through the matrigel were fixed in 4% paraformaldehyde and stained with 0.1% crystal violet for counting. Three separate fields were counted for each filter with microscope (UOP-DSZ2000X; Chongqing, China). The assay was also repeated three times.

A brief flowchart of our study is shown in Figure 1 . Clinical data and transcriptome expression profiles of 80 patients were downloaded from the TCGA-UVM cohort. Similarly, information about the 63 patients in the validation cohort (GEO-GSE22138) was obtained from the GEO database. The clinical characteristics of these 143 patients are displayed in Table 1 .

After assessment using an algorithm, the immune scores, stromal scores, and ESTIMATE scores for 80 samples in the TCGA-UVM cohort were calculated. For each scoring system, 40 samples each were classified into the high and low score groups based on the median scores, and OS was compared in each pair of samples. OS was found to be significantly lower in patients with higher immune scores (p <0.001, log-rank test; Figure 2A ), stromal score (p = 0.015, log-rank test; Figure 2B ) and ESTIMATE scores (p = 0.001, log-rank test; Figure 2C ), indicating that immune and stromal components in the TME are associated with poor prognosis in patients with UM.

To identify genes that might be responsible for these differences in OS, gene expression data were analyzed in 80 patients from the TCGA-UVM dataset. These analyses identified 1,031 DEGs (834 upregulated and 197 downregulated) between the high and low immune score groups ( Figure 3A ) and 1,367 DEGs (1,311 upregulated and 56 downregulated) between the high and low stromal score groups ( Figure 3B ). Comparison of these DEGs found that 739 overlapped, comprised of 48 co-down regulated genes and 691 co-upregulated genes ( Figure 3C ). These findings suggested that these 739 genes were likely to be key factors that affect both the prognosis and TME of UM. Furthermore, GO enrichment analysis of these 739 genes found that several GO terms were enriched in immune activity. The ten most significant terms in biological processes (BP), molecular functions (MF), and cellular components (CC) are shown in Figure 3D . KEGG analysis also showed that these DEGs mainly involved immune-related pathways ( Figure 3E ), further indicating that the immune composition within the TME is significantly prognostic of survival in patients with UM.

Univariate Cox regression analysis showed 450 genes potentially prognostic of survival in patients with UM ( Table S3 ). These genes were input into the LASSO Cox regression model ( Figure 4A ), with optimal performance attained by a combination of eight genes (CTF1, LFNG, ZNF667, ARC, CCL24, PRRX1, PARP8, and ISG20) ( Figure 4B ). Multivariate Cox analysis further showed that five of these genes (LFNG, ZNF667, CCL24, PRRX1, and PARP8) were independently prognostic of survival in patients with UM ( Figure 4C ). Previous studies have demonstrated that PRRX1, serve as a EMT-TF, is involved in the EMT process and tumor progression in several cancers (30–32), however, the role of PRRX1 in UM was not reported. We therefore chose to focus our analysis on PRRX1.

Based on median PRRX1 expression as the cut-off value, OS was found to be poorer in patients with high than low PRRX1 expression ( Figure 5A ). The negative prognostic impact of higher PRRX1 expression was validated in an independent GEO cohort, with MFS being significantly poorer with high than low PRRX1 expression ( Figure 5B ). The distribution of PRRX1 expression profiles is shown in Figures 5C, D . Figures 5E, F show the variations in patient survival and metastasis as a function of PRRX1 expression. The percentages of deceased patients ( Figure 5E ) and patients with metastasis ( Figure 5F ) were greater and the survival time shorter in the high than in the low PRRX1 expression group.

To evaluate the prognostic predictive power of PRRX1 expression, univariate and multivariate Cox analyses were performed to assess the ability of PRRX1 expression to predict OS in the TCGA-UVM cohort and to predict MFS in the GEO-GSE22138 cohort. Both univariate (hazard ratio [HR]: 2.324, 95% confidence interval [CI]: 1.514–3.566, p <0.001) and multivariate (HR: 2.571 95% CI: 1.465–4.513, p <0.001) analyses of 75 patients in the TCGA-UVM cohort with complete information on age, gender, stage, T classification, tumor diameter, tumor thickness, and extra scleral extension showed that high PRRX1 expression was a powerful and independent prognostic predictor of reduced OS ( Figure 6A ). The areas under the time-dependent ROC curves for 1-, 2-, and 3-year OS were 0.830, 0.738, and 0.870, respectively ( Figure 6B ). The area under the ROC curve (AUC) of PRRX1 expression for predicting OS was 0.794 ( Figure 6C ), which was less than that for the stage but ranked second among evaluation metrics. Similarly, both univariate and multivariate analyses of the 42 patients in the GEO cohort with complete information, including, age, gender, tumor diameter, tumor thickness, retinal detachment, and extra scleral extension, showed that high expression of PRRX1 was significantly predictive of reduced MFS ( Figure 6D ). The AUCs of the time-dependent ROC curves for 1-, 2-, and 3-year MFS were 0.677, 0.748, and 0.696, respectively ( Figure 6E ), and the AUC of PRRX1 expression for predicting MFS was 0.718 ( Figure 6F ), greater than that for all other factors. These results indicate that PRRX1 expression was an independent prognostic factor that could effectively predict the prognosis of UM patients.

Univariable Cox regression analysis was performed to assess the association of PRRX1 expression levels with OS in patients with various types of cancer in 33 TCGA cohorts. These pan-cancer survival analyses indicated that PRRX1 expression was an adverse prognostic factor in patients with seven types of cancer ( Figure 7A ): adrenocortical carcinoma (ACC), kidney renal clear cell carcinoma (KIRC), kidney renal papillary cell carcinoma (KIRP), lung adenocarcinoma (LUAD), stomach adenocarcinoma (STAD), thyroid carcinoma (THCA), and uveal melanoma (UVM). Assessment of the correlation between PRRX1 expression and chromosome copy number using Spearman correlation analysis showed that PRRX1 expression correlated negatively with chromosome 3 copy number but positively with chromosome 8q copy number ( Figure 7B ).

GSEA revealed that 26 Hallmark pathways were significantly enriched in patients with high PRRX1 expression ( Table S5 ). The top 20 ranked Hallmark pathways included the angiogenesis (p = 0.010, FDR = 0.005), epithelial–mesenchymal transition (p <0.001, FDR = 0.002), and TGF-β signaling (p <0.001, FDR = 0.002) pathways, and also various immune- and inflammation-related sets of genes, namely, inflammatory response, IFN-γ response, IL6-Jak-Stat3 signaling, IL2-Stat5 signaling, TNF-α signaling via NF-kB, and IFN-α response ( Figure 7C ). These results provide evidence for the possible phenotypes that may involve PRRX1 expression.

The differences in expression of five immune checkpoint genes in the high and low PRRX1 groups are shown in Figure 8A . The expression levels of all five genes in the high PRRX1 expression group were significantly higher than that of the low expression group. TMB was found to decrease with increasing PRRX1 expression (R = −0.31, p = 0.0052) ( Figure 8B ). Subsequent evaluation of the association between PRRX1 expression and 10 TME-related cell populations from TCGA transcriptomic data using the Microenvironment Cell Populations-counter (MCPcounter) method (33, 34) showed that PRRX1 expression was positively associated with eight types of cells: CD8 T cells, cytotoxic lymphocytes, NK cells, cells of the monocytic lineage, myeloid dendritic cells, neutrophils, endothelial cells, and fibroblasts ( Figure 8C ).

The previously stated GSEA results revealed that PRRX1 expression significantly and positively correlated with the EMT signatures ( Figure 9A ). To investigate the biological function of PRRX1 in UM progression, we then performed loss-of-function experiments that silence PRRX1 in human UM cell line MuM-2B. The efficiency of the siRNA knockdown was confirmed by WB ( Figure 9B ), si-PRRX1-1 was chosen for subsequent experiments, which showed the most significant gene silencing efficiency. Compared to the control, downregulation of PRRX1 upregulates E-cadherin expression and downregulates Snail and N-cadherin expression in UM cells ( Figure 9C ). Furthermore, results of transwell invasion and scratch assays also demonstrated that the knockdown of PRRX1 significantly reduced UM cell migration and invasion ( Figures 9D, E ). Together, we initially confirmed PRRX1 promotes UM progression through affecting the EMT process.