Paired comparisons of mutational profiles before and after brachytherapy in asian uveal melanoma patients

Uveal melanoma(UM) is the most common primary intraocular malignancy in adults. However, the incidence of UM in Asia is 10 to 20 times less than in Western populations. Therefore, for the first time, we report our whole exome sequencing (WES) data analysis to discover differences in the molecular features of Asian and Western UM, and to determine the disparities between the primary tumor before brachytherapy and enucleated samples after brachytherapy. WES of 19 samples (13 primary tumors, 5 enucleation samples after brachytherapy, and 1 liver metastasis) from 13 patients diagnosed with UM and treated between 2007 and 2019 at the Yonsei University Health System (YUHS) were analyzed using bioinformatics pipelines. We identified significantly altered genes in Asian UM and changes in mutational profiles before and after brachytherapy using various algorithms. GNAQ, BAP1, GNA11, SF3B1 and CYSLTR2 were significantly mutated in Asian UM, which is similar that reported frequently in previous Western-based UM studies. There were also similar copy number alterations (M3, 1p loss, 6p gain, 8q gain) in both groups. In paired comparisons of the same patients, DICER1 and LRP1B were distinctly mutated only in tumor samples obtained after brachytherapy using rare-variant association tests (P = 0.01, 0.01, respectively). The mutational profiles of Asian UM were generally similar to the data from previous Western-based studies. DICER1 and LRP1B were newly mutated genes with statistical significance in the regrowth samples after brachytherapy compared to the primary tumors, which may be related to resistance to brachytherapy.

Uveal melanoma(UM) is the most common primary intraocular malignancy in adults. However, the incidence of UM in Asia is 10 to 20 times less than in Western populations. Therefore, for the first time, we report our whole exome sequencing (WES) data analysis to discover differences in the molecular features of Asian and Western UM, and to determine the disparities between the primary tumor before brachytherapy and enucleated samples after brachytherapy. WES of 19 samples (13 primary tumors, 5 enucleation samples after brachytherapy, and 1 liver metastasis) from 13 patients diagnosed with UM and treated between 2007 and 2019 at the Yonsei University Health System (YUHS) were analyzed using bioinformatics pipelines. We identified significantly altered genes in Asian UM and changes in mutational profiles before and after brachytherapy using various algorithms. GNAQ, BAP1, GNA11, SF3B1 and CYSLTR2 were significantly mutated in Asian UM, which is similar that reported frequently in previous Western-based UM studies. There were also similar copy number alterations (M3, 1p loss, 6p gain, 8q gain) in both groups. In paired comparisons of the same patients, DICER1 and LRP1B were distinctly mutated only in tumor samples obtained after brachytherapy using rare-variant association tests (P = 0.01, 0.01, respectively). The mutational profiles of Asian UM were generally similar to the data from previous Western-based studies. DICER1 and LRP1B were newly mutated genes with statistical significance in the regrowth samples after brachytherapy compared to the primary tumors, which may be related to resistance to brachytherapy.
Uveal melanoma (UM) is the most common primary intraocular malignancy in adults. The incidence of UM varies widely between races. In Western populations, the annual incidence is 5 to 10 per million population per year 1,2 , whereas in Asia, it is reportedly much lower at 0.4 to 0.6 3,4 . Various radiotherapeutic and local therapeutic options have been applied to UM. The Collaborative Ocular Melanoma Study (COMS) trial showed that UMrelated mortality rates were not significantly different between enucleation and plaque brachytherapy in the patients with medium-sized UM 5,6 . The results justified the use of plaque radiotherapy rather than enucleation for most medium-sized UMs. Since then, this treatment has been expanded to apply to small as well as large UMs. Despite these treatment efforts, in over one-quarter to one-third of patients, metastasis develops within 10 years, usually involving the liver, and death typically occurs 1-3 years after treatment 7,8 . Until now, the mutation profiles of Asia UM were little known; for the first time, we performed whole exome sequencing (WES) in Asian UM, and compared the results with the Cancer Genome Atlas (TCGA) 9 data, which is the largest cohort composed mainly of Western patients. In addition, we compared the mutation profiles by WES before and after brachytherapy in eyes enucleated due to post-brachytherapy tumor regrowth. Schematic diagram of workflow is written in Supplementary Fig. 1

Results
Patient characteristics. This study included 13 UM patients that have been diagnosed and treated at the Yonsei university health system (YUHS) from August 2007 to December 2019. The study was approved by the Institutional Review Board (IRB) at YUHS. and written informed consent was obtained All study protocols adhered to the tenets of the Declaration of Helsinki. All patients were treated with brachytherapy after local resection. Enucleation was offered as the principal treatment if the patient had a large tumor (height > 10.0 mm), however, the patients who strongly refused primary enucleation received the brachytherapy. Local resection was performed for tumor debulking or diagnostic confirmation by endoresection for choroidal tumors and exoresection for iris or ciliary body tumors. Afterwards, brachytherapy with 106 Ru plaques (Eckert & Ziegler BEBIG, Berlin, Germany) was performed, with target tumor apex radiation doses ranging 85-100 Gy. When local recurrence (tumor regrowth) was noted during follow up, enucleation was performed.
Demographical and clinical data are summarized in Table 1. The analyzed samples were composed of triple types of samples (primary, enucleation after brachytherapy, and liver metastasis). Overall, 13 samples were primary tumors and 5 were paired enucleation samples after brachytherapy. Triple -paired samples were obtained from one patient.

