Integrative Analysis of DCE-MRI and Gene Expression Profiles in Construction of a Gene Classifier for Assessment of Hypoxia-Related Risk of Chemoradiotherapy Failure in Cervical Cancer

This article is featured in Highlights of This Issue, p. 3985

Improvements in the therapy of locally advanced cervical cancer are highly needed, as many patients experience local or systemic relapse and severe side effects after treatment (1, 2). Tumor hypoxia is an adverse factor associated with radiation resistance and metastasis (3–5). Thus, hypoxia is both a potential marker for predicting risk of recurrence and a target in new treatment strategies combined with radiotherapy. Adjuvant hypoxia-modifying treatment has been shown to improve the efficacy of radiotherapy in cervical cancer (6), and combination trials to find the optimal drug have been performed (7, 8) and are underway (ClinicalTrials.gov identifier NCT02394652). In many trials, difficulties in detecting drug effects and toxicity are major problems, as patients are randomized to the adjuvant treatment without knowledge of the hypoxia status of the tumor (9). Development of methods for classifying patients according to the hypoxia status is therefore an important requirement for reliable drug evaluation and to avoid added toxicity to patients with no expected benefit.

Gene expression signatures can provide an indirect measure of the hypoxia level in a tumor through the cellular transcriptional response to low oxygen concentrations and therefore be of value for classifier development. Several hypoxia signatures have been constructed in different tissues (10–18), but a common cancer signature seems difficult to derive (19). The only hypoxia signature presented for cervical cancer was discovered in our previous work by combining global gene expression profiles and dynamic contrast enhanced (DCE)-MRI of the same patients (20). The amplitude A_Brix in the Brix pharmacokinetic model (21) extracted from DCE-MR images was found by us and others to reflect hypoxia in cervical tumors (20, 22). The signature contained 31 hypoxia-responsive genes, for which expression correlated with the patient-derived A_Brix data (20). This approach increased the likelihood of finding a signature that described a clinically relevant hypoxic phenotype. Both DCE-MRI and the 31-gene expression signature may have potential as a tool for assessing hypoxia in cervical cancer (20, 22). A current major limitation of MRI-based classifiers is, however, the lack of image standardization across MR machines (23, 24). Gene expression assays are easier to standardize, and multigene signatures have shown to be useful for guiding treatment decisions in many cancer types, like breast and prostate cancer (25, 26). Furthermore, gene signatures provide information that may increase the understanding of resistance mechanisms at play in individual tumors and thereby means for choosing the optimal targeted therapy in an individualized, combined radiotherapy regime (27).

Today, imaging and genomics appear as distinct facets in the treatment decision support system. Development of a direct link between the two disciplines would facilitate implementation of a multifactorial tool for a more accurate response prediction (24, 28). In this work, we aimed to prepare our signature for clinical use by constructing an image-associated, dichotomous hypoxia gene classifier with a predefined threshold for classification that could be utilized in a prospective setting. By integrative analysis of A_Brix and gene expression profiles, we identified genes in the signature that combined showed the optimal discrimination between more and less hypoxic tumors at the time of diagnosis. The power of this novel approach was demonstrated by the presentation of a hypoxia gene classifier that was reflected in DCE-MR images and showed prognostic impact across cohorts and test assays.

