Transcriptional Changes of DNA Replication and Repair Factors Over Uveal Melanoma Subtypes

Background Uncontrolled replication is a process common to all cancers facilitated by the summation of changes accumulated as tumors progress. The aim of this study was to examine small groups of genes with known biology in replication and repair at the transcriptional and genomic levels, correlating alterations with survival in Uveal Melanoma tumor progression. Selected components of Pre-Replication, Pre-Initiation, and Replisome Complexes, DNA Damage Response and Mismatch Repair have been observed. Methods We have generated two groups for each gene examined above and below the average alteration level, and compared relative expression and survival across TCGA UVM subtypes based on somatic copy number alteration supported by DNA methylation and mRNA/miRNA/lncRNA expression. Significance between subtypes monosomic or disomic for chromosome 3 was determined by Fisher’s exact test. Kaplan Meier survival distribution based on disease specific survival was compared by log-rank test. Results Specific genes with significant alteration include MCM2 MCM4 and MCM5 of the Minichromosome Maintenance helicase complex, CDC45, MCM10, CIZ1, PCNA, FEN1, LIG1, POLD1, POLE, HUS1, CHECK1, ATRIP, MLH3, and MSH6. We found evidence of Exon 4 skipping in CIZ1 previously identified as a cancer variant and reportedly used as an early serum biomarker in lung cancer, accompanied by evidence of instability of a mononucleotide repeat in Intron 3. Mismatch Repair protein MLH3 was found to have splicing variations with deletions to both Exon 5 and Exon 7 simultaneously. PCNA, FEN1, and LIG1 had increased relative expression levels not due to their mutation or to copy number variation. Conclusion We have observed differences in relative and differential expression that support the concept that selected replication and repair genes and their products are causally involved in the origin and progression of uveal melanoma, suggesting specific avenues for early biomarker identification and also therapeutic approach.


Background
Comparatively uncontrolled replication carried out by highly evolutionarily conserved multiprotein complexes, is a process shared by all cancers. Although many other processes including immortality, epithelial to mesenchyme transition, telomere metabolism, metastasis etc. contribute to tumorigenesis, the summation of genomic alterations in a tumor must facilitate replication. Duplication in the transformed cell is achieved at the expense of decreased fidelity, making replication a focal point where heterogeneity is created that upon clonal selection leads to tumor expansion and survival. The status of several individual critical replication genes has been examined in The Cancer Genome Atlas (TCGA) efforts and auxiliary studies, with computational methods placing alterations into known pathways to identify potential targets for precision medicine. However, the behavior of replication factors as a group did not receive analogous systematic investigation, likely because they are thought of as being part of a process rather than a pathway.
In this study we have examined the status of small groups of genes with known biology in replication and repair at the genomic and transcriptional levels in one tumor type, Uveal Melanoma (UVM). TCGA has recently conducted an integrative analysis of this cancer, and the data providing its molecular underpinnings are now available to the public as part of the "Rare Tumor Project" [1]. While having a low incidence of 5.1 per million [2], UVM is the most common intraocular malignancy. It is highly lethal, with 50% of patients developing metastatic disease followed by 6-12 month survival from metastatic diagnosis [3]. There are no effective treatments for metastatic UVM, early diagnosis could lower lethality.
TCGA has examined this tumor type and generated datasets through whole exome sequencing (WES), whole genome sequencing (WGS), mRNA miRNA and lncRNA expression, DNA methylation, identification of immune infiltration, and detailed pathology with clinical outcome. Based on the combination of somatic copy number alterations (SCNA), DNA methylation patterns, mRNA miRNA and lncRNA expression, tumors were classified into four unsupervised clusters designated 1 through 4. Pathology and clinical data shows increasing malignancy across the clusters, with subtypes 3 and 4 having poorer outcomes that included shorter time to metastasis.
Selected components from the Pre-replication, Pre-initiation, Replisome, DNA Damage Repair (DDR), and Mismatch Repair (MMR) complexes have been investigated in this study for genomic and relative transcriptional changes across the UVM clusters. We have found alterations indicative of replication stress correlating with aneuploidy, increased malignancy, and decreased survival. Replication stress is thought to be an early and strong driving force in tumorigenesis. It is defined as impediment to the replication fork that causes slowing or stalling of the replication machinery and is brought on by a variety of factors broadly classified as "exogenous" and "endogenous" [4]. Examples of exogenous causes are radiation, therapeutic treatment, and diet, where endogenous causes are DNA structures, protein-DNA complexes, nucleotide pool availability, reactive oxidation species, transcription and replication complex collision, and mutation and expression alteration in tumor suppressor and oncogenes. Replication stress is common in cancer and has been previously suggested as therapeutically targetable [5,6].
Common components integral to replication involved in error prone bypass of damage are also likely effective points of intervention [7,8].
In a normal cell replication stress will activate DDR which prevents DNA damage from becoming fixed during replication and passed on in mitosis [9]. It does this in part by coordinating cell cycle control and providing a pause for repair, and also in some cases by triggering apoptosis. Dysregulated DDR in a tumor cell creates genomic instability that is cancer enabling, and makes the tumor become dependent on alternative pathways such as base excision repair (BER) that may be targetable (e.g. PARP inhibitors). In this study we find evidence of dysregulated DDR by relative expression differences, as well as possible dysregulation to Mismatch Repair pathway. Several replication and repair genes were observed to have transcripts with aberrant splicing. This included CIZ1 which may be useful as an early biomarker, ATRIP a component of DDR, and MLH3 which may indicate abnormalities in mismatch repair. While there is some evidence that MMR and DDR are connected [10], we have examined and discuss the repair pathways separately. Interestingly, components of the replication machinery that lack mutation and SCNA correlate none-the less with increased expression, suggesting transcriptional mechanisms were used to overcome replication stress and fork collapse in UVM progression.

