A Combined Strategy of SAGE and Quantitative PCR Provides a 13-Gene Signature that Predicts Preoperative Chemoradiotherapy Response and Outcome in Rectal Cancer Free

Preoperative chemoradiotherapy (CRT) is the standard treatment for stage II/III rectal cancer (RC). Compared with postoperative therapy, it improved local control and tolerability in the CAO/AROAIO 94 trial (1). More recently, the NSABP R-03 trial has reported a significant benefit in 5-year disease-free survival (DFS; ref. 2). Unfortunately, 5-fluorouracil (5FU)-based CRT produces a complete pathologic response in only 8% to 14% of the patients, and distant metastases develop in one-third of the cases after 5 years (1, 3, 4). Furthermore, patients who do not respond to preoperative CRT are exposed to unnecessary toxicities and resection of their primary tumor is delayed. Thus, the ability to predict treatment response and outcome before CRT is of the utmost importance to clinicians.

A number of postsurgical prognostic factors have been proposed (5–7), but CRT response cannot be clinically predicted. The identification of basal markers of resistance biomarkers could offer great help on this regard. Directed strategies that explore individual markers have not yielded clinically validated assays (7–9). Instead, gene expression arrays are appealing to predict the multifactorial event of treatment response (10–13). TaqMan Low Density Arrays (TLDA) can be done from routine formalin-fixed paraffin-embedded (FFPE) biopsies, making this platform of great interest for clinical application, as shown for the Oncotype Dx Assay (Genomic Health Inc.) signature in breast cancer (14), endorsed by ASCO (15) and NCCN guidelines (16). Selecting appropriate genes for analysis seems critical.

In this study, we sought to define a clinically exploitable gene signature that could predict treatment response and outcome in preoperatively treated RC patients. First, we carried out an undirected screening of candidate genes by evaluating dynamic changes in the RC transcriptome, before and after CRT in selected treatment-refractory patients. Serial analysis of gene expression (SAGE) was ideal because it does not require a priori knowledge of the transcripts that are present in cells and it accurately quantitates mRNA levels of the transcriptome, allowing the discovery of new molecular targets (17). It provided a set of 24 genes whose expression was altered after CRT. In the second stage, we evaluated 29 genes selected from the literature and the SAGE-emerging candidates (for a total of 53 genes) by using TLDA cards in FFPE initial biopsies from 94 locally advanced preoperatively treated RC patients. Basal changes in gene expression were correlated with tumor regression grade (TRG) and patient outcome (defined as 3-year disease free survival). A set of 13 genes was identified that could predict both endpoints.

The experimental design is depicted in Figure 1.

For the initial SAGE experiment, 25 patients from a RC neoadjuvant CRT trial were prospectively recruited (18). Chemotherapy consisted of a preoperative combination of oxaliplatin and raltitrexed (130 mg/m² and 3 mg/m², days 1, 21, and 42) on 3 cycles, together with preoperative radiotherapy (RT; 50.4 Gy in 28 fractions). We sought to evaluate dynamic changes in tumor transcriptome after CRT in refractory patients, as defined by lack of downstaging and minimal or null pathologic response according to Dworak's classification of TRG (TRG 0–1; ref. 19), with similar initial clinical stage (cT3N0, as diagnosed by using endorectal ultrasounds and MRI) and adequate fresh tissue samples. Successful libraries were obtained in 3 cases (patients A, B, and C) fulfilling the criteria. Patients A and C developed lung metastases 16 and 54 months after diagnosis, respectively, supporting the scenario of resistant disease.

For the quantitative PCR (qPCR) study, 125 additional patients were retrospectively recruited. The patients had clinical stages II and III adenocarcinoma (defined by endoscopic ultrasound and/or MRI) of the rectum within 10 cms from the anal verge and were treated with 2 or more preoperative cycles of chemotherapy, concomitant with 50.4 Gy in 28 fractions. The specified clinical and pathologic variables were recorded (Table 1). Thirty-one patients were excluded because of unsuitable pathologic samples or incomplete gene expression or clinical data. Ninety-four patients were selected for the response prediction statistical analysis. Patients received 4 different chemotherapy regimens as follows: (a) 29 UFT (300 mg/m² on days 1–14 concurrent with radiotherapy) and leucovorin (250 mg/m² i.v. on day 1 and 7.5 mg/m²/12 h p.o. on days 2–14; ref. 20; b) 17 5FU 225 mg/m²/d in continuous infusion during RT; 48 patients were treated with concurrent oxaliplatin-based CRT, including 27 (c) receiving raltitrexed and oxaliplatin (3 mg/m² and 130 mg/m², respectively, on days 1, 21, and 42), and 21 (d) 2 cycles of standard XELOX followed capecitabine (825 mg/m²/12 h) on radiotherapy days. Follow-up occurred at 3-month intervals for 2 years, 6-month intervals the third year, and yearly thereafter. Evaluations consisted of anamnesis, physical examination, and blood tests. Complete colonoscopy, chest radiography, abdominal ultrasound, and abdominopelvic CT were also scheduled according to local guidelines. Local and distance recurrence was histologically confirmed whenever possible. Alternate criteria included sequential enlargement of a mass or distant nodules in radiologic studies.