A total of 281 patients with locally advanced squamous cell carcinomas of the uterine cervix, prospectively recruited to our chemoradiotherapy protocol at the Norwegian Radium Hospital (Oslo, Norway) from 2001 to 2012, were included. Fifty-eight patients enrolled during the same period were not included because the tumor biopsies contained less than 50% malignant cells in hematoxylin and eosin–stained sections from the central part of the specimen. The patients constituted three independent cohorts; training cohort (n = 42), cohort 1 (n = 108), and cohort 2 (n = 131; Fig. 1; Supplementary Table S1). The patients were divided into the different cohorts based on the Illumina BeadArray version they were assayed by, that is, WG-6 v3 (training cohort, cohort 1) or HT-12 v4 (cohort 2). The patients assayed by the WG-6 arrays were further divided into the training cohort consisting of the patients with available A_Brix data and cohort 1 for the remaining patients. The patients were diagnosed in 2001 to 2004 for the training cohort, 2001 to 2009 (mainly 2004–2008) for cohort 1, and 2001 to 2012 (mainly 2005–2012) for cohort 2. The number of patients in cohort 1 and cohort 2 was sufficient for detecting an expected increase in the survival probability from 0.55 to 0.85 between more and less hypoxic tumors, respectively (20), with a power of 0.9 or larger (α = 0.05). Five additional patients with adenocarcinoma or adenosquamous carcinoma were included in the DCE-MRI analyses only. The study was approved by the Regional Committee for Medical and Health Research Ethics in southern Norway (S-01129). Written informed consent was achieved from all patients.

Tumor volume was determined from diagnostic MR images. Pathologic lymph nodes at diagnosis were also detected by MRI, or in a few cases by CT, according to the RECIST version 1.1 (29). One to four tumor biopsies (approximately 5 × 5 × 5 mm) were taken at diagnosis, immediately snap frozen, and stored at −80°C. All patients received external radiation of 50 Gy to the primary tumor and elective areas and 45 Gy to the remaining pelvis in 25 fractions five times per week. Pathologic lymph nodes received additional 14 Gy. This was followed by brachytherapy of 25 Gy in 5 fractions. Concurrent cisplatin (40 mg/m²) was given weekly in maximum 6 courses according to tolerance. Follow-up consisted of clinical examinations every third month for the first two years, twice a year for the next 3 years, and thereafter once a year. When symptoms of relapse were detected, MRI of pelvis and retroperitoneum as well as X-ray of thorax were performed. The study adheres to the REMARK criteria (30), and gene expression assays were performed blinded to the DCE-MRI results and clinical outcome.

Diagnostic DCE-MR images of 64 patients (42 patients in the training cohort, 17 patients in cohort 2, and 5 patients with adenocarcinoma or adenosquamous carcinoma) were applied for construction and validation of the classifier (Fig. 1). Imaging was performed using a 1.5 T Signa Horizon LX tomography (GE Medical Systems) with a pelvic-phased array coil. The MR examination included standard T1 and T2-weighted imaging, followed by DCE-MRI with 0.1 mmol/kg body weight Gd-DTPA (Magnevist, Schering) administered as a fast bolus injection, as described previously (31). The uptake of contrast agent, that is, the relative signal intensity increase as a function of time after Gd-DTPA injection, was analyzed with the Brix pharmacokinetic model (21), and the amplitude A_Brix was derived for each tumor voxel (31). An A_Brix histogram was generated for each tumor, and the mean A_Brix value over the 20^th to 30^th percentile interval, which showed the strongest association to outcome in previous work (20), was calculated and used as an image parameter reflecting hypoxia.

Hypoxia-induced genes were assessed in eight human cervical cancer cell lines from the ATCC: HeLa, SW756, C-33 A, C-4I, ME-180, HT-3, SiHa, and CaSki, and used for the first refinement of the signature. The cell lines were authenticated by STR/DNA profiling using PowerPlex 21 (Promega) by Eurofins Genomics. The cells were cultured at 37°C in a humidified 5% CO₂ incubator in DMEM with GlutaMAX (Gibco, Invitrogen) supplemented with 100 U/mL penicillin/streptomycin (Gibco) and 10% FCS. Cells were seeded in 10-cm plastic dishes (1–1.8 × 10⁶ cells to obtain ∼4 × 10⁶ cells after 48 hours of incubation) 24 hours before exposure to hypoxic (0.2% O₂, 5% CO₂) or normoxic (95% air, 5% CO₂) conditions for 24 hours at 37°C. The hypoxia treatment was performed in an Invivo2200 chamber (Ruskinn Technology).

