This experimental study was approved by the Ethics Committee of our hospital, and three cases of patients with limited prostate adenocarcinoma were collected from the Department of Urology of our hospital in the year of 2022, under the guidance of the physicians of the Department of Pathology. The PRAD tumour tissue specimens were used as the tumour group, and the paracancerous normal prostate tissue specimens were used as the normal group. Tissue specimens were collected and stored in liquid nitrogen to send for eccDNA sequencing.

PRAD and paracancerous normal prostate tissue specimens were subjected to eccDNA sequencing assisted by CloudSeq Biotech Inc. (Shanghai, China) using the circle-seq (Møller et al., 2018) method. Briefly, cell deposits were resuspended in L1 buffer (Plasmid Mini AX; A&A Biotechnology) supplemented with protease K (ThermoFisher) prior to digestion at 50°C overnight. Digested samples were alkali-treated and column-purified by following the instructions of the Plasmid Mini AX kit. Column-purified DNA samples were digested by FastDigest MssI (ThermoFisher) at 37°C for 16 h to remove mitochondrial circular DNA. Then, the samples were incubated with Plasmid-Safe ATP-dependent DNase (Epicentre) at 37°Cfor 1 week to remove the remaining linear DNA. The samples were then supplemented with 30U of DNase and a proportional amount of ATP every 24 h. The treated samples were used as templates for eccDNA amplification by using the RCA DNA Amplification Kit (GenSeq Inc.), followed by purification with the MinElute Reaction Cleanup Kit (Qiagen). Library preparation of purified DNA was performed with the GenSeq^® Rapid DNA Lib Prep Kit (GenSeq Inc.). High-throughput sequencing was performed on an Illumina NovaSeq 6000 sequencer in 150 bp double-ended mode to obtain the raw data. Quality control was performed with Q30 as following sequence, low-quality reads were removed firstly by using cutadapt software (v1.9.1), and high-quality clean reads were aligned to the reference genome by using bwa software (v0.7.12). Next, all eccDNAs were identified with circle-map software (v1.1.4) and then raw soft-clipped read counts of the break point were obtained by using SAMtools (v1.9) software. Normalisation and differential analysis were performed by using DESeq2 [7] (v1.38.3) software. Annotation of eccDNA was performed by using bedtools software (v2.27.1) and enrichment analyses were performed by using the eDEGs. eccDNA visualisation was performed by using IGV (v2.4.10) software.

