Integrated Analysis of Cell-Free DNA and Novel Protein Biomarkers for Stratification and Therapy Monitoring in Stage IV Pancreatic Cancer

1. Introduction

Clinical management of patients diagnosed with pancreatic adenocarcinoma (PDAC) still remains very challenging due to limited therapeutic options and a lack of biomarkers to better guide systemic treatment and to predict clinical course of disease. Even today, most patients present with stage IV disease and succumb to it within one year of diagnosis [1,2]. Although the disease is perceived as uniformly aggressive, there is a great heterogeneity among patients in terms of response to treatment. While genomic studies have shown that pancreatic cancer is predominantly characterized by mutations in the four genes KRAS, CDKN2A, TP53 and SMAD4, it remains unclear how clinical heterogeneity arises when most tumors are caused by the same mutation pathway [3,4,5]. Lack of predictive and prognostic biomarkers and impossibility of detecting intrinsic resistance of the disease to molecularly targeted therapy make it difficult to develop personalized treatment strategies.

Even though classical protein-based tumor markers (CA 19-9, CEA) have been used for years to monitor clinical course of disease in patients with known cancers, they are not suitable for screening purposes due to their sometimes low sensitivity or specificity [6,7,8,9,10]. A promising novel tool that has become the focus of the biomarker field in recent years is the detection of circulating nucleic acid fragments, in particular circulating tumor DNA (ctDNA), in the bloodstream of cancer patients [11,12,13,14]. This non-invasive blood-based technology can not only reflect tumor changes under therapy in real time [15], but also allows statements to be made regarding therapy response through longitudinal monitoring during the course of disease [15,16,17,18,19,20,21,22,23,24]. Hereby, highly sensitive detection methods such as digital droplet PCR allow sensitivities of 0.01 to 0.001% [25,26].

An alternative concept is the establishment of protein expression patterns from serum or plasma that may have the potential to map tumor-specific signatures from a few milliliters of blood [27,28]. A general observation, however, is that each individual method, whether protein-, DNA- or RNA-based, has its specific advantages and limitations, and possibly only the combination and integration of different methodological approaches can meet all clinical requirements. The biomarker field is thus moving towards a combination of complementary approaches [27,29]. Therefore, the aim of this exploratory biomarker study was to investigate the potential of integrated biomarker analysis for clinical application in PDAC through the combined highly sensitive detection of cell-free mutated KRAS (KRAS^mut cfDNA) and a protein biomarker signature that we have already been able to establish in a systems biology approach [30].

2. Materials and Methods

2.1. Human Samples and Patient Cohort

Blood samples were collected from patients with clinically and histologically confirmed metastatic PDAC at time of diagnosis and during routine clinical follow-up, in addition to clinic-pathologic, treatment and outcome data. Diagnosis of metastatic disease was based on a radiological and cyto-/histological confirmation. A total of 47 patients undergoing first-line systemic treatment or best supportive care following diagnosis of metastatic PDAC were recruited at the University Medical Center in Freiburg. All patients gave written informed consent to collection and analysis of blood samples. Local institutional review board (IRB) approved all relevant procedures and analyses (EK-Freiburg project number 46/18). Treatment was performed as per standard of care and blinded to KRASmut cfDNA results. Follow up included 3-monthly clinical and radiological examination using computed tomography and/or magnetic resonance imaging. Progressive disease (PD) was determined based on routine evaluation of radiological imaging.

2.2. Extraction of Cell-Free DNA from Plasma Samples

Venous blood samples from mPDAC patients were collected using commercially available EDTA tubes. Plasma was extracted through two subsequent centrifugation steps at 3000 rpm and 14.000 rpm for each 10 min at 4°C within one hour of collection, as previously described [31], and frozen at -80°C without further treatment until extraction of cell-free DNA. CfDNA was extracted from each 4ml plasma sample following the SEP/SBS protocol of the PME-free circulating DNA extraction kit (Analytik Jena, cat. no. 845-IR-0003050) according to manufacturer’s instructions. Two subsequent elution steps with each 30 µl Elution Buffer were performed. DNA was stored at -20°C until ctDNA quantification. CfDNA was evaluated with fragment analyzer and quantified using Qubit 2.0 fluorometer. DNA yield from 4ml of plasma typically ranged from 1-20 ng/µl up to 80 ng/µl for metastatic pancreatic cancer patients.

