Identification of genes and pathways associated with cytotoxic T lymphocyte infiltration of serous ovarian cancer

Background: Tumour-infiltrating lymphocytes (TILs) are predictors of disease-specific survival (DSS) in ovarian cancer. It is largely unknown what factors contribute to lymphocyte recruitment. Our aim was to evaluate genes and pathways contributing to infiltration of cytotoxic T lymphocytes (CTLs) in advanced-stage serous ovarian cancer. Methods: For this study global gene expression was compared between low TIL (n=25) and high TIL tumours (n=24). The differences in gene expression were evaluated using parametric T-testing. Selectively enriched biological pathways were identified with gene set enrichment analysis. Prognostic influence was validated in 157 late-stage serous ovarian cancer patients. Using immunohistochemistry, association of selected genes from identified pathways with CTL was validated. Results: The presence of CTL was associated with 320 genes and 23 pathways (P<0.05). In addition, 54 genes and 8 pathways were also associated with DSS in our validation cohort. Immunohistochemical evaluation showed strong correlations between MHC class I and II membrane expression, parts of the antigen processing and presentation pathway, and CTL recruitment. Conclusion: Gene expression profiling and pathway analyses are valuable tools to obtain more understanding of tumour characteristics influencing lymphocyte recruitment in advanced-stage serous ovarian cancer. Identified genes and pathways need to be further investigated for suitability as therapeutic targets.

Epithelial ovarian cancer is the most common cause of death from gynaecological malignancies (Parkin et al, 2005). The 5-year survival rates for ovarian cancer patients do not exceed 40%. This high mortality is best attributed to the absence of specific symptoms combined with the lack of reliable screening methods, which prevent diagnosis in the early stages of the disease in a majority of patients. Treatment generally consists of cytoreductive surgery followed by platinum-and taxane-containing chemotherapy. Classic prognostic factors are stage of disease at diagnosis, histological tumour type and grade, residual disease after primary surgery, and response to chemotherapy (Crijns et al, 2006). An increasing body of evidence suggests that next to these established prognostic factors, the presence of tumour-infiltrating lymphocytes (TILs) also independently contributes to prognosis (Zhang et al, 2003;Sato et al, 2005;Leffers et al, 2009a). Although the presence of TILs is generally considered a reflection of antitumour immunity, it is largely unknown why TILs are present in high numbers in some tumours and largely absent in others. It has been shown that endothelial factors and chemokines secreted by the tumour may have an important role (Curiel et al, 2004;Buckanovich et al, 2008). The existence of an antitumour immune repertoire in a selection of patients forms the rationale for the development of cancer immunotherapy. Although immunotherapy strategies generally induce potent peripheral immune responses in ovarian cancer patients, clinical responses have so far been disappointing (Hung et al, 2008). The combination of targeted agents that enhance lymphocyte recruitment to tumour sites with these immunotherapy strategies might be a lucrative approach to obtain clinical responses to immunotherapy. For instance, the in vivo addition of BQ-788, an endothelin B receptor antagonist, to previously immunogenic, but clinically ineffective, immunisation strategies resulted in enhanced homing of lymphocytes to tumours as well as improved clinical responses (Buckanovich et al, 2008).
To investigate what tumour factors contribute to the recruitment of lymphocytes, we analysed which genes and pathways were associated with the presence of tumour-infiltrating CTLs in a homogeneous group of 49 advanced-stage serous ovarian cancer patients previously profiled at our institute as part of a larger study (Crijns et al, 2009). The prognostic value of identified genes and pathways was subsequently validated on all 157 previously profiled late-stage serous ovarian cancer patients (Crijns et al, 2009). Furthermore, immunohistochemical staining of tissue microarrays was performed to validate findings.