Transcriptional profiles, clinical features, tumour mutation burden (TMB) and microsatellite instability (MSI) scores of PRAD were downloaded from The Cancer Genome Atlas Program (TCGA, https://portal.gdc.cancer.gov/) database. Validation set data were obtained from cBioPortal-SU2C/PCF (https://www.cbioportal.org/) and GEO70770 (https://www.ncbi.nlm.nih.gov/geo/). Data preprocessing and DEGs analysis were perfoemed by the “limma” and “affay” packages in the R environment. The coding genes amplified by differential eccDNA were taken to intersect with the DEGs of TCGA-PRADt to obtain eccDNA-amplified differentially expressed coding genes (eDEGs). Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) functional enrichment analyses were performed using the “clusterProfiler” package. GO enrichment analysis described the potential functions of genes in terms of Molecular Function (MF), Cellular Component (CC) and Biological Process (BP). KEGG analysed the major metabolic and signal transduction pathways in which the genes were involved through pathway annotation. Cox regression analysis and Lasso analysis were used to further identify prognostic genes and construct risk models. Cox regression analysis was performed using the “survival” package and Lasso analysis was performed using the “glmnet” package. Kaplan-Meier (KM) survival analysis plots and risk factor association plots based on the risk models were constructed using the “survival” and “ggplot2” packages. Nomograms were used to visualise the 1-year, 3-year, and 5-year survival predictions of the risk model. Time-dependent ROC curves were used to verify the model accuracy. Gene Set Enrichment Analysis (GSEA) was performed using the “Cluster Profiler,” “org. Hs.eg.db”, and “enrichplot” package to identify biological processes and enrichment pathways of the key gene. Quantification of the immune microenvironment was performed using XCELL, MCPCOUNTER, CIBERSORT, TIMER, EPIC, and QUANTISEQ algorithms. Identify the characteristics and differences in immune infiltration (performed with the “CIBERSORT” and “reshape2” packages), immune function (performed with the “RColorBrewer” package), and immune subtypes (performed with the “RColorBrewer” package) of risk models based on quantitative immune microenvironment results. ESTIMATE analysis identified specific signals associated with stromal and immune cell infiltration in tumour tissue and predicted the level of infiltrating stromal and immune cells by calculating stromal and immune scores, which were performed with the “utils” package. Gene mutation frequency and mutation burden of the risk models were analyzed with the“maftools” package. Immunophenotype score (IPS) was obtained from The Cancer Immunome Database (TCIA) (Charoentong et al., 2017), which was used to show the response of PRAD patients to immunotherapy. Also, Tumor Immune Dysfunction and Exclusion (TIDE) (Jiang et al., 2018) algorithm was used to assess patients’ immunotherapy response. Data for drug sensitivity analysis were obtained from Genomics of Drug Sensitivity in Cancer (GDSC) (Yang et al., 2013) by using the “oncoPredict” package. All the above data visualisations relied on R language implementation.

The eccDNA expression profiles in PRAD tumour tissues and paracancerous normal prostate tissues were obtained by eccDNA sequencing. The results showed that the tumour group contained 76,636 eccDNAsand 11,967 eGenes (Figure 1A; Supplementary Table 1), and 58,781 eccDNAs and 10,694 eGenes were identified in the normal group (Figure 1A; Supplementary Table 2). The eccDNA can encode some or all exons of a gene to affect protein expression, and different eccDNA can encode the same gene. In the normal group, we found 9,460 eGenes derived from 1-5 eccDNAs, 885 eGenes derived from 6–10 eccDNAs, 216 eGenes derived from 11–15 eccDNAs, 74 eGenes derived from 16–20 eccDNAs, 33 eGenes derived from 21–25 eccDNA derivatives, 17 eGenes were derived from 26–30 eccDNA types, 5 eGenes were derived from 31–35 eccDNA types, 3 eGenes were derived from 36–40 eccDNA types. The number of eccDNA types amplifying the CNTNAP2 gene are even more than 45, reaching up to 47 (Figure 1B; Supplementary Table 3). In the PRAD tumour group, 10,174 eGenes were detected to be derived from 1-5 eccDNA types, 1,151 eGenes were derived from 6–10 eccDNA types, 359 eGenes were derived from 11–15 eccDNA types, 141 eGenes were derived from 16–20 eccDNA types, 65 eGenes were derived from 21–25 eccDNAs, 37 eGenes were derived from 26–30 eccDNAs, 17 eGenes were derived from 31–35 eccDNAs, 8 eGenes were derived from 36–40 eccDNAs, 10 eGenes were derived from 40–45 eccDNAs, and five eGenes were derived from more than 45 eccDNA types, namely CNTNAP2, TRAPPC9, DAB1, RBFOX1 and CAMTA1 (Figure 1C; Supplementary Table 4).

At the same time, multiple coding genes could be amplified from the same eccDNA. In the normal group, there were 27,150 eccDNAs amplified 1 eGene, 1,983 eccDNAs amplified 2 eGenes, 99 eccDNAs amplified 3 eGenes, 12 eccDNAs amplified 4 eGenes, 2 eccDNAs amplified 5 eGenes, and 4 eccDNAs amplified 6–10 eGenes, 7 eccDNAs amplified 11–15 eGenes, 1 eccDNA amplified 16–20 eGenes, and 4 eccDNAs amplified more than 20 eGenes (Figure 1D; Supplementary Table 5). In the PRAD tumour group, there were 35,812 eccDNAs amplified 1 eGene, 1,958 eccDNAs amplified 2 eGenes, 156 eccDNAs amplified 3 eGenes, 17 eccDNAs amplified 4 eGenes, 1 eccDNA amplified 5 eGenes, and 8 eccDNAs amplified 6 to 10 eGenes, 7 eccDNAs amplified 11–15 eGenes, 4 eccDNAs amplified 16–20 eGenes, and 2 eccDNAs amplified more than 20 eGenes (Figure 1E; Supplementary Table 5).

In addition, we found that the size distribution of eccDNA in the normal group ranged from 10 bp to 14,000 kb (Figure 1F), and that in the tumour group ranged from 10 to 6,000 bp (Figure 1G). Both groups had an emergent peak at 300 bp (Supplementary Figures S1A–C), and there was no significant difference in the size distribution of eccDNA between the tumour and normal groups (Supplementary Figure S1D). The GC content enrichment of eccDNA in normal and tumour tissue were both significantly higher than other genomic regions (Figures 1H, I).

We further analyzed the genomic distribution of eccDNA on different chromosomes, including intact eccDNA (Figures 2A, B), eccDNA amplifying coding genes (Figures 2C, D), and eccDNA with unamplifying coding genes (Figures 2E, F). The results show that gene-rich chromosome contributed to a higher average frequency of eccDNAs per Mb than other chromosomes, such as chromosome 1, while gene-poor chromosome contributed to a lower average frequency of eccDNAs per Mb, such as chromosome Y. It suggests that regions with gene-rich are more preferentially producing eccDNA. eccDNA distribution on chromosomes between the normal group and PRAD tumour group was not significant differentiation (Figure 2G).

Finally, we explored the possible origins of eccDNAs by mapping the eccDNAs to different genomic elements (Figure 2H) and repetitive elements (Figure 2I). Notably, eccDNA was more significantly enriched in both 5′ UTR genomic region and repetitive elements, such as long interspersed elements (LINEs) and short interspersed elements (SINE), suggesting that these regions are more preferentially producing eccDNA in PRAD.

Based on the eccDNA sequencing results, a total of 4,290 differentially expressed eccDNA were screened in PRAD tissues compared with normal tissues (Figure 3A; Supplementary Tables 6, 7). Among them, 1,667 eccDNAs were higher expressed in the tumour tissues, and these eccDNAs amplified 798 eGenes. 2,623 eccDNAs were lowly expressed and amplified 1,183 eGenes (|FC(fold change)| ≥ 2, P < 0.05) (Supplementary Figures S1E, F). The transcriptome data of PRAD and normal samples were obtained from the TCGA database, and 5,960 DEGs were obtained by screening (Supplementary Figure S1G). The coding genes amplified by differential eccDNA and the DEGs of TCGA-PRAD were taken to be intersected (Figure 3B), and 499 eDEGs were obtained. Further analysis of the distribution of eDEGs on chromosomes (Figure 3C) showed that eDEGs were enriched on chromosomes 1 to 22 without appear on the Y chromosome. Functional enrichment analysis of eDEGs (Figure 3D) showed that their roles were mainly focused on post-translational modification, signal transduction, and cell intercellular communication pathways.

To further explore the eKDEGs of PRAD, this study combined the TCGA-PRAD transcriptome dataset and further performed one-way Cox regression analyses (Figure 3E), LASSO analyses (Figure 3F), and multifactorial Cox regression analyses (Figure 3G) on the eDEGs. ZNF330 and PITPNM3 were finally identified as eKDEGs and independent risk factors for PRAD. ZNF330 was highly expressed and PITPNM3 was lowly expressed in PRAD (Supplementary Figure S1H), which corresponded to eccDNA sources of ZNF330^{circle142141735-142142329}, PITPNM3^{circle6458635-6459156}. In addition, the differential expression was validated in the GSE70770 dataset set (Supplementary Figure S1I) with consistent results.

Survival analysis (Figure 3H) showed that PRAD patients with high expression of ZNF330, PITPNM3 had lower overall survival (OS), which suggested that ZNF330, PITPNM3 may be a risk factor for poor prognosis. A novel PRAD risk model was constructed by basing on these two genes, and it was found that the risk score of patients increased with higher expression of ZNF330 and PITPNM3 (Figure 3I). Patients were classified into high and low risk groups based on the risk scores, and the higher risk group had lower OS (Figure 3J), which may lead to poor prognosis. Nomograms were plotted to predict patient survival based on risk scores (Figure 3K), and the 1-year, 3-year and 5-year survival rates of patients gradually decreased with increasing risk scores. The area under the curve (AUC) of the time-dependent ROC at 1-year, 3-year and 5- years were 1.000, 0.912, and 0.900 respectively (Figure 3L), suggesting that the risk model had good predictive performance. cBioPortal-SU2C/PCF data set validated the model, and survival analysis still showed that the high-risk group was associated with poor prognosis (Figure 3M). GSEA analysis (Figure 3N) showed that the high-risk group was significantly enriched in the cytokine signalling pathway (KEGG_CYTOKINE_CYTOKINE_RECEPTOR_INTERACTION, KEGG_HEMATOPOIETIC_CELL_LINEAGE) and cancer-related pathways (KEGG_JAK_STAT_SIGNALING_PATHWAY). Meanwhile, clinical correlation analysis (Figure 3O) showed that high risk scores were positively correlated with later clinical T stage and higher Gleason scores, suggesting that PRAD in the high-risk group were more malignant. KEGG_HEMATOPOIETIC_CELL_LINEAGE) and cancer-related pathways (KEGG_JAK_STAT_SIGNALING_PATHWAY). Meanwhile, clinical correlation analysis (Figure 3O) showed that high risk scores were positively correlated with later clinical T stage and higher Gleason scores, suggesting that PRAD in the high-risk group were more malignant.

