VEGF-Related Germinal Polymorphisms May Identify a Subgroup of Breast Cancer Patients with Favorable Outcome under Bevacizumab-Based Therapy—A Message from COMET, a French Unicancer Multicentric Study

The prospective multicenter COMET trial followed a cohort of 306 consecutive metastatic breast cancer patients receiving bevacizumab and paclitaxel as first-line chemotherapy. This study was intended to identify and validate reliable biomarkers to better predict bevacizumab treatment outcomes and allow for a more personalized use of this antiangiogenic agent. To that end, we aimed to establish risk scores for survival prognosis dichotomization based on classic clinico-pathological criteria combined or not with single nucleotide polymorphisms (SNPs). The genomic DNA of 306 patients was extracted and a panel of 13 SNPs, covering seven genes previously documented to be potentially involved in drug response, were analyzed by means of high-throughput genotyping. In receiver operating characteristic (ROC) analyses, the hazard model based on a triple-negative cancer phenotype variable, combined with specific SNPs in VEGFA (rs833061), VEGFR1 (rs9582036) and VEGFR2 (rs1870377), had the highest predictive value. The overall survival hazard ratio of patients assigned to the poor prognosis group based on this model was 3.21 (95% CI (2.33–4.42); p < 0.001). We propose that combining this pharmacogenetic approach with classical clinico-pathological characteristics could markedly improve clinical decision-making for breast cancer patients receiving bevacizumab-based therapy.


Introduction
The vascular endothelial-derived growth factor A (VEGFA) is a positive regulator of angiogenesis and a coordinator of vascular homeostasis. In addition, the capacity of VEGFA to act as an immunosuppressor is well-established [1]. Therefore, targeting VEGFA and/or its receptors (VEGFR1, 2 and 3) constitutes a direct means to influence the tumoral environment. Bevacizumab, a VEGFA targeting monoclonal antibody, entered clinical practice more than 15 years ago [2] and is currently the most prescribed antiangiogenic treatment. Indeed, bevacizumab is indicated for colorectal cancer, breast cancer, non-small cell lung carcinoma (NSCLC), glioblastoma, renal cell carcinoma, ovarian cancer and cervical cancer [3]. The immunosuppressive properties of VEGFA have recently opened new perspectives regarding bevacizumab's immunomodulatory potential. Indeed, bevacizumab has recently demonstrated clinical benefits for NSCLC and hepatocellular carcinoma in combined cancer immunotherapy [2,4]. However, despite intense investigations, reliable biomarker signatures that would allow us to pinpoint the individuals most susceptible to benefit from bevacizumab-based treatments remain to be identified and validated. Germline genetic polymorphisms may influence the outcome of cancer treatments in several aspects, such as individual variability in drug metabolism, immune response [5] and target availability [6,7].
The French multicentric cohort study COMET aimed to identify biological factors that could predict the clinical benefit of a bevacizumab-paclitaxel combination therapy as first-line treatment in metastatic breast cancer. The COMET study was specifically designed to survey key candidate genes potentially linked to the pharmacogenetics of both bevacizumab and paclitaxel. Notably, the elimination of paclitaxel is known to be influenced by CYP2C8 and ABCB1 polymorphisms [8,9]. In addition, by following a small cohort of metastatic breast cancer patients receiving another bevacizumab-based combination (bevacizumab-taxane), we previously established the prognostic value of the single nucleotide polymorphism (SNP) 936 C > T (rs3025039) in VEGFA [10]. The ECOG 2100 trial demonstrated the importance of polymorphisms in the VEGFA and VEGFR2 genes for predicting the outcome of patients treated with the bevacizumab-paclitaxel combination compared with paclitaxel alone [11].
In the present study, 13 different SNPs across six genes were analyzed in an effort to identify a specific SNP signature significantly associated with bevacizumab-paclitaxel treatment outcome in terms of survival among metastatic breast cancer patients. The risk score analyses evaluating their predictive value as stand-alone variables versus when combined with common pathophysiologic criteria used in clinical practice are presented.