Total RNA was isolated from frozen samples using TRIzol reagent (Invitrogen) (training cohort, cohort 1), miRNeasy Mini Kit (Qiagen; cohort 2), or RNeasy Mini Kit (Qiagen; cell lines). RNA from different biopsies of the same tumor was pooled. Gene expression profiling was performed using two different Illumina array generations; HumanWG-6 v3 (training cohort, cohort 1, cell lines) and HumanHT-12 v4 (cohort 2, adeno/adenosquamous carcinomas, cell lines; Illumina Inc.), both including approximately 48,000 transcripts. cRNA was synthesized, labeled, and hybridized to the arrays. Signal extraction and quantile normalization were performed using the software provided by the manufacturer (Illumina Inc.), and log2-transformed data were used in the analyses. The data were deposited in the gene expression omnibus database (GSE72723).

The transferability of the classifier across assay platforms was tested in a reverse transcription quantitative PCR (RT-qPCR) assay on 74 patients from the training cohort and cohort 1, using the 7900HT Fast Real-Time PCR System (Applied Biosystems). Genomic DNA was removed from the RNA samples by using the DNA-free Kit (Ambion) according to the manufacturer's instructions. Reverse transcription was performed by the High-Capacity cDNA Reverse Transcription Kit (Life Technologies) with 2,000 ng DNase-treated RNA in a 20 μl reaction. 1 μl cDNA was amplified using TaqMan Fast Universal PCR Master Mix (2×) and Custom TaqMan Array 96-Well Fast Plates with dried-down TaqMan Assays (Life Technologies) for three reference genes: that is, CHCHD1, SRSF9, and TMBIM6, which were identified as good endogenous controls for studying hypoxia in cervical cancer, and six classifier genes (Supplementary Table S2). Minus reverse transcriptase controls were tested for residual genomic DNA for all TaqMan assays. Nontemplate controls without RNA into the cDNA synthesis were included. All negative controls had undetermined quantification cycle (C_q). The data were analyzed using the comparative C_q method (32) by normalizing the C_q value of each classifier gene to the average C_q value of the three reference genes: dC_q,gene = C_q,gene − (C_q,CHCHD1 + C_q,SRSF9 + C_q,TMBIM6)/3. All samples were run in duplicate. Samples with mean C_q > 37 or SD > 0.4 were excluded.

To construct a gene classifier for separating the more and less hypoxic tumors, a previously presented algorithm (11) was utilized with small modifications, as described in Supplementary Document S1. A_Brix was used as a measure of hypoxia, and the construction was performed blinded to the clinical outcome. In brief, the classifier consisted of three parameters for each gene included, that is, the mean expression of the gene in a group of more hypoxic tumors (⁠|${\mu _{more}}$|⁠) and a group of less hypoxic tumors (⁠|${\mu _{less}}$|⁠) and its common variance (⁠|$W$|⁠) across both groups. The classifier parameters were determined as described below based on the training cohort, which was assumed to represent a typical patient cohort, and tabulated for further use in the classification of new tumors.

For a new tumor to be classified, the expression of each classifier gene in this tumor was compared with the tabulated classifier parameters |${\mu _{more}}$| and |${\mu _{less}}$| of the same gene, using the third classifier parameter, W, as weighting factor. This was performed for all classifier genes, and two expression distances, D_more and D_less, were calculated for the tumor:

where |$i$| is the more or less hypoxic group, |$m$| is the |$m$|th gene in the classifier, |${\mu _{im}}$| is the tabulated mean expression of gene |$m$| in group |$i$|⁠, |${y_m}$| is the expression of gene |$m$| in the tumor being classified, and |${W_m}$| is the tabulated common variance of gene |$m$|⁠. The tumor was classified into the group for which the D value was smallest, and hence, the weighted gene expressions were closest. Thus, by considering the difference D_less-D_more, tumors with a negative difference were classified as less hypoxic, whereas tumors with a positive difference were classified as more hypoxic. The classification threshold is therefore zero for the D_less-D_more difference and predefined by the method.