Institutional approval from our ethical committee was obtained for the conduct of the study. Patients provided written consent so that their samples and clinical data could be used for investigational purposes.

Two independent expert surgical pathologists who were blinded to patient outcome evaluated TRG according to the classification of Dworak, as described by Rodel and colleagues (19, 21). Response to CRT was considered as TRG greater than 0–1. The kappa index was used to evaluate concordance between the 2 observers and the pathologists discussed discordant results.

To identify genes involved in RC treatment response, a previous exploratory SAGE study was done. SAGE libraries were generated from initial tumor biopsies before CRT and postoperative specimens from patients A, B, and C. Tumor specimens were snap frozen in liquid nitrogen at −70°C and stained with hematoxylin and eosin (H&E). Eligible samples were composed of at least 90% of tumor cells. SAGE was done as described in the Supplementary Material. The pairwise comparisons (before and after treatment for each patient) were determined by using Monte Carlo simulation to identify transcripts in which difference was statistically different at P ≤ 0.01 and P ≤ 0.05 (Supplementary Tables S1 and S2). Gene matches for significant tags were manually verified by using SAGE Genie (22). Eligibility of SAGE genes for qPCR analysis was based on statistically significant altered expression, similar regulation trend (up/downregulation) across the 3 patients, and biological plausibility according to the published literature.

FFPE tumor sections from the initial biopsies of the 125 patients were reviewed by an expert pathologist. Ninety-four patients were selected after excluding cases with insufficient tumor cells, incomplete clinical data, or unsuccessful amplification results. More than 75% tumor cells enrichment was ensured, when necessary, by subsequent macrodissection with the use of a safety blade and new confirmatory H&E staining in postoperative biopsies. RNA was extracted from 10 to 15 five-micrometer sections by using the Masterpure RNA Purification Kit (Epicentre Biotechnologies) according to the manufacturer's instructions, detailed in Supplementary Material.

qPCR was done by using the ABI Prims 7900 HT Sequence Detection System (Applied Biosystems), according to manufacturer's instructions. TLDAs design was done as described in the Supplementary Material.

Testing for associations between response and clinical variables was done with χ² test (Yates corrected) or Fisher's exact tests for categorical variables.

The expression of each gene was measured in triplicate, and the average intensity in the microarray was assessed. To reduce variation between microarrays, the intensity values for each sample in each microarray were rescaled by means of a quantile normalization method. Each intensity value was then log-transformed to a base-2 scale. qPCR data results were normalized against the 3 most stable genes (GAPDH, B2M, and PMSB4) before analysis. The stability of the expression of the candidate control genes was determined with genNorm (23).

The total sample size was randomly assigned to either a training cohort (60%) or a testing cohort (40%; Supplementary Fig. S1). Patients were classified into 3 groups (TRG 0–1, TRG2–3, and TRG4). We used a supervised class prediction by stepwise forwards and backwards discriminant analysis method to define a subset of genes predictive of response. This model was applied to the testing cohort to validate its usefulness. Accuracy, sensitivity, specificity, positive predictive values (PPV), and negative predictive values (NPV) were calculated (more detailed in Supplementary Materials).

Two-factor analyses were calculated (F₁ and F₂, the discriminant functions) for the genes that were selected by discriminant analysis. A score was obtained for every patient with the formula:

To determine differences in DFS and overall survival (OS; measured from the time of diagnosis), univariate analyses were done according to the Kaplan–Meier method. Comparisons between curves were evaluated with the log-rank test. Cox regression was used for the multivariate analysis. A value of P < 0.05 was considered statistically significant, and all tests were 2-tailed. The proportional hazards assumption for covariates was investigated by examining the scaled Schoenfeld residuals. SPSS 17.0 and SLPUS v. 6. were used for statistical analysis.

Six SAGE-libraries were generated from patients A, B, and C. Every patient was studied before and after treatment. Pairwise comparisons revealed statistically significant dynamic changes in gene expression when libraries were compared before and after treatment for each patient (P < 0.01). The numbers of tags, unique tags for each library, and dynamic changes are shown in Supplementary Tables S1 and S2. Certain tag-corresponding genes significantly altered after treatment included tissue specific genes such as the trefoil-family factor 3 (TFF3) and the gastrointestinal glutathione peroxidase 2 (GPX2) or transcripts associated with treatment response as glutathione S-transferase pi or with fibrotic response after injury, such as collagen type I alpha 1 (COL1A1) and secreted protein acidic cysteine-rich (SPARC). For qPCR analysis, we focused on 24 genes that were significantly regulated in pairwise comparisons (Supplementary Table S2) in addition to 29 genes previously reported to be involved in colorectal carcinogenesis or treatment response (Supplementary Table S3).