PATIENTS AND METHODS Patients
The presence of tumour-infiltrating T lymphocytes was previously evaluated by our group in 306 ovarian cancer patients (Leffers et al, 2009a). In short, we used a tissue microarray containing core biopsies (0.283 mm 2 each) of tumour tissue obtained from each of 306 ovarian cancer patients at primary surgery. Sections were stained with anti-CD8 (1 : 20; Dako Cytomation, Glostrup, Denmark). The number of intraepithelial CD8 þ CTLs were separately counted for each core and an average per core was calculated when at least two cores were present for a single tumour sample. Patients with 48 CD8 þ T lymphocytes per 0.283 mm 2 were found to have a better prognosis. For this study, we selected from this heterogeneous population only advanced-stage serous ovarian cancer patients with low (p5 per 0.283 mm 2 of tumour) or high (X8 per 0.283 mm 2 of tumour) CTL numbers who were also included in the previously published microarray study (Crijns et al, 2009). This large microarray study contained 157 advancedstage ovarian cancer patients for whom fresh frozen tumour tissue was available. Only 63 of these patients were also represented on our tissue microarray, of whom only 49 met the requirements stipulated above.
Patients were treated at the University Medical Center Groningen by a gynaecological oncologist and staged according to FIGO classification (Cancer Committee of the International Federation of Gynaecology and Obstetrics, 1986). Tumours were graded and classified according to the WHO criteria by a gynaecological pathologist (Scully, 1999). Adjuvant chemotherapy generally consisted of different platinum-based treatment regimens. Response to chemotherapy was evaluated according to the WHO criteria (World Health Organization, 1979). After treatment, patients were followed-up for at least 10 years with gradually increasing intervals. Informed consent was obtained for the collection and storage of tumour samples in a tissue bank for future research. Information on clinicopathological characteristics and follow-up of patients was obtained from a computerised database in which information of all patients with epithelial ovarian cancer treated at our institute is prospectively recorded. For this study, relevant data were retrieved into a separate anonymous database. In this separate database, patient identity was protected by study-specific, unique patient codes. In case of uncertainties with respect to clinicopathological and follow-up data, the larger databases could only be checked through two data managers who have daily responsibilities for the larger database, thereby ascertaining the protection of patients' identity. According to Dutch law no approval from our IRB was needed.

Microarray analysis
As mentioned above, we selected 49 advanced-stage serous ovarian cancer patients based on the relative absence or presence of CD8 þ TILs, who were previously profiled as part of a larger study (Crijns et al, 2009). In brief, after RNA extraction and amplification, samples were hybridised to two-colour 70-mer oligonucleotide microarrays (B35 000 Operon v3.0 probes, Ebersberg, Germany). All samples were hybridised at least twice and samples were loaded using a random design to prevent systematic biases (Churchill, 2002;Kerr, 2003;Hsu et al, 2007). Arrays were scanned using the Affymetrix GMS428 (Affymetrix, Santa Clara, CA, USA). Expression values were calculated using Bluefuse software (BlueGnome, Cambridge, UK). Operon v3.0 probe identifiers were converted to official gene symbols using probe annotations provided by The Netherlands Cancer Institute (NKI, Amsterdam, The netherlands). Only oligonucleotides specifically responding with a single hit on a gene during a BLAST search were used. Expression values of multiple probes targeting a single gene were averaged, resulting in a total of 15 909 distinct genes. Subsequently, expression data of the multiple hybridisations per tumour sample were averaged. Microarray data of the previous, larger study from which our patients were selected are available at http://www.ncbi.nlm.nih.gov/geo/ under number GSE13876.

Class comparison between low TIL and high TIL samples
The BRB Array Tools 3.6.0 software package, developed by the Biometric Research Branch of the US National Cancer Institute, was used for class comparison between low-TIL and high-TIL samples. Differentially expressed genes were identified using a paired T-test (threshold Po0.05).