Mutation analysis
Mutational findings in this report are based upon data generated by the TCGA Research network and can be found at: http://cancergenome.nih.gov/ [11,12].

Relative Expression Analysis
RNA-seq-derived exon expression levels were first visualized in heat maps. Integrative Genome Viewer (IGV) [13,14] using RNAseq data (https://portal.gdc.cancer.gov/ ) from the same case, to confirm the authenticity of the exon. The "Sashimi Plot" function in IGV was used to identify alternative splicing and isoforms from RNAseq data, with "Minimum junction coverage" routinely set at "4".
In the "R" environment (https://cran.cnr.berkeley.edu/) Z scores were calculated for each exon of each gene by mean-centering, using the tumor cohort as average, the log2 transformed RPKM values and dividing by the standard deviation, visualizing high (red), no change/no expression (white), and low (blue) and arranging data by UVM cluster assignments (1)(2)(3)(4) in heat maps.

Placement of Cases into "High" and "Low" Expression Groups
An output file containing Z scores for each exon was created and used to calculate an average Z score for each gene. This was regarded reasonable as structural variations to genes were found only with low frequency. UVM cases were sorted into two groups, those with average Z scores "Above" and those "Below" zero. After group designation, the cluster (1)(2)(3)(4) each case belonged to was identified and the numbers of cases "High" and "Low" for each cluster counted and displayed graphically using GraphPad Prism 6. Significant differences between "1 and 2" versus "3 and 4" were determined by Fisher's exact (two tailed), using GraphPad Quick Calcs. Because very few genes and few cases (total n =80, in each cluster n=15 through 23) were examined, "q" values were not calculated, with the rationale that doing so might increase "type II" errors. Placement of cases into "high" and "low" grouping was made relative to the total tumor cohort average. Conjecturally, all cases "Below" the tumor cohort average could also be "Above" an average made from appropriate adjacent normal tissues, which was not available for RNA analysis in UVM samples (normal tissues for each UVM case were available for DNA analysis). The procedure was evaluated qualitatively using BAP 1 and RPS19 as test genes, and compared to mRNA differential expression Z scores generated by cBioPortal whose outcomes were made from RNAseq Version 2 using RSEM to perform quantitation, with the expression distribution of each gene compared only to tumors diploid for that gene when a "normal" cohort was not available.

Clinical Data and Survival Analysis
The TCGA UVM cohort was made up of eighty matched tumor and normal specimens.
Tumors were obtained from patients that did not have previous systemic chemotherapy or focal radiation, with appropriate consent obtained from institutional review boards. A panel of five histopathologists with expertise in ocular pathology and melanoma, examined hematoxylin and eosin stained sections from paraffin embedded tumors defining tumor extent, cell morphology, pigmentation, mitotic index, and the presence of tumor-infiltrating lymphocytes and macrophages. The following information was also curated: "Tumor status" (date of last contact), "Vital status" (dead/alive), "Date of last contact", "Date of Death", "Cause of Death", "Other Cause of Death", "New Tumor Event after Initial Treatment", "Histology of New Tumor Event", "Site of New Tumor", "Other Site/new event", "Date of New Event", "Additional surgery", "Additional Treatment/Radiation", "Additional Treatment/Pharmaceutical".
In principal four types of survival analysis can be made with TCGA clinical data.
"Overall Survival" (OS) defined as the period from date of diagnosis until death from any cause, "Progression-Free Interval" (PFI) from date of diagnosis until the occurrence of an event in which the patient with or without the tumor does not get worse, "Disease-Free Interval" (DFI) date of diagnosis until first recurrence, and "Disease Specific Survival" (DSS) diagnosis date until death from the specific cancer type. All UVM Survival Curves constructed for this study were "Disease Specific Survival" curves, as recommended by TCGA Pan Cancer Guidelines https://wiki.nci.nih.gov/display/TCGAM/2017-03-23+PCA+Network. Kaplan Meier survival plots were constructed using GraphPad Prism 6.0 software. The Log-rank (Mantel-Cox) and Hazard Ratio tests were used to determine significance.

BAP1 and RPS19
A total of thirty-seven genes were studied (Table 1). Mutations and differential expression for the gene set have been depicted in an "Oncoprint" from cBioPortal (http://www.cbioportal.org/public-portal) ( Figure 1). In this study we examine expression by two algorithms, and use the terms "differential" versus "relative" to distinguish between them.
"Differential" expression is defined as being made with respect to a non-neoplastic or, as is the case for UVM, an estimated normal sample (see Materials and Methods). "Relative" expression is specified as comparison to the total UVM tumor cohort average and is made for both mRNA and individual exons.
BRCA 1 Associated Protein 1 (BAP 1) and RPS19 were used to compare the methods of analysis ( Table 2). Examination of BAP 1 differential expression Z scores calculated by cBioPortal, showed a two-fold greater inclusion of cases "Below" average expression. For this particular tumor type due to the relationship of monosomy 3 to subtypes and survival, the approach using an estimated reference alters comparison to relative expression for genes located on chromosome 3. SCNA subtypes 3 and 4 are both monosomic for chromosome 3 and constitute approximately half the total tumor cohort. In contrast using RPS19 as a test gene, a gene which codes for a 40S Ribosome complex protein with cytological location at chromosome 19q13.2, showed no significant difference in relative expression found between the study method and cBioPortal values. UVM do not have significant SCNA for chromosome 19, the four additional cases found in the "Below" group of cluster 3 are due to the use of RSEM verses RPKM. These results show incongruity between differential expression about an estimated normal value and relative expression about a tumor cohort average, when high numbers of cases are not diploid. We would like to note explicitly that presentation of the discrepancy is not meant to claim one set of calculations superior to the other, but to explain why additional calculations were made for relative expression.

Pre-Replication and Pre-Initiation Complexes
The  (Table 3), correlating with increased malignancy and decreased disease specific survival (P= 0.0001). Four of these cases can be seen using differential expression ( Figure 1). Half of the MCM2-7 complex components are located on chromosomes that the TCGA finds by SNP-based Copy Number Analysis to have copy number alterations that include monosomy chromosome 3, 8q and 6p gains. Comparing relative expression for all genes in this study found on chromosome 3 indicates relative expression levels do not always simply correlate with SCNA ( Figure 6A-B), reflecting the TCGA finding that expression subtypes are only partially concordant with SCNA subtypes [1].

MCM4 (8q11.21) has increased relative expression in clusters 3 and 4, with increased
expression having worse survival that correlate with CNV gains to chromosome 8q in cluster 4 (P= 0.0018) ( Figure 5, Table 3). MCM5 (22q12.3) had lower relative expression in UVM clusters 3 and 4, with cases having worse survival. For MCM5, the P values are significant but less convincing (P = 0.036), (see supplemental Figure 1 for data covering all MCM helicase components).
We include CDC45, GINS1-4, MCM10, and CIZ1 in the Pre-replication complex (see discussion). While CDC45 and MCM10 appeared highly expressed in the higher risk subtypes, alterations did not correlate with survivability (Table 3)  showed significant difference in the number of cases "above" and "below" the mean (P = 0.0042 by Fisher's exact test, two tailed), indicating clusters 3 and 4 had more cases with higher than average POLE expression than POLD1 (Figure 8 C).
Upon DNA damage the 9-1-1 response is activated to stop cell cycle progression ( Figure   3). HUS1 is a member of this pathway that acts as part of the replication "sliding clamp". Its relative expression was significantly increased across the UVM clusters (Table 4). RAD1, RAD9A, and RAD17, also members of the DDR complex, did not have significant expression or survival differences. RAD1 showed evidence of alternative splicing involving exon 3 previously reported as a natural isoform ( Figure 6D).
MLH3 and MSH6 transcripts had significant expression differences between clusters 1 and 2 versus 3 and 4 however there were no survival advantages or disadvantages ( Table 5). The MLH1 relative expression pattern across UVM SCNA subtypes did not completely behave as anticipated from simple correlation with CNV ( Figure 6). Because MSI was reported in CIZ1 Examination of UVM RNA-seq data for MLH3 led to the identification of an isoform by Sashimi plots with deletions in both exons 7 and exon 5 ( Figure 6C) simultaneously. MSH6 was also examined in IGV for alternative splicing, which was not observed.

Discussion
In this report, selected components of replication and repair processes of UVM were replication rate while possible, is not verifiable from current data. Another putative function for MCM2 due to its histone binding and chaperone capacity may be to orchestrate histone dynamics during replication [32], and it has also been linked conceptually to transcription [33]. UVM appear to have MCM concentration alterations including decreased MCM2, increased MCM4 and decreased MCM5 that correlate with increased malignancy and decreased survival.
Overexpression of MCM4 has been found previously in multiple transformed cell lines and numerous tumor types [34][35][36]. The unbalanced expression of MCM2-7 individual components likely affects the quantitative amount of functional helicase available, lowering the ability to load in excess and contributing to unregulated replication.
We discuss CIZ1 in the context of the pre-initiation complex. genuinely appear have alterations in the germline that also appear in the tumor ( Figure 7F).
Interestingly, examination of differential expression across tumor types ( Figure 7A) shows CIZ1 not much different between tumors and their normal tissue quantitatively, unlike MCM2 for example ( Figure 5A).
Reviewing relative expression for PCNA, FEN1, and LIG1 we found increase that implied their involvement in overcoming the replication stress evidenced by factors in the prereplication and pre-initiation complexes. PCNA is a ring shaped homo-trimer that encircles DNA [37]. It interacts competitively with many other factors at the PIP motif [38] [39] and is an essential co-factor for DNA polymerases during replication. It is involved in repair processes and can also be modified post-translationally by phosphorylation, ubiquitination, SUMO proteins, ISGylation, Acetylation, and S-nitrosylation. Each modification has corresponding biological functions that include proliferation, MMR inhibition, translesion synthesis, homologous recombination, genomic stability, and apoptosis, respectively. One of its major roles is to promote tolerance to DNA damage during replication [40]. Because of its role in proliferation, PCNA is a target for cancer therapy. Several small molecules that either block PCNA-Pol or PCNA trimer formation have been identified with proliferation-inhibitory effects in vitro [37] [41-43].
Polymerase processivity is activated by PCNA. Two polymerases, POLD1 and POLE, important in the replication of B form DNA were selected for observation, as examination of polymerase behavior it its entirety was beyond the scope of this study. POLD1 relative expression was decreased in subtypes 3 and 4, POLE was increased. The expression levels of POLD1 and POLE support a recently proposed model [44] in which a switch to POL  and away from POL occurs upon DNA damage. We have not yet addressed the behavior of POL , or of the specialized polymerases that recently have been hypothesized to assist the replication machinery in the prevention of replication stress [4].
In our examination of DDR components we saw significant decrease in relative expression of ATRIP (ATR Interacting Protein) in clusters 3 and 4 [45][46][47][48][49]. Nine of these cases (11% of UVM) were also found by differential expression (Figure 1). Exon 3 deletion was identified in ATRIP transcripts across the UVM subtypes ( Figure 4). PARADIGM and MARINa algorithms previously used to examine the DDR pathway did not find DDR impairment. ATRIP is located at 3p21.31 (very close to BAP1, 3p21.1), suggesting the disparity might be a result of using an estimated normal reference value. Interestingly as mentioned in the TCGA study [1], recent reports suggest that BAP1 itself may function in DDR [50] [51].

Conclusions
In summary, we have examined selected replication and repair factors in UVM with known biology for expression differences with and without mutation across TCGA SCNA subtypes. We have observed differences in expression that support the concept that these genes and their products are causally involved in the origin and progression of uveal melanoma. These data also suggest further avenues of research for biomarker identification and therapeutic approach. The few genes examined in this study are prototypical for the kind of findings that can be made examining RNAseq data using relative expression. We are currently modifying the methodology to be more inclusive, routine, and high throughput to rapidly screen all replication and repair genes across all TCGA tumor types.

Abbreviations
The