The tissue samples were fixed in 10% formalin, embedded in paraffin, cut into 4 µm thick sections, stained with hematoxylin and eosin and examined under a light microscope. Fresh tumor and adjacent normal tissue samples were fixed in 10% neutral-buffered formalin at room temperature (RT) for 24 h. The paraffin-embedded tissue blocks were cut into 4 µm thick continuous sections. Hematoxylin-eosin staining was performed at RT (hematoxylin staining for 5 min and eosin staining for 2 min). Hematoxylin-eosin staining images were captured under a light microscope; Subsequently, the FFPE tissues were analyzed using IHC. All steps were performed on FFPE tumor tissue sections (4 µm thick) and adjacent normal renal tissue sections. Tissue sections were dewaxed in xylene, followed by graded ethanol hydration. Hydrated sections were incubated in sodium citrate buffer and heated at 95°C for 20 min to retrieve masked epitopes of target proteins (CD56, Syn, insulinoma associated protein 1, CgA and neuron specific enolase). Sections were treated with 0.5% hydrogen peroxide for 20 min to quench endogenous peroxidase activity and eliminate non-specific signal interference. Sections were permeabilized with 0.1% (v/v) Triton X-100 to enhance intracellular and membrane protein accessibility for antibodies targeting internal epitopes (CD56, Syn, insulinoma associated protein 1, CgA and neuron specific enolase). Permeabilized sections were incubated with 5% (w/v) bovine serum albumin (Sigma-Aldrich, Merck KGaA, Darmstadt, Germany; Cat. No. A9647) at RT for 60 min to block non-specific antibody binding sites. Sections were incubated with target-specific primary antibodies (Table SI) overnight at 4°C to allow specific antigen-antibody binding. Sections were incubated with biotin-labeled goat anti-rabbit/mouse IgG secondary antibodies (1:200; Catalogue No. BA-1000/BA-9200; Vector Laboratories, Burlingame, CA, USA) at RT for 30 min, per the manufacturer's standard protocol. Sections were incubated with HRP-labeled avidin (per manufacturer's protocol) and developed with DAB. Sections were stained with Mayer's hematoxylin at RT for 2 min to label cell nuclei. IHC-stained sections were visualized and images captured using an Olympus BX53 light microscope (Olympus Corporation, Tokyo, Japan).

Using a Covaris instrument, genomic DNA was fragmented into 180–280 bp segments for WES. The whole exome was captured according to the manufacturer's guidelines using an Agilent SureSelect Human All Exon Kit (Agilent Technologies, Inc.). Using an AMPure XP system (Beckman Coulter, Inc.), the products were purified and subsequently quantified using an Agilent high-sensitivity DNA assay on an Agilent Bioanalyzer 2100 system (Agilent Technologies, Inc.). DNA libraries, ~150 bp in insert size, were sequenced on an Illumina HiSeq platform (Illumina, Inc.) at a sequencing facility in Novogene Co., Ltd. Paired-end sequencing was performed using the Illumina HiSeq 4000 SBS kit (Catalog No. FC-404-1003). The final libraries were normalized to a loading concentration of 2 nM. The concentration was accurately measured using the Qubit^® dsDNA HS Assay Kit (Thermo Fisher Scientific) and confirmed by quantitative PCR (qPCR) using the KAPA Library Quantification Kit (Roche, Catalog No. KK4824) to ensure precise molarity calculation for cluster generation.

Raw sequencing data typically contained a few adapter reads, previously unexplained nucleotides and low-quality nucleotides. Ensuring the quality of the analysis required filtering raw reads to obtain clean reads, which served as the foundation for subsequent analysis. Paired reads that had adapter contamination, >10% uncertain bases or >50% low-quality bases were eliminated if they had a Phred quality score <5. The quality and integrity of the processed clean reads were verified using FastQC (v0.11.9, bioinformatics.babraham.ac.uk/projects/fastqc/), with key metrics including Phred quality scores (Q20/Q30), GC content distribution and sequence length uniformity validated; the average Q30 ratio of clean reads was >97% for both tumor and normal tissue, indicating high sequencing accuracy. Clean reads were aligned to the UCSC hg19 reference genome using BWA-MEM, and the aligned BAM files were further quality-checked using SAMtools and Picard to assess sequencing depth, mapping quality (MQ ≥30) and PCR duplication rate (<5%). Only clean, high-quality data that passed all the above quality verification steps were collected and used to enable meaningful downstream analysis.