Study Population and Clinical Characteristics of Tumors
Out of the 342 patients included, 306 genomic DNA samples were available for gene polymorphism analysis. The patients' disposition is presented in Figure 1. The patients' characteristics, the disease Pharmaceuticals 2020, 13, 414 3 of 12 and the treatment baselines are summarized in Table 1. The mean age at inclusion was 55 years (range:  28-80). At initial diagnosis, the main histological type was invasive ductal carcinoma (79.5%). Most patients had received prior neoadjuvant/adjuvant chemotherapy (n = 210, 68%). Patients presented a small number of metastatic sites (90% had less than three metastatic sites); 62.5% of patients had a metastasis-free interval greater than 24 months. One hundred and thirty-five patients (48.5%) had histological Grade II tumors and 43.5% had Grade III tumors. The vast majority of patients had a positive hormone receptor status (n = 206, 76.5%) and 23.5% had triple-negative subtype tumors. Twenty-eight patients (9.35%) showed a complete clinical response (CR), 139 (46.5%) a partial response (PR), 101 (33.8%) a stable disease (SD) and 131 (10.35%) disease progression (PD).

Study Population and Clinical Characteristics of Tumors
Out of the 342 patients included, 306 genomic DNA samples were available for gene polymorphism analysis. The patients' disposition is presented in Figure 1. The patients' characteristics, the disease and the treatment baselines are summarized in Table 1. The mean age at  inclusion was 55 years (range: 28-80). At initial diagnosis, the main histological type was invasive ductal carcinoma (79.5%). Most patients had received prior neoadjuvant/adjuvant chemotherapy (n = 210, 68%). Patients presented a small number of metastatic sites (90% had less than three metastatic sites); 62.5% of patients had a metastasis-free interval greater than 24 months. One hundred and thirty-five patients (48.5%) had histological Grade II tumors and 43.5% had Grade III tumors. The vast majority of patients had a positive hormone receptor status (n = 206, 76.5%) and 23.5% had triplenegative subtype tumors. Twenty-eight patients (9.35%) showed a complete clinical response (CR), 139 (46.5%) a partial response (PR), 101 (33.8%) a stable disease (SD) and 131 (10.35%) disease progression (PD).
The relationship between PFS and patient characteristics and 13 SNPs was assessed by univariate analyses (Tables S1 and S2). PFS was significantly associated (Table 2 and Figure S1 with the triple-negative subtype and histological Grade III tumors had a median PFS of 5.9 months (95% CI (4.50-7.85)), while patients with hormonal receptor-positive and histological Grade I or II tumors had a median PFS of 12.94 months (95% CI (11.83-14.82); p < 0.001). Table 2. Univariate analysis for the polymorphisms and clinical characteristics according to progression-free survival (PFS) and overall survival (OS).

Significant clinical characteristics
Histological grade

Associations among Overall Survival, SNPs and Patient Characteristics
Univariate analyses of OS were performed (Table S3 and Table S4). Among all tested variables, OS was significantly correlated (Table 3 and Figure S3 1-1.9); p = 0.01). As shown in Figure S4, triple-negative subtype and hormonal receptorpositive status were predictive of the worst and best prognosis, respectively. More precisely, triplenegative subtype patients had a median OS of 14.4 months (95% CI (12-21.7)) when having a histological Grade III tumor and a median OS of 15.7 months (95% CI (11.6-NA)) with histological Grade I or II tumors. In contrast, hormonal receptor-positive patients had a median OS of 38 months (95% CI (33.9-45.6)) when having a histological Grade I or II tumor and a median OS of 31.8 months (95% CI (25.3-41.7)) for histological Grade III patients.