Quantitative analysis of the immune microenvironment was performed on the PRAD risk model based on multiple immunological algorithms. Different immune infiltration patterns were observed in patients with the high and low-risk groups. The immune infiltration analysis (Figure 4A) showed that the immune microenvironment in the high-risk group had increased levels of T cells CD4 memory resting, Macrophages M0, Macrophages M2 and Tregs, and decreased levels of T cells follicular helper and NK cells activated. Meanwhile, the infiltration abundance of some immune cells correlated with prognosis (B cells naive and Tregs infiltration were associated with poor prognosis, and Macrophages M1 and M2 infiltration were associated with better prognosis) (Supplementary Figure S2A). Single gene immune infiltration analysis of ZNF330 and PITPNM3 (Supplementary Figures S2B, C) also revealed multiple immune cell content changes. Interestingly, differential analysis of immune cell function showed that the majority of immune cells were functionally active in the high-risk group (Figure 4B), and altered immune function was associated with patient prognosis (DCs functionally active was associated with a poorer prognosis, and APC_co_inhibition, Mast_cells and Tfh functionally active were associated with a better prognosis) (Supplementary Figure S2D). In addition, there were differences of immune subtypes distribution in the risk model (Figure 4C). ESTIMATE analysis found a negative correlation between the risk model and stromal score, immune score and estimated scores (Figure 4D), suggesting that tumour purity was higher in the high-risk group. IPS scores were calculated to predict response of PRAD patients to two ICIs, anti-CTLA-4 and anti-PD-1 (Figure 4E). We found a positive correlation between the high-risk group andips-CTLA4 (−)/PD1 (−), while a negative correlation was existed between the high-risk group and ips-CTLA4 (−)/PD1 (+) and ips-CTLA4 (+)/PD1 (+) in, indicating that the PRAD model risk can influence the immunotherapy response and the high-risk group of PRAD patients had poor responses to ICIs. These results suggested that tumour immunity may play a key role in PRAD.