A total of 94 patients were analyzed (77 patients from Hospital La Paz, Madrid and 17 from Hospital Clinic of Barcelona). A flow diagram of the patients is depicted in Figure 1. Clinicopathologic data are shown in Table 1. After a median follow-up of 46 months (range 6–127 months), local relapse occurred in 6 patients (7%) and distant metastases occurred in 22 (23%). The 3-year DFS was 73% (95% CI: 64–82) and 3-year OS was 90% (95% CI: 84–96). TRG 0–1, 2, 3, and 4 were seen in 18, 53, 11, and 12 patients, respectively, with good concordance between pathologists (Kappa = 0.7).

From the 53 genes, 13 were selected by discriminant analysis to build a predictor of response (Table 2). In the training cohort, the model correctly classified 86% of the cases (49 of 57). Two nonresponsive patients (TRG 0/1) and 6 responsive cases (TRG 2, 3, and 4) were misclassified. The test correctly classified 87% of responders (sensitivity) and 82% of nonresponders (specificity). The PPV was 95% and the NPV was 60%. In the testing cohort, the model correctly classified 6 of 7 nonresponders (specificity = 86%) and 22 of 30 responders (sensitivity = 73%), with an overall accuracy of 76% (Supplementary Fig. S1) and a PPV of 95%.

The 13-gene score was applied to the entire cohort of patients and was inversely correlated with DFS. Prognostic groups were defined as low (>2) and high risk (≤2) (Figure 2). Clinical variables, except for basal carcinoembryonic antigen level (CEA), did not differ between risks groups. On univariate analysis, the risk score was significantly correlated with DFS along with pathologic variables [AJCC Stage, ypT (post–CRT-yatrogenic pathologic T) and ypN (post–CRT-yatrogenic pathologic N; Table 3)]. Given the heterogeneity of the sample, we investigated whether the type of chemotherapy received impacted the scoring results. When stratifying by treatment, the score remained statistically significant (P = 0.03 and P = 0.01 for fluoropyrimidines alone and oxaliplatin-based CRT, respectively; Supplementary Fig. S2). Univariate stratified analysis for other factors, such as ypT and ypN, confirmed the independent prognostic value of the score for DFS (Supplementary Figs. S3 and S4). Interestingly, the score was such a strong predictor of DFS that it remained the only statistically significant factor when Cox multivariate analysis was used, adjusted for age, ypT, ypN, CEA, and tumor response (HR = 4.35 95% CI: 1.20–15.75; Table 4).

Treatment and tumor-related factors determined after surgery, primarily CRM and nodal status, have been associated with RC recurrence (5). However, individualized treatment selection for patients with RC relies on pretreatment clinical variables that do not predict the efficacy of neoadjuvant therapy. This study used a sequential genomic strategy to identify a basal 13-gene signature that correctly predicted RC preoperative CRT response in 86% of patients in a training set and 6 of 7 nonresponders in a confirmatory testing cohort. Also, it predicted 2 groups of patients with clear differences in DFS.

This predictor model offers some advantages over other published tumor gene signatures that evaluate response to preoperative RT (24) or combined CRT (10, 12, 13, 25) in terms of clinical feasibility, relevance of the surrogate marker, and prognostic ability. This is the only signature obtained from small amounts of routinely prepared FFPE endoscopic samples by using a qPCR platform. In addition, the endpoint of tumor downstaging used by Ghadimi and colleagues (10) is subjected to numerous evaluation pitfalls. In contrast, TRG provides direct biological information about cytotoxicity, clearly beyond downstaging (26), and it seems to correlate with prognosis in both RC (21, 27–30) and other tumors (31–35). Furthermore, the 13-gene signature is the only one that predicts Dworak's TRG 0–1, which is the group containing the poorest responders. These patients likely have a bad prognosis and would be good candidates for alternative strategies to standard CRT, such as new drugs, doses, and time strategies. TRG 0–1 accounted for 25% of the patients in the largest data set from a randomized trial (21) and 19.1% of the patients in our series. Conversely, identifying patients with some grade of response (95% PPV in the training and testing cohorts) would be very important to further explore the promising strategy of induction chemotherapy that has been recently supported by the GEMCAD GCR3 randomized phase II trial (36). It would exclude those patients who do not benefit from CRT and prevent them from delaying surgical definitive treatment. Because CRT is standard treatment for stage II/III RC, the clinical interest of predicting Dworak's TRG 1–2–3 versus TRG 4 or patients with less than 10% tumor viable cells, as reported by Kim (12) and Rimkus (13), respectively, is less evident. In contrast, a new kinase profile has shown correlation between basal kinase activity and response to CRT in good, intermediate, and poor responders, although it was only validated among good and intermediate responders (25). Besides, these signatures have not been associated with prognosis. Our signature prognostic score remained the only statistically significant factor, even above postsurgical ypT and ypN, following multivariate analysis (HR = 4.35 95% CI: 1.2–15.75, P = 0.02). In our series, the TRG predictive signature-derived score, but not the TRG classification, correlated with outcome. A conceivable explanation for this phenomenon could be that the molecular portrait retrieves deeper informative biological aspects beyond TRG. In fact, a prognostic versus merely predictive capability cannot be excluded.

