The K14E7 and FVB mouse models have been previously described (4). Briefly, mice were obtained from the mouse repository located at National Cancer Institute at Frederick campus (Frederick, MD, USA). The K14E7 transgenic mice were generated and maintained in the FVB/N inbred strain and used as transgenic heterozygotes harbouring overlapping HPV16 E6 and E7 open reading frames (ORFs) spanning nucleotides 79-883 with a TTL in the E6 gene to preclude E6 expression. Persistent oestrogen administration and E7 expression were sufficient to produce high-grade cervical dysplasia and invasive cervical malignancies (17). For this study, four groups were created: FVB (not transgenic, not treated, control group), FVB+E₂ (not transgenic, treated group), K14E7 (transgenic, untreated group) and K14E7+E₂ (transgenic, treated group). The mice were housed and treated according to guidelines from the Association for Assessment and Accreditation of Laboratory Animal Care International and the National Institutes of Health Guide for the Care and Use of Laboratory Animals (publication no. 8023, revised 1978). The euthanasia procedure via use of a CO₂ chamber was followed according to the American Veterinary Medical Association's Guidelines for the Euthanasia of Animals (2013 edition).

All experiments and procedures were approved by the Research Unit for Laboratory Animal Care Committee (NOM-062-ZOO-1999; Unit for the Production and Experimentation of Laboratory Animals, The Centre for Research and Advanced Studies of the National Polytechnic Institute, México City, México). The environmental conditions for the mice included free access to food/water with a standard diet of LabDiet^®, a 12:12 light:dark cycle, a room temperature of 25°C and 30% humidity, with housing in Super Mouse 750™ Micro-Isolator^® cages in groups of 4–5 mice. Female, 1-month-old, virgin transgenic and non-transgenic mice, weighing 12–13 g, were dorsally implanted with continuous release pellets under the skin, which delivered 0.05 mg E₂ (Innovative Research of America, Sarasota, FL, USA) over 60 days. After 1 month of treatment, the mice, weighing 13-15 g, were sacrificed at 2 months of age. Control and treatment groups were formed using 20 to 24 mice in each.

Untreated or treated K14E7 hemizygous and non-transgenic FVB control virgin female mice were sacrificed by CO₂ chamber. All specimens were fixed in 4% paraformaldehyde at 4°C overnight and embedded in paraffin. Sections were deparaffinized and rehydrated as previously described (18). Serial sections were cut (at 5-μm thickness) and stained with haematoxylin and eosin (H&E) at room temperature (10 min haematoxylin and 2 min eosin).

Untreated or treated K14E7 hemizygous and non-transgenic FVB control virgin female mice were sacrificed by CO₂ chamber. Cervical biopsy samples were immediately stored in RNAlater^® solution (Thermo Fisher Scientific, Inc., Waltham, MA, USA) at 4°C overnight. Tissues were recovered from RNAlater solution with sterile forceps, quickly blotted to remove excess solution and immediately snap-frozen in liquid nitrogen. Total RNA was extracted from snap-frozen tissues using standard procedures with TRIzol reagent (Invitrogen; Thermo Fisher Scientific, Inc.). In total, 2 pooled samples from 3 mice were used to avoid variability and to identify biomarkers for the early stages of carcinogenesis in a population of K14E7 mice treated with E₂, as recommended by previous studies (19–21). Total RNA collected from 6 female mice from each group was pooled. RNA quantity and quality were assessed on an Agilent 2100 Bioanalyzer (Agilent Technologies, Inc., Santa Clara, CA, USA). RNA samples with an RNA integrity score of >8.0 were further processed for microarray analysis. cDNA synthesis, amplification and gene expression profiling were performed according to the manufacturer's protocols (Affymetrix^® WT Sense Target Labelling assay; Thermo Fisher Scientific, Inc.). Affymetrix microarray experiments were performed in triplicate for each group in the GeneChip Mouse Gene 1.0 ST Array platform.