TMB and MSI are molecular markers for determining the suitability of immunotherapy for tumour patients, which also suggest genomic instability. In the risk model, we observed that the mutation rate was lower in the high-risk group (Figure 4F), and the same five most mutated genes in two risk groups were SPOP, TTN, TP53, KMT2D, and FOXA1. What’s more, the TMB and MSI scores were negatively correlated with the risk scores (TMB: R = −0.15; MSI: R = −0.2) (Figure 4G). In addition, Correlation analysis of mutations in the risk model gene set and Gleason scores (Figure 4H) showed that advanced Gleason scores were higher ratioin the mutation group in PRAD patients, which may be associated with a poor prognosis.

Subsequently, the TIDE score and immune exclusion was calculated for each PRAD patient based on the TIDE analysis (Figure 5A), and the risk scores were positively correlated with them. It is suggested that patients in the high-risk group are more likely to experience immune escape and poor immunotherapy. Based on the previous analysis, we found that both immune infiltration and mutation have important roles in the risk model. Therefore, we further explored the correlation between mutation and immune phenotype through the mutation profile of the gene set in the risk model. Both gene amplification (Amp) and deletion (Dele) were types of copy number variation (CNV) (Zhang et al., 2009), which were belongs to mutations. The results showed that the infiltration score was higher in the gene amplification and deletion groups compared with the wild type (WT) (Figure 5B), and the deletion group had lower exhausted abundance, whereas there was no significant difference in exhausted abundance between the wild type and amplification group (Figure 5C). The abundance of iTreg (Figure 5D) in the amplification group, and B cell (Figure 5E) and DC cell (Figure 5F) in the deletion group were higher compared to the wild type.

Drug sensitivity analysis based on the PRAD risk model gene set (Figure 5G) revealed that the high- Drug sensitivity analysis based on the PRAD risk model gene set (Figure 5G) revealed that the high-risk group decrease the sensibility of Ipatasertib, AZD5363, Oxaliplatin, MK-2206, AZD8055, and AZD8186in PRAD, and increase the sensibility of Sinularin, Cyclophosphamide, AZD4547, and Osimertinib in PRAD.