Copy number variation and clinical response to chemotherapy and bevacizumab in the Czech metastatic colorectal cancer patients

Summary

Background: Despite bevacizumab being the first biological agent approved for the treatment of metastatic colorectal cancer (mCRC), there is not any established DNA biomarker to improve its efficacy and personalize the treatment. Materials and methods: Thirty patients with mCRC on bevacizumab therapy (15 with a good response and 15 with a poor response) from the University Hospital Olomouc were followed. Formalin-fixed paraffin-embedded (FFPE) samples were used for copy number variation (CNV) analysis using the OncoScan FFPE Assay Kit in order to capture approx. 900 tumor genes. Results: In the group of good responding patients, 102 genes (classified as ATPases, type AAA, neuronal signal transmission, regulation of transcription, and superior domain PH type), potentially significant positive predictive tumor biomarkers of bevacizumab treatment, were found. In the poorly responding group, 74 potentially negative predictive genes (classified as galectines, Jak-STAT signalling pathway, MAPK cascade, differentiation, and F-box associated domain) were identified. Conclusion: In the pilot study, we found promising copy number variation biomarkers of bevacizumab response in FFPE samples of mCRC patients. The validation phase should be focused especially on the genes associated with angiogenesis (AGRN, MAPK8, ARHGAP22, LGALS13, LGALS4, ZFP36, and MYC), tumorigenesis (DVL1), and tumor proliferation (IFNL1, IFNL2, IFNL3, MAP3K10, and MAP4K1).

Keywords:

bevacizumab – colorectal carcinoma – structural genetic variation

Introduction

Bevacizumab (Avastin®, F. Hoffman-La Roche AG, Basel, Switzerland) is a recombinant humanized monoclonal antibody that binds extracellularly to vascular endothelial growth factor A (VEGF-A), preventing its interaction with VEGF receptors (VEGFR) on the surface of endothelial cells. Bevacizumab inhibits the angiogenic activity of VEGF-A, thus limiting the formation of new blood vessels (antiangiogenic effect) and reducing the density of existing vasculature (antivascular effect) [1,2]. Bevacizumab used to be the mainstay of targeted biological colorectal cancer (CRC) therapy for patients with mutations in RAS genes (40–50% cases) [3] in whom cetuximab or panitumumab cannot be used [4]. According to Institute of Biostatistics and Analyses, Brno, Czech Republic, bevacizumab is usually part of the FOLFOX or XELOX combination treatment for CRC in Czech cancer centers. Compared to cetuximab, bevacizumab prolongs overall survival and progression-free survival (PFS) in right-sided RAS wildtype and BRAF wildtype and BRAF mutant tumors [5,6].

Despite bevacizumab being the first biological agent approved for metastatic colorectal cancer (mCRC), there is not any established DNA biomarker to improve its efficacy [7]. It does not mean that such a marker was not searched for: somatic mutations in NRAS, BRAF, and/or PIK3CA gene were suggested as a negative prognostic biomarker [8]. Even mutation in KRAS that disqualifies patients for cetuximab and panitumumab treatment and originally steered patient to bevacizumab treatment is suggested as a negative prognostic factor of bevacizumab [8,9]. Also, many single nucleotide polymorphisms (germ-line variants) in genes are involved in inflammation, immune system, and RAS signalling [10–13]. However, none of these markers was confirmed in independent studies or used clinically.

We hypothesize that somatic DNA variants, biomarkers of bevacizumab treatment, exist but were for a long time overlooked because of technical limitations – because they pose copy number variants (CNVs). Such type of variation now became accessible for testing using technology of molecular inversion probes.

 

Material and methods

Clinical records of 142 patients with metastatic colorectal cancer treated at University Hospital Olomouc were used to select 15 + 15 patients according to their clinical response to bevacizumab treatment. The tumor tissue was examined by an experienced pathologist and the percentage of tumor cells was determined. Formalin-fixed paraffin-embedded (FFPE) sections were used to isolate genomic DNA using the Cobas DNA Sample Preparation Kit (Roche). DNA quantification was performed by qPCR and was related to the concentration of the operational gene GAPDH in the sample compared to the standard [14].

Analysis of DNA segment copy number variation (CNV) was performed using the OncoScan FFPE Assay Kit 1.0 (Thermo Fisher Scientific), which uses molecular inversion probe (MIP) technology and is designed primarily for analysis in limited amounts of degraded FFPE-derived DNA. An input amount of 80 ng of DNA is sufficient. The probes are designed for more than 220,000 single nucleotide polymorphisms (SNPs) located every 50–120 kb in approx. 900 tumor genes. The OncoScan FFPE Assay protocol was performed according to the manufacturer’s instructions.

