|Home | About | Journals | Submit | Contact Us | Français|
We used oligonucleotide microarrays to investigate gene expression changes associated with multi-step human papillomavirus type 16 (HPV16)-mediated carcinogenesis in vitro. Gene expression profiles in 4 early passage HPV16-immortalized human keratinocyte (HKc) lines derived from different donors were compared with their corresponding 4 late-passage, differentiation-resistant cell lines, and to 4 pools of normal HKc, each composed of 3 individual HKc strains, on Agilent 22 k human oligonucleotide microarrays. The resulting data were analyzed using a modified T-test coded in R to obtain lists of differentially expressed genes. Gene expression changes identified in this model system were then compared with gene expression changes described in published studies of cervical intraepithelial neoplasia (CIN) and cervical cancer. Common genes in these lists were further studied by cluster analysis. Genes whose expression changed in the same direction as in CIN or cervical cancer (concordant) at late stages of HPV16-mediated transformation in vitro formed one major cluster, while those that changed in the opposite direction (discordant) formed a second major cluster. Further annotation found that many discordant expression changes involved gene products with an extracellular localization. Two novel genes were selected for further study: overexpression of SIX1 and GDF15, observed during in vitro progression in our model system, was confirmed in tissue arrays of cervical cancer. These micro-array-based studies show that our in vitro model system reflects many cellular and molecular alterations characteristic of cervical cancer, and identified SIX1 and GDF15 as 2 novel potential bio-markers of cervical cancer progression.
Persistent high-risk human papillomavirus (HPV) infection is the major cause of cervical cancer. Once integrated into the host genome, high-risk HPV exert their oncogenic effects primarily through the continuous expression of the oncoproteins E6 and E7 (reviewed in Ref. 1). Many activities have been described for both E6 and E7, among which the following are best characterized and critical for transformation: E6 binds to E6-associated protein (E6-AP) resulting in the ubiquitination and degradation of p532; E7 binds to pocket protein family members, in particular, the retinoblastoma protein (Rb) causing inactivation and degradation of Rb (reviewed in Ref. 3). In the cervix, HPV-infected epithelial cells may appear as low-grade squamous intraepithelial lesions (LSIL), also referred to as cervical intraepithelial neoplasia (CIN) I. In the majority of cases, LSIL regresses spontaneously; however, 10 to 15% of LSIL may progress into high-grade squamous intraepithelial lesions (HSIL), also referred to as CIN II/III. Ultimately, only 0.3% of LSIL will progress into invasive cervical cancer.4–6
To study the molecular mechanisms of HPV-mediated transformation, we established an in vitro model system in which normal human keratinocytes (HKc), immortalized by transfection with HPV16 DNA (HKc/HPV16), progress toward malignancy through growth factor-independent (HKc/GFI) and differentiation-resistant (HKc/DR) stages.7–9 HKc/HPV16 reproduce the morphological features of dysplastic epithelium in vivo,10 and HKc/DR produce squamous cell carcinomas in nude mice, upon transfection with activated ras or Herpes Simplex Virus 2 DNA.11 Since growth factor-independent proliferation and loss of response to differentiation signals are acquired phenotypes of cancer cells, our in vitro model system may be useful to elucidate the molecular changes that underlie these phenotypes. Our previous work demonstrated that the HKc/GFI and HKc/DR phenotypes are at least partially due to increased expression and constitutive activation of the epidermal growth factor (EGF) receptor12,13 and loss of sensitivity to growth inhibition by transforming growth factor (TGF)-beta,9,14 respectively. Both EGF receptor overexpression and loss of sensitivity to growth inhibition by TGF-beta also occur in cervical cancer in vivo.15 To further explore changes that mark the progression of HPV16-transformed cells toward malignancy, we compared the gene expression profiles of HKc/HPV16 to normal HKc and of HKc/DR to HKc/HPV16. The resulting gene lists were then compared with literature-derived lists of gene expression changes described in CIN, cervical cancer,16–28 as well as in other HPV-transformed cell lines.16,29–33
Normal HKc from neonatal foreskin were maintained in complete serum-free MCDB153-LB medium as previously described.7,8,12 Normal HKc transfected with a plasmid containing a head-to-tail dimer of HPV16 DNA gave rise to immortalized HKc/HPV16, whose growth still requires EGF and bovine pituitary extract (BPE). HKc/GFI were selected by culturing HKc/HPV16 in MCDB153-LB medium lacking EGF and BPE. HKc/DR were then selected by maintaining HKc/GFI in complete medium containing calcium chloride (1 mM) and fetal bovine serum (5%).7,8,12 All cells were washed with Dulbecco’s phosphate buffered saline (PBS) and fed with complete serum-free MCDB153-LB medium 24 hr prior to RNA extraction.
Total RNA was extracted from normal HKc, HKc/HPV16 and HKc/DR using the Total RNA Isolation Mini Kit from Agilent Technologies (Santa Clara, CA) according to the manufacturer’s instructions. RNA quality and quantity were assessed on an Agilent 2100 Bioanalyzer. Only RNA samples with a RNA Integrity Number greater than 9.0 were used for subsequent amplification and labeling. mRNA (500 ng) was amplified and labeled with Cyanine 3 (Cy3) or Cyanine 5 (Cy5) using the Low RNA Input Linear Amplification Kit from Agilent Technologies. To minimize amplification bias, 2 aliquots of total RNA from each sample were individually amplified and labeled, then combined, fragmented and hybridized to the microarrays as described later.
Labeled amplified RNA was hybridized to the Human 1A (version 2) Oligoarray from Agilent Technologies (Santa Clara, CA) which contains 22,575 60-mer probes. Hybridization (17 hr at 60°C) was conducted in an Agilent DNA Microarray Hybridization Oven using Gene Expression Hybridization and Wash Buffer Kits from Agilent Technologies. Arrays were placed in Agilent Stabilization and Drying Solution prior to scanning on either an Agilent G2565AA scanner or a Perkin Elmer ScanArray Express HT scanner. Images were quantified using ImaGene software from BioDiscovery (El Segundo, CA).
A schematic representation of the experimental design is provided in Supplementary Figure 1a. RNA samples were obtained from 4 pools of normal HKc; each pool was composed of cells from 3 different and unique normal HKc strains; thus, 12 normal HKc strains were used in total. Amplified and labeled RNA samples from the normal HKc pools (referred to as N1, N2, N4 and N5, to match the numeric designation of the HKc/HPV16 and HKc/DR lines) were cohybridized to a set of arrays with amplified and labeled RNA from 4 independently-derived HKc/HPV16 lines, each established from a different donor: HKc/HPV16d-1 (referred to as IM1), HKc/HPV16d-2 (IM2), HKc/HPV16d-4 (IM4) and HKc/HPV16d-5 (IM5) at passage numbers ranging from 13 to 19 since primary culture. These lines have been previously described and extensively studied in our laboratory. 7–9,12–14,34 A second set of arrays was hybridized with amplified and labeled RNA from the 4 HKc/HPV16 lines and from their respective HKc/DR lines (referred to as DR1, DR2, DR4 and DR5) at passage levels near or above 150 postselection in medium containing serum and high calcium, with the exception of HKc/DR5 which was at passage 77 postselection. To test for dye bias and improve measurement accuracy, dye-swap replicates were performed for the N1/IM1, IM1/DR1 and the N5/IM5, IM5/DR5 series (supplementary Fig. 1a, double arrows).
Data normalization was performed by the global LOESS method prior to further analysis.35 The significance (Y) of change of signal intensity for a gene (g) is estimated by a modified t test, as follows:
where xgR and xgG are the mean log intensities of technical replicates of sample R and G, respectively. σg is the standard deviation of xgR − xgG and c is a significance cutoff (c = 2). σg is estimated from a pool of genes with similar expression level. We applied this method to the analysis of the normalized data from the microarrays to generate a list of significantly changed genes. The significance cutoff of the algorithm is shown in the MA plot in supplementary Figure 1b. The common changes among biological replicates were then identified. The implementation of this method is in R, based on Bioconductor (http://www.bioconductor.org).
To assess how well data from different stages of transformation, different individuals, and dye swap replicates correlated with each other, we devised an adaptive correlation method able to handle the intensity-dependent variance that affects the results of typical correlation methods. To obtain an adaptive correlation curve for data from 2 microarrays, intensity values from 4 measurements (two Cy3 and 2 Cy5) for each spot were averaged; the resulting average intensities were then sorted and divided into 50 bins containing equal numbers of genes. The correlation of the log2 of the Cy5/Cy3 ratios between the 2 microarrays within each bin was plotted in relation to the average intensities, and the points connected into a trend line smoothed by LOESS. supplementary Figure 2a shows the adaptive correlation curves of several paired comparisons: two dye-swap pairs (IM5/N5 versus dsIM5/N5 and DR5/IM5 versus dsDR5/IM5); paired comparisons of expression profiles at different stages of in vitro progression (IM5/N5, versus DR5/IM5); and paired comparisons between different individuals (IM1/N1 versus IM5/N5). With the exclusion of the low-intensity spots (log2 intensity < 6), the dye-swap replicates correlated with each other very well (adaptive correlation of about 0.9), indicating that dye bias was minimal. As expected, data from different stages of transformation or different individuals did not correlate well (adaptive correlations between 0.3 and 0.5, at average log2 intensities of 6 or greater). Increasing the bin number to 200 did not change the correlation trend lines (data not shown). Overall, the adaptive correlation curves showed that the differences in expression profiles are greatest between different transformation stages and minimal between technical replicates (dye-swap).
All cluster analysis was unsupervised using 1- correlation to generate the distance matrix, and clustered by average linkage.36 One-way clusters of cell line comparisons were obtained using all gene expression ratio measurements, while two-way clusters were obtained using subsets of gene expression ratio measurements selected based on specific criteria, i.e. early and late changes, which are also detected in cervical cancer.
Unsupervised one-way hierarchical cluster analysis was performed on the microarray data using all ratio measurements (supplementary Fig. 2b). The 4 comparisons between HKc/HPV16 (IM1, IM2, IM4 and IM5) and the normal HKc pools (N1, N2, N4 and N5) including their respective dye swaps, clustered together, and the 4 HKc/DR (DR1, DR2, DR4 and DR5) and their respective HKc/HPV16 lines comparisons and dye swaps formed a separate cluster.
To confirm the validity of our microarray data and methods of analysis, the mRNA levels for 9 genes (LTBP1, PLAU, CTGF, APC, TGFBI, FN1, INHBA, SIX1 and GDF15) were compared between the 4 HKc/HV16 lines and the 4 normal HKc pools and between the HKc/DR lines and their respective HKc/HPV16 lines by real-time RT-PCR. Selection of these genes for testing was guided by the following criteria: we chose examples of genes that were highly expressed, and some that had low levels of expression; we included genes whose expression changed in all 4 cell lines (i.e., SIX1, GDF15); some that changed in less than 4 cell lines (i.e., LTBP1); and 2 which did not change (CTG, APC). We included genes that were consistently up-regulated during progression (i.e., SIX1, GDF15) and others that were down-regulated, at least in some cell lines (i.e., PLAU.) cDNA synthesis was performed using the iScript™ cDNA Synthesis Kit from BioRad (Hercules, CA) following the manufacturer’s protocol. An aliquot (2.5%) of the cDNA synthesized from 0.5 μg of total RNA was used in the real-time PCR reaction, using iQ™ SYBR® Green Supermix (BioRad) per manufacturer’s instructions. Each PCR reaction was performed in triplicate. Serial dilutions of the cDNA product were subjected to real-time PCR amplification to estimate primer efficiency. Quantification of real-time PCR results was performed using the Pfaffl model.37 RT-PCR reactions for the gene GUSB were performed in parallel, and the results used for normalization of the real-time PCR data for the genes of interest.
The correlation between the real-time RT-PCR and the microarray data for all 9 genes, measured in all cell lines (R2 = 0.787), is shown in supplementary Figure 2c. All the significant changes (whether up or down) identified by microarray analysis for these 9 genes were confirmed by real-time RT-PCR. No correlation between microarray data and RT-PCR results was found when the change detected by microarray analysis fell outside of the significant range (data not shown).
Cervical cancer tissue arrays (CR802 for SIX1 staining, CR801 for GDF15 staining) were purchased from US Biomax (Ijamsville, MD). Each slide contains 80 cores consisting of cervical squamous cell carcinoma (40 cores) and adjacent normal cervical tissue (40 cores). Some of the normal adjacent cervical tissue cores (8 cores in the CR802 array and 4 cores in the CR801 array) did not contain epithelial tissue and therefore were excluded from the scoring of staining intensity. Immunohistochemistry was conducted in the Histology Core of the South Carolina Cancer Center. Briefly, tissue arrays were subjected to antigen retrieval (EDTA pH 9.0 at 98°C for 40 min for SIX1, or sodium citrate pH 6.0 at 98°C for 40 min for GDF15) followed by incubation with either rabbit anti-SIX1 antibody at a 1:100 dilution (HPA001893, Atlas Antibodies) for 60 min or goat anti-GDF15 antibody at a 1:100 dilution (sc-46250, Santa Cruz) for 45 min at room temperature. Slides were rinsed twice with Tris-Buffered Saline (TBS), incubated with either goat anti-rabbit antibody at a 1:500 dilution (sc-2040, Santa Cruz) or donkey anti-goat antibody at a 1:500 dilution (sc-2042, Santa Cruz) for 45 min at room temperature, followed by 2 TBS washes. Immunohistochemical staining was achieved using 3,3′ diaminobenzidine and hydrogen peroxide.
Three evaluators independently scored each core on the tissue array based on a staining intensity scale of 0 to 3, classifying each specimen’s staining as negative, weak, medium or strong. In most cases, at least 2 evaluators gave identical scores. In those few cases where there was disagreement over the score, samples were rescored, again independently, and the scores reconciled.
Direct, paired comparisons of gene expression were made between 4 independently-derived HKc/HPV16 and 4 pools of normal HKc strains, each derived from 3 separate and unique donors, and between 4 HKc/DR and their respective HKc/HPV16 lines, to identify genes whose expression was altered during in vitro progression of these cells toward a premalignant phenotype. A detailed description of the experimental design, the methods of data analysis and the results of assays performed to assess data quality are presented in the Materials and Methods section and in Supplementary Figures 1 and 2.
As a first approach to the analysis of the gene expression data, we annotated and analyzed the gene expression changes that occurred in all 4 individual cell lines, since these gene expression changes characterize in vitro progression in our model system (Fig. 1a). We divided changes into early stage, between HKc/HPV16 and normal HKc (supplementary Table I), and late stage, between HKc/DR and HKc/HPV16 (supplementary Table II). In all, we detected expression changes in 138 genes at the early stage, and 60 genes at the late stage. The expression of 4 genes, CRCT1, CRYAB, KLK5, KLK7, changed at the early stage and continued to change, in the same direction, at the late stage hence, these genes are listed with both the early and late stage changes in supplementary Tables I and II. An indirect comparison between HKc/DR and normal HKc revealed 34 additional genes changed in all 4 HKc/DR lines when compared with normal HKc (supplementary Table III). These genes were not identified by the direct comparisons because significant change was detected in less than 4 cell lines in the individual comparisons. However, the combined dataset allowed them to be identified as significantly changed in all 4 HKc/DR lines when compared with normal HKc. Genes were grouped based on their functional category, obtained from SOURCE, an online annotation tool (http://source.stanford.edu/cgi-bin/source/sourceSearch). Genes with multiple functions are classified into all appropriate groups, and therefore they appear multiple times in supplementary Tables I, II and III.
As expected, we detected mostly down-regulated expression of a group of immune defense-related genes (15 down-regulated and 1 up-regulated) among the early changes, and mostly down-regulated expression of epithelial structural proteins (20 down-regulated and 2 up-regulated) such as keratins 4, 13 and 24, among both early and late stage changes. Down-regulation of 7 keratin genes in HKc/DR correlates well with the differentiation-resistant phenotype of these cells. We find gene expression changes known to occur in cervical cancer, such as increased p16 and ICAM1, among the early changes.38,39 We also detected gene expression changes that have not been reported in the literature, but could be confirmed at the protein level by searching a high throughput immunohistochemistry staining database constructed on a panel of cancers (http//:www.hpr.se). These changes included increased TNFRSF8, increased UCHL1, reduced KRT4 (in HKc/HPV16 versus normal HKc), reduced KRT6B and increased SIX1 (in HKc/DR versus HKc/HPV16).
We also analyzed gene expression changes that occurred in at least 2 out of 4 individual cell lines and compiled a list of 1963 changes (Fig. 1b). To assess whether or not gene expression changes in our in vitro model system correspond to cervical carcinogenesis in vivo, we compared our list of 1963 genes whose expression changed in 2 out of 4 cell lines (2 out of 4 list, not shown) to a list of 1300 unique gene expression changes detected in CIN or in cervical cancer, compiled from 12 published microarray studies16–27 and 1 differential display study28 (Fig. 1b). In the resulting common set, we further divided gene expression changes into 2 groups, based on whether the direction of change was the same (concordant) or opposite (discordant) in our cell lines, as in CIN or cancer (Fig. 1c). Lists of concordant and discordant gene expression changes detected in at least 2 out of 4 individual HKc/HPV16 lines versus normal HKc (early stage changes), and at least 2 out of 4 HKc/DR versus HKc/HPV16 (late stage changes) are provided in supplementary Table IV. Within these lists, genes that change in expression levels in the same direction at both stages are classified as early plus late changes. Altogether, we detected 118 concordant gene expression changes and 65 discordant gene expression changes (Fig. 1c). In a group of 45 genes, expression changed at the early stage, but levels returned to normal or near normal at the late stage. These genes (Fig. 1c, IM-specific) are listed as a separate group in supplementary Table V, and were not included in further analysis for this paper. However, some of these genes may be of considerable interest because they appear to be specific to the earlier stages of HPV16-mediated transformation: among these, for example, is the EGF receptor, the RNA and protein levels of which (as we had described previously) increase in HKc/HPV16 and remain elevated in HKc/GFI, at the intermediate stage between HKc/HPV16 and HKc/DR12,13 but return to normal or near normal in HKc/DR.
Unsupervised 2-way cluster analysis was performed on the concordant and discordant changes of gene expression at early and late stages, and the 2 resulting heatmaps are shown in Figure 2. Gene expression changes detected at both early and late stages are included in the early-stage cluster (Fig. 2, early). To the left of each heatmap, a color column shows the concordance (blue) or discordance (yellow) of the observed changes with the gene expression changes reported in cervical cancer. An interesting pattern emerged from this analysis: In the 2 major clusters of gene expression changes detected at the late stage (Fig. 2, late) most of the discordant changes were in the upper cluster, while most concordant changes were in the lower cluster. There are no obvious concordance patterns in the cluster of gene expression changes detected at early stage (Fig. 2, early).
The clustering of late changes according to concordance with cervical cancer prompted us to further explore the concordant and discordant changes of gene expression at early and late stages, in the context of cellular localization of gene products, molecular processes, and biological processes, by performing a Gene Ontology analysis using the software package Ermine J. This analysis revealed that extracellular localization was overrepresented in the discordant late-stage expression changes, which clustered together in our unsupervised cluster analysis. Utilizing cellular localization information gathered through SOURCE (http://source.stanford.edu/cgi-bin/source/sourceSearch) from Stanford University, and Pathway Studio from Ariadne Genomics (Rockville, MD), we found that, while only about 10% (6 out of 62) of the concordant changes affected genes with extracellular products (Fig. 3a), about 28% (13 out of 46) of the discordant late-stage expression changes affected genes that encoded proteins with extracellular localization (Fig. 3b). Gene expression changes that occurred in both early and late stages are included in both graphs but colored differently. No such pattern of subcellular localization was detected at early stages (data not shown). The fact that a significant portion of discordant late-stage changes localized to the extracellular space indicated that some of the differences between HKc/DR and cervical cancer cells in vivo may be due to the absence of epithelial-stroma interactions in monolayer culture.
To assess the degree of similarity of our in vitro model system to other HPV-transformed cell lines, we compared our list of 1963 genes changed in 2 out of 4 cell lines, to a list of 413 unique gene expression changes described in other HPV-transformed epithelial cell lines, compiled from 6 published microarray studies16,29–33 (Figs. 1b and 1d). Except for the study by Pett et al.,33 which used Affymetrix arrays, all other previous studies used different array platforms with fewer genes than the Agilent arrays used in this study. Therefore, many concordant changes may have gone undetected. Among the 413 unique gene expression changes reported by others in HPV-transformed cell lines, 120 could be detected in at least 2 out of 4 of our cell lines, with 87 of these changing in the same direction (Fig. 1d). The gene list resulting from the comparison between our cell lines and other HPV-transformed cell lines is shown in supplementary Table IV.
Among the genes significantly over-expressed in HKc/DR compared to HKc/HPV16, the homeobox gene SIX1 is of special interest because it has been reported to be over-expressed in breast cancer and rhabdomyosarcoma, in association with a metastatic phenotype.40–42 To determine whether SIX1 was overexpressed in cervical cancer specimens, we performed immunohistochemical staining of cervical cancer tissue arrays and obtained scores of the staining intensity from 3 independent evaluators. Staining for SIX1, which is a transcription factor, was mainly localized in the nucleus in both cancer and adjacent normal cells. While none of the adjacent normal tissue cores exhibited strong SIX1 staining (Fig. 4a) strong SIX1 staining was detected in 23% of the cervical cancer tissue cores (Fig. 4b).
In addition, we observed an increase in GDF15 mRNA in all 4 HKc/DR lines compared with their respective HKc/HPV16 lines and the normal HKc pools. GDF15, also called MIC 1, or placental TGF-beta, is up-regulated in breast, prostate and colon cancer.43 We performed immunohistochemical staining for GDF15 in a cervical cancer tissue array. GDF15 was detectable in all specimens. In adjacent normal tissue, GFD15 was localized to the cytoplasm, in basal and immediately supra-basal layers of the epithelium (Fig. 4c). In cancer tissue, uniform cytoplasmic staining was detected throughout the epithelium (Fig. 4d). We detected more often medium intensity staining and less often weak staining in cervical cancer tissue, compared to adjacent normal tissue. Strong GDF15 staining was only detected in cervical cancer specimens (13%) (Fig. 4d).
One of the most important questions this work begins to address is: How well does our in vitro model system for HPV16-mediated carcinogenesis recapitulate the molecular changes that underlie the development of cervical cancer? Our model cell lines clearly possess many of the phenotypic characteristics of cancer cells, including a progressive acquisition of autonomous growth, loss of differentiation, overexpression of the EGF receptor, resistance to growth inhibition by TGF-beta, and other transformed cell phenotypes, the molecular basis of which we have explored and documented over the years.7–9,12–14,34 However, HPV-transformed cells in vivo are exposed to selective pressures that are different from those which we applied in vitro in our system, and therefore may very well take different paths to progression. Despite this, we believed that it may still be possible to demonstrate common pathways through which progression takes place, and therefore establish and define the suitability of our model system as a research tool to study cancer progression.
We utilized DNA microarray technology to perform a systematic comparison of the gene expression changes that characterize 2 key stages of in vitro progression of HPV16-immortalized HKc, to gene expression changes that have been documented in cervical cancer. It was encouraging to observe that unsupervised clustering of the data neatly groups together cell lines according to their status (immortalized or DR). Also, the finding that many expected gene expression changes were confirmed in our dataset strengthened our conclusion that we had a robust and valid dataset to rely upon. As expected, there was more variability in gene expression profiles among the HKc/DR than among HKc/HPV16 at an earlier stage of progression. This too was encouraging, as HKc/DR have likely accumulated more gene expression changes, including both changes underlying their common differentiation resistant phenotype, and also individual changes, specific to each cell line, thus diverging more extensively from the initial relatively uniform gene expression profiles.
For a first round of analysis, we compiled a list of all the genes that change in expression level during progression in all 4 HKc/DR lines, when compared with HKc/HPV16, and in all 4 HKc/HPV16 lines, when compared with normal HKc. This gene list contains changes that are common with 4 entirely independent transformation and progression events, on the different genetic background of 4 individual HKc strains. Therefore, these lists are likely to include a subset of genes representing the common denominator of in vitro transformation and progression in this model system. Next, we compiled a list of the early and late changes in gene expression detectable in at least 2 out of 4 cell lines, and compared this second gene list to a list of changes described in CIN and cervical cancer.
A first glance at the functional information on genes included in these lists revealed several aspects known to be associated with HPV-mediated transformation, including down-regulation of interferon-inducible genes. The general picture, suggesting a suppression of interferon responses and immunomodulatory activities in HPV16-transformed cells, is in line with a body of work showing down-regulation of immune responses by oncogenic HPV (reviewed in ref. 44). A recent microarray study on epithelial cells and fibroblasts collected by laser capture microdissection of normal cervical tissue, CIN 1 to 3 tissue specimens and invasive cervical cancer proposed that CIN1 is characterized by a “pro-proliferative/immunosuppressive gene expression signature,”27 which is also present in our in vitro model system. For example, we observed the over-expression of several proliferation related genes (including MCM2, MCM4, CDC7, TOP2A, and the EGFR) among the early stage changes present in our model system (supplementary Table IV).
We have also observed over-expression of proliferation-related genes in HKc/DR, such as over-expression of centromere protein F, cyclin E1, H2A, and histone family member X. In the case of HKc/DR, a high proliferating rate, a characteristic selected for by growth conditions in vitro, may recapitulate the molecular features of a sub-population of cancer cells. In fact, when we compared our 2 out of 4 list of gene expression changes also found during cervical carcinogenesis in vivo (supplementary Table IV) with a “proliferation cluster” identified in a subgroup of cervical cancers characterized by a poor prognosis,24 we found that 34 out of the 35 common genes changed in a concordant fashion. This high concordance rate indicated that our cell line model might recapitulate, at the late stage, the high proliferating feature of this sub group of cervical cancers.
A major feature of HKc/HPV16 and HKc/DR is a disruption of the differentiation program. Accordingly, a series of differentiation-related genes are modulated at the early and late stages of HPV16-mediated transformation, including several keratin genes and other structural molecules expressed in differentiating keratinocytes. In addition, we observed down-regulation of a group of kallikreins (kallikrein 5, 6, 7, 10 and 11)45 among which kallikrein5 is specific to skin and plays a role in epithelial cell differentiation,46 and kallikrein10 is found profoundly down-regulated in breast cancer cells47: these changes are potentially of considerable interest under a functional point of view. While a general decrease in expression of kallikreins is an early event, down-regulation of the expression of selected SERPIN molecules (SERPINB1, B2, B5, E2) is mainly a late event. SERPINB5 was classified as a tumor suppressor molecule, expressed in normal mammary epithelium but lost in breast cancer.48 These and other leads warrant further investigation.
Our microarray results also identified several individual genes that may serve as potential biomarkers for cervical cancer: for example, we discovered increased expression of SIX1 and GDF15 in HKc/DR and confirmed over-expression of the products of these 2 genes in cervical cancer tissue arrays.
We also observed a 73% concordance between the 2 of 4 list and the list of genes changed in HPV transformed cell lines in vitro.16,29–33 Several studies comparing gene expression profiles of tumor tissue and cultured cell lines from various sites, including breast, colon and lung cancer, found that cell lines clustered together, and their profiles were distinct from those of their corresponding tumors.49–52 Interestingly, the mRNA levels of 5 (IFI27, S100A9, MMP9, DSP, and FGFR2) out of the 15 genes that are common among the cervical cancer list, the other HPV-transformed cell line list, and the 2 out of 4 list changed in opposite directions in cervical cancer and in HPV-transformed cell lines. For example, IFI27 had been found to be reduced in HPV transfected cells29,30 and in our study, while it is reported as over-expressed in cervical cancer.18 Thus, reduced expression of IFI27 is probably a feature of HPV-transformed cell lines compared to normal cells in vitro, but not in vivo.
The lack of stroma in the in vitro model system may contribute to the discordance of gene expression profiles between cell lines and cancer tissue. Stroma interactions have been found to be essential for carcinogenesis in vivo. A recent microarray study on laser capture microdissected cervical tissue found that the transition to CIN2 coincides with the activation of pro-angiogenesis pathways in fibroblasts, and the transition to CIN3 and squamous cell carcinoma is characterized by a proinvasive gene expression signature, both attributed to the demands of an environment in which nutrients become limiting.27 We have not observed the same proangiogenesis or a proinvasive signature in our HKc/DR cells, probably because of the abundance of nutrients in the culture medium. Instead of pro-invasive features, we observe a down regulation of MMPs (MMP1, MMP9, MMP10). Previous studies found that cancer cell lines in culture in general down-regulate cell adhesion molecules and membrane signaling proteins, compared to tumor and normal tissue.53 The observation that many of the discordant gene expression changes in HKc/DR involve genes with extracellular products has motivated us to begin studies exploring gene expression profiles of our cells in organotypic culture.
Finally, we must observe that the list of genes whose expression is changed in cervical cancer is far from being complete, and that the search and comparison task has just begun: in this first phase of study, we constructed a literature-derived database of gene expression changes detected in cervical cancer, and relied on the resulting gene list for our comparison. This is both a strength (we relied on a host of well documented information) and a limitation (we do not have a direct comparison, made with the same tools) of this study. In addition, while all of our cell lines were produced with HPV16, we have not taken into account the specific HPV types represented in the cervical cancer gene expression database, and their potential differential effects on gene expression.
We conclude that HKc/HPV16 and HKc/DR in monolayer culture recapitulate important molecular features of HPV-transformed cells in vivo, and that some of the discordant gene expression changes detected in HKc/DR may be attributed to cell-micro-environment interactions, rather than to intrinsic properties of the transformed cells. The details of these influences and interactions warrant further investigation.
This work was performed with the support of the DNA Microarray Facility and the Histology Core Facility of the South Carolina Cancer Center. Grant support: National Institutes of Health grants R01CA89502, P20CA096427 and 5P20MD001770-030002 (Dr. K.E. Creek) and 5P20RR016461 (Dr. L. Pirisi).
Grant sponsor: National Institutes of Health; Grant numbers: R01CA89502, P20CA096427, 5P20MD001770-030002, 5P20RR016461.
This article contains supplementary material available via the Internet at http://www.interscience.wiley.com/jpages/0020-7136/suppmat