Discussion
The present pharmacogenetic study was based on a homogeneous group of 306 patients with respect to cancer type and advanced disease status with an identical treatment protocol. An optimal predictive model for OS was obtained with the tumor subtype (triple-negative) and polymorphisms within the gene encoding the target of bevacizumab itself: VEGFA (rs833061 C/C), as well as within the genes encoding its receptors (VEGFR1; rs9582036 C/A or C/C and VEGFR2; rs1870377 T/T). The germline genome is sometimes considered to have limited applications regarding cancer therapy optimization [7]. However, the present study supports the idea that gene polymorphisms may help predict bevacizumab-based treatment outcomes. Bevacizumab is increasingly used outside its classical clinical area of prescription [12]. Moreover, the immunomodulatory properties of bevacizumab have recently led to very promising associations in combined immunotherapies for NSCLC and Pharmaceuticals 2020, 13, 414 9 of 12 hepatocellular carcinomas [2]. Thus, the present data may further encourage the identification and validation of reliable and much needed bevacizumab biomarkers. The clinical impact of VEGFA gene polymorphisms was recently examined in patients with advanced colorectal cancer treated with bevacizumab in association with chemotherapy. So far, the prognostic value of VEGFA rs699947 is unclear [13,14]. In contrast, VEGFA rs833061 demonstrated strong prognostic value in advanced colorectal cancer patients and, to a lesser extent, in breast cancer patients [13]. Notably, our data further support the prognostic value of VEGFA rs833061 for breast cancer patients.
Considering that rs833061 is located in the 5'-UTR of VEGFA, we then sought to investigate whether this specific polymorphism might have an impact on VEGFA expression, especially considering that this SNP is associated with nine SNPs (r 2 > 0.8) all located in the promoter region or in the first introns of VEGFA. Together, these SNPs were associated with the high expression of VEGFA in thyroid tissue or different alternative splicing events in the left heart ventricle vs. muscle (according to GTEX portal query results). These observations warrant further investigation to establish whether this SNP is functionally linked to the phenotype (poor outcomes) observed in the context of the present study.
rs9582036 is located on the FLT1 gene (also known as VEGFR1) and is part of a cluster of 15 SNPs with close linkage disequilibrium (0.6 < r 2 < 0.8). Fourteen of these SNPs are intronic, including the synonymous variant rs9582036 associated with an alternative splicing event in adipose tissue and the thyroid (according to the GTEX portal). All SNPs are located between the 25th and the 28th intron (30 exons in total). An in silico analysis did not show evidence of reliable variations. Nevertheless, the presence of a shorter isoform (ENST0000061 7835.4) that starts at Exon 28 and codes for the last 85 amino acids of the protein could be regulated differently due to the surrounding SNPs. Further phenotypic analyses would be required to understand the relationship between the SNP and the isoform (ENST00000617835.4). The rs1870377 T/A variant located on the KDR gene (also known as VEGFR2) introduces a missense glutamine (Gln) to histidine (His) substitution at the amino acid position 472, potentially impacting VEGFR2 degradation.
The present study contains a certain number of limitations, as we aimed to focus on a predefined set of SNPs. GWAS would be required in order to gain a broader insight into the SNPs involved in bevacizumab-based treatment outcomes. However, GWAS necessitates a larger number of patients. The impact in terms of OS conferred by the specific SNPs touching the targets of bevacizumab and its cellular receptors merits further investigation in the case where bevacizumab is associated with immunotherapy. In this respect, recent experimental data [15] highlighted a direct effect of VEGFA on endothelial receptors that modifies T cell diapedesis and allows TRegs rather than T cytotoxic cells to migrate to the tumor bed, thus decreasing the tumor's immune capacity. The SNPs described herein that have an impact on the outcome of bevacizumab-based treatment may thus be of potential prognostic value in the context of bevacizumab-checkpoint inhibitor combinations, a particularly promising domain of application for this widely used antiangiogenic and immunomodulatory agent.