After scanning the arrays, raw data were obtained in the form of CELL files, which were analyzed using OncoScan Console 1.3 (Thermo Fisher Scientific). Quality control of MAPD and ndSNPQC parameters was also performed. Subsequently, the data were analyzed in R software [15] using the rCGH package [16]. The resulting segmentation data obtained by the segmentCGH function were normalized by the EMnormalize function with default parameter settings, except for the mergeVal parameter set to 0. This normalized segmentation file was divided into two subsets according to the experimental groups. Both these segmentation subsets were further analyzed with the GISTIC 2 tool [17] to identify significant tumor targets in the genome. The threshold log2 ratio value for the section with loss of DNA segments (deletion section) was set to −0.3 (19% decrease) and for the section with increased DNA segment copy number (amplified section) to 0.3 (23% increase). This resulted in the identification of sites in the genome with a demonstrably higher or lower copy number in each experimental group compared to the normal population and the genes located in these sections. The DAVID database version 6.8 [18] and its Functional Annotation Chart tool were used to annotate genes.

Results

The combination of chemotherapy and biological therapy (first-line FOLFOX and bevacizumab), mutated KRAS gene, wildtype BRAF gene, and PFS were taken into account in the selection process, which divided the 30 patients into a poor responding group (PFS ≤ 9 months) and a well responding group (PFS ≥10 months). The experimental groups consisted of 15 men and 15 women aged 42 to 76 years at the time of the initiation of bevacizumab treatment (with a median age of 62 years), from whom 13 metastases (meta) and 17 primary tumors (prim) were collected (Tab. 1, 2). Tissue samples from study patients were obtained with informed consent and the study was approved by the ethical committee of University Hospital Olomouc (NV15-31230A).

In the group of good responding patients, 102 genes, potentially significant positive predictive tumor biomarkers of bevacizumab treatment, were found. In the poor responding group, 74 potentially negative predictive genes were identified.

In the group of good responders, only the 18p11.32 region, where 14 genes were located, was significantly amplified, while multiple regions were significantly deleted: 1p36.33 (61 genes), 8p11.22 (2 genes), 10q11.23 (21 genes), 14q32.33 (2 genes), 16p13.3 (1 gene), and 20p12.1 (1 gene) (Fig. 1 A, B, on the left). In the poor responding group, regions with increased copy number of segments 8q24.21 (1 gene), 14q12 (1 gene) and 19q13.2 (72 genes) were found (Fig. 1A, on the right), but no deletion exceeded the threshold of 0.25 (Fig. 1B, on the right).

Using the Functional Annotation Chart function, the DAVID gene ontology tool version 6.8 divided the genes identified in the group of good responders into four groups at the 5% significance level. Almost all genes were lying within regions with a loss of DNA segments, except for the THOC1 gene, and accounted for no more than 5% of the total number of genes in this group (identified in this group of patients). The first group consisted of ATPase enzymes (ATAD3A, ATAD3B, and ATAD3C) occurring in the 1p36.33 region. The next group consisted of genes involved in neuronal signal transduction (AGRN, MAPK8, DVL1, CHAT, and SLC18A3) from 1p36.33 and 10q11.23. ERCC6 from 10q11.23 and THOC1 from 18p11.32 were included in the group involved in transcription regulation. A group of genes from the 1p36.33 and 10q11.23 regions encode proteins with the PH domain involved in signal transduction processes. Some of them have GTPase activity, while others bind to phospholipids (AGAP4, ACAP3, ARHGAP22, WDFY4, and PLEKHN1).

In the group of patients with poor response to bevacizumab treatment, five groups of amplified genes were generated at a 5% significance level after analysis in the DAVID tool. The genes accounted for a maximum of 11% of the total number of genes in this group and are almost exclusively found in the 19q13.2 region. The CLC, LGALS13, LGALS17A, LGALS14, LGALS16, LGALS4, LGALS7, and LGALS7B genes belong to the carbohydrate-binding galectins involved in apoptotic and hydrolytic processes.

The IFNL1, IFNL2, and IFNL3 genes are involved in the positive regulation of the immune response and contribute to the Jak-STAT signalling pathway. The MAP3K10, MAP4K1, ZFP36, PSMC4, PSMD8, MYC, and RASGRP4 genes are part of the MAPK signalling pathway. The EID2, EID2B, SIRT2, CATSPERG, DLL3, and GGN genes are associated with differentiation. The FBXO17, FBXO27, and NCCRP1 genes encode F-box proteins involved in protein degradation (Tab. 3).

