Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Int J Cancer. Author manuscript; available in PMC 2010 May 18.
Published in final edited form as:
PMCID: PMC2872618

Gene expression changes during HPV-mediated carcinogenesis: A comparison between an in vitro cell model and cervical cancer


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.

Keywords: HPV, papillomavirus, cervical cancer, microarray, in vitro model

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.46

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.79 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,1628 as well as in other HPV-transformed cell lines.16,2933

Material and methods

Cell lines and cell culture conditions

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.

RNA isolation, amplification and labeling

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.

Microarray hybridization, scanning and spot quantitation

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).

Overall design of the microarray studies

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. 79,1214,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 and identification of significantly changed genes

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 xgRxgG 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 (

Analysis of microarray data

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).

Cluster analysis

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.

Cluster analysis of all gene expression changes

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.

Real-time RT-PCR analysis of mRNA levels for selected genes

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).

Immunohistochemistry for SIX1 and GDF15

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.

Scoring of staining intensity

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.


Alterations of gene expression during in vitro progression of HPV16-immortalized HKc

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 ( Genes with multiple functions are classified into all appropriate groups, and therefore they appear multiple times in supplementary Tables I, II and III.

Figure 1
(a) composition of 232 gene expression changes detected in 4 out of 4 individuals in our in vitro model system. Early: HKc/HPV16 compared with normal HKc; Late: HKc/DR compared with HKc/HPV16; DR_N, changes between HKc/DR and normal HKc, which were derived ...

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// 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).

Comparison of gene expression changes in our model system with changes detected in cervical cancer

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 studies1627 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).

Figure 2
Unsupervised two-way cluster analyses of early and late gene expression changes that overlap with changes detected in cervical cancer or CIN. Cluster to the left of the heat map, gene expression ratios; cluster above the heat map, comparison among cell ...