Study Design and Patients
This prospective, multicenter COMET trial covered a cohort of 342 consecutive patients included between 2012 and 2015. The inclusion criteria were: age > 18 years, histologically confirmed metastatic HER2-breast cancer, ECOG ≤ 2, life expectancy ≥ 12 weeks, patients receiving first-line chemotherapy with bevacizumab (10 mg/kg every 2 weeks) and paclitaxel (90 mg/m 2 at Day 1, Day 8 and Day 15). Treatment was repeated every 4 weeks, according to routine practice, until disease progression or unacceptable toxicity. The exclusion criteria were prior chemotherapy for metastatic breast cancer, concomitant endocrine therapy or radiation therapy with curative intent for oligometastatic disease. The tumor response was defined according to Response Evaluation Criteria in Solid Tumours (RECIST 1.1) criteria as complete response (CR), partial response (PR), stable disease (SD) or progressive disease (PD). The objective response rate (ORR) was defined as the proportion of patients whose best overall Pharmaceuticals 2020, 13, 414 10 of 12 response over the entire treatment period was CR or PR. Informed written consent was obtained from each patient. The French Institutional Ethics Committee approved the study in June 2012 and it was registered (NCT01745757).

In Silico Analysis
The potential functional impact of SNPs was investigated using several in silico applications. HaploReg (http://www.broadinstitute.org/mammals/haploreg/haploreg.php) was used to collect all linkage disequilibrium and to examine the haplotype blocks within the genes. Regulome DB (http://www.regulomedb.org/) was accessed to explore the chromatin status, conservation and regulatory motif alterations within sets of genetically-linked variants. The GTEX Portal (https://gtexportal.org/home/) served to identify all cis-eQTL SNPs affecting the expression of genes of interest and as a microRNA binding-site prediction tool. SNPMIR (https://www.genomique.info/joomla/) was used to predict whether a SNP within the 3'-UTR of the genes of interest would disrupt/eliminate or enhance/create a microRNA binding site.

Statistical Analysis
An evaluation of the missing data rate was performed on 13 SNPs. All SNPs included in the analysis had a missing data rate of less than 10%. Missing genotypes were assigned using multiple imputations by chained equations (MICE) [18]. Dominant and recessive models were investigated to test possible associations between SNPs and progression-free survival (PFS) and overall survival (OS). For each SNP, risk alleles were coded as 1 and non-risk alleles as 0. For each SNP, the hazard ratio (HR) and 95% confidence intervals were calculated by Cox regression for associations between the genotype and PFS or OS. PFS was defined as the time between inclusion and the date of disease progression or death due to any cause. OS was defined as the time between inclusion and death due to any cause. Patients showing no event (death or progression) or lost to follow-up were censored at the date of last contact. PFS and OS were performed using the Kaplan-Meier method. Median follow-up with a 95% confidence interval was calculated with the reverse Kaplan-Meier method. For PFS and OS, we constructed a genetic and a clinico-pathological predictive model using a multivariable Cox regression method with backward elimination (p-value ≤ 0.05). The variables selected in each of the previous two predictive models were included in a third multivariable Cox regression with backward elimination to create a combined predictive model. For each model, the proportional hazard assumption was checked using statistical tests and graphical diagnostics based on the scaled Schoenfeld residuals [19]. For each multivariable model, a classifier predicting the risk of progression or death was Pharmaceuticals 2020, 13, 414 11 of 12 based on the linear predictor given by the model. The predictive ability of each model was compared using the area under the ROC curve (AUC). For PFS and OS, the clinical endpoint was fixed at t = 24 months and t = 60 months, respectively. Risk groups for predictive models were obtained using the best threshold value to obtain the highest sensitivity with a specificity of at least 80%. The sensitivity was defined as the proportion of patients experiencing progression or death before the time t attributed to the poor prognosis group. The specificity was defined as the proportion of event-free patients beyond the time t attributed to the good prognosis group. A p-value of <0.05 was considered statistically significant and all tests were two-sided. All statistical analyses were performed with R.3.5.2 software on Windows ® and the survMisc [20], timeROC [21], Survminer [22] and LDcorSV packages [23].

Conclusions
This study identified genetic polymorphisms related to a bevacizumab-paclitaxel treatment combination and established a genetic risk score in association with the clinico-pathological characteristics of patients. The highlighted risk groups identified patients with a poor or good prognosis and could be used to improve clinical decision-making.