To determine the classifier parameters, we defined the more and less hypoxic groups in the training cohort by separating the tumors into a low and high A_Brix group, respectively. A leave-one-out cross-validation approach was used to determine both the optimal separation of the two groups based on A_Brix and which genes to include in the classifier to achieve correct classification according to A_Brix in the same analysis (Supplementary Document S1). The classifier parameters |${\mu _{more}},\ {\mu _{less}}$|⁠, and |$W$| were then determined for each classifier gene in the training cohort and tabulated (Table 1).

The classifier parameters are specific for the Illumina WG-6 assay. To use the classifier on Illumina HT-12 and RT-qPCR data, transformation of the gene expression levels was performed before calculation of the D values (Supplementary Document S2).

Clinical endpoints were progression-free survival (PFS), locoregional control (LRC), and disease-specific survival (DSS) for follow-up until 5 years. For patients with less than 5 years of follow-up, status was assessed on October 3, 2014. Patients were censored at their last appointment or at 5 years. For PFS, time from diagnosis to disease-related death or first event of relapse was used, LRC was defined as control within the irradiated pelvis, and DSS was defined as not dead from cervical cancer. Deaths within 2 years after inclusion were regarded as disease related unless another cause was documented.

Associations were estimated by Spearman rank correlation or Pearson correlation, and differences between groups were assessed with Fisher exact test. A competing risk model was used for calculating cumulative incidence, with death from other causes than cervical cancer as competing event. Fourteen of the 281 patients died of other causes within the follow-up period of 5 years. Gray test was used to assess differences between cumulative incidence curves. Cox uni- and multivariate proportional hazard (PH) analyses were performed, and assumptions of PHs were tested graphically using log-minus-log plots. Existing clinical markers, that is, tumor volume, Federation International de Gynecologie et d'Obstetrique (FIGO) stage, and lymph node involvement, were used as covariates to evaluate the independent prognostic value of the classifier. Significance level was 5%, and all P values were two sided. Cox PH analyses were performed in SPSS Statistics 21, whereas all other analyses were performed using R (33) version 3.1.1.

The classifier was based on a 31-gene signature identified in previous work (Fig. 1; ref. 20). To enrich for genes specific for hypoxia in cervical cancer, the 31-gene signature was first reduced to include only the 17 genes with at least a 2-fold upregulation by hypoxia in at least one of eight cervical cancer cell lines (Supplementary Table S3). These 17 genes served as a basis for classifier construction in the training cohort. We found an optimal separation of the cohort into 11 more hypoxic (low A_Brix) and 31 less hypoxic (high A_Brix) tumors (Fig. 2A; Supplementary Document S1). This separation led to correct classification in 33 to 34 of the 42 leave-one-out classifications when including the 5 to 12 top-ranked genes in the classifier (Fig. 2A). To determine how many genes to include in the classifier, the ranking of the genes was recorded for the 42 leave-one-out classifications (Supplementary Document S1). For the 6-gene classifier, the same genes were present in at least 41 of the classifications. Results for the other candidates were slightly less stable, and the final classifier was chosen to consist of these six genes (ERO1A, DDIT3, KCTD11, P4HA2, STC2, UPK1A; Fig. 2B; Supplementary Document S1).