Gene set enrichment analysis (GSEA)
As it is unclear whether large differences in the expression of a single gene are biologically more relevant than more subtle, although coordinated, differences in a set of genes belonging to a single biological pathway, we performed gene set enrichment analysis. Expression data of all 15 909 genes were compared against functional gene sets to determine whether any of these sets were enriched in samples containing high or low numbers of CD8 þ TILs. The comparison was performed using 340 gene sets reported in two databases (174 sets from BIOCARTA: http://www.biocarta. com; 166 sets from Kyoto Encyclopedia of Genes and Genomes database (KEGG): http://www.genome.jp/kegg/). Statistical significance of enrichment was determined using an empirical genebased permutation test using 1000 permutations. Gene sets with an enrichment P-value of o0.05 are reported. We also calculated false discovery rates (FDRs) for each functional gene set, which represent the estimated probability that a given enrichment score represents a false-positive finding. We only report gene sets with an FDR of o0.25. With such an FDR, the results are likely to be valid at least three out of four times, which is considered a suitable cutoff for the generation of interesting hypotheses (Gu et al, 2007). The GSEA was executed using the GSEA 2.0 software package (Broad Institute of MIT and Harvard, Cambridge, MA, USA).

Leading-edge subset analysis
The leading-edge subset is defined as the subset of genes in a functional gene set that appears in the ranked list of 15 909 genes at, or before, the point in which the running enrichment score reaches its maximum deviation from zero. The genes within this subset can be interpreted as the most important in the enrichment of the functional gene set. Leading-edge subsets were defined for all statistically enriched function gene sets (Po0.05). Subsequently, overlap between leading-edge subsets from significantly enriched functional gene sets identified in the different databases was determined to discover genes belonging to more than one leading-edge subset, that is, possible key genes. patients previously profiled at our institution (Crijns et al, 2009). The log expression levels of individual genes were entered into a univariate Cox proportional hazards regression model. Genes with a Po0.05 were considered to be associated with DSS. Furthermore, GSEA was performed on the cohort of 157 ovarian cancers to evaluate associations between DSS and identified pathways.

Immunohistochemistry
Protein expression of selected enriched genes from identified pathways was evaluated by immunohistochemistry using tissue microarray sections. Tissue microarrays were constructed from paraffin-embedded tumour tissue obtained at primary debulking surgery, performed by the gynaecological oncologists from the University Medical Center Groningen between May 1985 and April 2003. The tissue microarrays contain four 0.6 mm core biopsies from each of 361 patients. For this study, tissue samples obtained at primary debulking surgery from 108 advanced-stage serous ovarian cancer patients were analysed, for whom staining of TILs and HLA-A and HLA-B/C has previously been performed and published (Leffers et al, 2009a, b). In addition, staining was performed for HLA-DP/DQ/DR (clone CR3/43, DAKO, Heverlee, Belgium). In brief, after dewaxing and rehydration, 4 mm sections were microwaved in 10 mM citrate buffer pH 6.0 for antigen retrieval. Subsequently, sections were incubated overnight with the primary antibody (dilution 1 : 100). Sections were subsequently incubated with DAKO Envision þ (DAKO). Antigen -antibody reactions were visualised with 3,3-diaminobenzidine. Tissue was counterstained with haematoxylin.
Sections were scored independently by two observers (MG/NL) who were unaware of clinicopathological characteristics and the TIL status of patients. A semiquantitative quality control system was used taking into account both the intensity of staining and the percentage of positive tumour cells as previously described (Leffers et al, 2009b). The sum of both scores was used to identify three categories of expression: no expression (total score 0 -2), positive expression in a proportion of cells or weak expression in all cells (total score 3 -6), and positive expression in a majority of cells (total score 7 and 8). Only patients for whom at least two evaluable cores were available were included for further analysis.
Associations between CD8 þ TIL and protein expression of selected enriched genes in identified pathways were evaluated using the Jonckheere -Terpstra test (a form of the Kruskal -Wallis test that also tests for linearity). The Mann -Whitney U-test was used to evaluate associations between CD8 þ TIL and well-known prognostic factors (i.e., age, FIGO stage, histological grade, and residual tumour after primary debulking surgery). Disease-specific survival was defined as the date of surgery until the date of death because of ovarian cancer (including fatal complications of treatment) or the date of last follow-up. Differences in DSS based on protein expression levels were plotted using Kaplan -Meier survival curves and evaluated using log-rank tests. SPSS software package for Windows, version 16.0 (SPSS Inc., Chicago, IL, USA) was used. P-values of o0.05 were considered statistically significant.

Patient characteristics
Microarray data of 49 advanced-stage serous ovarian cancer patients with either low (n ¼ 25) or high numbers of CD8 þ TILs (n ¼ 24) were evaluated (Figure 1). Clinical and pathological characteristics as well as tumour percentage of the samples used for microarray did not differ between patients with low or high numbers of CD8 þ TILs (Table 1). Median DSS was higher for patients with high CD8 þ TILs (log-rank test P ¼ 0.025).

Differential expression of genes and biological pathway analysis
A comparison of the expression levels of all 15 909 genes showed differential expression of 320 genes between tumours containing low and high numbers of CD8 þ TILs (Po0.05). The differences between low and high CD8 þ TIL samples were small (Supplementary Table 1). In view of these small differences, we subsequently evaluated whether coordinated differences in genes belonging to a single biological pathway existed rather than differences in the expression levels of single genes. In all, 14 pathways in BIOCARTA and 8 pathways in KEGG (Po0.05, FDR o0.25) were enriched in tumour samples with high numbers of CD8 þ TILs, whereas only one pathway was enriched in tumour samples with low CD8 þ TIL numbers (Table 2). Interestingly, and conveniently serving as an internal control, one of the pathways enriched in high CD8 þ TIL tumours was the CTL-mediated Study on tumour-infiltrating T lymphocytes (n = 306) Late-stage serous ovarian cancer patients with known CD8 + T lymphocyte status in primary ovarian tissue (n = 108) Gene profiling study Late-stage serous ovarian cancer patients (n = 157) 320 differentially expressed genes 23 differentially activated pathways Validation 2: Immunohistochemical evaluation of enriched genes from antigen processing and presentation pathway (n = 108) Figure 1 Flowchart illustrating patient selection and validation. Index and validation patients were selected from previous studies at our institute investigating the prognostic effect of gene expression (Crijns et al, 2009) and tumour-infiltrating lymphocytes on ovarian cancer patients (Leffers et al, 2009a).
immune response pathway. We subsequently performed leadingedge subset analysis to identify key regulatory genes common to the enriched pathways (Table 3), which among others identified several genes encoding for HLA molecules.

Effect on DSS
Univariate survival analysis was performed to assess the prognostic value of the 320 genes identified as differentially expressed between high and low CD8 þ TIL tumours, using 157 late-stage serous ovarian cancer patients previously profiled at our institute (Crijns et al, 2009). A significant association with DSS was observed for 54 genes (Table 4). Of the 23 genes associated with high CD8 þ TIL tumours, 21 were associated with improved survival. Conversely, 27 of 31 genes associated with low CD8 þ TIL tumours were associated with decreased survival (Figure 2). We next performed GSEA to identify which pathways were associated with DSS. Eight pathways associated with the presence of CD8 þ TIL, such as the antigen processing and presentation pathway, were also positively associated with DSS ( Table 2). The interdependency of these pathways has been visualised in Figure 3A and B.

Immunohistochemical validation
On the basis of GSEA and leading-edge subset analysis, which showed that the presence of TILs was associated with the expression of MHC class I and II genes as part of the antigen processing and presentation pathway, we evaluated immunohistochemical staining of HLA-A, HLA-B/C, and HLA-DP/DQ/DR in 108 advanced-stage serous ovarian cancer patients for whom information on CD8 þ TILs was available, part of which was previously published for a larger patient cohort (Leffers et al, 2009a, b). Partial or total loss of HLA-A and HLA-B/C was observed in 70.4 and 62.0% of patients, respectively, whereas HLA-DP/DQ/DR upregulation was observed in 68.5% of patients. Increasing levels of MHC class I and II protein expression strongly correlated with increased numbers of CD8 þ TILs (Table 5). However, no association was observed between the expression of HLA-A, HLA-B/C, and HLA-DP/DQ/DR with DSS (P ¼ 0.114, P ¼ 0.599, and P ¼ 0.692, respectively).

DISCUSSION
Although it has been well established that TILs are predictors of DSS in ovarian cancer, it is largely unknown what factors contribute to the presence of these lymphocytes. By comparing gene expression profiles of 25 tumours containing low and 24 tumours containing high numbers of CD8 þ TILs, we identified 320 genes differentially expressed by primary tumours of late-stage serous ovarian cancer patients. In addition, for 54 of these genes, an association with survival was observed in a large validation cohort containing 157 advanced-stage serous ovarian cancer  patients. Genes connected to high CD8 þ TIL tumours were associated with improved survival. With GSEA, next to pathways merely reflecting the presence of lymphocytes, several pathways were identified to be associated with the (lack of) CD8 þ TILs, some of which were also associated with survival in our validation set. Lastly, the association of a number of genes from the enriched antigen processing and presentation pathway with the presence of CD8 þ TIL was confirmed using immunohistochemistry. Our study illustrates that gene expression profiling is a valuable approach to elucidate what tumour cell characteristics contribute to or impede recruitment of lymphocytes into serous ovarian cancer. However, a problem inherent to the design of our study is the impossibility to discern what genes were expressed by tumour cells and what signal was derived from TILs themselves. This is substantiated by the fact that some of the pathways found to be enriched in high CD8 þ TIL tumours are lymphocyte-specific pathways, for example, T cytotoxic pathway. To our advantage, this observation can also be regarded to validate the immunohistochemical evaluation of CD8 þ T-cell count used for patient selection. To avoid the signal from lymphocytes, one should profile a limited number of tumour cells, which could be accomplished by microdissection of tumour cells (Glanzer and Eberwine, 2004). However, as with decreasing cell numbers, RNA yield also reduces such assays heavily depend on mRNA amplification. Especially genes with low numbers of transcripts may be under-represented after amplification and thus not identified in subsequent profiling studies (Nygaard et al, 2005). For this study, we therefore decided that, as the percentage of tumour cells did not differ between samples with low and high numbers of CD8 þ TILs, it was safe to assume that differences in non-lymphocyte-restricted genes and pathways reflected differences in gene expression by tumour cells.
Ultimately, our study was intended not only to establish what tumour factors contribute to lymphocyte recruitment, but also to discover putative factors that might enhance clinical results of immunotherapy for ovarian malignancies by improving lymphocyte recruitment when targeted. In this respect, one of the interesting genes identified as differentially expressed between high CD8 þ TIL and low CD8 þ TIL tumours and associated with DSS of late-stage serous ovarian cancer patients is interferon regulatory factor 1 (IRF-1). Recently, IRF-1 was reported to be a positive prognostic factor in ovarian and colorectal cancer (Galon et al, 2006;Zeimet et al, 2009) and was found to be associated with infiltration of CD8 þ T lymphocytes in ovarian cancer (Callahan et al, 2008). Binding of IFN-g to the IFN-g receptor leads to upregulation of IRF-1, which in turn results in: (1) induction of IFN-g-inducible genes, such as the TAP1, LMP2, and b 2 -microglobulin genes of the MHC class I-dependent pathway, as well as (2) activation of CIITA, a critical transcription factor for MHC Abbreviation: KEGG ¼ Kyoto Encyclopedia of Genes and Genomes. class II gene expression. Thus, IRF-1 facilitates recognition of tumour cells by immune cells, ultimately resulting in an IFN-gdependent positive feedback loop. Correspondingly, in high CD8 þ TIL tumours we observed differential expression of several MHC class I and II genes, all part of the selectively activated antigen processing and presentation pathway. Immunohistochemical evaluation confirmed the positive association of intra-tumoural cytotoxic T cells with surface expression of MHC class I and II molecules HLA-A, HLA-B/C, and HLA-DP/DQ/DR. Although several MHC class II alleles were differentially expressed in high CD8 þ TIL tumours, no association was found between HLA-DP/ DQ/DR membrane expression and survival. Moreover, only HLA-DQB2, a virtually non-polymorphic gene, was associated with DSS in our large validation cohort (Gaur et al, 1992). As it encodes six thus far unknown putative proteins, the association of HLA-DQB2 with CD8 þ TILs and survival is intriguing and deserves further investigation. Although we previously reported decreased survival for ovarian cancer patients in association with HLA-B/C downregulation, an association with survival was observed neither at mRNA nor protein level in this study (Leffers et al, 2008). The difference in prognostic effect observed in this study could be explained by differences in study population and/or size. In this study a homogeneous population of 108 late-stage serous ovarian cancer patients was used for immunohistochemical validation. Renewed analysis of the previously published data using only late-stage serous ovarian cancer patients (n ¼ 151) did not yield a correlation of HLA-B/C expression with DSS either (data not shown).
To our knowledge, only one study attempting a better understanding of mechanisms underlying the infiltration by CTLs of serous ovarian cancer by gene profiling of tumour samples has previously been published (Callahan et al, 2008), although the genetic signature of ovarian cancer containing high numbers of regulatory T cells was recently published (Barnett et al, 2010). It seems that here too genes of the antigen presentation pathway have a major role. Three genes seem to be important for both recruitment of regulatory FoxP3 þ and cytotoxic CD8 þ lymphocytes, that is, CCL5, APOL6, and IRF-1. In the study of Callahan et al (2008) that profiled 38 high-grade advanced-stage ovarian carcinomas, 81 genes were associated with CD8 þ T-cell infiltrate. Only two of these genes, IRF-1 and CXCR6, were also found to be differentially expressed in our cohort of 49 serous advanced-stage ovarian cancer patients (P ¼ 0.018578 resp. P ¼ 0.035424). Several possible explanations for the lack of concordance in identified genes exist. Technique-specific issues, including choice of microarray platform and randomisation of samples on arrays, may impede overlap in results (Draghici et al, 2006). Furthermore, differences in patient population exist between the two studies, for example, only high-grade tumours vs low and high-grade tumours in our study. An additional difficulty, inherent to microarray studies, is the use of small patient cohorts to evaluate large numbers of potential predictors of lymphocyte recruitment. This raises the likelihood of finding distinctive patterns based on chance rather than on biology, a phenomenon called overfitting (Fehrmann et al, 2007). Overfitting reduces the chances of finding overlap on the level of individual genes between studies, which can be countered by paying more attention to overlap in functional gene sets rather than individual genes (Crijns et al, 2009).
Despite these challenges, the expression of IRF-1, described above, and of CXCR6 was positively associated with the presence of CD8 þ TILs in both studies. Recently, it was established that both the chemokine receptor CXCR6 and its ligand CXCL16 are not only expressed by immune cells, but also by carcinomas (Meijer et al, 2008). Moreover, radiation was shown to recruit lymphocytes to carcinomas through the release of CXCL16 by tumour cells (Matsumura et al, 2008). Whether the expression of CXCR6 is similarly induced by radiation and/or chemotherapy and also influences lymphocyte attraction remains to be investigated.
Not only was there a discrepancy between identified genes between our study and that of Callahan et al (2008), we also found a substantial difference in genes and pathways associated with survival between the present and our previous study (Crijns et al, 2009). An important reason for this divergence is a difference in approach. Whereas the association with survival was the primary focus of our previous study, the current study was designed to discover genes and pathways that might be linked to lymphocyte recruitment. Only identified genes and pathways were    subsequently evaluated for an association with survival. A second, although related, reason is the fact that for the present study we were less stringent in the selection of genes (Po0.05), as we were primarily interested in associations with lymphocyte recruitment rather than survival, although for the 86-gene profile only genes meeting a Po0.001 were selected (Crijns et al, 2009). Thus, only two genes associated with survival in the current study were also part of the 86-gene profile, that is, BRSK1 and C1orf151. Neither gene has so far been further investigated to explain its prognostic effect and/or role in lymphocyte recruitment.
In summary, this study shows that gene expression profiling and pathway analysis are valuable strategies to obtain more insight into what tumour characteristics contribute to lymphocyte recruitment to advanced-stage serous ovarian carcinomas. Identified genes and pathways need to be further validated and evaluated for their value as a therapeutic target.