To identify inherited genetic variants that may predispose individuals to PMMC, germline variants were called from the WES data of the three adjacent non-cancerous tissue samples. This analysis aimed to pinpoint rare, potentially pathogenic variants in known cancer predisposing genes (CPGs).

Variant Calling: Germline single-nucleotide polymorphisms (SNPs) and short InDels were identified using the SAMtools mpileup and BCFtools call pipeline. For each sample’s analysis-ready BAM file, genotype likelihoods were calculated at each position using mpileup. These likelihoods were then piped to bcftools call with the -mv options to perform variant calling using the multiallelic-aware model and to output only sites containing a variant allele.

Variant Filtering for Susceptibility Genes: The resulting raw variant calls were subjected to a rigorous, multi-stage filtering cascade designed to remove common polymorphisms and technical artifacts, thereby enriching for rare variants with a higher likelihood of being pathogenic:

1. Coverage Filter: Variants at sites with a read depth below 10× were excluded, as calls with insufficient read support are unreliable.

2. Population Frequency Filter: To isolate rare variants, any SNP or InDel with a minor allele frequency (MAF) greater than 0.001 (0.1%) in large-scale public population databases was removed. The databases used for this filtering step were the 1000 Genomes Project, the Exome Aggregation Consortium (ExAC) (22), and the NHLBI GO Exome Sequencing Project (ESP6500) (23).

3. Database Annotation Filter: Variants cataloged in the dbSNP database were filtered out to remove known common variants. However, an exception was made for variants that, despite being in dbSNP, are also annotated in the COSMIC database, as this suggests a potential role in somatic carcinogenesis.

4. Functional Consequence Filter: To focus the analysis on variants most likely to have a functional impact, variants located in intergenic regions, introns, or other non-coding regions were excluded. Additionally, synonymous mutations within coding sequences that do not alter the resulting amino acid were removed.

5. Genomic Context Filter: Variants situated within known repetitive elements or low-complexity regions of the genome were discarded, as these regions are notoriously prone to read misalignment and can lead to false-positive variant calls.

6. In Silico Pathogenicity Prediction Filter: The remaining variants were annotated with predictions of their functional impact using a suite of in silico tools. A variant was retained for final consideration only if it met stringent criteria for predicted pathogenicity: classified as “damaging” or “harmful” by at least one of the following algorithms—SIFT, PolyPhen-2 HumVar, PolyPhen-2 HumDiv, or MutationTaster—or classified as “possibly damaging” by two or more of these tools.

Somatic mutations—genetic alterations unique to the tumor cells—were identified by comparing the WES data from each tumor sample against its matched adjacent non-cancerous tissue sample. This paired analysis is essential for distinguishing true somatic events from germline variants. For enhanced accuracy and sensitivity, somatic single-nucleotide variants (SNVs) and InDels were called using the GATK4 MuTect2 (v4.2.0.0) (24) pipeline, following the GATK Best Practices for somatic short variant discovery.

For each of the three tumor-normal pairs, MuTect2 was run in paired mode, providing both the tumor and normal BAM files as input. This allows the algorithm to construct a model of the normal sample’s germline variation and identify variants present only in the tumor with high confidence. The raw VCF output from MuTect2 was then processed with GATK’s FilterMutectCalls tool, which applies a series of learned filters to flag and remove common sequencing artifacts and low-confidence calls. Only variants that passed all filters (i.e., had “PASS” in the FILTER column of the VCF file) were retained. Finally, these high-confidence somatic variants were annotated to determine their location and to predict their effect on protein function, and only those occurring within coding regions were carried forward for downstream analyses.

Somatic CNVs, representing gains (amplifications) and losses (deletions) of genomic segments in the tumor cells, were detected from the WES data using Control-FREEC (v11.6) (32). This tool analyzes read counts in predefined genomic windows to infer copy-number status. For each tumor sample, the matched normal sample was used as a control to establish a baseline diploid copy-number, enabling the specific detection of somatic events.

The software was configured with parameters optimized for targeted sequencing data, like WES, which is characterized by non-uniform coverage. Key parameters included: degree = 1 (for control-read-count-based normalization), forceGCcontentNormalization = 1 (to correct for GC-content bias), minCNAlength = 3 (requiring a CNA to span at least three consecutive windows), minimalSubclonePresence = 30 (to detect subclones present in at least 30% of cells), noisyData = TRUE (to account for non-uniform capture efficiency), and minMappabilityPerWindow = 0.85 (to exclude windows with low mappability).