Somatic variant detection. Genome Analysis
Tool kit (GATK) 10 4.1.0.0. Mutect2 11 without matched normal pipeline was used to call somatic variant. The average target region sequencing depth and the average of ontarget rate of preprocessed bam were 87.21x (SD 11.28) and 96% (SD 0.5%), respectively. The sequencing coverage and quality statistics are provided in Supplementary data.
Since we have called variant of FFPE samples without matched blood samples, the optimized filtration pipeline for tumor only sequencing 12 was applied applied on somatic variant candidates from MuTect2. After filtration to distinguish somatic alterations from germline and sequencing artifacts, the average number of somatic alterations of primary samples (n = 13) was 21.46 (SD 12.38). This was not statistically different to the number of somatic mutations in TCGA UM (Wilcoxon Rank Sum test, P = 0.203, Supplementary Fig. 2).
Significant mutated gene analysis. In MutSigCV, a total of 10 genes were found to be significant in the YUHS data (Fig. 1A, Supplementary Fig. 3). GNAQ, BAP1, GNA11, SF3B1, EIF1AX, PTPRD and CYSLTR2 were similarly found to be significant genes in the TCGA data using the same algorithm, whereas the other 3 www.nature.com/scientificreports/ genes (SLFN11, KTN1 and FANCL) were not. Using the OncodriveCLUST algorithm, 6 genes (GNA11, BAP1, GNAQ, SF3B1, CYSLTR2 and SLFN11) were noted to be significant in the YUHS data (Fig. 1A, Supplementary  Fig. 4); and except for BAP1, CYSLTR2 and SLFN11, these were also significant in the TCGA data. Altogether, GNAQ, BAP1, SF3B1, GNA11, and CYSLTR2 were significant in both algorithms, and have been reported as frequently mutated genes in a previous UM study 13 . We could find common recurrent alterations (p. Molecular changes between before and after brachytherapy. In paired comparisons using the SKAT-O test, DICER1 and LRP1B were newly mutated genes with statistical significance in enucleation samples when compared to the primary tumor of the same patients (P = 0.01 and 0.01, respectively; Fig. 1A). The 4 types of somatic alterations of DICER1 and 5 somatic alterations of LRP1B were exclusively present in the enucleation samples (Table 2). In a triple comparison of one patient (YUHS 12), comprising of the primary tumor, an enucleated sample after brachytherapy, and a liver metastasis sample, when the latter samples were compared to the primary tumor, one gene (METRK) mutations in the enucleation sample and six gene (FANCA, KMT2D, PCLO, PGR, TET3, and ZBTB7A) mutations in the metastasis were found to be newly occurred, and were each exclusive without common mutations (Fig. 1B).
Mutational signature analysis. Mutational signature analysis with primary variants data indicated mainly 5 signatures (Cosine similarity: 0.919), which were single base substitution (SBS) 5 (34.44%), SBS32 (23.1%), SBS1 (17.46%), SBS7b (17.28%), and SBS7a (7.76%) in primary samples (Fig. 2). SBS7b and SBS7a have been found in skin cancers from sun exposed areas and thus likely to be due to exposure to ultraviolet (UV) light. However, none of the signatures from TCGA UM was related to exposure to UV (Supplementary Fig. 13). www.nature.com/scientificreports/ Additionally, SBS5 was also found in TCGA UM signatures and SBS1 has been known to be related to various malignancies and aging. Therefore, we found both known signatures (SBS1, SBS5) from previous Western-based UM studies. Similar to a previous study showing ultraviolet radiation signatures were observed in iris tumors 14 , SBS7b and SBS7a were predominantly (≥ 50%) observed in YUHS iris melanoma samples (YUHS7, YUHS9) ( Supplementary Fig. 14).

Discussion
We analyzed Asian UM and confirmed that significantly mutated genes were similar to those of TCGA, the largest western cohort.
In primary sample analysis, we found similar mutational profiles with western data. The 5 genes (GNAQ, BAP1, SF3B1, GNA11 and CYSTLR2) that have been reported as recurrent mutated or driver genes in other UM study were significant in two algorithms (MutsigCV and OncodriveCLUST). In addition, SLFN11 that never have been reported as UM related gene was also significant in both of them. This could be novel UM related cancer gene. The previously reported recurrent CNAs of UM were similarly observed in Asian UMs. Although the proportion of some CNAs (6p gain, M3) was different between the two cohorts, it is difficult to determine  www.nature.com/scientificreports/ whether the results are from a difference between Asian and Western populations due to the small sample size of our cohort. In the cases where the tumor regrows post-brachytherapy, enucleation should be considered. We performed paired comparison of WES data before and after brachytherapy and confirmed that two novel genes, DICER1 and LRP1B, were significantly mutated in the regrowth samples compared to the primary tumors. These two genes have already been associated with malignancy in non-ocular tumors 15 . It is not clear whether these newly developed mutations are the result of cancer evolution or the effect of brachytherapy; however, they could be related to resistance to brachytherapy.
In a triple comparison of one patient (Fig. 1B), the enucleated regrowth sample after brachytherapy and liver metastasis sample did not share any newly developed mutations. This may indicate that the liver metastasis event occurred earlier and independently of the tumor regrowth after brachytherapy.

Conclusion
This report, albeit a small sample size, is the first WES analysis in Asian UM, and through paired comparisons, novel and insightful results can be drawn.

Methods
Sample preparation and data generation. DNA was extracted from 19 formalin fixed paraffin embedded (FFPE) tissue from 13 UM patients that have been diagnosed and treated at the Yonsei university health system (YUHS) from August 2007 to December 2019 and captured using SureSelectHuman All Exon V6 for whole exome sequencing. After library QC of DNA, the libraries were sequenced on the Illumina HiSeq X ten platform. This study was approved by the Institutional Review Board (IRB) at YUHS and written informed consent was obtained. All study protocols adhered to the tenets of the Declaration of Helsinki.
Sequencing data processing and somatic variant analysis. The paired end reads of FASTA files from sequencing were aligned using BWA-mem 16 on human genome (hg) 19. After duplicated reads of the aligned bam files were marked and removed using Picard, the base quality of reads in bam files was recalibrated using Genome Analysis Tool kit (GATK) 10

Variant filtration for only tumor sequencing and Annotation.
To remove sequencing artifact, the variants with "bad_haplotype", "chimeric_original_alignment","base_quality","duplicate_evidence","fragment_ length","low_avg_alt_quality","mapping_quality","multiallelic","n_ratio","read_orientation_artifact", "read_ position", "str_contraction","strand_artifact" and "strict_strand_bias" were filtered out. We selected genes that have been reported as cancer genes in both OncoKB 18 , cBioportal 19 , for uveal melanoma. We also added MAP-KAPK5 which was not listed in OncoKB but it has been reported to be a frequently mutated gene in uveal melanoma data from TCGA. If the population allele frequency of the variants is more than 1% in any subpopulation of 1000 genomes project 20 , ExAC 21 , KOVA 22 , gNomad and Korean Genome Project 23 data, the variants were excluded as germline mutations. The second exclusion criteria is that the variant were present in the Korean 1,000 depression exome data 24 , which is a panel of normal. It removes not only germline mutation but also potential platform specific artifacts. Then we reviewed the all variants using Integrative Genomics Viewer 25 . If there is a missed variant that is found in another tissue from a same patient, we used blastn 26 to align all sequence on neighboring region(± 50 bp) of the variant. We consider the variant exist if more than each 2 paired aligned reads that encoding same nucleotide of variants. The variants which passed the above criteria were selected as analysis ready somatic variants. All passed variants were annotated with SIFT 27 , Polylphen2 28 and CADD 29 score using ANNOVAR 30 .
Significant mutated genes. The filtered non-synonymous variants were analyzed for significantly mutated genes using two algorithms, MutSigCV 31 and OncodriveCLUST 32 . We defined the genes that have P less than 0.05 as significant in MutSicCV. on Genepattern 33 . with default parameter. OncodriveCLUST is used to identify genes with a significant bias of mutation clustering within the protein sequence. If Q value is less than 0.05 in OncodriveCLUST, the gene is significant. Maftools 34 was used to run OncodriveCLUST. We also ran the same algorithms on uveal melanoma somatic variants data from TCGA. The somatic variant data were produced from tumor bam files from GDC data portal with same calling pipeline for YUHS.
Identifying copy number alterations. CNVkit 35 was used to identify copy number alterations in UMs.
Because construction of a normal reference from pooled normal samples is necessary for the somatic copy number calling pipeline, 1,000 WES of normal blood from the Korean 1000 database were used. Fused Lasso (flasso) algorithm was used as a segmentation method. The thresholds of log2 transformed relative ratio to reference ploidy for copy number 0,1,2,3 were -1.1, -0.4, 0.3, 0.7, respectively. The same pipeline was also utilized on TCGA tumor-normal paired data.
Rare variant association test. Since single variant association test has little power for testing of rare variants in relatively small samples, an alternative approach is needed. Sequence kernel association optimal test (SKAT-O) 36 is an optimal unified approach for rare variant association testing in case control sequencing studies. The test was applied on non-synonymous variants of genes between primary and matched enucleation samples. R library SKAT was used to run the statistical test. www.nature.com/scientificreports/ Mutational signature analysis. The 96 different contexts of single base substitutions from filtered variants were analyzed for generating mutational signatures. The signatures are identified as causes of mutational process, due to their unique mutational pattern and specific activity on the genome. The python packages, Sig-ProfilerMatrixGenerator and SigProfilerExtractor, were used for this analysis. The iteration time was set 100 times to be performed to extract each signature.

Data availability
All data generated or analyzed during this study are included in this article. Raw