Discussion

We performed the analysis on a cohort of 30 patients that was divided into two subcohorts depending on their good or poor response to bevacizumab treatment. We discovered several groups of genes that may be related to cancer initiation, development, and spread.

In the good responding patients’ group, we found three genes belonging to the AAA-ATPase family (namely ATAD3A, ATAD3B, and ATAD3C) in the DNA segment loss region. Of this group, the ATAD3A gene is most frequently found in patients with lung adenocarcinoma. Overproduction of ATAD3A is associated with increased resistance to therapy and poorer prognosis [19,20]. Similar findings also apply to astrocytomas, in which the production of ATAD3B may be involved in chemoresistance. In contrast, oligodendrogliomas do not produce ATAD3 protein at all and are among the treatment-sensitive gliomas [21]. ATAD3B somatic mutation is associated with shorter overall survival and appeared only in metastatic patients [22]. Treatment resistance has also been confirmed for these genes in patients with breast and prostate cancer [23, 24].

Furthermore, a group of genes involved in neuronal signal transmission was created in the annotation tool for this group of patients, including AGRN, MAPK8, DVL1, CHAT, and SLC18A3 genes, which were also located in the region with loss of DNA segments. Agrin (AGRN), a heparan sulphate polysaccharide found both on the cell surface and intracellularly, is important in angiogenesis [25,26]. This protein has been produced in increased amounts by bile duct and liver cancer cells as well as squamous cell carcinoma cells, where it has been involved in cell migration, cell adhesion and resistance to treatment [27,28].

Mitogen-activated protein kinases (MAPK8, formerly JNK1) are involved in many functions in the body – cell proliferation, differentiation, survival, stress response, apoptosis, and cell transformation by activated oncogenes in many cell types [29]. The overproduction of MAPK8 protein in healthy organisms promotes invasiveness and angiogenesis and is involved in tumor progression and resistance to cytostatic drugs [30–32]).

The DVL1 gene is part of the Wnt signalling pathway, which has a function in embryogenesis and tumorigenesis. Increased DVL1 production accelerates breast tumor growth and in CRC, it is associated with the perineural spread of tumor and liver metastasis [33–35]. Increased choline acetyltransferase (CHAT) production has been observed in patients with squamous cell carcinoma of the lung [36]. Increased SLC18A3 gene activity has been reported in CNS metastases of papillary thyroid carcinoma [37].

Another group of genes regulating transcription consisted of the ERCC6 gene from the region with a loss of DNA segments and the THOC1 gene from the region with an increased number of DNA segments. The ERCC6 protein, a DNA protection factor, is involved in cell hypertrophy, but when it is knocked out of function, cell proliferation is reduced, and apoptosis is triggered. B cells are also more sensitive to commonly used chemotherapeutic drugs [38]. THOC1, a subunit of the THO complex involved in the formation of mRNA ribonucleoprotein complexes, is produced to a greater extent in lung and ovarian tumors, but less so in the skin and testicular tumors [39]. In general, more THOC1 is formed in tumor cells, and in CRC patients, the expression level can distinguish patients with an aggressive phenotype and associated poor prognosis [40,41].

The last group of genes in patients with a good response to treatment consisted of genes functionally superior to the PH-type domain, namely ACAP3, PLEKHN1, AGAP4, ARHGAP22, and WDFY4. Of this group of genes, according to the literature, only the product of the ARHGAP22 gene is related to angiogenesis and has been linked to the development of diabetic retinopathy in patients with diabetes [42].

In patients with a poor response to treatment, the first identified group was a family of abundantly represented galectins, which find their applications in diverse cellular processes such as embryonic development, wound healing, apoptosis, intercellular junction formation, cell migration, immune response and, last but not least, malignant proliferation [43–45]. Representatives of this group were specifically the CLC gene, LGALS13, LGALS17A, LGALS14, LGALS16, LGALS4, LGALS7, and LGALS7B. The CLC gene was more highly expressed in patients with early CRC compared to those with late-onset disease [46]. The amount of protein corresponds to the number of eosinophils at the site of inflammation [47]. A group of genes on chromosome 19 encoding placental galectins (LGALS13, LGALS14, and LGALS16) trigger T cell apoptosis and are involved in immune tolerance [48]. During pregnancy, the product (PP13) of the LGALS13 gene is involved in the vasodilatation of maternal blood vessels required to increase blood flow to the fetus [49]. The LGALS4 gene has been reported to be less expressed in the tissue of CRC patients and was essentially absent in patients with the invasive form of the disease. The product of this gene is involved in cell cycle arrest and slowing cell migration and motility. LGALS4 expression reduces the production of proteins of the Wnt signalling pathway. The Wnt signalling pathway is important in many biological processes, such as cell differentiation, migration, and polarity. Dysregulation of this signalling pathway is a hallmark of CRC [50]. However, another research group found elevated levels of galectin 4 in the serum of CRC patients examined, and levels were even higher in metastatic disease [51,52]. Elevated levels of circulating galectins 2, 4, and 8 in the blood induce the secretion of cytokines and chemokines from the vascular endothelium, which promotes the formation of endothelial adhesion molecules and endothelial ducts as a part of angiogenesis [53].