To place the genomic alterations found in PMMC into a broader context of skin malignancies, the mutational landscape of PMMC was compared with that of SKCM, a well-characterized cancer type.

Data Acquisition: Publicly available somatic mutation data for SKCM were acquired from two major cancer genomics databases. First, data for 1,366 SKCM cases were extracted from the COSMIC database (v94) (https://cancer.sanger.ac.uk/cosmic) using the search terms “skin” for the primary site and “malignant melanoma” for the disease type. Secondly, Masked Somatic Mutation (SNV) data for 467 patients with SKCM and 289 patients with cervical cancer (CC) were downloaded from The Cancer Genome Atlas (TCGA) portal (https://portal.gdc.cancer.gov/).

Intersection Analysis: A list of gene mutation frequency was generated from the TCGA SKCM cohort. To identify somatic mutations potentially unique to PMMC, the list of SMGs identified in this study was intersected with two lists derived from the SKCM data (1): the top 300 most frequently mutated genes in the TCGA cohort, and (2) genes mutated in 10 or more samples in the COSMIC SKCM cohort. Genes present in the PMMC SMG list but absent from both SKCM lists were highlighted as candidate genes specific to the pathobiology of PMMC.

CNV Comparison: For a comparative analysis of large-scale genomic alterations, Masked Copy Number Segment data for 472 TCGA SKCM tumor samples were also downloaded from the TCGA portal (data accessed on May 31, 2021) to compare CNV profiles between PMMC and SKCM.

HMV-II human melanoma cells were purchased from Sigma-Aldrich (92042701, USA). HMV-II cells were removed from liquid nitrogen and quickly placed in a 37 °C water bath. The cryovial was gently shaken to thaw. After thawing, the cells were transferred to a centrifuge tube containing 5 mL of culture medium and centrifuged at 1000 rpm for 5 min at room temperature. The supernatant was discarded, and the cells were resuspended in RPMI 1640 Medium (Gibco, USA) supplemented with fetal bovine serum, plated onto a culture dish, and gently pipetted to mix. The cells were incubated at 37 °C in a humidified atmosphere of 5% CO2. Cells were passaged when the cell density reached 80% and carefully handled to prevent contamination and authenticated by short tandem repeat sequence analysis to ensure quality.

HMV-II cells were seeded at a density of 5 × 10⁵ cells per well in a 6-well plate. Next, 15 nM siRNA-1 or siRNA-2 and siRNA NC (GenePharma, Shanghai, China) were transfected into HMV-II cells using Lipofectamine 3000 (Thermo, New York, USA). The sequence of siRNA was: P53-siRNA-1: AGACCUAUGGAAACUACUU; P53-siRNA-2: CCAGAUGAAGCUCCCAGAA; AKT3-siRNA-1: UCAUCUUUCUCCUUCAUUA; AKT3-siRNA-2: AGAGGUGUUAGAAGAUAAU; SMYD4-siRNA-1: AGACUCCAUCUCAAAGAAA; SMYD4-siRNA-2: AAAGGUUACUUGGUGGGAA; ATAD3A-siRNA-1: ACUUCUCAAUGAGGAGAAU; ATAD3A-siRNA-2: AGGACAAAUGGAGCAACUU; siRNA-NC: UUCUCCGAACGUGUCACGUTT.

To overexpress SIRPB1, AHNAK2, and XKR6, plasmids containing the coding sequences of these three genes or empty vectors were transfected into HMV-II cells (Supplementary Table 2). See 2.11 for the transfection procedure.

After washing with phosphate-buffered saline (PBS), cells were trypsinized and resuspended in complete culture medium. Cells were seeded in 96-well plates (5,000 cells per well). CCK-8 reagent (Promega, USA) was added to each well after 1, 2, 3, and 4 days. After incubation at 37 °C for 1 h, the optical density (OD) at 450 nm was measured. Experiments were repeated three times.

Transwell assays were performed using Corning Transwell chambers (Corning, USA) featuring 8-μm membrane pores. First, 15 μL of diluted Matrigel (diluted at a ratio of 1:8) was added to the upper chamber, followed by seeding 5 × 10⁵ HMV-II cells. After 24 h of incubation, the upper chamber was removed, and the excess liquid in it was discarded. The cells attached to the lower surface of the upper chamber were fixed in methanol for 30 min at room temperature and then stained with 0.1% crystal violet for 10–20 min at room temperature. The cell images were captured using a microscope. This experiment was repeated three times.

HMV-II cells were lysed using TRIzol reagent (Invitrogen, New York, USA) to extract total RNA. qRT-PCR analysis was carried out with the HiScript II One Step qRT-PCR SYBR Green kit (Vazyme, Jiangsu, China) on a Bio-Rad CFX96 PCR system (Bio-Rad, Hercules, California, USA). Primers were designed and synthesized by RuiBiotech (Beijing, China). Using GAPDH as the internal reference, the expression levels of target genes were calculated via the 2^−ΔΔCt method. Primer sequences are presented in Supplementary Table 3. The experiment was repeated three times.

We extracted total proteins from cells using a lysis buffer (AR0105, Boster, China). Subsequently, the proteins were separated by 10% sodium dodecyl sulfate-polyacrylamide gel electrophoresis and then transferred onto a polyvinylidene fluoride (PVDF) membrane. Prior to adding specific primary antibodies, the membrane was blocked with 5% nonfat milk at room temperature for 1.5 h. After that, the membrane was incubated overnight at 4 °C with the following primary antibodies: anti-P53 (345567, 1:1000) and anti-ATAD3A (124203, 1:1000) from Zenbio, and anti-AKT3 (21641, 1:1000), anti-SMYD4 (17594, 1:1000), anti-NF-κB (66535, 1:1000), anti-p-NF-κB (82335, 1:2000), anti-AKT (10176, 1:2000), and anti-p-AKT (66444, 1:2000) from Proteintech. To detect the protein bands, the membrane was incubated with a horseradish peroxidase-conjugated secondary antibody and an enhanced chemiluminescence reagent (AR1191S, Boster, China). Finally, band density was quantified using ImageJ software and normalized to the loading control.

For Co-IP analysis, the treated cell lysates were incubated with target antibodies overnight at 4 °C. Then, immunoprecipitation was carried out using Protein A+G agarose (P2012, Beyotime, China) for 4 h at 4 °C. Subsequently, the samples were subjected to western blot analysis to detect co-precipitated proteins. The total cell lysates were also analyzed by western blotting as input controls. Finally, the protein interactions were quantitatively analyzed using ImageJ software.

We stained paraffin sections using an immunohistochemical staining kit (Elabscience, China, E-IR-R215). Tissue sections were first dewaxed, hydrated, subjected to antigen retrieval, inactivated for endogenous enzymes, and blocked with serum. Anti-SIRPB1 (11811, 1:200, Proteintech), anti-AHNAK2 (17682, 1:200, Proteintech), and anti-XKR6 (23968, 1:200, Proteintech) were then applied to the sections and incubated at 37 °C for 2 h. After three washes with PBS, secondary antibodies (Elabscience, China, E-IR-R215B) were applied, incubated at 37 °C for 30 min, and then washed thrice with PBS. Finally, DAB staining solution was applied to the sections; a brownish-yellow color indicated a positive result. Pathological images were quantified and analyzed using ImageJ.

All data in this study were based on at least three biological replicates, and all statistical analyses were conducted in R (version 4.4.1). The Mann-Whitney U test was used for intergroup comparisons. PathScan was used for pathway enrichment analysis to analyze KEGG-related signaling pathways associated with high-frequency mutation genes; the screening condition was P.adjust < 0.05. SKCM in the TCGA database uses the Chi-squared (χ2) test to analyze differences between the tumor and non-tumor groups CNV. All P-values were calculated using two-tailed tests, and a P-value < 0.05 was considered statistically significant.

The surgically resected tissues from three patients with PMMC underwent WES (tumor sample: 200×depth, normal tissue: 100×depth) (Supplementary Table 4). We used the muTect software (33) to search for Somatic SNV sites and showed that MM01, MM02, and MM03 have Somatic SNVs of 282, 871, and 381, respectively. The vast majority of SNVs are distributed in the exon coding region (CDS) and intronic. The SNVs in the CDS are mainly synonymous mutations and missense mutations (Table 1).

InDels in the coding region or at a splice site may alter protein translation. Frameshift mutation occur when the length of an inserted or deleted base string is a non-integer multiple of 3, causing the entire reading frame to shift. Compared with non-frameshift mutations, frameshift mutation exert a greater impact on gene function, and are also subjected to stronger screening pressure (34). Using Strelka to detect Somatic InDel information, this study showed that MM01, MM02, and MM03 have 16, 30, and 99, Somatic InDels respectively. The vast majority of InDel occurs in the CDS and intronic regions, among which exons, the coding regions are mainly frameshift deletion, frameshift insertion, and non-frameshift deletion (Table 2).

Samtools software was used to detect Germline Mutations (SNPs, InDels) in the patient’s normal tissues, and then the detected mutant genes were screened for possible cancer-predisposing genes. Finally, using the Cancer Gene Census (CGC) database (http://cancer.sanger.ac.uk/cancergenome/projects/census/), the FACD (Familial Cancer Database, http://www.familialcancerdatabase.nl/) database, intoGen (Rubio-Perez C, 2015) database, the above screening results were annotated. In this study, the susceptibility genes in the adjacent tissues of PMMC were mainly NCOA4, ATRX, BPTF, NTRK1, RASGRP1, BAZ2B, NCOR2, KMT2D, and the main mutation types were Missense Mutation and In Frame Ins (Figure 1, Supplementary Table 5).

In our PMMC study, the most common mutation type was C>T/G>A, followed by T>C/A>G, C>A/G>T, C>G/G>C, T>A/A>T, and T>G/A>C (Figure 2A). Based on the frequency of 96 mutation types in PMMC samples, we performed NMF clustering and identified three mutation features (signature A, signature B, and signature C) (Figure 2B). MM01 and MM03 corresponded to signatures C and A, respectively, and MM02 corresponded to signature B (Figure 2C). After that, unsupervised hierarchical clustering of cosine similarity values further identified signature A (cosine similarity = 0.77). Moreover, signature C (cosine similarity = 0.75) showed a similarity with SBS5, and signature B (cosine similarity = 0.79) with SBS40 (Supplementary Figure 1). SBS5 showed similarity with the signature of a clock, with number of mutations correlating with aging, and SBS40 was also found to be associated with the age of cancer patients (Supplementary Table 6).

First, we compared the somatic mutations of the tumor sample with the known driver genes in the database and screened out the known driver genes in the tumor sample. It was found that the driver gene shared by MM01 and MM03 was TP53. The driver genes of MM01 are CACNA1D and ATRX, MM03 are PATZ1, EPHA3, MN1, GATA3, BCL9, NSD1, and EWSR1. The driver genes of MM02 are different from MM01 and MM03, and are NRAS, NEB, GRIN3A, SMYD4, SPEG, KEAP1, TLR1, and APC (Figure 3A). Driver gene mutations in PMMC tumor samples are similar to those in para-cancerous tissues, and remain dominated by Missense and In Frame Ins (Supplementary Table 7).

After comprehensively considering the somatic SNV and InDel mutations, the high-frequency mutation genes were screened out as AHNAK2, SIRPB1, XKR6, ATAD3A, TP53, AKT3, SMYD4, GBF1, and DSCAML1 (Figure 3B). The main types of mutations were missense and in-frame (Supplementary Table 8).

Further analysis of 1366 SKCM patients in the COSMIC database, the top 30 genes were BRAF, NRAS, RAC1, CHCHD2, RPS27, PPP6C, OR4M2, AC134980.3, ISX, DNAH5, MAP2K1, HBE1, RGS7, OR51B2, BCL2L12, OR5B17, KNSTRN, IDH1, ARID2, POTEH, STK19, EPPK1, CNTNAP2, NDUFB9, ERGIC3, and ADAM2 (Figure 3C). At the same time, analysis of 467 SNP data of patients with SKCM in the TCGA database revealed the first 30 genes as TTN, MUC16, DNAH5, BRAF, PCLO, LRP1B, ADGRV1, ANK3, CSMD1, CSMD2, HYDIN, APOB, DNAH7, XIRP2, MGAM, RP1, DSCAM, FAT4, DNAH8, CSMD3, PKHD1L1, DNAH9, DNAH3, MXRA5, FAT3, PTPRT, MYH2, FLG, SCN10A, and USH2A (Figure 3D). At the same moment, when we analyzed the SNP data of 289 patients with CC in the TCGA database, we found that the first 30 genes were TTN, PIK3CA, KMT2C, MUC16, KMT2D, FLG, EP300, DMD, SYNE1, FBXW7, LRP1B, RYR2, DST, USH2A, ADGRV1, HUWE1, MUC17, LRP2, SYNE2, NAV3, NEB, PCLO, CREBBP, DNAH2, FAT1, MDN1, MUC5B, CSMD1, PRKDC, and PTEN (Supplementary Figure 2).

The high-frequency mutation genes in PMMC were compared with genes mutated ≥10 times in COSMIC SKCM and with the top 300 most frequently mutated genes of TCGA SKCM. The gene shared between COSMIC, CC, and PMMC was TP53; the gene shared by TCGA and PMMC was AHNAK2; and the gene shared by CC and PMMC was DSCAML1 (Supplementary Figure 3). In addition, CCDC6, SIRPB1, LIM2, VSTM1, WDR20, XKR6, ATAD3A, AR, AKT3, SMYD4, GBF1, and DSCAML1 were identified as high-frequency genes in PMMC, while their mutation frequencies were lower in COSMIC and TCGA (Figure 4A).

Finally, using the STRING website, WDR20, AR, TP53, AKT3, CCDC6, SMYD4, ATAD3A can form a protein-protein interaction network (Figure 4B); PathScan software analysis of high-frequency mutation gene: metabolic pathways mainly focus on mitogen-activated protein kinase (MAPK) signaling pathway, cGMP-PKG signaling pathway, Human T-cell leukemia virus 1 infection, adrenergic signaling in cardiomyocytes, cellular senescence, oxytocin signaling pathway, Wnt signaling pathway, breast cancer, Kaposi sarcoma-associated herpesvirus infection, B cell receptor signaling pathway, proteoglycans in cancer, estrogen signaling pathway, Human cytomegalovirus infection, gastric cancer, dilated cardiomyopathy, pathways in cancer, Cushing syndrome, mitophagy animal, salivary secretion, cortisol synthesis and secretion, hepatocellular carcinoma, Human papillomavirus infection, longevity regulating pathway, endometrial cancer, oocyte meiosis, sphingolipid signaling pathway, thyroid cancer, hypertrophic cardiomyopathy, central carbon metabolism in cancer, and cholinergic synapse (Figure 4C, Supplementary Table 9).

We selected P53-related AKT3, SMYD4, and ATAD3A as candidate genes. We then used siRNA to knock down P53, AKT3, SMYD4, and ATAD3A in HMV-II cells. qPCR results showed that both siRNA-1 and siRNA-2 reduced the mRNA expression of P53 (P < 0.0001), AKT3 (P < 0.0001), SMYD4 (P < 0.0001), and ATAD3A (P < 0.0001) compared with the siRNA NC group (Figure 5A). We then further evaluated the effects of knockdown of these four candidate genes on cell proliferation and invasion. We found that knockdown of P53 (P < 0.0001), AKT3 (P < 0.0001), SMYD4 (P < 0.0001), and ATAD3A (P < 0.0001) significantly reduced cell proliferation (Figure 5B) and invasion (Figures 5C, D) compared to their levels in the siRNA NC group. Based on this, we further explored the interactions between AKT3, SMYD4, ATAD3A, and P53. Co-IP experiments confirmed direct interactions between P53 and AKT3, SMYD4, and ATAD3A (Figure 5E).

Somatic CNVs in PMMC samples were assessed, revealing 114 gains and 40 losses in MM01, 93 gains and 10 losses in MM02, and 194 gains and 16 losses in MM03 (Figure 6A, Supplementary Table 10). The CNV burden in the three PMMC samples was further visualized using the Circos tool (Supplementary Figure 4, Supplementary Figure 5, Supplementary Figure 6). In TCGA SKCM (472 tumor samples and 467 non-tumor samples), 97 genes showed differential CNVs between tumor and normal groups (Figure 6B). When TCGA SKCM CNVs were compared with those from the three patients with PMMC (Figure 6C), 33 CNV genes were shared, and 5,712 genes exhibited CNVs only in 3 PMMCs. Intersecting these results with previously identified SMG showed that SIRPB1, AHNAK2, and XKR6 carried both CNVs and high-frequency mutations in all three PMMC samples (Figure 6D).

We found increased CNVs in SIRPB1, AHNAK2, and XKR6 in patients with PMMC and further evaluated their impact on gene expression and cellular pathways. Immunohistochemical staining of sections from patients with PMMC revealed statistically significant increases in expression of SIRPB1, AHNAK2, and XKR6 in tumor-associated regions compared with adjacent tissue (Figure 7A). Subsequently, we transfected HMV-II cells with SIRPB1, AHNAK2, and XKR6 genes. qPCR results demonstrated that, compared with the vector group, the expression levels of oe-SIRPB1 (P = 0.0007), oe-AHNAK2 (P < 0.0001), and oe-XKR6 (P = 0.0004) in HMV-II cells were significantly increased (Figure 7B). Finally, we found that when SIRPB1, AHNAK2, and XKR6 were overexpressed, p-NFκB and p-AKT were activated in the SIRPB1 (P < 0.0001), AHNAK2 (P = 0.0001), and XKR6 groups (P < 0.0001), respectively, compared with their levels in the vector group (Figures 7C, D). However, there were no statistically significant differences in NFκB and AKT levels between the SIRPB1, AHNAK2, and XKR6 groups and the vector group, respectively (Figures 7C, D).