On the basis of the expression levels of the six genes in a new tumor and the classifier parameters |${\mu _{more}}$|⁠, |${\mu _{less}}$|⁠, and |$W$| for each gene (Table 1), the expression distances D_more and D_less can be calculated for the independent tumor and its similarity to the more and less hypoxic groups determined. D_less-D_more represents degree of this similarity and is linearly correlated with the gene expression levels (P < 0.001; r = 0.99; n = 42; Supplementary Document S1). Classification of the tumors in the training cohort resulted in D_less-D_more values from −22.0 to 18.6 (Fig. 2C, left). Tumors with different D_less-D_more values showed differences in the A_Brix map (examples in Fig. 2C, right), and D_less-D_more showed a significant association to A_Brix (P < 0.001; ρ = −0.66; n = 42; Fig. 2C, left) that was just as strong as for the original 31-gene signature (20) represented by a gene score (P < 0.001; ρ = −0.68; n = 42). The association to A_Brix was thus retained after reducing the signature and constructing the classifier.

The prognostic significance of the classifier was evaluated in two independent datasets based on different Illumina array generations (WG-6 and HT-12; Fig. 1). For cohort 1, assayed using the WG-6 arrays, 40 of 108 patients (37%) were classified with a more hypoxic tumor. In this group, a significantly higher cumulative incidence of recurrence was found than in the less hypoxic group [40% vs. 15% at 5 years; P = 0.0035; HR, 3.06; 95% confidence interval (CI), 1.39–6.74; Fig. 3A]. Thus, the classifier could separate the patients into two groups with different outcome similarly to the original 31-gene hypoxia score (Fig. 3B), showing that the prognostic value of the signature was retained in the 6-gene classifier.

Of the 131 tumors in cohort 2, assayed by HT-12 arrays, 42 (32%) were classified as more hypoxic, and the patients in the more hypoxic group had a significantly higher risk of relapse than the others (36% vs. 20% at 5 years; P = 0.014; HR, 2.42; 95% CI, 1.20–4.90; Fig. 3C). In comparison, separation of patients based on the original 31-gene hypoxia score with the same threshold as in Fig. 3B gave a smaller, not statistically significant difference between the groups (31% vs. 18% at 5 years; P = 0.13; HR = 1.75; CI, 0.84–3.66). Moreover, a significant correlation between A_Brix and D_less-D_more was found in a subgroup of the patients with available A_Brix data (P = 0.027; ρ = −0.47; n = 22; Fig. 3D). These findings validated the prognostic significance and the A_Brix association of the classifier and demonstrated its robustness across array generations.

We evaluated the ability of the classifier to predict the hypoxia status defined by A_Brix using all 64 tumors with available A_Brix data. Forty-nine out of 64 tumors were correctly classified according to A_Brix, giving an accuracy of 0.77. Hence, there was a statistically significant association between classifier- and A_Brix-based classification (P < 0.001; Supplementary Table S4). Furthermore, the classifier identified most of the A_Brix-defined more hypoxic tumors correctly with a sensitivity of 0.93 but also classified some of the A_Brix-defined less hypoxic tumors as more hypoxic, resulting in a precision of 0.48 and a specificity of 0.72 for this group (Supplementary Table S5). The precision of less hypoxic classification was 0.97, showing that almost all of the tumors classified as less hypoxic were correctly classified according to A_Brix. Of importance, the D_less-D_more value and/or the A_Brix value of most of the incorrectly classified tumors were close to the borders between the more and less hypoxic groups, and as many as 6 of the 14 tumors incorrectly classified as more hypoxic experienced relapse (Supplementary Fig. S1).