2.3. Enzyme-linked Immunosorbent Assay (ELISA)

ELISA antibodies and preparation kits for SFN and TFF1 were manufactured by Wuhan USCN Business (USCN, Co. Ltd., Wuhan China; cat. no. #SEH179Hu, #SEB049Hu) and for CEMIP, COL10A1, HGF, LAMB3, POSTN, SERPINB5, SPP1 and TMPRSS4 by LifeSpan Biosciences (LSBio, Seattle, USA; cat. no. #LS-F7390, # LS-F13131, #LS-F2441, #LS-F11919, #LS-F3645, #LS-F13455, #LS-F13367, #LS-F7623). Plasma was frozen at -80°C before preparation according to manufacturers’ instructions. The colorimetric reactions were analyzed on a Tecan Infinite M200 Pro.

2.4. Droplet Digital PCR (ddPCR)

Locked nucleic acid (LNA)-based probes and associated primer pairs for detecting KRAS mutations in plasma cfDNA were developed using Beacon Designer v.8.20 soft-ware (Premier Biosoft, Palo Alto, CA, USA). These primers and probes were custom-synthesized by Integrated DNA Technologies (IDT, Inc., Coralville, IA, USA). Comprehensive details regarding the design process, including sequences, have been previously published [31]. The reaction mixture for droplet digital PCR (ddPCR) was prepared by combining ddPCR Supermix for Probes (Bio-Rad, catalog #186-3024) with primers, probes, 2 µl of template DNA per well, and nuclease-free water (Ambion). Droplets were created us-ing the QX100/200™ Droplet Generator (Bio-Rad, catalog #1863002) following the manufacturer’s protocol. Each sample was analyzed in four replicates. After generating the droplets, the mixture was transferred to a 96-well PCR plate (Bio-Rad, catalog #12001925) and PCR was run on a C1000 Touch™ Thermal Cycler (Bio-Rad, catalog #1851197). The amplified droplets were subsequently analyzed using the QX100/200™ Droplet Reader (Bio-Rad, catalog #1863003) with the QuantaSoft software v1.7.4.0917 (Bio-Rad, catalog #1864011). Details of the PCR protocols, assay controls, and data analysis methods for the KRAS assays have been previously published [31]. The number of copies per milliliter of plasma was determined using the following calculation:

copies/ml plasma = $-\ln ({\displaystyle \frac{N_{neg}}{N}})/V_{droplet}$ x 20 µL x 7,5 [31]

2.5. Statistical Analysis

The primary outcome measured was PFS based on routine evaluation of radiological imaging. PFS was defined from start of first-line therapy to the verified first radiologic progression based on standard restaging imaging or death due to any cause. For this purpose, patients were divided into two groups based on their radiological response: patients with complete or partial remission and patients with stable disease were assigned to the "non-PD" group and compared to patients with "PD". OS was the secondary end point and defined as time from the date of diagnosis till death. Kaplan–Meier survival analysis was performed to estimate progression-free and survival time. Univariate analyses were carried out using the log-rank test. Backward stepwise Cox regression modeling to estimate hazard ratio (HR) with 95% confidence interval (CI) was used to explore independent prognostic factors for PFS and OS. Fisher’s exact test and the Wilcoxon-Mann–Whitney test were carried out to compare independent variables. All statistical analyses were performed using GraphPad Prism Version 9.4.0 (GraphPad Software, Inc., La Jolla, California, USA) and SPSS 26 software Version 26.0.0.0 (IBM Corporation, New York, United States). P values < 0.05 were considered significant.

2.6. Risk Stratification and Classification

The protein levels were determined for CEMIP, COL10A1, HGF, LAMB3, POSTN, SERPINB5, SFN, SPP1, TFF1 and TMPRSS4. In order to illustrate the concentration level of each of the markers, the level of the markers were scaled, independently. In the case of COL10A1, TMPRSS4, CEMIP and KRAS^mut cfDNA, six levels with drastically higher values in comparison to the other outliers, were excluded to make the general trend of the biomarkers observable. To consider the effect of the therapy in respect to OS and the kinetics of protein levels, a binary variable was provided for each protein as follows: The binary variable represents whether the level of that protein increases or decreases after therapy start for the corresponding patient. For increasing protein level “1” is used and “-1” for decreasing one. For the patients with more than one follow up sample, we decided for the increasing or decreasing trend that most of the samples showed. In case of ties, the first follow up’s binary value was used. We performed the same procedure on the KRAS^mut cfDNA levels. Each univariate model stratified the patients based on the binary values, independently and was evaluated using log-rank test.

In order to identify a unified signature of the markers using a multivariate CoxPH model [32], the R package GLMnet was used using LASSO and 3-fold cross-validation based on the combination of the binary variables [33,34]. Risk labels were computed using the estimated model. The prognostic potential of the gained model was evaluated using log-rank test.

Using the labels of low and high risk of the patients of the PT-and-UT group provided by the multivariate CoxPH model, a Risk Classifier was trained to classify the patients to the corresponding risk groups. The learning process used the protein levels of the selected biomarkers of the prior to treatment condition as the feature space and the risk labels of the PT-and-UT patients as the labels. Binomial Logistic Regression [33] was applied for learning the classifier with alpha = 0.6 for the setting of the elastic net and 4-fold cross validation. Higher absolute values of the coefficients increases the impact of the corresponding variable on the Risk Classifier.

The Risk Classifier predicted then the risk label of each of PT-and-UT and PT-only patients. Three patients in the PT-and-UT group were misclassified (error rate 16%). In the case of PT-only group, the predicted risk labels were used to provide a Kaplan Meier plot, where the stratification significance was evaluated using log-rank tests.

3. Results

3.1. Patient Characteristics

47 patients with clinically and histologically confirmed mPDAC were included in the study. Clinico-pathological parameters of these stage IV PDAC patients have been published earlier [31]. Additional basic clinical characteristics of the study cohort are summarized in Supplemental Table 1 and 2. 42/47 patients (89.4%) received first-line chemotherapy, 5/47 (10.6%) underwent best supportive care (BSC). Median follow-up time among surviving patients was 9.5 months (95% CI 3.0 – 20.0 months). Median progression-free survival (PFS) was 3.0 months (95% CI 2.0 – 6.0 months), median overall survival (OS) was 9.0 months (5.0 – 13.0 months). At the time of final analysis 41/47 of the patients (87.2%) had died.

3.2. Analysis of Plasma KRAS^mut cfDNA

A total of 160 blood samples from 47 patients with mPDAC were analyzed in this study using ddPCR. First samples were taken at a median of 7.5 days (95% CI 4-11) prior to start of first-line treatment. Median number of samples collected was three samples per patient (95% CI 2-4). Median time interval between start of therapy and first follow-up was 55 days (95% CI 50 – 62). Mutant KRAS detection rate at baseline was 70.21% (33/47). Thereby, a correlation between KRAS^mut cfDNA positivity and the presence of liver metastases was shown (P = 0.0015, Supplemental Table 3). PDAC patients with liver metastases also presented significantly higher KRAS^mut cfDNA copies/ml plasma than patients with metastases located elsewhere (P = 0.0003; Supplemental Figure 1A). For the tumor marker CA 19-9 this difference was less pronounced (P = 0.0178; Supplemental Figure 1C). Regarding tumor grading, there was no difference between well and poorly differentiated tumors and the amount of KRAS^mut cfDNA or CA 19-9 levels (Supplemental Figures 1B, D). Although KRAS^mut cfDNA positive patients had higher levels of the tumor marker CA 19-9 (P = 0.0054; Supplemental Figure 1E), there was generally no correlation between the CA 19-9 level and the amount of KRAS^mut cfDNA in plasma (Supplemental Figure 1F).

3.3. Univariate and Multivariate Analyses of PFS and OS in mPDAC Patients

For the univariate analysis of OS in our cohort of mPDAC patients, 11 independent demographic and clinic-pathologic features were examined, which are presented in Table 1. Four independent variables – namely, tumor differentiation, the number of metastatic sites, administration of systemic treatment and change in levels of KRAS^mut cfDNA during chemotherapy were identified as prognostic factors. The multivariate Cox proportional hazards regression model showed that KRAS^mut cfDNA increase at time of the first follow-up (Hazard ratio (HR) = 10.9, 95% CI: 2.589-46.17, P = 0.001) and tumor differentiation (HR = 3.17, 95% CI: 1.175-8.535, P = 0.023) were the only significant factors for survival in our cohort of mPDAC patients (Table 1). Corresponding analyses were also performed for PFS and are presented in Supplemental Table 4. Again, the multivariate Cox proportional hazards regression model showed that KRAS^mut cfDNA increase at time of the first follow-up (HR = 10.9, 95% CI: 2.575-46.444, P = 0.001) was a significant factor for PFS. Furthermore, the number of metastatic sites (HR = 7.20, 95% CI: 1.149-45.080, P = 0.035), higher Age (HR = 4.39, 95% CI: 1.250-15.42, P = 0.021), tumor location (HR = 0.20, 95% CI: 0.053-0.721, P = 0.014) and CA 19-9 positivity in plasma before the start of systemic treatment (HR = 4.38, 95% CI: 1.180-16.285, P = 0.027) were also shown to be significant variables for PFS.

3.4. Predictive and Prognostic Value of KRAS^mut cfDNA and CA 19-9

To evaluate the predictive and prognostic relevance of KRAS^mut cfDNA and CA 19-9 levels, baseline values of these biomarkers before start of palliative first-line treatment were correlated with PFS and OS. Thereby, high KRAS^mut cfDNA levels (> 25 copies/mL) before start of systemic chemotherapy were inversely associated with PFS (HR = 2.07, 95% CI: 1.011-4.229, P = 0.0057; Figure 1A) and showed also a trend towards inferior OS (HR = 1.72, 95% CI: 0.850-3.489, P = 0.0775; Figure 1B). In addition, detection of a high amount of total cfDNA in plasma was associated with shorter PFS in our study cohort (HR 2.56, 95% CI: 1.217-5.367, P = 0.0002; P < 0.0001; Supplemental Figure 2A), but had no prognostic value (Supplemental Figure 2B). CtDNA positivity before palliative first-line chemotherapy was not associated with inferior PFS and OS in our study cohort (Supplemental Figure 2C, D). Analysis of the tumor marker CA 19-9 prior to therapy also had no predictive or prognostic value in our cohort (Figure 1C, D). A threshold value of 100 U/ml was chosen because previous data showed that this value is not found in healthy individuals without a clinical history of cancer disease [27,35].

3.5. Association of KRAS^mut cfDNA and CA 19-9 Dynamics with Survival

Protein tumor markers and KRAS^mut cfDNA are highly dynamic biomarkers. Therefore, in a next step, a second blood sample was analyzed at time of the first follow-up (median after 55 days) and biomarker changes in this time interval were correlated with PFS and OS. For 25/47 (53.2%) patients a sufficient number of follow-up samples was available for subsequent analysis. Increase of KRAS^mut cfDNA hereby was associated with significantly reduced PFS and OS (Figure 2A, B), while increase of the total amount of cell-free DNA did not correlate with PFS or OS (Supplemental Figure 2E, F). Increase of the tumor marker CA 19-9 was associated with a non-significant trend towards inferior OS but not PFS (Figure 2C, D). Furthermore, the radiological imaging based division into progressive disease (PD) and non-PD at time of first restaging did not show any difference regarding OS (Figure 2G). Interestingly, patients with PD at time of first follow-up already showed elevated KRAS^mut cfDNA levels in plasma prior to first-line therapy (Figure 1E). Also at time of first follow-up, they showed significantly higher KRAS^mut cfDNA levels in plasma (Figure 2E). On the contrary, no difference was found for the tumor marker CA 19-9 (Figure 1F, 2F). Corresponding to this, 22/23 patients (95.65 %) who presented with PD at time of first follow-up were positive for KRAS^mut cfDNA before the start of first-line therapy, whereas 7/11 patients (63.64 %) with stable disease were KRAS^mut cfDNA negative (Figure 1G). For the tumor marker CA 19-9, no correlation was found (Figure 1H).

3.6. Clinical Response Prediction by Kinetics of KRAS^mut cfDNA and CA 19-9

Figures 3A and 3B illustrate the correlation between KRAS^mut cfDNA and CA 19-9 dynamics during first-line therapy and radiologic response to treatment at time of first follow-up. 13/25 patients (52 %) showed PD at the first restaging. Increase of KRAS^mut cfDNA was hereby significantly associated with PD and outperformed CA 19-9 as dynamic marker (Figure 3C). Overall, kinetics of KRAS^mut cfDNA analyses nicely reflected individual patients’ course of disease, as illustrated in 6/6 exemplary shown cases (Supplemental Figure 3). Individual patient analyses showed that KRAS^mut cfDNA is a highly dynamic biomarker, which can be clinically especially relevant for patients with non-elevated CA 19-9 levels. However, prospective studies and larger cohorts are needed to better unravel the temporal relationship between biomarker dynamics and clinical disease course and to investigate the impact on therapeutic strategy.

3.7. Risk Stratification Based on the Kinetics of KRAS^mut cfDNA and Biomarker Proteins

Out of the total 47 patients in this study, only 19 cases had enough plasma left for ELISA measurements prior to and under therapy (PT-and-UT group). For a set of 17 patients plasma samples prior to therapy were available (PT-only group) and therefore used as the validation set. In order to improve the prognostic impact of KRAS^mut cfDNA in mPDAC, the role of ten proteins in plasma (CEMIP, COL10A1, HGF, LAMB3, POSTN, SERPINB5, SFN, SPP1, TFF1 and TMPRSS4) was investigated by considering the effect of the first-line therapy. The related genes were previously described as a subset of relevant biomarkers in distinguishing PDAC tissue from normal tissue by RNA sequencing and showed significantly elevated levels in plasma samples of PDAC patients [30]. To examine the prognostic role of these proteins in mPDAC, we referred to the protein levels measured via Enzyme-linked Immunosorbent Assay (ELISA) of two subgroups of patients (scaled values in Supplemental Figure 4D).

As a first step of our analysis, we examined the single biomarker impact on the OS of the prior to and under treatment mPDAC patients (PT-and-UT group) using univariate HR model. Most of the single biomarkers were not able to stratify patients significantly (Supplemental Table 5). However, the prognostic impact of the single proteins of LAMB3 and SERPINB5 was significant, although the corresponding significance was not high (P = 0.044; P = 0.037, labelled with a star, Supplemental Table 5, Supplemental Figure 4A and 4B). The prognostic impact of KRAS^mut cfDNA is the highest among all of the evaluated markers (P = 0.0042, Supplemental Figure 4C).

Based on the above results, we used a multivariate Cox Proportional Hazard model (CoxPH) model to identify the prognostic biomarkers together with KRAS^mut cfDNA based on the PT-and-UT group [32,33]. Similar to the previous univariate models, we use the binary variables to consider the kinetics of each protein and the effect of the therapy. However, the binary variables are used together to identify a unified signature of the biomarker proteins and KRAS^mut cfDNA (Figure 4A). Figure 4B shows the result of the multivariate CoxPH to stratify patients to high or low risk groups (HR = 10.2, 95% CI: 2.03-51.7, P = 0.0014). The significance of the multivariate CoxPH model was improved in comparison to each of the univariant ones. The selected set of prognostic biomarkers identified using LASSO together with related impacts were given in Supplemental Table 6. In detail, the selected set of biomarkers by CoxPH includes KRAS^mut cfDNA and six other proteins including CEMIP, TFF1, LAMB3, HGF, TMPRSS4 and SERPINB5. CEMIP was the only significant biomarker and the impact was higher as for KRAS^mut cfDNA.

To predict the therapy response in respect to OS, we used the selected set of markers to classify eligible patients. Patients predicted in the low risk group benefit more likely of further therapy. Based on the selected markers and by using the prior to therapy samples of the PT-and-UT group together with the identified label of low or high risk given by the multivariate CoxPH, a Risk Classifier was provided to predict the label of risk (Supplemental Table 7). The accuracy of the Risk Classifier for the PT-and-UT group is 84%. The actual risk label of each PT-and-UT patient and the predicted one by the Risk Classifier is given in the Supplemental Table 8.

The Risk Classifier was applied to recognize the risk label of the PT-only group to predict the therapy response in respect to OS. The OS Kaplan Meier plot of Fig 4C based on the predicted risk labels by the Risk Classifier presents an acceptable separation within the PT-only group, although not significant (P = 0.17). We believe that a higher number of patients especially with under therapy samples are necessary to validate the Risk Classifier. However, the Risk Classifier showed a significant prognostic performance for all available patients’ data of prior to therapy samples (P = 0.043, Figure 4D).

4. Discussion

Within this single-center exploratory study, we followed a cohort of 47 mPDAC patients undergoing palliative first-line treatment or BSC at our institution. We analyzed KRAS^mut cfDNA using previously established highly sensitive single-target and discriminatory multi-target ddPCR assays [31]. Results were integrated with CA19-9 levels and ten protein biomarkers markers (CEMIP, TFF1, LAMB3, HGF, TMPRSS4, SERPINB5, SPP1, POSTN, TFF1, COL10A1) for association with progression-free and overall survival endpoints. Over the past decades numerous studies have demonstrated the predictive and prognostic significance of ctDNA for closer disease monitoring in pancreatic cancer [19,20,21,23,24,36,37,38,39,40]. To our best knowledge, even today there are only limited data on the serial analysis of ctDNA in advanced or mPDAC available [20,21,22,24]. What takes our study apart is the combined longitudinal analysis of KRAS^mut cfDNA and novel protein biomarkers overcoming limitations of single-parameter based methods, enabling risk stratification, and molecular monitoring of patients with mPDAC. Our study highlights the importance of KRAS^mut cfDNA serial monitoring as a sensitive and highly specific biomarker for therapy response and disease progression. The use of discriminatory multi-target ddPCR assays herby allows for direct KRAS SNV detection without prior tumor NGS making it possible to assess tumor dynamics in a time relevant co-clinical setting analogous to CA19-9 levels [31,41]. Increase of KRAS^mut cfDNA during first-line palliative treatment in our study cohort was significantly associated with shorter progression-free and overall survival and outperformed CA 19-9 kinetics and standard imaging. In particular, for patients not expressing CA 19-9 or lacking dynamic changes of the tumor marker during systemic treatment analysis of ctDNA in plasma might provide additional diagnostic benefit. However, each individual method has its specific advantages and limitations and possibly only the combination and integration of different methodological approaches can meet all clinical requirements. For our cohort integrated analysis of KRAS^mut cfDNA combined with six evaluated circulating protein biomarkers, allowed basal risk stratification at time of first follow-up. Yet, routine clinical assessment of dynamic ctDNA changes and protein expression patterns during systemic treatment of mPDAC using a standardized approach is lacking.

Still today systemic treatment options for pancreatic cancer are limited to a small number of combination chemotherapy regimens [42,43] and new attempts in personalized treatment based on molecular profiling [44,45]. Given the fact that according to the current literature only approximately 50% of PDAC patients receive further line treatment, early switch of chemotherapy in case of lack of response is crucial [42,46]. Thus, we advocate serial ctDNA analysis to assess response dynamics within a time-relevant framework. KRAS^mut cfDNA guided therapy monitoring in pancreatic cancer appears feasible, but still requires extensive validation in prospective clinical trials to investigate whether treatment decisions based on serial ctDNA assessment results in improved outcome of PDAC patients.

Nevertheless, there are also some limitations of our study. Subgroup analyses of the present study were performed with a relatively small number of metastasized PDAC patients and survival outcomes were correlated retrospectively. Thus, results are of an exploratory nature and only hypothesis generating. In contrast to previously published data [22,24,47,48,49], KRAS^mut cfDNA positivity prior to palliative first-line therapy was not associated with shorter PFS and OS in our cohort. However, high KRAS^mut cfDNA levels (> 25 copies/mL) before start of systemic chemotherapy were associated with significantly reduced PFS and showed a trend towards inferior OS, highlighting the importance of identifying clinically validated cut-offs for ctDNA analysis [18,50,51]. Interestingly, Del Re et al. [20] also found no association between KRAS^mut cfDNA positivity prior to palliative first-line therapy and survival endpoints. Also in this study highly sensitive ddPCR was used to detect ctDNA in plasma and KRAS^mut cfDNA detection rate of 70.4% at baseline was identical to the detection rate in our study, whereas in previously published studies detectability ranged from 0% to 62.2% [48,49,52]. The significantly higher positivity rate could be another possible explanation for the lack of association with survival endpoints.

In summary, our pilot study proposed a clinically feasible way to assay KRAS^mut cfDNA in plasma together with novel protein biomarkers in mPDAC patients and outperformed CA 19-9. Through combination of those markers, patients could be better stratified in terms of treatment response and overall prognosis.