Limitations of our study include the small number of patients comprising the training and testing sets (57 and 37, respectively) and the heterogeneous treatment schedules. However, this study was larger than other RT/CRT predictive signatures published in RC. All patients received combined CRT with at least 2 cycles of chemotherapy. Therefore, it is not possible to determine whether their response was the result of either treatment component or the combined therapy. Nonetheless, combined CRT is the usual strategy, and the score identified significant DFS differences in both 5FU and oxaliplatin-based CRT groups (P = 0.03 and P = 0.015, respectively). Of note, recent data from randomized trials do not support the use of oxaliplatin during concurrent preoperative CRT in RC (37, 38).

In 2003, Crane and colleagues underlined an interest in identifying altered marker expression in response to preoperative treatment (39). We present the first global “dynamic” transcriptome study. Importantly, this study was done in CRT refractory patients by using SAGE, a powerful genomic tool. This labor-intensive method is very quantitative and is not restricted to a finite number of gene sequences localized to a chip. Although the exploratory study was restricted to 3 selected patients (cT3N0, TRG 0–1), the “intrapatient” changes in the Monte Carlo analysis provided plausible biological data. The majority of regulated genes were involved in oncogenic functions, and one-fourth of the observed changes (regulation of 24 genes) were shared among the 3 cases. These data highlighted interesting candidates for prognostic and functional evaluations. Moreover, regulated tags included an interesting 10% of expressed sequence tags.

The CRT-induced fibrosis in the surgical specimens was mirrored at the molecular level by upregulation of genes involved in the healing process. These genes included COL1A1, the key enzyme in the synthesis of collagens, PH-4 (40), and the TFF3 protein, which is involved in mucous restitution (41). The production of free radicals and hypoxia resulted in the regulation of GPX2 (42), HIG-2 (43), and angiopoietin-like 4 (44). Other genes that were consistently regulated have roles in cell-cycle proliferation and apoptosis and include CCNB1IP1 (45) and RCC1 (46). EEF2, an essential factor for protein synthesis that plays a role in gastrointestinal carcinogenesis (47), and WARS, a gene associated with poor prognosis in colorectal cancer (48), were also identified. Genes participating in the metastatic process were identified, including MMP14 (49), integrin B4 (50), and SPARC, a key controller of cell–matrix interactions that is also associated with tumor angiogenesis and metastasis (51, 52).

The SAGE screen of candidate genes added important value to the development of the gene predictor set. SAGE provided 5 of the 13 predictive genes (CKB, GPX2, HIG-2, PH-4, and TFF3) and a relatively small number of genes was required as compared with the other RC gene signatures [13 vs. 33 (ref. 11), 42 (ref. 13), 54 (ref. 10), 86 (ref. 25), and 95 (12)]. In general, genes in the signature suggest cell defense against CRT injury. This is exemplified by the downregulation of cell proliferation and apoptosis moderators, including the antiapoptotic BAK1 (53) gene and the antioxidant GPX2 (42), and the overexpression of the growth inhibitor CDKN1 (54) in responders. Consistent with most reports, thymidylate synthetase was overexpressed in nonresponders (55). Notably, only GPX2 has been previously reported as a component of other RC signatures (24). The lack of consistency among RC signatures is indicative of the complexity of resistance pathways in response to treatment injury.

In summary, by using a sequential genomic approach, we have developed a clinically exploitable gene signature that eventually may help select RC patients for customized treatment. Our results need to be confirmed and retrospective validation in well-controlled prospective clinical trials is currently ongoing by GEMCAD.

No potential conflicts of interest were disclosed.

This work has been funded by a FIS PI021094 grant. P. Cejas was supported by Asociación Española Contra el Cáncer (AECC) programa PAO 2010. V.M. Garcia was supported by Fundación Para la Investigación Biomédica del Hospital La Paz (FIBHULP) grant: REX 09.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.