To further evaluate its robustness, a technical transfer of the classifier to RT-qPCR was performed for 74 patients (Fig. 1). For 19 of the 74 samples, gene expression measurements for KCTD11 and/or UPKIA were unsuccessful. The robustness of the classifier when lacking UPK1A and/or KCTD11 as compared with the full 6-gene classifier was estimated from the Illumina data of the 281 patients included in the study. Prognostic significance was obtained for all four combinations of the 4–6 genes (Supplementary Table S6). All 74 patients could therefore be classified on the basis of the data of 4–6 genes. The expression levels of the six classifier genes showed a strong correlation between the RT-qPCR and Illumina WG-6–based measurements, with a Spearman rank correlation coefficient above 0.7 (Supplementary Fig. S2). Moreover, there was a significant correlation between the Illumina-based and RT-qPCR–based classifications (P < 0.001; ρ = 0.74; Fig. 3E). Of 74 patients, 56 (75.7%) were classified to the same group by the RT-qPCR and Illumina assay (Supplementary Table S7), showing reproducibility of the classifier using two different assays (P < 0.001). For the remaining 18 patients, 11 were classified with a more hypoxic tumor by Illumina and not by RT-qPCR, whereas 7 were classified with a more hypoxic tumor by RT-qPCR and not by Illumina. For many of the patients classified differently with the two methods, D_less-D_more was close to the border between the two groups (Fig. 3E). Moreover, the few patients with missing RT-qPCR data for one or two genes were not classified differently as compared with the others (Fig. 3E), illustrating the robustness of the classifier in that missing values can be tolerated.

The 26 patients (35%) classified to the more hypoxic group by RT-qPCR had a poor prognosis compared with the others (58% vs. 35% at 5 years; P = 0.035; HR = 2.12; CI, 1.06–4.26; Fig. 3F), and the separation between the groups were similar to that obtained with Illumina data for the same patients (58% vs. 34% at 5 years; Supplementary Fig. S3A). Also the A_Brix association could be validated across the platforms, although the correlation was somewhat stronger for Illumina WG-6 (P < 0.001; ρ = −0.70; n = 32) than RT-qPCR (P = 0.003, ρ = −0.51; n = 32; Supplementary Fig. S3B and S3C).

To compare classifier performance with existing clinical markers, all cohorts were merged (Fig. 1). Of 281 patients, 100 (36%) were classified with a more hypoxic tumor. The probabilities of having high tumor stage, large tumor volume, or lymph node involvement were significantly higher for the more hypoxic group compared with the less hypoxic group (Table 2). Among the patients in the more hypoxic group, 36% had a stage IIIA–IVA tumor, 63% had a tumor volume above the median value of 36.8 cm³, and 52% had lymph node involvement at the time of diagnosis as compared with 24% (P = 0.037), 43% (P = 0.0019), and 38% (P = 0.032) in the less hypoxic group.

In univariate Cox regression analyses, both the clinical markers and 6-gene classifier showed prognostic significance regardless of endpoint, with HRs above 2 (Table 3). Moreover, the significance of the classifier was retained in multivariate analyses, with HRs of 2.12 (PFS), 2.62 (LRC), and 2.12 (DSS; Table 3). Thus, even though there were associations between hypoxia classification and the other markers, the classifier emerged as an independent prognostic factor at all endpoints.

Our study is the first to construct a classifier for assessing hypoxia-related risk of chemoradiotherapy failure in cervical cancer. By focusing on an aggressive phenotype, rather than solely on the patient outcome, which is often done in development of prognostic tests (25, 34), our classifier may be more robust against possible gradual improvements in diagnostics and treatment over time. Furthermore, to our knowledge, our approach is the first to combine functional imaging and gene expression variations in the construction of a gene classifier for clinical use. Correlation analyses of image features and gene expression have proven powerful in exploration of the molecular background of the images by us (20) and others (35). Here, we demonstrate a novel way of combining imaging and gene expression to construct a genomic tool for utilization in the treatment decision making both alone and with a potential for future integration with DCE-MRI.

The gene classifier was constructed by using the DCE-MRI parameter A_Brix as an independent measure of tumor hypoxia. Previously identified A_Brix-associated hypoxia-responsive genes (20) were used as potential classifier genes, and the final classifier consisted of six genes that together showed the best separation of the tumors in the training cohort according to hypoxia status. A new tumor was evaluated by comparing the expression levels of these genes with the tabulated mean levels (classifier parameters), determined for the more and less hypoxic groups in the training cohort, using a predefined threshold of zero for the output, D_less-D_more, to classify the tumor as more or less hypoxic. Hence, our dichotomous classifier with a predefined threshold would be easy to interpret and applicable in a prospective setting in the clinic.