Signal intensities from each array were analysed using Partek Genomic Suite version 6.4 (Partek, Inc., Chesterfield, MO, USA). Raw intensity values for the probes were normalized using robust multiarray analysis background correction. Two-way analysis of variance (ANOVA) was performed to identify differentially expressed genes. Genes with statistically significant differences in expression levels (P<0.05) and a fold change of ≥1.5 and ≤−1.5 were included in the final set of differentially expressed genes. To identify the biological processes altered by E7 and/or E₂, Ingenuity Pathway Analysis (IPA; https://www.qiagenbioinformatics.com/products/ingenuity-pathway-analysis/; Ingenuity Systems, Redwood City, CA, USA), a bioinformatics tool for visualizing expression data in the context of Kyoto Encyclopaedia of Genes and Genomes-defined biological pathways (http://www.genome.jp/kegg/) was used.

RNA isolated by TRIzol reagent (Invitrogen; Thermo Fisher Scientific, Inc.) from 6 mice from each group (K14E7, FVB, K14E7+E₂ and FVB+E₂) was purified, and RNA quality determined by electrophoresis on a 2% agarose gel, with visualization of ribosomal RNA by ethidium bromide staining. RNA was quantified by spectrophotometric analysis at 260 and 280 nm. cDNA synthesis was performed using the Maxima First Strand cDNA Synthesis kit according to the manufacturer's protocols (Thermo Fisher Scientific, Inc.). RT-qPCR was performed using an Applied Biosystems 7300 instrument and a DNA Master SYBR Green I kit (both Thermo Fisher Scientific, Inc.). Templates were amplified using 45 cycles of a 3-step PCR protocol, which included 30 sec of denaturation at 95°C, 30 sec of primer-dependent annealing at 60°C and 30 sec of template-dependent elongation at 72°C. Each gene-specific RNA was quantified in triplicate by qPCR, and mRNA ratios relative to the house-keeping gene GAPDH were calculated for standardization of gene expression levels across samples using the 2^−ΔΔCq method (22). All primer sequences and product sizes are shown in Table I.

Human foreskin primary keratinocytes (HFKs) were obtained from the American Type Culture Collection (Manassas, VA, USA) and cultured at 37°C in a humidified atmosphere with 5% CO₂ in Dulbecco's modified Eagle's medium:F12 supplemented with 0.18 mM adenine (Sigma-Aldrich^®; Merck KGaA, Darmstadt, Germany), 0.1 μg/ml hydrocortisone (Sigma-Aldrich; Merck KGaA), 100X Human Keratinocyte Growth Supplement (Thermo Fisher Scientific, Inc.), 2% foetal bovine serum (Gibco; Thermo Fisher Scientific, Inc.) and an antibiotic-antimycotic mixture (Invitrogen; Thermo Fisher Scientific, Inc.). The cells were transfected with the HPV16-E7 oncogene using Attractene^® reagent (Qiagen GmbH, Hilden, Germany) according to the manufacturer's recommended protocols. The cells were seeded in 6-well plates (3×10⁵ cells/well), and the next day, the cells were transfected with a mixture of 4.5 μl Attractene/well and 1.2 μg plasmid (PCDNA3 HPV-16E7) per well. To obtain transient and stable transfected cell lines, the expression of β2-microglobulin was used as an internal control of gene expression in the qPCR assay. Cells were selected and maintained in growth media containing 200 mg/ml geneticin (G418; Invitrogen; Thermo Fisher Scientific, Inc.).

Cervical tissue sections were fixed in 4% paraformaldehyde overnight at 4°C for immunohistochemical staining. Tissues were cut into 5-μm sections with a microtome and placed on electro-charged slides (Van-Wessel, Querétaro, México) and preserved using GVA mount solution (Zymed^®; Thermo Fisher Scientific, Inc.). Protein detection for immunohistochemistry was conducted using the Mouse/Rabbit PolyDetector HRP/DAB Detection system^® (Bio SB, Santa Barbara, CA, USA), according to the manufacturer's protocols, and HRP-conjugated anti-rat antibodies. The samples were incubated at 4°C overnight with primary antibodies against Ki-67 (cat. no. BSB 5713; 1:100 dilution; Bio SB), p16^INK4a (cat. no. sc-468; 1:100 dilution; Santa Cruz Biotechnology, Inc., Dallas, TX, USA), Granzyme B (GZMB; cat. no. AB4059; 1:100 dilution; Abcam, Cambridge, MA, USA) or PI-9 (cat. no. sc-57531; 1:100 dilution; Santa Cruz Biotechnology, Inc.). Anti-rabbit HRP-conjugated antibodies (1:200 dilution; cat. no. A0545; Sigma-Aldrich; Merck KGaA) were used as secondary antibodies for the detection of proteins. Following immunohistochemical procedures, the tissues were counterstained with haematoxylin (10 min haematoxylin and 2 min eosin) at room temperature and mounted using GVA mount solution (Zymed; Thermo Fisher Scientific, Inc.). For immunofluorescence, cervical tissue sections were rinsed in 1X PBS and blocked for 2 h at 4°C with 1X PBS supplemented with 0.3% Triton X-100 and 1% bovine serum albumin (Sigma-Aldrich; Merck KGaA). Next, the sections were washed three times with 1X PBS and incubated for 1 h at 37°C with the aforementioned primary antibodies against GZMB or PI-9. The sections were then incubated with a fluorescein isothiocyanate (FITC)-labelled secondary antibody (rabbit; 1:50 dilution; cat. no. 65-6111; Zymed; Thermo Fisher Scientific, Inc.) for 30 min at room temperature. The sections were then rinsed as aforementioned, counterstained with propidium iodide (30 min at room temperature) and mounted using Vectashield^® (Vector Laboratories, Inc., Burlingame, CA, USA). The stained sections were examined by confocal microscopy using a Leica TCS SP2 (Leica Microsystems, Inc., Buffalo Grove, IL, USA). Captured images were imported into Adobe Photoshop CS6 (Adobe Systems, Inc., San Jose, CA, USA) to generate maximum projections.

UVC irradiation was performed with a 1300 Series Class II, Type A2 Biological Safety Cabinet (Thermo Fisher Scientific, Inc.). HFKs were grown at 37°C in a humidified atmosphere with 5% CO₂ in DMEM:F12 media supplemented with 0.18 mM adenine (Sigma-Aldrich; Merck KGaA), 0.1 μg/ml hydrocortisone (Sigma-Aldrich; Merck KGaA), 4 μg/ml insulin (Gibco; Thermo Fisher Scientific, Inc.), 20 ng/ml recombinant epidermal growth factor (PeproTech, Inc., Rocky Hill, NJ, USA), 5% foetal bovine serum (Gibco; Thermo Fisher Scientific, Inc.) and an antibiotic-antimycotic mixture (Invitrogen^®). HFKs were irradiated with UVC radiation for 30 sec in the absence of medium. Following the addition of DMEM:F12-HFK medium, HFKs were cultured for 24 h prior to harvesting.

For all microarray data comparisons, ANOVA with Tukey's honest significant difference post hoc test was performed using Partek Genomics Suite^® 6.6 Software (Affymetrix; Thermo Fisher Scientific, Inc). For all RT-qPCR data comparisons a one-way ANOVA with Bonferroni correction was used, performed in GraphPad Prism 5 (Graphpad Software, Inc., La Jolla, CA, USA). P<0.05 was considered to indicate a statistically significant difference.

H&E staining of cervical tissue sections from 2-month-old mice from the FVB, FVB+E₂, K14E7 and K14E7+E₂ groups showed that only the K14E7+E₂ group developed low-grade lesions corresponding to cervical intraepithelial neoplasia grade I (CIN-I) (H&E; Fig. 1). Histopathology of the exocervical tissues of FVB (non-transgenic) mice demonstrated normal stratified squamous epithelium (H&E-FVB; Fig. 1), with sporadic expression of p16^INK4a and basal PCNA expression in the epithelial basal layer (p16^INK4a-FVB and PCNA-FVB; Fig. 1). Histopathology of the exocervical tissues of the E₂-treated FVB mice demonstrated moderately hyperproliferative squamous epithelium (arrows, H&E-FVB+E₂; Fig. 1), with similar sporadic expression of p16^INK4a as the untreated FVB mice, but with a higher PCNA expression in the epithelial basal layer, as expected with E₂ treatment (arrows, p16^INK4a-FVB+E₂ and PCNA-FVB+E₂; Fig. 1). In K14E7 mice, wider stratified basal and middle epithelial layers were observed (H&E-K14E7; Fig. 1). As expected with high-risk HPV-associated lesions (23–25), p16^INK4a expression was increased in these samples compared with that in samples from FVB mice (arrows, p16^INK4a-K14E7; Fig. 1). Additionally, an even signal for PCNA expression along the basal layer and a hyperproliferative stratified squamous epithelium in the exocervical tissues was observed (arrows, PCNA-K14E7; Fig. 1). Notably, in the E₂-treated K14E7 mice, the development of mild dysplasia, indicated by the presence of epithelial cells with large and extended nuclei in the basal and middle layers of the exocervical tissues (arrows, Fig. 1), accompanied with infiltration of mixed inflammatory cells (neutrophils and lymphocytes) was observed (dotted arrows, H&E-K14E7+E₂; Fig. 1). Similar to that in K14E7 mice, the expression of p16^INK4a in the E₂-treated K14E7 mice was increased relative to that in the FVB mice (arrows, p16^INK4a-K14E7+E₂; Fig. 1). Finally, the expression of PCNA in the E₂-treated K14E7 mice was also detected in the basal and suprabasal layers in the hyperproliferative stratified squamous epithelium of the exocervix, consistent with the presence of proliferating cells (arrows, PCNA-K14E7+E₂; Fig. 1). All these traits are characteristic of low-grade lesions corresponding to CIN-I. Classification of cervical dysplastic lesions in 2-month-old E₂-treated K14E7 mice was determined according to the 'histopathological grading system for cervical squamous carcinogenesis in transgenic mouse' (4).

To evaluate the global gene expression profiles in the initial stages of cancer development, 2-month-old FVB and K14E7 mice treated for 1 month with E₂ pellets were used and compared against untreated FVB (control) or K14E7 mice. The 4 groups were examined with Whole Mice Genome Oligo Microarrays using ANOVA, with a fold-change criterion of ≥1.5 and ≤-1.5 and a P-value of <0.05. In summary, it was observed that the most marked changes in gene expression occurred in the E₂-treated groups (FVB or K14E7). Compared with those in the FVB group, 685 differentially expressed genes (409 upregulated and 276 downregulated) were detected in the FVB+E₂ group, 321 differentially expressed genes (181 upregulated and 140 downregulated) were detected in the K14E7 group and 933 differentially expressed genes (554 upregulated and 379 downregulated) were detected in the K14E7+E₂ group (Fig. 2A). A Venn diagram (Fig. 2B) was constructed to identify genes that were commonly and exclusively modulated by the E7 oncoprotein and/or E₂. The comparison between the K14E7+E₂ and FVB+E₂ groups yielded 248 overlapping genes, and the comparison between the K14E7+E₂ and K14E7 groups yielded 50 overlapping genes. In total, 108 genes were overlapping among all groups. Notably, 527 (K14E7+E₂), 290 (FVB+E₂) and 124 (K14E7) genes were exclusively regulated in each group (data not shown; dropbox data).

We hypothesize that the 527 exclusively downregulated or upregulated genes in the K14E7+E₂ group that lead to carcinogenesis are regulated by the cooperation between E₂ and E7 (Fig. 2B). To identify differentially expressed genes involved in diverse cellular processes, a Gene Ontology analysis was performed using the IPA software and it was observed that several processes associated with cancer, including 'cellular movement', 'cancer' and 'cellular growth and proliferation', were mainly modified in the K14E7+E₂ group, and to a slightly lesser degree, in the FVB+E₂ group (Fig. 3), whereas genes associated with 'inflammatory response', and 'dermatological diseases and conditions', were altered in the K14E7 group. Notably, even though the K14E7+E₂ and FVB+E₂ groups shared differentially expressed genes involved in several cellular processes, the gene profile for each group was quite different (Fig. 3).

To validate the microarray data, genes involved in different processes were selected (Fig. 3) and RT-qPCR analysis was performed. A total of 15 specific genes that showed differential expression in the microarray data obtained from the K14E7+E₂ group were validated. Among the upregulated genes, cytochrome P450 family 2 subfamily e polypeptide 1 (Cyp2e1; not synergistic), glycerophosphodiester phosphodiesterase domain containing 3 (Gdpd3), interleukin 1 receptor type II (Il1r2), lecithin-retinol acyltransferase (phosphatidylcholine-retinol-O-acyltransferase) (Lrat), MGAT4 family member C (Mgat4c), natriuretic peptide type C (Nppc) and serine (or cysteine) peptidase inhibitor clade B member 9 [Serpinb9b; (PI-9); synergistic] were validated, and among the downregulated genes, eomesodermin (Eomes; not synergistic), Gzmb, Gzmc, Gzmd, Gzme, Gzmf, Gzmg and glucoside xylosyltransferase 2 (synergistic) were validated. The results for all the genes analysed by RT-qPCR were consistent with the microarray data; the comparative results are shown in Table II. A synergistic effect at the transcriptional level was found in 13 of the validated genes, probably as a result of interactions between E₂ and E7. Another set of genes that showed differential expression in the K14E7 group (Cyp2e1, Il1r2, Gzmd and Gzmg) and in the FVB+E₂ group (Gdpd3, Gzmd, Gzmg and Eomes) was also validated. Validation of all groups was made using FVB untreated mice as the control. The RT-qPCR results are summarized in Figs. 4–6.

To validate whether GrB and PI-9 protein expression was associated with the mRNA expression, immunofluorescence assays were performed for GrB and its inhibitor, PI-9, in tissue sections of murine exocervix obtained from the FVB, FVB+E₂, K14E7 and K14E7+E₂ groups. Consistent with the mRNA expression results, immunofluorescence assays demonstrated cytoplasmic GrB protein expression in the suprabasal and granular layers of exocervical tissues from the FVB control group (GrB-FVB; Fig. 7). For the FVB+E₂ group, high GrB expression in the hyperplasic epithelial cells of the suprabasal and granular layers of the exocervix was also found (GrB-FVB+E₂; Fig. 7). A low signal for GrB in the K14E7 group was found (GrB-K14E7; Fig. 7), and as expected, a very weak signal was observed in the K14E7+E₂ mice, indicating the downregulation of GrB at the protein level (GrB-K14E7+E₂; Fig. 7). On the other hand, also in accordance with the microarray and RT-qPCR results, the expression of PI-9 was very low in the FVB group, and PI-9 expression was not observed in the FVB+E₂ group (PI-9-FVB and PI-9-FVB+E₂; Fig. 7). Moreover, a clear cytoplasmic signal for PI-9 was observed in the basal and suprabasal layers of the cervical epithelium and the stromal tissues in the K14E7 group (PI-9-K14E7; Fig. 7), and a very high signal was observed in the basal and suprabasal layers of the epithelial tissues in the K14E7+E₂ group (PI-9-K14E7+E₂; Fig. 7).

To verify the effect of the HPV16-E7 oncoprotein on human GrB and PI-9 expression, HFK cell lines transiently or stably transfected with plasmids carrying the ORF for the HPV16-E7 oncoprotein were generated. By RT-qPCR analysis, it was determined that compared with that in the HFKs transfected with empty vector, GrB expression in HFKs that were stably transfected with HPV16-E7 and exposed to UV was inhibited (Fig. 8A). PI-9 expression was upregulated by 2.0-fold when HFKs were transiently transfected with the HPV16-E7 oncoprotein ORF (Fig. 8B) and by 5.2-fold when HFKs were stably transfected (Fig. 8C). These results suggest that the HPV16-E7 oncoprotein induces PI-9 expression in vivo and in vitro. The expression of E7 was determined by RT-qPCR in stably transfected HFKs (see figure S1 available in the following Dropbox location: https://www.dropbox.com/s/h8kbqty7umkbg60/Figure%20S1.pdf?dl=0). The expression of E7 was determined by RT-qPCR in stably transfected HFKs (data not shown; dropbox data).