Subcellular localization of products of concordant and discordant genes

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 ( 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.

Figure 3
Cellular localization maps of gene products differentially expressed in HKc/DR versus HKc/HPV16. Plasma membranes are represented by a brown band separating the extracellular space from the intracellular space, and nuclei are depicted as light blue ellipses. ...

Comparison of gene expression profiles of our in vitro model system with other HPV-transformed cell lines

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,2933 (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.

Increased levels of SIX1 and GDF15 protein in cervical cancer tissue arrays corroborate microarray results showing overexpression of these 2 genes in HKc/DR

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.4042 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).

Figure 4
Representative immunohistochemical staining for SIX1 and GDF15 in cervical cancer and adjacent normal cervical epithelial tissue, from tissue arrays. (a) SIX1 staining (brown) in a normal cervical specimen (adjacent to cancer tissue). The pie chart presents ...

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.79,1214,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,2933 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.4952 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.

Supplementary Material

Supplementary Figure 2

Supplementary Table 1

Supplementary Table 2

Supplementary Table 3

Supplementary Table 4

Supplementary Table 5

Supplementary Table 6

Supplemetary Figure 1


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.



1. Snijders PJF, Steenbergen RDM, Heideman DAM, Meijer CJLM. HPV-mediated cervical carcinogenesis: concepts and clinical implications. J Pathol. 2006;208:152–64. [PubMed]
2. Werness BA, Levine AJ, Howley PM. Association of human papillomavirus types 16 and 18 E6 proteins with p53. Science. 1990;248:76–9. [PubMed]
3. Dyson N, Howley PM, Munger K, Harlow E. The human papilloma virus-16 E7 oncoprotein is able to bind to the retinoblastoma gene product. Science. 1989;243:934–7. [PubMed]
4. Melnikow J, Nuovo J, Willan AR, Chan BK, Howell LP. Natural history of cervical squamous intraepithelial lesions: a meta-analysis. Obstet Gynecol. 1998;92:727–35. [PubMed]
5. Ostor AG. Natural history of cervical intraepithelial neoplasia: a critical review. Int J Gynecol Pathol. 1993;12:186–92. [PubMed]
6. Cantor SB, Atkinson EN, Cardenas-Turanzas M, Benedet JL, Follen M, MacAulay C. Natural history of cervical intraepithelial neoplasia: a meta-analysis. Acta Cytol. 2005;49:405–15. [PubMed]
7. Pirisi L, Yasumoto S, Feller M, Doniger J, DiPaolo JA. Transformation of human fibroblasts and keratinocytes with human papillomavirus type 16 DNA. J Virol. 1987;61:1061–6. [PMC free article] [PubMed]
8. Pirisi L, Creek KE, Doniger J, DiPaolo JA. Continuous cell lines with altered growth and differentiation properties originate after transfection of human keratinocytes with human papillomavirus type 16 DNA. Carcinogenesis. 1988;9:1573–9. [PubMed]
9. Creek KE, Geslani G, Batova A, Pirisi L. Progressive loss of sensitivity to growth control by retinoic acid and transforming growth factor-β at late stages of human papillomavirus type 16-initiated transformation of human keratinocytes. Adv Exp Med Biol. 1995;375:117–35. [PubMed]
10. Woodworth CD, Waggoner S, Barnes W, Stoler MH, DiPaolo JA. Human cervical and foreskin epithelial cells immortalized by human papillomavirus DNAs exhibit dysplastic differentiation in vivo. Cancer Res. 1990;50:3709–15. [PubMed]
11. DiPaolo JA, Woodworth CD, Popescu NC, Koval DL, Lopez JV, Doniger J. HSV-2-induced tumorigenicity in HPV16-immortalized human genital keratinocytes. Virology. 1990;177:777–9. [PubMed]
12. Zyzak LL, MacDonald LM, Batova A, Forand R, Creek KE, Pirisi L. Increased levels and constitutive tyrosine phosphorylation of the epidermal growth factor receptor contribute to autonomous growth of human papillomavirus type 16 immortalized human keratinocytes. Cell Growth Differ. 1994;5:537–47. [PubMed]
13. Akerman GS, Tolleson WH, Brown KL, Zyzak LL, Mourateva E, Engin TS, Basaraba A, Coker AL, Creek KE, Pirisi L. Human papillomavirus type 16 E6 and E7 cooperate to increase epidermal growth factor receptor (EGFR) mRNA levels, overcoming mechanisms by which excessive EGFR signaling shortens the life span of normal human keratinocytes. Cancer Res. 2001;61:3837–43. [PubMed]
14. Mi Y, Borger DR, Fernandes PR, Pirisi L, Creek KE. Loss of transforming growth factor-β (TGF-β) receptor type I mediates TGF-β resistance in human papillomavirus type 16-transformed human keratinocytes at late stages of in vitro progression. Virology. 2000;270:408–16. [PubMed]
15. Kim JW, Kim YT, Kim DK, Song CH, Lee JW. Expression of epidermal growth factor receptor in carcinoma of the cervix. Gynecol Oncol. 1996;60:283–7. [PubMed]
16. Vazquez-Ortiz G, Pina-Sanchez P, Vazquez K, Duenas A, Taja L, Mendoza P, Garcia JA, Salcedo M. Overexpression of cathepsin F, matrix metalloproteinases 11 and 12 in cervical cancer. BMC Cancer. 2005;5:68. [PMC free article] [PubMed]
17. Cheng Q, Lau WM, Tay SK, Chew SH, Ho TH, Hui KM. Identification and characterization of genes involved in the carcinogenesis of human squamous cell cervical carcinoma. Int J Cancer. 2002;98:419–26. [PubMed]
18. Chen Y, Miller C, Mosher R, Zhao X, Deeds J, Morrissey M, Bryant B, Yang D, Meyer R, Cronin F, Gostout BS, Smith-McCune K, et al. Identification of cervical cancer markers by cDNA and tissue microarrays. Cancer Res. 2003;63:1927–35. [PubMed]
19. Wong YF, Selvanayagam ZE, Wei N, Porter J, Vittal R, Hu R, Lin Y, Liao J, Shih JW, Cheung TH, Lo KW, Yim SF, et al. Expression genomics of cervical cancer: molecular classification and prediction of radiotherapy response by DNA microarray. Clin Cancer Res. 2003;9:5486–92. [PubMed]
20. Ahn WS, Bae SM, Lee JM, Namkoong SE, Han SJ, Cho YL, Nam GH, Seo JS, Kim CK, Kim YW. Searching for pathogenic gene functions to cervical cancer. Gynecol Oncol. 2004;93:41–8. [PubMed]
21. Sopov I, Sörensen T, Magbagbeolu M, Jansen L, Beer K, Kühne-Heid R, Kirchmayr R, Schneider A, Dürst M. Detection of cancer-related gene expression profiles in severe cervical neoplasia. Int J Cancer. 2004;112:33–43. [PubMed]
22. Hudelist G, Czerwenka K, Singer C, Pischinger K, Kubista E, Manavi M. cDNA array analysis of cytobrush-collected normal and malignant cervical epithelial cells: a feasibility study. Cancer Genet Cytogenet. 2005;158:35–42. [PubMed]
23. Steinau M, Lee DR, Rajeevan MS, Vernon SD, Ruffin MT, Unger ER. Gene expression profile of cervical tissue compared to exfoliated cells: impact on biomarker discovery. BMC Genomics. 2005;6:64. [PMC free article] [PubMed]
24. Rosty C, Sheffer M, Tsafrir D, Stransky N, Tsafrir I, Peter M, de Crémoux P, de La Rochefordière A, Salmon R, Dorval T, Thiery JP, Couturier J, et al. Identification of a proliferation gene cluster associated with HPV E6/E7 expression level and viral DNA load in invasive cervical carcinoma. Oncogene. 2005;24:7094–104. [PubMed]
25. Wong YF, Cheung TH, Tsao GS, Lo KW, Yim SF, Wang VW, Heung MM, Chan SC, Chan LK, Ho TW, Wong KW, Li C, et al. Genome-wide gene expression profiling of cervical cancer in Hong Kong women by oligonucleotide microarray. Int J Cancer. 2006;118:2461–9. [PubMed]
26. Chao A, Wang TH, Lee YS, Hsueh S, Chao AS, Chang TC, Kung WH, Huang SL, Chao FY, Wei ML, Lai CH. Molecular characterization of adenocarcinoma and squamous carcinoma of the uterine cervix using microarray analysis of gene expression. Int J Cancer. 2006;119:91–8. [PubMed]
27. Gius D, Funk MC, Chuang EY, Feng S, Huettner PC, Nguyen L, Bradbury CM, Mishra M, Gao S, Buttin BM, Cohn DE, Powell MA, et al. Profiling microdissected epithelium and stroma to model genomic signatures for cervical carcinogenesis accommodating for covariates. Cancer Res. 2007;67:7113–23. [PubMed]
28. Seo MJ, Bae SM, Kim YW, Kim YW, Hur SY, Ro DY, Lee JM, Namkoong SE, Kim CK, Ahn WS. New approaches to pathogenic gene function discovery with human squamous cell cervical carcinoma by gene ontology. Gynecol Oncol. 2005;96:621–9. [PubMed]
29. Chang YE, Laimins LA. Microarray analysis identifies interferon-inducible genes and Stat-1 as major transcriptional targets of human papillomavirus type 31. J Virol. 2000;74:4174–82. [PMC free article] [PubMed]
30. Nees M, Geoghegan JM, Hyman T, Frank S, Miller L, Woodworth CD. Papillomavirus type 16 oncogenes downregulate expression of interferon-responsive genes and upregulate proliferation-associated and NF-κB-responsive genes in cervical keratinocytes. J Virol. 2001;75:4283–96. [PMC free article] [PubMed]
31. Ruutu M, Peitsaro P, Johansson B, Syrjanen S. Transcriptional profiling of a human papillomavirus 33-positive squamous epithelial cell line which acquired a selective growth advantage after viral integration. Int J Cancer. 2002;100:318–26. [PubMed]
32. Duffy CL, Phillips SL, Klingelhutz AJ. Microarray analysis identifies differentiation-associated genes regulated by human papillomavirus type 16 E6. Virology. 2003;314:196–205. [PubMed]
33. Pett MR, Herdman MT, Palmer RD, Yeo GSH, Shivji MK, Stanley MA, Coleman N. Selection of cervical keratinocytes containing integrated HPV16 associates with episome loss and an endogenous antiviral response. Proc Natl Acad Sci USA. 2006;103:3822–7. [PubMed]
34. Borger DR, Mi Y, Geslani G, Zyzak LL, Batova A, Engin TS, Pirisi L, Creek KE. Retinoic acid resistance at late stages of human papillomavirus type 16-mediated transformation of human keratinocytes arises despite intact retinoid signaling and is due to a loss of sensitivity to transforming growth factor-β Virology. 2000;270:397–407. [PubMed]
35. Yang YH, Dudoit S, Luu P, Lin DM, Peng V, Ngai J, Speed TP. Normalization for cDNA microarray data: a robust composite method addressing single and multiple slide systematic variation. Nucleic Acids Res. 2002;30:e15. [PMC free article] [PubMed]
36. Datta S, Datta S. Evaluation of clustering algorithms for gene expression data. BMC Bioinformatics. 2006;7(Suppl 4):S17. [PMC free article] [PubMed]
37. Pfaffl MW. A new mathematical model for relative quantification in real-time RT-PCR. Nucleic Acids Res. 2001;29:e45. [PMC free article] [PubMed]
38. Coleman N, Stanley MA. Characterization and functional analysis of the expression of vascular adhesion molecules in human papillomavirus-related disease of the cervix. Cancer. 1994;74:884–92. [PubMed]
39. Nasu K, Narahara H, Etoh Y, Kawano Y, Hirota Y, Miyakawa I. Serum levels of soluble intercellular adhesion molecule-1 (ICAM-1) and the expression of ICAM-1 mRNA in uterine cervical cancer. Gynecol Oncol. 1997;65:304–8. [PubMed]
40. Ford HL, Kabingu EN, Bump EA, Mutter GL, Pardee AB. Abrogation of the G2 cell cycle checkpoint associated with overexpression of HSIX1: a possible mechanism of breast carcinogenesis. Proc Natl Acad Sci USA. 1998;95:12608–13. [PubMed]
41. Reichenberger KJ, Coletta RD, Schulte AP, Varella-Garcia M, Ford HL. Gene amplification is a mechanism of Six1 overexpression in breast cancer. Cancer Res. 2005;65:2668–75. [PubMed]
42. Yu Y, Davicioni E, Triche TJ, Merlino G. The homeoprotein six1 transcriptionally activates multiple protumorigenic genes but requires ezrin to promote metastasis. Cancer Res. 2006;66:1982–9. [PubMed]
43. Welsh JB, Sapinoso LM, Kern SG, Brown DA, Liu T, Bauskin AR, Ward RL, Hawkins NJ, Quinn DI, Russell PJ, Sutherland RL, Breit SN, et al. Large-scale delineation of secreted protein biomarkers overexpressed in cancer tissue and serum. Proc Natl Acad Sci USA. 2003;100:3410–5. [PubMed]
44. Stanley M. Immune responses to human papillomavirus. Vaccine. 2006;24 (Suppl 1):S16–S22. [PubMed]
45. Yousef GM, Diamandis EP. The expanded human kallikrein gene family: locus characterization and molecular cloning of a new member. KLK-L3 (KLK9) Genomics. 2000;65:184–94. [PubMed]
46. Hansson L, Stromqvist M, Backman A, Wallbrandt P, Carlstein A, Egelrud T. Cloning, expression, and characterization of stratum corneum chymotryptic enzyme. A skin-specific human serine proteinase. J Biol Chem. 1994;269:19420–6. [PubMed]
47. Liu XL, Wazer DE, Watanabe K, Band V. Identification of a novel serine protease-like gene, the expression of which is down-regulated during breast cancer progression. Cancer Res. 1996;56:3371–9. [PubMed]
48. Sager R, Sheng S, Pemberton P, Hendrix MJ, Maspin A tumor suppressing serpin. Adv Exp Med Biol. 1997;425:77–88. [PubMed]
49. Alon U, Barkai N, Notterman DA, Gish K, Ybarra S, Mack D, Levine AJ. Broad patterns of gene expression revealed by clustering analysis of tumor and normal colon tissues probed by oligonucleotide arrays. Proc Natl Acad Sci USA. 1999;96:6745–50. [PubMed]
50. Perou CM, Jeffrey SS, van de Rijn M, Rees CA, Eisen MB, Ross DT, Pergamenschikov A, Williams CF, Zhu SX, Lee JC, Lashkari D, Shalon D, et al. Distinctive gene expression patterns in human mammary epithelial cells and breast cancers. Proc Natl Acad Sci USA. 1999;96:9212–17. [PubMed]
51. Alizadeh AA, Eisen MB, Davis RE, Ma C, Lossos IS, Rosenwald A, Boldrick JC, Sabet H, Tran T, Yu X, Powell JI, Yang L, et al. Distinct types of diffuse large B-cell lymphoma identified by gene expression profiling. Nature. 2000;403:503–11. [PubMed]
52. Virtanen C, Ishikawa Y, Honjoh D, Kimura M, Shimane M, Miyoshi T, Nomura H, Jones MH. Integrated classification of lung tumors and cell lines by expression profiling. Proc Natl Acad Sci USA. 2002;99:12357–62. [PubMed]
53. Birgersdotter A, Sandberg R, Ernberg I. Gene expression perturbation in vitro-a growing case for three-dimensional (3D) culture systems. Semin Cancer Biol. 2005;15:405–12. [PubMed]