Application of the classifier in different cohorts and test assays suggested that the method provides a robust assessment of hypoxia-related aggressiveness. Associations with outcome were demonstrated for two cohorts recruited over a time period of 11 years, during which MR image-guided brachytherapy has been implemented in our hospital and replaced brachytherapy based on 2D X-ray imaging (36). Moreover, the second cohort was analyzed with another Illumina array generation, which required data transformation prior to classification. Successful utilization of the classifier was also demonstrated in a RT-qPCR assay, with samples collected over a period of 8 years. In all tests of the classifier, not only its prognostic impact but also its association to A_Brix was retained. Notably, the classifier seemed to be more robust than the original 31-gene score, for which association to clinical outcome failed to reach statistical significance in the second cohort.

The classifier was associated with outcome independent of conventional markers for PFS, LRC, and DSS and may provide information about disease progression that is not covered by current diagnostics. Even though FIGO is a strong prognostic factor in cervical cancer, the classifier emerged as an independent prognostic factor. Importantly, through its biologic fundament, the classifier could be a well-needed tool for selecting patients to hypoxia-modifying adjuvant therapy and to avoid added side effects to those who are not expected to respond. The latter scenario is important because the risk of damage to critical pelvic organs is already a significant challenge in the conventional chemoradiotherapy (1, 2, 37). In addition, many patients classified with hypoxia in our work had a high risk of both failure and side effects by conventional markers, presented with large tumor, lymph node involvement, and/or stage III or IV disease. These patients belong to a group where the gain of treatment improvements would be particularly high (1, 2). A higher incidence of lymph node–positive patients in the more hypoxic group and a weak correlation between hypoxia and tumor size have been found in earlier cervical cancer studies using oxygen electrodes to assess hypoxia (4, 38, 39), in accordance with our findings.

Hypoxia gene signatures have been constructed in other tissue types by different approaches, and some signature genes are in common. Among 32 signatures reviewed by Harris and colleagues (19), two of our six classifier genes, P4HA2 and ERO1A, are among the 20 most frequently appearing ones. The presence of common genes indicates that it might be possible to identify a signature that is applicable in different tumor sites, which has indeed been demonstrated in some cases (10, 15, 17, 40, 41). However, although the head and neck cancer signature by Eustace and colleagues identified patients with laryngeal cancer who seemed to benefit from adjuvant hypoxia-modifying treatment, this was not the case when the signature was used on patients with bladder cancer (12). Moreover, when comparing six of the previously published hypoxia gene signatures (10–15) with the classifier genes in this work, we found that our genes had the strongest prognostic impact in the cervical cancer dataset (Supplementary Document S3). In particular, only the signature developed in breast cancer by Hu and colleagues (15) showed a significant association to outcome (P = 0.026) in Cox PH univariate analysis. Out of our classifier genes, only KCTD11 and P4HA2 were present in three (10–12) or one (11) of the six signatures, respectively. These results support the view that there are molecular differences in the aggressive hypoxic phenotype across cancers (19) and that cancer-type specific hypoxia gene classifiers are more effective (18). Furthermore, although some of the other signatures were statistically significantly correlated with A_Brix (Supplementary Document S3), the signature of our classifier genes showed a stronger A_Brix-association and was thus the best one for describing the aggressive phenotype reflected in the DCE-MR images of cervical cancer.

The gene classifier appeared to classify a higher number of tumors as more hypoxic compared with A_Brix-based classification, and there might be an intermediate group that should be handled separately in a clinical trial. However, as 43% of the tumors classified as more hypoxic only with the classifier experienced recurrence, we speculate that the classifier might be a better method for identifying the aggressive hypoxic tumors than A_Brix. Furthermore, the classifier was shown to have a high precision in classifying the less hypoxic tumors, thus correctly identifying the patients that may not benefit from hypoxia-modifying treatment and potentially sparing many patients for unnecessary toxicity. A combination of the gene classifier with DCE-MRI–based hypoxia evaluation for improved prediction of hypoxia-related failure would be of interest to evaluate in future work.