The second group consisted of lambda interferons (IFNL1, IFNL2, and IFNL3), functionally classified as a group of genes associated with the Jak-STAT signalling pathway that triggers an immune response and antiviral, antiproliferative, and antitumor responses [54,55].

The third group was represented by the MAP3K10, MAP4K1, ZFP36, PSMC4, PSMD8, MYC, and RASGRP4 genes and functionally aligned to the MAPK cascade. In ductal adenocarcinoma of the pancreas, MAP3K10 and MAP4K1 genes have opposing effects. Loss of MAP3K10 expression decreases cell proliferation, whereas the absence of MAP4K1 is associated with the progression and development of invasive disease [56,57]. The opposite situation occurs when MAP4K1 is knocked out of function in colorectal cancer, and when tumor cell invasiveness is suppressed [58]. The ZFP36 gene has a tumor suppressor function. It suppresses the growth of colorectal cancer cells by regulating the expression of VEGF and COX-2. Its production is reduced in cancer cells, but it is produced in high amounts in healthy mucosa [59–61]. The PSMC4 gene promotes cell death caused by stress factors and increases intracellular protein ubiquitination [62]. The MYC (c-Myc) gene influences the regulation of angiogenesis, and an increased copy number of this gene indicates a worse disease prognosis for patients in clinical stage II and III disease [63,64]. Other genes are not well described in the literature.

The fourth group of genes was associated with differentiation and includes the following genes: EID2, EID2B, SIRT2, CATSPERG, DLL3, and GGN. The SIRT2 gene plays an important role in the cellular response to hypoxia through the regulation of the hypoxia-inducible factor (HIF-1a) [65]. The remaining genes have not been described at all or only marginally.

The last group consisted of genes from the F-box associated domain. These were the FBXO17, FBXO27, and NCCRP1 genes. None of these genes has been described in the literature concerning angiogenesis or cancer.

Increased expression of the DVL1 gene may be important in breast cancer carcinogenesis due to disruption of the Wnt signalling pathway [33,34]. It may be also involved in the development of cervical cancer, again through disruption of the Wnt signalling pathway [66]. The protein is formed to a greater extent in prostate cancer and this may be related to progression through the Wnt/beta-catenin pathway [67].

The hERG1 and aHIF-2a genes were found to be independent prognostic factors for a positive response to bevacizumab [68]. It appears that hERG1 is upstream of HIF-2a and the entire proangiogenic signalling pathway, confirming what has been previously shown in CRC cells. Indeed, in the latter model, hERG1 activity positively regulates HIF--2a expression and subsequently VEGF-A secretion, suggesting that hERG1-positive patients who have more aHIF-2a and greater VEGF-A secretion would benefit from VEGF-A blockade via bevacizumab treatment, as shown in the reported survival analysis. The importance of hERG1 and its downstream pathway is also suggested by the finding that hERG1-positive patients with KRAS mutations have longer disease PFS than patients with the same mutation but negative for hERG1 [68, 69].

 

Conclusion

In the pilot study of FFPE samples from mCRC patients (15 good responders and 15 poor responders to bevacizumab), using OncoScan method, we found promising copy number variation biomarkers. This small cohort of 30 patients can serve as a learning set, ready for validation using method amenable for a low number of targets but a higher throughput (such as PCR methods) in a larger cohort of patients. The validation phase should be focused especially on genes associated with angiogenesis (AGRN, MAPK8, ARHGAP22, LGALS13, LGALS4, ZFP36, and MYC), tumorigenesis (DVL1), and tumor proliferation (IFNL1, IFNL2, IFNL3, MAP3K10, and MAP4K1).

 

Contributors and supporting agencies

This study was supported by the following research programs: NV16-32198, LM2023033, LM2023053, EF16_013/0001674, LX22NPO5102, TN02000109, and IGA LF UP 2023_006.

 

Availability of data and materials

The data used in this study are available from the corresponding author upon request.