Using the Burrows-Wheeler Aligner software (v0.7.17, http://bio-bwa.sourceforge.net/), valid sequencing reads were aligned to the reference genome (University of California Santa Cruz human genome version 19) to generate BAM files. Subsequently, Picard software (v2.26.10, Broad Institute, http://broadinstitute.github.io/picard/) was used to mark duplicate reads after sorting the BAM files with SAMtools (v1.15.1, http://www.htslib.org/). The final BAM files were used to calculate the sequence coverage and depth. Carcinogenic signaling pathways from The Cancer Genome Atlas (TCGA) database cohorts were collected and their enrichment in the sample was assessed.

To ensure meaningful analysis, variants [single nucleotide polymorphisms (SNPs) and INDELs] were identified and filtered using bcftools v1.23, SAMtools/HTSlib Team; htslib.org/doc/bcftools.html and mpileup integrated in SAMtools v1.23, htslib.org/doc/samtools-mpileup.html in SAMtools. Using the Genome Analysis Toolkit (GATK, version 4.6.2.0, gatk.broadinstitute.org/), all potential polymorphic SNP/INDEL mutation sites in the exons were identified. ANNOVAR software (v20210220, annovar.openbioinformatics.org/en/latest/) was used to describe and annotate variant call formats, with reference to the 1000 Genomes Project, Database of Single Nucleotide Polymorphisms (dbSNP) and other comparable databases. SNP/INDEL mutation data were imported into R software (version 4.5.0, R Foundation for Statistical Computing, Vienna, Austria) for mutational analysis. Oncoplot (waterfall) graphs were generated using the maftools package (version 2.24.0, bioconductor.org/packages/maftools/). Somatic single nucleotide variations (SNVs) were detected using MuTect (v2.2.7, gatk.broadinstitute.org/hc/en-us/articles/360037593851-MuTect2), and somatic INDELs were identified using Strelka (v2.9.10, Illumina, Inc.; http://github.com/Illumina/strelka). Furthermore, VarScan2 (version 2.4.6, github.com/dkoboldt/varscan) and Control-FREEC software (v11.6, github.com/BoevaLab/FREEC) were applied to identify variations in somatic copy number variations (CNVs).

A precipitation plot was used to visualize hypermutated genomic regions, known as kataegis loci, by mapping mutation frequency across chromosomes. Using the kataegis criterion (21), genomic sites with hypermutation were identified, characterized by an average of six consecutive mutations within 1,000 bp. Using the baseline site sequence as a starting point, six consecutive mutations were added and the average mutation distance was calculated. When the average interchange distance was >1,000, an element was appended to the end of the queue and one was removed from the front. If the average mutation distance was ≤1,000, the number of mutations was increased until it was >1,000. Subsequently, all mutations in the sequence were documented and created as a kataegis. The rainfall plot and kataegis cluster identification were performed using R software (version 4.2.2; R Development Core Team) with the GenVisR package (version 1.34.0; Bioconductor).

Predisposing genes can confer susceptibility to disease or increase vulnerability when exposed to specific environmental stimuli. SAMtools was used to identify germline mutations (SNVs and INDELs), and the findings were filtered using the Cancer Gene Census (CGC, v98, cancer.sanger.ac.uk/census) database to identify potential cancer-predisposing genes. The detailed methods in filtering mutation sites were as follows: i) Variant sites with a depth of <10× were filtered out; ii) the variant loci in the 1000 Genomes Project database (frequency >0.01 in the population) were filtered, the diversity loci among individuals were removed and rare mutations that were likely to cause disease were obtained. The mutation sites with frequencies <0.01 in the 1000 Genomes Project database and Catalogue Of Somatic Mutations In Cancer (COSMIC, v98; cancer.sanger.ac.uk/cosmic) database were retained; iii) the variations in the exonic region or the splicing region (2 bp upstream of the splicing site) were retained. Specifically, the variation sites located in the intron, non-coding and intergenic regions were filtered out; iv) synonymous mutations (mutations that did not cause changes in amino acid coding) were removed to obtain mutations that had an impact on the protein product (retaining mutations including frameshift and non-frameshift mutations in INDEL); and v) based on the scores of Sorting Intolerant From Tolerant (SIFT), Polymorphism Phenotyping v2-Human Variation (PolyPhen2_HVAR) and PolyPhen2_ Human Diversity (HDIV) for filtering, if a mutation was determined as a harmful mutation in at least one database or as a moderately harmful mutation in ≥2 or more databases, it was retained. The dbNSFP database was annotated for non-synonymous mutations, primarily assessing amino acid conservation and pathogenicity using corresponding computational scores. These scores included SIFT (v6.2.1, sift.bii.a-star.edu.sg/), PolyPhen-2 (v2.2.2, genetics.bwh.harvard.edu/pph2/) (including both HVAR and HDIV models), Functional Analysis Through Hidden Markov Models (v2.3, University of Bristol; http://fathmm.biocompute.org.uk/), MutationTaster (v2.0, Charité-Universitätsmedizin Berlin; http://www.mutationtaster.org/), Genomic Evolutionary Rate Profiling++ (v2.1, Stanford University; http://github.com/sidowlab/gerp++), MutationAssessor (r3, Memorial Sloan Kettering Cancer Center; http://mutationassessor.org/), Phylogenetic P-values (v1.5, Cold Spring Harbor Laboratory; http://compgen.cshl.edu/phast/), Likelihood Ratio Test (v1.5, Cold Spring Harbor Laboratory; http://compgen.cshl.edu/phast/phyloP-tutorial.php), and Site-specific Phylogenetic (v1.5, Cold Spring Harbor Laboratory; http://compgen.cshl.edu/phast/). Different scoring systems are based on different algorithms to evaluate the conservation and pathogenicity of mutations. Each score corresponds to the original score (score), the converted score (score_converted) and the predicted classification (pred). The present study primarily screened for mutations classified as predicted by SIFT (score ≤0.05) and PolyPhen2_HVAR (score ≥0.909). Scores were interpreted as follows: Variants with SIFT scores 0.00–0.05 were considered harmful and a PolyPhen score of 0.00 indicated a benign variant. Based on the American College of Medical Genetics (ACMG) scoring annotation results, the loci annotated as harmful, potentially harmful or ambiguous were retained. The final screening was conducted in accordance with ACMG guidelines (22). Furthermore, SAMtools was used to explore germline mutations, including SNVs and INDELs, in the normal tissues from patients. Potential predisposing genes based on germline mutations were identified and compared against CGC.

High-frequency mutant genes detected in a large number of samples, which may serve as candidate genes for disease occurrence and progression, were identified. Driver genes were those genes considered notable in the occurrence and development of diseases. Once driver genes were identified, personalized and precise medical treatment may be provided to patients in a targeted manner. SNP and INDEL data from all previous sample analyses were used to identify potential driver genes using the OncodriveCLUSTL algorithm (23). Furthermore, MutSigCV (24) was used to identify high-frequency mutant genes and these results were combined with analyses from the OncodriveCLUSTL algorithm to identify driver mutations. Additionally, a secondary verification was performed by amplifying the DNA fragment using PCR and the products were sequenced using Sanger sequencing to identify the mutant base. Genomic DNA extracted from the original FFPE tumor tissue using the QIAamp DNA FFPE Advanced Kit (Qiagen) was used as the template for PCR. Amplification was performed using the KAPA2G Fast HS PCR Kit (Kapa Biosystems, Inc.) according to the manufacturer's instructions. Prior to sequencing, PCR products were electrophoresed on a 2.0% (w/v) agarose gel containing GelRed Nucleic Acid Gel Stain (Biotium, Inc.). Target amplicons were visualized using a Gel Doc XR + Imaging System (Bio-Rad Laboratories, Inc.), excised, and purified using the QIAquick Gel Extraction Kit (Qiagen). Purified products were then subjected to Sanger sequencing. Primer sequences and detailed PCR conditions are presented in Tables SII and SIII.

Functional enrichment analyses of the identified predisposing and driver genes were performed to annotate their biological functions and associated signaling pathways using the Gene Ontology (GO) database (geneontology.org/), Kyoto Encyclopedia of Genes and Genomes (KEGG) database (genome.jp/kegg/), and Eukaryotic Orthologous Groups (KOG) database (ncbi.nlm.nih.gov/COG/KOG/). For GO analysis, the genes were categorized into three main ontologies: biological process (BP), cellular component (CC), and molecular function (MF); significantly enriched GO terms were screened based on the biological relevance to PRSCC pathogenesis. For KEGG analysis, gene enrichment in canonical signaling pathways and biological pathways was identified to reveal the molecular mechanisms underlying PRSCC development associated with the candidate genes. For KOG analysis, the genes were classified into orthologous functional categories to further annotate their conserved cellular and molecular functions across eukaryotes.

Literature published before October 2025 was searched for studies on patients with PRSCC in PubMed (https://pubmed.ncbi.nlm.nih.gov/) and EMBASE (https://www.elsevier.com/products/embase) databases using the keywords ‘renal’, ‘small cell carcinoma’ and ‘kidney’ with the appropriate Boolean modifiers. The analysis included only patients with complete follow-up information. As detailed in Table SIV, 60 patients were included from the literature. Their clinicopathological features, such as age, sex, symptoms, tumor size, laterality, primary site, pathological tumor (T) stage, lymph node metastasis, distant metastasis, surgery, chemotherapy, follow-up time and survival status, were collected.

Univariate Cox regression analyses were performed to identify variables significantly associated with overall survival (OS). Variables with a P<0.05 from the univariate analysis were subsequently included in the multivariate analysis. A multivariate analysis was performed using a Cox proportional-hazards model to identify independent prognostic factors for PRSCC. According to the univariate or multivariate analyses, hazard ratios (HR) and their corresponding 95% confidence intervals (CI) were calculated. P<0.05 was considered to indicate a statistically significant difference. All statistical analyses were performed using R software (version 4.2.2; R Development Core Team).

The histopathological manifestation of PRSCC was confirmed by protein expression profiling using IHC. IHC analysis revealed that CD56 and Syn were diffusely and strongly positive, whereas INSM1, CgA and NSE were focally and weakly positive, aligning with the tumor characteristics of PRSCC (Fig. 1). Proliferative activity was high, as indicated by a Ki-67 proliferation index of ~80% (Fig. 1). Furthermore, moderate expression level of CK and CK8/18 suggested that the tissue originated from epithelial cells (Fig. 1). Negative staining of S100, WT1, CD45, CD99 and p63 helped differentiate PRSCC from paraganglioma, Wilms' tumor, malignant lymphoma, Ewing sarcoma/peripheral primitive neuroectodermal tumor and high-grade urothelial carcinoma, respectively (Fig. 2); negative staining of TTF1 helped differentiate PRSCC from metastatic small cell carcinoma; and negative staining of CK20, PAX-2, PAX-8 and vimentin helped differentiate PRSCC from urothelial and renal cell carcinoma (including papillary renal cell carcinoma and ccRCC) (Fig. 2).

Typically, the genome of an individual includes ~3.6 million SNPs. High-frequency SNPs are predominantly recorded in the dbSNP database. The distribution of SNPs and INDELs across different genomes and coding regions is presented in Table II. A total of 83,934 SNPs were identified in PRSCC specimens and 28,228 SNPs in adjacent normal specimens, with most being identified in coding sequence (CDS) areas. Deletions occurring at splicing sites or coding regions may impact protein translation. Additionally, the genome of each individual generally harbors ~350,000 INDELs. A total of 2,489 INDELs were identified in PRSCC specimens and 2,206 INDELs in adjacent normal specimens, with majority of INDELs identified in CDS areas.

The mutation sites were filtered and screened based on repeatability, depth and quality value to obtain a high-confidence mutation dataset annotated using ANNOVAR. A total of 63,677 missense mutations and 3,294 nonsense mutations were identified. Furthermore, somatic mutations were detected in genes not previously known to be associated with PRSCC, and the following were the top 10 types of mutant genes: TTN, ZNF717, MUC3A, MUC16, MUC6, OBSCN, CDC27, MUC4, HYDIN and ANKRD36C (Fig. 4). Subsequently, the mutation analysis is presented in the form of oncoplot (waterfall) graphs (Fig. 5). An SNV is a variation that arises from the replacement of a single nucleotide in the genome. ‘MuTect’ was primarily used to identify somatic SNV mutation sites. A total of 113 somatic SNVs were identified, the majority of which were located in the intronic and CDS regions. As presented in Table III, the somatic SNVs were distributed across different genomic regions. Additionally, using ‘Strelka’, a total of 26 somatic INDELs were identified. Similarly, the majority of the somatic INDEL sites were located in the intronic and CDS regions. Table III presents the somatic INDEL detection results across the various genomic regions.

CNVs, defined as changes in the number of copies of genomic segments, are the primary source of SVs and are classified into duplications and deletions. The gain and loss counts of CNVs were 76 and 193 in tumor specimens, respectively. There were 356 frameshift deletions (95.7%), 8 non-frameshift deletions (2.2%), and 8 unknown (2.2%) mutation types (Fig. 6A). Regarding the region of mutation, 365 were exonic (98.1%), 5 ncRNA exonic (1.3%), 1 upstream (0.3%) and 1 downstream (0.3%) (Fig. 6B). By contrast, the gain counts and loss counts of CNVs were 65 and 186 in adjacent normal specimens, respectively. There were 253 frameshift deletions (98.1%), 2 non-frameshift deletions (0.8%) and 3 unknown (1.2%) mutations (Fig. 6A). Furthermore, 257 were exonic (99.6%) and 1 ncRNA_exonic (0.4%) (Fig. 6B).

Genomic SV is a general term for sequence variants involving >50 base pairs. Delly was used for SV detection and to count the various SV types. Delly can identify five types of SV mutations: Deletion (DEL), insertion (INS), duplication (DUP), inversion (INV) and translocation (BND). For tumor specimens, 5,894 unknown (87.6%), 595 start-loss (8.8%), 141 frameshift deletions (2.0%), 68 non-frameshift deletions (1.0%), 24 stop-loss (0.4%) and 7 stop-gain (0.1%) functional types were identified (Fig. 7A). Regarding the distribution of the mutations by region, 3,727 were exonic (58.9%), 1,659 intronic (26.2%) and 945 intergenic (14.9%) mutations (Fig. 7B). Furthermore, the present study identified 2,832 BND (42.1%), 1,936 INV (28.8%), 1,063 DEL (15.8%) and 898 DUP (13.3%) (Fig. 7C). For adjacent normal specimens, a total of 3,831 mutations (91.5%) were classed as unknown, 329 start-loss (7.9%), 15 frameshift deletions (0.4%), 8 non-frameshift deletions (0.2%), 3 stop-loss (0.07%) and 3 stop-gain (0.07%) (Fig. 7A). Regarding the regions of mutations, 1,752 were exonic (44.6%), 1,362 intronic (34.7%) and 812 intergenic (20.7%) (Fig. 7B). Furthermore, 2,692 BND (64.3%), 763 INV (18.2%), 386 DEL (9.2%) and 348 DUP (8.3%) mutations were identified (Fig. 7C).

The correctness of the SNP dataset was evaluated using the transition/transversion (Ts/Tv) ratio. Based on permutation and combination principles, SNPs can undergo six types of substitutions: C ↔ G, C ↔ T, A ↔ G, A ↔ T, A ↔ C, and G ↔ T. Structural reasons make the probability of transition higher than that of transversion. The Ts/Tv ratio across the whole gene was ~2.2, whereas it was ~3.2 in the coding region; in the present study, 83,890 C>T substitutions were identified (~75.2%). The distribution of chromosome-based variants was displayed as a Circos plot to depict the somatic cell variation in the PRSCC sample (Fig. 8).

A total of 75 genes were associated with the ‘receptor tyrosine kinase (RTK)-Ras’ pathway, followed by 60 genes for ‘NOTCH’ and 54 genes for ‘Wnt’. Additional carcinogenic signaling pathways are listed in Table SV.

Clusters of mutations were identified by examining entire cancer genome mutation catalogs. Kataegis foci, characterized by localized substitution hypermutation, are defined by clusters of C>T and/or C>G mutations significantly enriched at TpCpN trinucleotides and on the same DNA strand. Kataegis foci range from a few to thousands of mutations and are generally located near chromosomal rearrangements. Different tumor types influence various genomic regions. Starting from the reference sequence, the six base substitutions were analyzed in classes and the average intermutation distance was determined. The rainfall plots for the tumor tissues are depicted in Fig. 9.

Predisposing genes can confer susceptibility to disease or increase vulnerability when exposed to specific environmental stimuli. SAMtools was used to identify germline mutations (SNVs and INDELs) and the results were filtered against the CGC database to identify potential cancer-predisposing genes. The mutation annotation sites were screened after retaining the scoring factors for harmfulness and uncertainty of significance, and applying the filtering conditions for susceptible genes. A total of 350 genes were identified with SIFT scores between 0–0.05. Additionally, the SIFT score was set to 0 and the PolyPhen2_HVAR score to 1, and the genes that deviated the most from the relative score were used to identify genes most confidently predicted as deleterious, yielding eight genes (DST, OR10H3, PTK2B, APOBR, ZNF606, CCN4, ADCK1 and MYH2). Genes with scores that were predicted to be deleterious were primarily screened based on mutation classifications predicted by SIFT (score ≤0.05) and PolyPhen2_HVAR (score ≥0.909) (Fig. 10A). The majority of functionally enriched pathways for predisposing genes in PRSCC were involved in ‘developmental processes' and ‘cellular component organization or ‘biogenesis’ (Fig. 10B). According to the KEGG pathway analysis, PRSCC primarily indicated enrichment in ‘membrane transport’, ‘signal transduction’ and ‘cellular community’ (Fig. 10C). Enrichment analysis using KOG pathway analysis indicated enrichment of ‘signal transduction mechanisms’, ‘transcription’ and ‘general function prediction only’ (Fig. 10D).

Previously established driver genes that exhibited somatic variances in the tumor sample were filtered. Subsequently, potential driver genes were identified using SNP and INDEL information from the previous analysis of all samples. Mutations in 10 driver genes, KRTAP10-9, HYDIN, ZNF665, KRTAP10-2, GPAM, MUC12, KRT9, CCDC168, DUSP27 and MDC1, were identified (Table SVI; Fig. 11A). The present study indicated that the majority of functionally enriched pathways for driver genes in PRSCC were ‘multicellular organismal process’ and ‘cellular component organization or ‘biogenesis’ (Fig. 11B). Based on KEGG pathway analysis, PRSCC primarily indicated enrichment in ‘signal transduction’ (Fig. 11C). Enrichment analysis using KOG demonstrated ‘signal transduction mechanisms’ as the primarily enriched pathways (Fig. 11D). Secondary verification was performed using Sanger sequencing to identify the mutant base. Sanger sequencing of germline DNA identified a HYDIN A/G variant (Chr16:70908736), the first such report of this mutation in PRSCC, to the best of our knowledge (Fig. 12).

According to the univariate analyses (Table SVIII), tumor size, T and M stages were all significantly associated with a worse OS (all P<0.05). Furthermore, Kaplan-Meier curves revealed that the patients with T3/T4 tumors had a worse OS compared with those patients with T1/T2 tumors (P<0.0001; Fig. S3). Multivariate analysis using a Cox proportional-hazards model was used to identify independent prognostic factors for PRSCC. Variables with a P-value <0.05 from the univariate analysis were included in the multivariate analysis. The multivariate analysis of OS demonstrated that patients with T3/T4-stage tumors had worse survival outcomes compared with those patients with T1/T2 tumors (HR=6.266; 95% CI, 1.779–22.080; P=0.004). The multivariate results demonstrated that T stage was an independent prognostic factor for OS (Table SVIII).