The Brix model used to calculate A_Brix represents a simplistic description of the uptake of contrast agent in tumors (21), and this weakness in the DCE-MRI analysis might explain why A_Brix failed to identify some of the aggressive tumors detected only with the classifier. This limitation of our work needs to be addressed before a multifactorial tool based on the combined use of our classifier and imaging can be developed for clinical use. Furthermore, the predictive value of the classifier in identifying the tumors that will benefit from hypoxia-modifying treatment remains to be shown. Two head and neck cancer signatures have shown promise as predictive hypoxia biomarkers (12, 42, 43), illustrating the potential of gene signatures for hypoxia stratification. For cervical cancer, the optimal hypoxia-targeting drug has yet to be defined (7, 8), and careful considerations in the study design have to be done due to toxicity problems. Future work should involve multicenter studies for verifying the classifier's prognostic value and testing the predictive impact in a clinical trial.

The classifier genes provide biologic insight into the disease that may aid the choice of drugs for combination trials. STC2, P4HA2, and ERO1A have been identified as direct target genes of the hypoxia-inducible factor 1 (HIF1; refs. 44–47). In addition, DDIT3 (48), ERO1A (49), and STC2 (50) may be activated through the unfolded protein response (UPR). Thus, both HIF1 and UPR pathways, which are major oxygen-sensitive pathways promoting hypoxia tolerance (51), might be involved in mediating the aggressive phenotype reflected by our classifier. Several HIF1 and UPR-targeting drugs are currently under development (48, 52–54) and could be relevant to evaluate in the selected patients.

In conclusion, we have developed a prognostic gene classifier that is associated with tumor hypoxia as evaluated by DCE-MRI and can characterize a tumor as either more or less hypoxic based on the expression level of six genes in a biopsy. The classifier may give an early indication of a patient's risk of hypoxia-related recurrence and has potential to serve as a biomarker for patients that would be eligible for clinical trials testing hypoxia-modifying drugs. This could help implementing hypoxia-guided treatment strategies for patients with cervical cancer. Furthermore, the approach of combining functional imaging and gene expression variations in tumors may also be beneficial for the development of other phenotype-specific gene classifiers for clinical use.

C.H. Julin and H. Lyng are listed as inventors on a patent application covering the clinical use of the hypoxia gene signature that is owned by Oslo University Hospital (WO2013/124738).

No potential conflicts of interest were disclosed by the other authors.

The funding sources of the study have no role in the study design; in the collection, analysis, or interpretation of the data; or in the writing of the article or the decision to submit the article for publication.

Conception and design: C.S. Fjeldbo, C.H. Julin, H. Lyng

Development of methodology: C.S. Fjeldbo, C.H. Julin, J. Alsner

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): C.S. Fjeldbo, C.H. Julin, M.F. Forsberg, E.-K. Aarnes, G.B. Kristensen, E. Malinen

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): C.S. Fjeldbo, C.H. Julin, M. Lando, M.F. Forsberg, E.-K. Aarnes, J. Alsner, E. Malinen, H. Lyng

Writing, review, and/or revision of the manuscript: C.S. Fjeldbo, C.H. Julin, M. Lando, M.F. Forsberg, E.-K. Aarnes, J. Alsner, G.B. Kristensen, E. Malinen, H. Lyng

Study supervision: H. Lyng

This work was supported by The Research Council of Norway (grant no. 226120/O30), The Norwegian Cancer Society (grant no. 107438 - PR-2007-0179), and The South-Eastern Norway Regional Health Authority (grant no. 2012015).

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.