Next-generation sequencing shows the genomic features of ovarian clear cell cancer and compares the genetic architectures of high-grade serous ovarian cancer and clear cell carcinoma in ovarian and endometrial tissues

Sample collection

The ECCC sequencing data (MSK panel) were downloaded from the literature (DeLair et al., 2017) and the HGSOC sequencing data were downloaded from the cBioportal database (http://www.cbioportal.org/). Moreover, the OCCC data (1,021 panel) were collected from the Geneplus genomic data bank from May 2015 to August 2022. The detailed clinico-pathological and sample information were shown in Table S1. The retrospective study was designed and conducted in accordance with the Declaration of Helsinki. Written informed consent was granted in sample collection, gene sequencing, and data analysis, with the obtained information authorized for publishing. It was approved by the Institutional Review Board (IRB) of Taizhou Hospital of Zhejiang Province (No. K20220844).

Targeted capture sequencing

A tissue kit was used to extract genomic DNA (gDNA) from OCCC samples and matched normal tissues (Qiagen, Hilden, Germany). Following the manufacturer’s instructions, sequencing libraries were created using the KAPA DNA Library Preparation Kit (Kapa Biosystems, MA, USA). For both biopsy samples, barcoded libraries were hybridized to a panel of 1,021 genes with full exons, chosen introns from 288 genes, and selected regions from 733 genes. The comprehensive gene list was described in our earlier research (Wang et al., 2020). The DNBSEQ-T7RS (BGI, Shenzhen, China) with 100 bp paired-end reads was used to sequence the DNA.

Mutation, somatic interactions, and somatic copy number variation calling

MuTect was utilized to identify minor insertions and deletions (Indels), as well as single nucleotide variations (SNVs) (version 1.1.4). For quality control, somatic mutations were only found if they met the following criteria: (i) they were present in less than 1% of the population in the 1000 Genomes Project, the Exome Aggregation Consortium, and the Genome Aggregation Database (gnomAD) (https://gnomad.broadinstitute.org); (ii) they weren’t present in paired germline DNA from normal tissues; and (iii) they were found in at least five high-quality reads and without paired-end reads bias. The maftools was used to delineate the mutational landscape (Mayakonda et al., 2018).

The somatic interactions landscape was drawn by maftools (Mayakonda et al., 2018) to identify gene sets mutated in a mutually exclusive or co-occurring manner. A Fisher’s exact test on a 2*2 contingency table including frequencies of mutant and nonmutated samples is used to assess the pattern of exclusivity or co-occurrence for a pair of genes.

Focal somatic copy number variation (SCNV) was identified by CONTRA (Li et al., 2012) and the frequency of larger fragmental CNV was calculated. To compare the mutational frequency across each tumor subtype, the overlapped region was obtained from each panel. The top10 CNV genes were extracted from the HGSOC segment file and OCCC 1021 panel results.

Analysis of mutational signature and clonal architecture

The mutational signature analysis was performed with unfiltered somatic mutations using the R package YAPSA (Hubschmann et al., 2021) and matched to the COSMIC signature database (https://cancer.sanger.ac.uk/cosmic/signatures, Mutational Signatures V3.3).

The variant allele frequency (VAF) ratio of mutations was used to infer the clonality of mutational events (Gerlinger et al., 2012). Specifically, the VAF ratio was calculated by dividing the VAF of each mutation by the maximum VAF observed in the same sample. A higher VAF ratio suggested the respective event occurred at an earlier stage during tumor progression. Mutations with VAF ratios >0.6 were determined as clonal mutations while the rest were considered subclonal mutations.

Statistics

Two-sided Mann–Whitney and Fisher’s exact tests were performed on GraphPad Prism (version 8.02) or R (version 3.6.1) to generate the P value. For all tests, a P value <0.05 was considered statistically significant.

Somatic mutations and indels in OCCC, HGSOC, and ECCC

We examined the mutational burdens in 34 samples with OCCC and observed 538 somatic single nucleotide variations (SNVs) and small insertions and deletions (Indels) in the tumor samples. In OCCC, mutations of top10 genes were all relatively high (12%–56%), the main mutations were ARID1A (56%) and PIK3CA (53%), following TP53 (38%), KRAS (21%), KMT2C (18%). Low-density lipoprotein receptor-related protein 1B (LRP1B) mutation occurred only in samples without ARID1A truncating mutations in OCCC, while the ARID1B/CREBBP/ MET/RBM10 mutation occurred in the ARID1A in truncating group (Fig. 1A). After taking the intersection of the OCCC sequencing panel, 1,653 and 750 mutations were identified in HGSOC and ECCC cohorts respectively. According to a prior study, HGSOC had the highest TP53 mutation frequency (92%) and the low mutation rates for other genes (NF1, KMT2C, CDK12, LRP1B, and TOP2A were the only mutations greater than or equal to 5%) (Fig. S1A). In the ECCC cohort, the genes TP53 (43%), PIK3CA (40%), PPP2R1A (30%), ARID1A (23%), PIK3R1 (23%), FBXW7 (23%), and KMT2D (20%) had high-frequency mutations (Fig. S1B). We then compared the frequencies of mutations (top 10 genes) in known cancer-associated genes in OCCC, HGSOC, and ECCC. In result, we observed similar mutational frequencies in some canonical cancer driver genes between OCCC and ECCC cohorts, such as PIK3CA (OCCC 53% vs. ECCC 40%), TP53 (OCCC 38% vs. ECCC 43%), KRAS (OCCC 21% vs. ECCC 13%), APC (OCCC 9% vs. ECCC 13%), and KMT2C (OCCC 18% vs. ECCC 13%), while only the mutation rate of ARID1A in OCCC was significantly higher than that in ECCC (p < 0.01; Fig. 1B). In addition, only the mutation rate of TP53 was significantly higher in HGSOC than that in OCCC, while 17/27 the mutation rates of genes in OCCC were significantly higher than that in HGSOC (p < 0.001; Fig. 1B).

Further gene interaction analysis showed that DNMT3A (not shown in top 10 genes) and PIK3CA were mutually exclusive occurrences, APC and TP53/KMT2C/ACIN1/FBXW7/EGFR were co-occurrence relationships in OCCC samples (p < 0.05; Fig. 2A). In HGSOC, TP53 and FAT1 were mutually exclusive occurrence relationships, APC co-occurred with KMT2C/BRAC1, and MTOR co-occurred with KMT2A (p < 0.05, Fig. 2B). In ECCC, ARID1A co-occurred with PTEN/ KMT2C/ EPHA5/DICER1/ECOR/PI3KR1, and TP53 co-occurred with MEL, and PIK3CA co-occurred with PTEN/ EPHA5/BCOR/SMARCA4 (p < 0.05, Fig. 2C). These results suggested that there were differences in the nature of interacting oncogenes in these three tumor types.

Clonal architectures in OCCC, HGSOC, and ECCC

To explore the probable timing order of the mutation events arose in OCCC, HGSOC, and ECCC, we gauged the variant allele frequency (VAF) ratio of somatic mutations and Indels in these tumors with the method mentioned in previous literature (Gerlinger et al., 2012). The median of high/low VAF represents an early or late event, respectively. As shown in Fig. 3, driver genes like TP53 and RB1 were determined as clonal mutations with high VAF ratios (>0.6), indicating they might occur at an earlier stage compared with other subclonal mutations with low VAF ratios (<0.5), and might play crucial roles in HGSOC tumorigenesis. By contrast, the median VAF value of all genes was less than 0.5 in OCCC and ECCC, suggesting that all alterations were the later events in these two malignancies as subclonal mutations. Although the mutational frequency of ARID1A and PIK3CA were higher than other genes, the median VAF value of KRAS was highest and the TP53 mutation arose later in OCCC. In ECCC, the median VAF values of TP53, ARID1A, KRAS and PIK3CA mutation were low as later events. These results suggested that the driver genes triggered cancer in OCCC, ECCC, and HGSOC were different, and different mutations of genes were presented in these tumors at different stages, some mutations were newly obtained during OCCC and ECCC malignant evolution.

Mutational signatures in OCCC, HGSOC, and ECCC

The mutational signatures can reflect the potential carcinogens exposed in tumorigenesis and progress. We identified somatic SNVs and Indels corresponding Catalogue of Somatic Mutations in Cancer (COSMIC) single base signatures (SBS) and insertion and deletion signatures (ID) in OCCC. Somatic SNVs were mainly attributed to two SBSs in the COSMIC database, SBS7b, and SBS31 specifically, while Indels were mainly attributed to ID5, ID11, ID16, and ID17 (Fig. 4A). We then identified the SNVs mutational signatures patterns of OCCC, HGSOC, and ECCC to gain insight into the different etiology of those tumors. As shown in Fig. 4B, mutation signatures of OCCC were shared with ECCC, including SBS31 and SBS7b. C >T (cytosine>thymine transitions) mutational pattern of SBS7b was frequently found in cancers of the skin and associated with exposure to ultraviolet light (Alexandrov et al., 2020). SBS31, characterized by C >T mutations, may be due to platinum drug treatment. The proposed etiology of ECCC-specific SBS23 of the C >T pattern is still unknown. SBS39, featured with predominant C >G alterations, was shared in HGSOC and ECCC. It was found with high distribution in head and squamous cell carcinoma, breast cancer, and ovary-adenoma and without known proposed etiology (Alexandrov et al., 2020). SBS1 and SBS38 were private signatures in HGSOC. C >T change of SBS1 is an endogenous mutational process initiated by spontaneous or enzymatic deamination of 5-methylcytosine to thymine which generates G: T mismatches in double-strand DNA, which is considered as a cell division/mitotic clock. SBS38, a pattern of C>A mutation, was found only in melanomas with potential indirect damage from ultraviolet light. These results suggested that the potential mutational signature patterns in OCCC were more similar to ECCC than that in HGSOC.

Somatic copy number variation in OCCC, HGSOC, and ECCC

The focal somatic copy number variation (CNV) amplification and deletion were extracted from OCCC 1021 panel by CONTRA (Li et al., 2012). The results as shown in Fig. 5A and Table S2, we exhibited the frequencies of gain or loss of CNVs located on chromosomes (Chr) and focus the high frequency of CNV gain events in genes such as GNAS (0.97) on Chr20q13.32, NF1 (0.94) on Chr17q11.2, ASXL1 (0.91) on Chr20q11.21, TCF7L2 (0.91) on Chr10q25.2, CREBBP (0.81) on chr16p13.3, LRP1B (0.59) on 2q22.2 and loss events in ATM (0.97) on 11q22.3, BRCA2 (0.97) on 13q13.1, FANCM (0.97) on 14q21.2, RB1 (0.97) on 13q14.2, FLT3 (0.94) on 13q12.2, ROS1 (0.94) on 6q22.1 and JAK2 (0.91) on 9p24.1. Further comparing the top10 genes with a high frequency of CNV amplification or deletion between OCCC and HGSOC, we obtained significantly different cancer-associated genes at the CNV level in these two subtypes (Fig. 5B). In result, canonical oncogenes such as ARID1B, ARID1A1, ASXL1, CREBBP, EPHA2, GNAS, NBN, LRP1B, NF1, STAT3, TCF7L2, TOP1, TOP2A, KMT2C were gains in OCCC, whereas CAMTA1, CSF1, ETS1, OPCML were gained in HGSOC (p < 0.001, Fig. 5B). Tumor suppressor genes such as NKX3 −1, CBFB and WWOX were loss in HGSOC and AMT, BRCA2, FLT3, JAK2, KMT2C, PDGFRA, PIK3R1, RB1, ROS1 were loss in OCCC (p < 0.001, Fig. 5B).