Next Article in Journal
In Vitro Evaluation of the Therapeutic Potential of Phage VA7 against Enterotoxigenic Bacteroides fragilis Infection
Next Article in Special Issue
The Reservoir of Persistent Human Papillomavirus Infection; Strategies for Elimination Using Anti-Viral Therapies
Previous Article in Journal
Inactivation of HIV-1 in Polarized Infant Tonsil Epithelial Cells by Human Beta-Defensins 2 and 3 Tagged with the Protein Transduction Domain of HIV-1 Tat
Previous Article in Special Issue
The Not-So-Good, the Bad and the Ugly: HPV E5, E6 and E7 Oncoproteins in the Orchestration of Carcinogenesis
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
Viruses
Volume 13
Issue 10
10.3390/v13102045
viruses-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

The Mouse Papillomavirus Epigenetic Signature Is Characterised by DNA Hypermethylation after Lesion Regression

by
Allison M. Tschirley
1,
Peter A. Stockwell
1,
Euan J. Rodger
1,
Oliver Eltherington
1,
Ian M. Morison
1,
Neil Christensen
2,
Aniruddha Chatterjee
1,† and
Merilyn Hibma
1,*,†
1
Department of Pathology, Dunedin School of Medicine, University of Otago, Dunedin 9054, New Zealand
2
Department of Pathology, Pennsylvania State University, State College, PA 16802, USA
*
Author to whom correspondence should be addressed.
†
Joint senior authors.
Viruses 2021, 13(10), 2045; https://doi.org/10.3390/v13102045
Submission received: 22 September 2021 / Accepted: 7 October 2021 / Published: 11 October 2021
(This article belongs to the Special Issue HPV and HPV Vaccine)
Browse Figures
Versions Notes

Abstract

:
Papillomaviruses (PVs) are double-stranded DNA tumour viruses that can infect cutaneous and mucosal epidermis. Human papillomavirus (HPV) types have been linked to the causality of cutaneous squamous cell carcinoma (cSCC); however, HPV DNA is not always detected in the resultant tumour. DNA methylation is an epigenetic change that can contribute to carcinogenesis. We hypothesise that the DNA methylation pattern in cells is altered following PV infection. We tested if DNA methylation was altered by PV infection in the mouse papillomavirus (MmuPV1) model. Immunosuppressed mice were infected with MmuPV1 on cutaneous tail skin. Immunosuppression was withdrawn for some mice, causing lesions to spontaneously regress. Reduced representation bisulphite sequencing was carried out on DNA from the actively infected lesions, visibly regressed lesions, and mock-infected control mice. DNA methylation libraries were generated and analysed for differentially methylated regions throughout the genome. The presence of MmuPV1 sequences was also assessed. We identified 834 predominantly differentially hypermethylated fragments in regressed lesions, and no methylation differences in actively infected lesions. The promoter regions of genes associated with tumorigenicity, including the tumour suppressor protein DAPK1 and mismatch repair proteins MSH6 and PAPD7, were hypermethylated. Viral DNA was detected in active lesions and in some lesions that had regressed. This is the first description of the genome-wide DNA methylation landscape for active and regressed MmuPV1 lesions. We propose that the DNA hypermethylation in the regressed lesions that we report here may increase the susceptibility of cells to ultraviolet-induced cSCC.

1. Introduction

Human papillomaviruses (HPV) are a heterogenous, epitheliotropic, and ubiquitous group of non-enveloped viruses that range from harmless to cancer-causing. High-risk mucosal HPVs such as type 16 and 18 cause almost 100% of cervical cancers, as well as at least 30% of head and neck cancers [1]. Some HPV types infect cutaneous sites, causing hyperproliferation of skin cells that results in visible warts (papillomata), whereas infection with β-HPV types can occur in the absence of any visible lesion [2]. β-HPV types have been linked to the causality of cutaneous squamous cell carcinoma (cSCC), especially in the context of immune suppression [3,4]. However, this proposition is confounded by the observation that, unlike cervical cancer, β-HPV DNA is not always detected in the resultant tumour [5,6,7,8,9]. The lack of a need for the retention of viral DNA suggests that changes that support tumorigenicity occur in these cells following resolution of infection.
DNA methylation is an epigenetic change that can permanently alter expression of genes. It has been linked to many cancers, including cSCC [10,11,12], and some of the most frequently hypermethylated genes in cSCC include CADM1, CDH1, and DAPK1 [10]. Several DNA tumour viruses, such as Epstein–Barr virus, Hepatitis B virus, and HPV dysregulate host gene DNA methylation [13]. Additionally, studies have described regulation of immune-related and tumour suppressor host genes by viral methylation following expression of selected HPV proteins in cells or have reported methylation patterns in cancer cells containing high-risk HPV E6 and E7 [14]. However, the extent of methylation of host genes that occurs during infection and after resolution of lesions is currently unknown.
Mouse papillomavirus (MmuPV1) is a pi papillomavirus that was first identified as a cutaneous infection in immune suppressed mice [15,16]. The MmuPV1 genome encodes seven open reading frames (ORFs) found in many HPV types (E1, E2, E4, E6, E7, L1, and L2). As with other papillomaviruses, each of the MmuPV1 ORFs is differentially expressed in the viral life cycle [16]. MmuPV1 does not contain an E5 ORF, similar to the cutaneous β HPV types [17]. Importantly, MmuPV1 encodes genes that have oncogenic potential. MmuPV1 E6 and E7 viral proteins increase proliferation of cells and tumorigenicity [18], and mice infected with MmuPV1 following ultraviolet (UV) irradiation can develop lesions that progress to squamous cell carcinoma [19,20]. Interestingly, although viral DNA was abundant in the initial warts in those lesions, very few cells contained amplified viral DNA once lesions had progressed to cSCC. MmuPV1 has recently been reported to frequently integrate into the mouse genome [21].
We hypothesise that viral infection alters methylation of the host DNA, with the premise that this may increase the tumorigenic potential of those cells. We tested this hypothesis in MmuPV1 infected cutaneous skin in immunosuppressed BALB/c mice, and methylation changes that occur during active lesion growth and following regression of lesions were measured. Interestingly, we found a high level of host DNA hypermethylation only after resolution of the visible lesion, and not in active lesions following MmuPV1 infection, compared to mock-infected control skin. We propose that the hypermethylation of host DNA we observe here contributes to an increased susceptibility to tumorigenesis following resolution of papillomavirus lesions.

2. Materials and Methods

2.1. Animals

Female specific pathogen-free (SPF) BALB/c mice, aged 8–12 weeks, were obtained from the Hercus Taieri Research Unit at the University of Otago and housed in a SPF environment. A total of 63 mice were initially in the study. This strain of mouse is considered moderately susceptible to MmuPV1 and only develop lesions when immune-suppressed.

2.2. Immune Suppression and Group Descriptions

For systemic immune suppression, 100 mg/mL Cyclosporin (CsA; Novartis Neoral®, Basel, Switzerland) diluted in phosphate buffered saline (PBS) was administered to 30 control and 30 infected mice by gavage five times per week (75 mg/kg, 11.25 mg/mL), starting one week prior to MmuPV1 infection. Three mice were non-CsA treated control mice. Groups consisted of infected and actively increasing lesion mice (A) and MmuPV1 infected then regressed (R). Their mock-infected or mock-infected then regressed paired samples were denoted (C)A and (C)R, respectively. A comparison between CsA treated and untreated mice (C) was performed to assess any methylation differences due to immune suppression. Methylation effects of scarification were also controlled for by comparing the non-scarified region of mock-infected (C)A and the CsA treated groups. Commencing at day −3, Baytril (Bayer Comporation, Leverkusen, Germany) was added to drinking water for all immune suppressed mice to reduce the risk of bacterial infection.

2.3. Ethics

All animal experiments were approved by the Animal Ethics Committee at the University of Otago (ethical approval 12 Nov 2014; AEC56/14), were performed in accordance with the New Zealand Animal Welfare Act (1999), and adhered to the National Animal Ethics Committee (Ministry of Primary Industries, New Zealand) guidelines.

2.4. MmuPV1 Virus and Infection Protocol

For the initial infection of mice, mouse papillomavirus was previously isolated from lesions on the tails of Foxn1nu/Foxn1nu mice and diluted 1:10 in PBS, as previously described, was used [22]. Briefly, lesions were scraped from the tail and homogenised, and the supernatant was collected. Viral supernatant was then diluted with PBS and passed through a sterile syringe filter before storage at –20 °C. For infection of the mice, tail skin was scarified using a hand-held rotary tool with a felt wheel attachment on day −3 [23]. Animals were anaesthetised with Isoflurane (Bayer Corporation, Leverkusen, Germany), a 1.5 cm long mark was placed at the base of the tail using a felt tip pen, and the rotary tool was moved 10 times along the mark at a speed of approximately 15,000 rpm to remove the upper epidermal layer. Marcain (Pfizer, New York City, NY) was applied dropwise over the wound area for pain relief. At day 0, mice were placed in restrainers, and the scarified region was scored lightly with a needle tip. Inoculation was performed by placing 20 µL of prepared virus onto the scored region of the scarified area (Figure 1A). Mock-infected mice received PBS. At day 49, animals who had yet to show signs of a lesion were re-infected with 20 µL of native tail virus isolated from lesions on the tails of BALB/c mice, diluted 1:5 in PBS and sterile filtered. At day 102, CsA was discontinued for the regressed group, to allow the lesions to heal naturally.

2.5. Lesion Monitoring and Harvesting

Mice were monitored weekly for signs of lesions. Measurements of lesions were taken either immediately prior to harvesting or, for lesions in the regressed group, immediately prior to withdrawing CsA. A lesion was “actively increasing” if its largest volume measurement was its last one. A lesion was considered “partially regressed” if it had reached its greatest volume and was decreasing in volume. A lesion that was “fully regressed” had appeared during the study but was no longer overtly present. Beginning at day 102 for the regressed group, the lesions were outlined with marker pen for identification of tissue to be harvested once lesions regressed. The actively increasing lesion group mice and their mock-infected controls were culled at day 122 post-infection. Tail skin was removed, and epidermal dissociation of the epithelium was performed following the protocol of Lichti et al. [24]. DNA was extracted from epidermal tissue (approximately 3–5 mm2 per lesion) using a DNAeasy blood and tissue kit (Qiagen; Hilden, Germany). Regressed and CsA group mice and their controls were culled at day 133 post-infection once lesions were no longer palpable (31 days post withdrawal of CsA for regressed group). Lesion tissue was harvested as described for the actively infected group above.

2.6. Reduced Representation Bisulphite Sequencing Library Preparation

Reduced representation bisulphite sequencing (RRBS) was performed on the 16 samples of genomic DNA extracted from tail epidermal tissue as previously described [25,26,27]. Briefly, DNA was digested into fragments overnight using the MspI restriction endonuclease. The DNA was then purified following QIAquick kit (Qiagen, Hilden, Germany) instructions. Subsequently, DNA was end-repaired by incubation for 30 min with Nano End Repair Mix 2 (Illumina; San Diego, CA, USA) and again purified using the MinElute PCR purification kit (Qiagen; Hilden, Germany). Adenylation of 3′ ends was followed by ligation of sequencing adaptors, with a specific index to identify each sample. The fragments were size selected (150–325 bp) by excising the appropriate bands from 3% Nusieve agarose gels. These fragments were then bisulphite converted using EZ DNA methylation kit (Zymo Research; Irvine, CA, USA) prior to a PCR amplification step. The quality of the 16 libraries was confirmed using bioanalyser traces before being sequenced on an Illumina HiSeq2500 machine (Table 1).

2.7. DNA Methylation Analysis

Differential methylation analysis was performed using the previously reported differential methylation analysis package (DMAP) [28,29,30]. Briefly, we used a fragment-based approach, using the MspI fragments that had been digested in the RRBS library preparation as the unit of analysis for identifying variable methylation. Only fragments with at least two CpG sites and covered by at least ten sequenced reads per CpG were included in the data (hereon called high coverage fragments). High coverage fragments were mapped to the mouse genome and the closest protein coding genes were annotated.

2.8. Group Comparisons

Five different group comparisons were performed on the data produced by DMAP; regressed vs. control, active vs. control, active vs. regressed, CsA vs. no CsA, scarified vs. non-scarified (Figure 1B). Active and regressed tissues were the experimental groups. Non-scarified vs. scarified was a methylation comparison only and compared animal samples from the CsA and actively increasing lesion control groups. These control tissues were tested to determine whether scarification impacted on methylation. CsA vs. no CsA tissues was assessed to determine whether CsA administration impacted methylation.

2.9. Analysis of Viral Reads and Pathway Analysis

Methylation data were aligned to the MmuPV1 genome using bowtie in the Bismark Bisulphite Read Mapper [31]. Sequences that mapped to the MmuPV1 genome were also aligned to the Mus musculus genome to ascertain whether mapping was unique to the viral genome. The mapper BSMAP(z) was then used to verify the alignments shown by Bismark. This mapping software confirmed that the reads were not indicating artefactual behaviour of Bismark alone. The Integrative Genomics Viewer (IGV) was used to visualise reads within the MmuPV1 genome [32]. Lanes were auto-scaled to show proportionate representation between the samples.

2.10. Statistical Analyses

All data are presented as the mean ± standard error of the mean and all statistical analyses of methylation data were done in GraphPad Prism v7® (Graphpad Software, San Diego, CA, USA). For methylation analysis, an ANOVA test was applied to the fragments of the experimental and control groups, and regions showing methylation differences with a fold change of ≥ 1.5 and significant p-values were identified. Mann–Whitney U was used to test for differences between groups. The Benjamini-Hochberg false-discovery rate (FDR) method under the multiple t-tests function was applied to increase the power of the data and identify significant p values with an FDR of 5%, as indicated by the volcano plots. Unsupervised hierarchical clustering was performed in the R environment using the Euclidean distance metric of all analysed high coverage autosomal fragments (n = 126,681) from the 16 RRBS methylomes. This shows samples with bigger differences in total amount of methylation farther away from each other in the hierarchy. Pathway analysis was performed on a gene list containing the 221 genes differentially expressed in the core promoter region using the online platform Enrichr [33]. The enrichment p values calculated by Enrichr are from a modified Fisher’s exact test, which is a proportion test that assumes a binomial distribution and independence for probability of any gene belonging to any set.

3. Results

3.1. Immune-Suppressed BALB/c Mice Develop and Maintain Active MmuPV1 Infection

Thirty mice that were CsA treated were infected with MmuPV1 and were monitored for the development of visible lesions and another thirty CsA treated control mice were mock-infected with PBS (Figure 1A). Infected mice that had not developed visible lesions were re-infected at day 49. Out of 30 total mice, six developed lesions visible to the naked eye between days 49 and 63 post-infection, with a mean lesion length of 2.4 ± 0.3 mm (Figure 1B). The remainder of the infected mice did not develop visible lesions or developed lesions that spontaneously regressed while the animal was still immune-suppressed and were therefore excluded from the study. The final groups consisted of, firstly, mice (n = 3) that were harvested at day 122 with “actively increasing” lesions (A1, A2, A3) and three matched mock-infected controls ((C)A1, (C)A2, (C)A3), and secondly, mice that produced active lesions but were harvested at day 133 with “regressed” lesions (R1, R2, R3) following the withdrawal of CsA at day 102, and three matched CsA treated, mock-infected controls ((C)R1, (C)R2, (C)R3). All mice in the regressed group had healed their lesions (assessed by visual and tactile tests) by 31 days post CsA withdrawal. Lesion length measurements taken immediately prior to CsA withdrawal for the regressed group were 2 mm, 1.5 mm, and 3 mm for R1, R2, and R3, respectively. Lesion length measurements for the active group were 3 mm, 2mm, and 3 mm for A1, A2, and A3, respectively. Groups of two CsA treated mice (CsA1, CsA2) and two controls ((C)1, (C)2) were included in the study to test whether any changes in DNA methylation occurred as a result of CsA treatment alone.
RRBS was carried out on all 16 tissue samples. On average, we generated 9.3 × 106 ± 0.5 × 106 uniquely aligned reads for each of the samples tested (Table 1). The mean mapping efficiency was 47 ± 1.3%, which is consistent with the expected level of efficiency for a bisulphite converted genome.
To assess the relationship of global RRBS methylation between the samples that were tested, we performed hierarchical clustering of all analysed high coverage autosomal fragments from the 16 RRBS methylomes in this study (Figure 1C). We found that all the infected samples (A1-3, R1-3) grouped into two clusters, irrespective of whether they were actively increasing lesions or if the lesion had regressed (Figure 1C). Samples from uninfected mice ((C)A1-3, (C)R1-3, CsA1-2, (C)1-2) generally grouped into two other clusters, separated from the infected mice.
An analysis of differentially methylated fragments (DMFs) was carried out between test and control groups. The criteria for the inclusion of a DMF was a fold change difference of ≥ 1.5, a significant p-value (p < 0.05), and a false discovery rate (FDR) of 5% or less using the stringent Benjamini–Hochberg multiple test correction. We were surprised to see that no significantly DMFs were identified for actively increasing lesions (A1-3 c.f. (C)A1-3; Table 1). As might be predicted, treatment of mice with CsA (CsA1-2 c.f. (C)1-2) or the scarification of tissue alone ((C)A1-3 c.f. CsA1-2) had minimal effects on differential methylation. In contrast, 834 DMFs were identified when we compared tissues from lesions that had regressed compared to matched control tissues.

3.2. Genome-Scale DNA Methylation Identifies Extensive Hypermethylation in Regressed Lesions

We carried out further analysis to determine if the DMFs were hyper or hypomethylated. The striking finding from this analysis was that 98% of the DMFs identified in the regressed group were significantly hypermethylated (Mann–Whitney U, p < 0.0001) (Figure 2A), whereas none of the methylated fragments in the actively increasing lesion group were significantly hyper or hypomethylated.
To further probe the DMFs in the regressed MmuPV1 lesions relative to their location in the genome, we generated DNA methylation maps for the specific genomic regions. We found that methylation was increased overall in the DMFs from all genomic regions that we assessed in regressed skin when compared with controls (Figure 2B). Significantly increased methylation was detected in DMFs genome-wide (Mann–Whitney U, p < 0.0001), proximal to gene promoters (defined as ±1 kb of transcription start site) (Mann–Whitney U, p < 0.0001), in the core promoter region (defined as ±500 bp of transcription start site) (Mann–Whitney U, p < 0.0001), in gene introns (Mann–Whitney U, p < 0.0001) and gene exons (Mann–Whitney U, p < 0.0001), and in intergenic regions (Mann–Whitney U, p < 0.0001).
To further investigate the methylation landscape, we performed an analysis of the wider CpG island topography (10 kb upstream to 1 kb downstream of transcription start sites). We found that most of the hypermethylation was in the CpG island cores, being defined as stretches of DNA 500–1500 bp long with a CG:GC ratio of more than 0.6 [34]. Of 328 hypermethylated fragments detected within the region of interest (Figure 2C), 66% were found in the core, 4% were in the CpG island shore (up to 2 kb upstream and 2 kb downstream of the CpG island core), 1% were in the interface between core and shore (termed core-shore), and <1% were in the shelf (up to 2 kb upstream and 2 kb downstream of the shore). Additionally, 27% of the regions were outside of these CpG island features. Of the nine hypomethylated fragments that were detected (Figure 2D), five were in the CpG island core.

3.3. Genes in the Core Promoter Region Are Skewed Extensively towards Hypermethylation in Regressed Lesions

As promoter DNA methylation strongly influences expression from the corresponding gene, we specifically investigated differential methylation in the core promoter region. Of the 221 DMFs in the core promoter region, 214 were hypermethylated, and seven were hypomethylated. When comparing individual samples within the groups, all 100 genes with >20% methylation difference were hypermethylated in the regressed group (Figure 3A). This was also the case at > 40% methylation difference (Figure 3B), where all 27 genes were hypermethylated in regressed tissue compared to control tissue.

3.4. Hypermethylated Genes in Regressed Lesions Are Enriched for Cellular Senescence and Cancer Related Pathways and for CTCF and RNA Polymerase II Binding

To document the function of the genes that were hypermethylated in the core promoter regions of regressed skin, we performed functional enrichment analysis using two different sources (KEGG and Wiki Pathways for Mouse). We found genes involved in cellular senescence/cell cycle, mRNA processing, and p53 signalling were significantly enriched (p < 0.05) in both the analyses (Figure 4A). Disease associated pathways, particularly cancer pathways, were also enriched in our analysis. To place our findings in a broader epigenomic context, we utilised histone modification ChIP-Seq data (consisting of active and repressive histone marks) from the ENCODE project as well as from the Epigenomics roadmap (Figure 4B). The genes harbouring hypermethylation in the core promoters of regressed skin were predominantly enriched for the repressive chromatin mark H3K9me3. We also found significant enrichment for several active histone marks, especially H3K27ac and H3K9ac (Figure 4B). Overlap analysis of ENCODE transcription factor ChIP-seq data identified enrichment of several transcription factors for the hypermethylated genes (Figure 4C). The most frequently identified transcription factors were CTCF (which encodes a zinc-finger binding protein that is a transcriptional and chromatin regulator) and POLR2A (encoding the RNA Polymerase II Subunit).

3.5. MmuPV1 Viral Reads Are Highly Concentrated in the E4 Portion of the Genome

We probed the data for MmuPV1 viral reads to determine whether viral DNA remained after lesion regression (Figure 5A). Traces of the viral genome were present in all actively increasing tissue (5.7%, 0.2% and 5.4% of reads were viral in origin for samples A1 to A3, respectively). Interestingly, viral reads were also detected in two of three regressed tissues (7.2%, 0.0% and 5.3% for samples R1 to R3, respectively). The small number of background reads mapping to the MmuPV1 genome in uninfected groups (Table 2) most likely occurred through the misalignment of bisulphite reads, as previously reported [35]. The most intriguing finding of viral read analysis was the high proportion of viral reads within the E4 region of the MmuPV1 genome relative to other regions of the genome (Figure 5B,C). In two of the actively increasing lesion samples (A1 and A3), more than half of all reads aligned to this region. In one of the regressed lesions, sample R1, 74% of the viral reads mapped to E4. The L2 region also had high clustering, contributing up to 25% of all viral reads in the actively increasing lesion samples.

4. Discussion

Here, we show the first reported analysis of global and gene-specific DNA methylation during MmuPV1 infection and following regression of the lesion. Interestingly, differential methylation was only detected following lesion regression, and not during actively increasing lesions. This differential methylation was highly skewed towards hypermethylation.
Hypermethylation has been observed in the pre-cursor lesions of cervical cancer, with the increasing methylation of CCNA1 from low grade to high grade cervical lesions being suggested as a diagnostic marker for progression [36]. DNA methylation is variable during cervical carcinogenesis, with maximal variation in DNA methylation in “risk” CpG loci immediately prior to the onset of cervical cancer [37]. Additionally, HPV-associated head and neck cancer tumour growth was found to be suppressed by methylation inhibitors [38]. Thus, hypermethylation has been implicated in HPV-related cancers.
In our analysis, DAPK1 was hypermethylated in core promoters in regressed skin compared to control tissue. DAPK1 is a tumour suppressor protein down-regulated in many cancer types [39]. Li and colleagues found that hypermethylation of DAPK1 was associated with cSCC [10]. Although Li et al. found no correlation between hypermethylation of DAPK1 and the presence of cutaneous HPV, DAPK1 hypermethylation is implicated in other HPV associated human cancers, such as cervical cancer [40]. The hypermethylation of DAPK1 that we observed here suggests that this gene may be “switched off” in MmuPV1-regressed lesions; however, changes in protein levels as a result of hypermethylation needs to be confirmed for this and the other hypermethylated genes that we have identified.
The identification of hypermethylation of POLR2A shown here is consistent with the observation that one of the enriched pathways that we identified was mRNA processing, suggesting a role of hypermethylation in RNA processing in regressed skin. We also identified hypermethylation of the CTCF promoter. Variable CTCF binding has been strongly linked with differential DNA methylation in multiple human cell types, and it has been shown that CTCF binding patterns are markedly different in normal versus immortal human cells. In immortal cells, disruption of CTCF binding was strongly associated with hypermethylation [41]. Our analysis provides evidence for hypermethylation of CTCF which may lead to disruption of CTCF binding in regressed skin compared to normal skin controls.
Although our study did not investigate the effects of UV on MmuPV1 infection or lesion progression, UV is known to be a mutagen and an immunosuppressant and important in the causality of cSCC [42]. Uberoi et al. [20] found that immune-competent FVB mice infected with MmuPV1 and exposed to UV radiation were more likely to develop lesions, and that some of those lesions progressed to cSCC. Around 58% of FVB mice that were exposed to UVB either 24 h prior to, or 24h post-infection developed papillomas by 12 weeks; however, this was strain specific and could not be replicated in BALB/c mice. Some of the lesions (between 37% and 58%) harvested at 6 months showed some evidence of malignant progression. They concluded that UVB increases susceptibility to MmuPV1 disease by inducing immunosuppression, but their data also support a relationship between UVB, MmuPV1, and the development of cSCC.
Defective DNA mismatch repair following UV exposure can contribute to cSCC. We found that the DNA mismatch repair genes MSH6 and PAPD7 were hypermethylated within the core promoter region of regressed skin (5% and 6.5% methylation difference, respectively). MSH6 is important for UVB-induced apoptosis [43,44]. Its hypermethylation indicates that the damage response no longer occurs. Key genes such as MSH6 and PAPD7, whose promoters are hypermethylated in regressed MmuPV1 lesions, may render the skin more sensitive to the damaging effects of UV.
It is interesting to note that the methylation changes that we observe here only occurred after regression of the visible lesion, not during the active infection. Although it seems most likely that these DNA methylation changes occur in keratinocytes and as a result of the previous infection, there are other possible explanations. For example, epigenetic changes that include DNA methylation are associated with the healing response [45]. CsA withdrawal did not result in was not the source of the methylation changes we observed here, as this was controlled by the CsA also being withdrawn from the matched control mice. Although the predominant source of the DNA in this study is keratinocytes, some infiltration into the epidermis of immune cells that are specifically associated with immune-mediated regression of the lesion may contribute to changes in DNA hypermethylation. Mitra et al. (2020) reported distinct immune methylation signatures in metastatic melanoma that related to patient survival, and it is feasible that the DNA hypermethylation that we measured in this study may reflect changes in infiltrated immune cells [46].
A limitation of this study is that the experimental group sizes were smaller than planned (n = 3), because the number of mice that developed visible papillomas in our study was about a third of what had been reported by others. It has previously been reported that BALB/c mice demonstrate intermediate susceptibility to MmuPV1 infection when immunosuppressed with CsA [47]. In that study, five out of eight (63%) BALB/c mice produced MmuPV1 lesions following infection. In our study, the proportion of mice that developed lesions was much lower (20 %), and some mice required a second dose of virus for a visible lesion to appear. The differences between these studies may be attributable to differences in the amount of virus administered or in the administration of the CsA treatment (oral gavage c.f. subcutaneous injection) between the two studies. Regarding the group sizes, we acknowledge that other studies reporting methylation in mouse skin typically use slightly larger group sizes (n = 5) [48,49]. In MmuPV1 infectivity studies, Xue et al. (2017) used n = 3 for their analyses, but with three sites per mouse, and Wei et al. (2020) used group sizes of 3–8 mice for tongue infections using MmuPV1 [50]. We note that hypermethylation in the regressed mice was pronounced and consistent in all three animals in our study and was not detected in any of the actively infected mice, giving us confidence in the validity of these data. Furthermore, statistical rigor was applied in the methylation analysis, with an FDR of 5% on Benjamini–Hochberg-corrected data. However, we also note that in HPV16-E5 transgenic mouse exophytic papillomas developed following dimethylbenz[a]anthracene/12-O-tetradecanoyl/phorbol-13-acetate (DMBA/TPA) treatment but that the carcinomas that developed were flattened out endophytic lesions. Although these models differ in several regards, it is possible that the regressing lesions were flattened out endophytic cSCC lesions [51]. This should be confirmed in future studies.
Our analysis of viral reads found evidence of viral DNA in regressed tissues with no remaining visible lesion. Viral gene expression was detected by Xue et al. [52] well before lesions were visible. These investigators observed that the viral copy number did not directly correlate with lesion size. This was similarly found in our analysis of viral reads and the relationship with lesion size, where one visible lesion in the actively increasing lesion group (sample A2) yielded only 0.2% viral reads, while skin with no visible lesion in the regressed group (sample R1) showed the highest level of viral DNA at 7.2%. A limitation of using RRBS data to look at viral reads is that we cannot visualise the whole genome, only fragments, and consequently, we cannot comment on copy numbers or integration of the viral genome.
The low viral reads of sample A2 could be partly explained by the size of the lesion, which was smallest of the three actively increasing lesions (height: 0 mm, length: 2 mm, width: 1 mm) and could contribute to its methylation profile clustering with two of the samples from the regressed lesion group (Figure 1C). Interestingly, R1, which had the highest level of viral reads, clustered with A2 (0.2%), R2 (0%), and (C)1 (did not receive virus). Except for (C)1, this shows that the largest difference in methylation is between four mice that were mock infected and three mice that were infected with MmuPV1, suggesting that infection status, not level of viral reads, is more important in methylation differences. Rabbit oral papillomavirus DNA and RNA has been found to persist for at least a year post lesion regression [53]. The high viral reads in two out of three regressed lesions suggests viral latency in MmuPV1 is possible, or that viral integration has occurred. The diversity in the quantity of viral DNA found could be explained by the stage of epithelial differentiation that the cells were in when collected, as productive papillomavirus infection is dependent on the differentiation of the host epithelial cell. A progression of this research could involve re-establishment of immunosuppression in regressed animals to see if the latent infection re-emerges at the site of the original lesion.
One of the striking findings within the analysis of viral reads was the high clustering of reads in the E4 region. E4 protein is found to be abundant in upper layers of the epithelium of productive lesions [54]. E4 expression is followed by the late structural proteins, explaining the concurrent, though considerably lower, expression of L2. This pattern was also found by Xue et al. [52], who looked solely at active MmuPV1 infections in the tail, muzzle, and ear of mice.
There is considerable variability in the pattern of expression of different papillomavirus types, depending on the animal model and tissue tropism [55]. High levels of E4 transcripts in the upper layers of the epithelium, similarly to the DNA viral reads found here, follow the pattern found in low-grade squamous intraepithelial lesions of high-risk cervical HPV infections [56]. In an analysis of 92 different HPV genomes, the E4 region was found to have the highest proportion of CpG sites, sites where methylation predominantly occurs [57]. This is even though CpG sites tend to be under-represented in viruses, which may be a way to avoid methylation by host methyltransferases or CpG-mediated immune responses.
In summary, our study demonstrates the notable hypermethylation and the retention of some viral DNA in previously MmuPV1-infected tissue following resolution of the visible lesion. We highlight several genes hypermethylated in the core promoter regions of regressed skin that we speculate could influence susceptibility of those cells to UV-induced cSCC.

Author Contributions

Conceptualisation of research: M.H., A.M.T., A.C., N.C., I.M.M., and O.E.; data curation: A.M.T., P.A.S., and A.C.; formal analysis: P.A.S., A.M.T., M.H., A.C., and E.J.R.; methodology: A.M.T., N.C., M.H., and A.C.; software: A.C., P.A.S., and E.J.R.; funding acquisition: M.H., A.M.T., A.C., and I.M.M.; writing: A.M.T., M.H., O.E., A.C., I.M.M., and N.C. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by University of Otago Research Grant (2016) and a University of Otago Dean’s Bequest Grant (2016).

Institutional Review Board Statement

The study was approved by the Ethics Committee of University of Otago (protocol code AEC56/14).

Informed Consent Statement

Not applicable.

Data Availability Statement

The RRBS data generated as part of this study were submitted to Gene Expression Omnibus (GEO) with the accession number GSE160868.

Acknowledgments

We thank Jiafen Hu (Pennsylvania State University) for supplying additional RNAseq data from MmuPV1 tissues. A.C. would like to thank the Royal Society of NZ (Rutherford Discovery Fellowship) for their support.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

References

  1. de Martel, C.; Plummer, M.; Vignat, J.; Franceschi, S. Worldwide burden of cancer attributable to HPV by site, country and HPV type. Int. J. Cancer 2017, 141, 664–670. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  2. Antonsson, A.; Forslund, O.; Ekberg, H.; Sterner, G.; Hansson, B.G. The ubiquity and impressive genomic diversity of human skin papillomaviruses suggest a commensalic nature of these viruses. J. Virol. 2000, 74, 11636–11641. [Google Scholar] [CrossRef] [Green Version]
  3. De Jong, S.J.; Imahorn, E.; Itin, P.; Uitto, J.; Orth, G.; Jouanguy, E.; Casanova, J.-L.; Burger, B. Epidermodysplasia verruciformis: Inborn errors of immunity to human beta-papillomaviruses. Front. Microbiol. 2018, 9, 1222. [Google Scholar] [CrossRef] [PubMed]
  4. Mittal, A.; Colegio, O.R. Skin Cancers in Organ Transplant Recipients. Am. J. Transplant. 2017, 17, 2509–2530. [Google Scholar] [CrossRef] [PubMed]
  5. Bavinck, J.N.B.; Plasmeijer, E.I.; Feltkamp, M.C.W. β-Papillomavirus Infection and Skin Cancer. J. Investig. Dermatol. 2008, 128, 1355–1358. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  6. Hufbauer, M.; Akgül, B. Molecular Mechanisms of Human Papillomavirus Induced Skin Carcinogenesis. Viruses 2017, 9, 187. [Google Scholar] [CrossRef] [Green Version]
  7. Karagas, M.R.; Nelson, H.; Sehr, P.; Waterboer, T.; Stukel, T.; Andrew, A.; Green, A.C.; Bavinck, J.N.B.; Perry, A.; Spencer, S.; et al. Human papillomavirus infection and incidence of squamous cell and basal cell carcinomas of the skin. J. Natl. Cancer Inst. 2006, 98, 389–395. [Google Scholar] [CrossRef]
  8. Tommasino, M. HPV and skin carcinogenesis. Papillomavirus Res. 2019, 7, 129–131. [Google Scholar] [CrossRef]
  9. Underbrink, M.P.; Howie, H.L.; Bedard, K.M.; Koop, J.I.; Galloway, D.A. E6 proteins from multiple human betapapillomavirus types degrade Bak and protect keratinocytes from apoptosis after UVB irradiation. J. Virol. 2008, 82, 10408–10417. [Google Scholar] [CrossRef] [Green Version]
  10. Li, L.; Jiang, M.; Feng, Q.; Kiviat, N.B.; Stern, J.E.; Hawes, S.; Cherne, S.; Lu, H. Aberrant methylation changes detected in cutaneous squamous cell carcinoma of immunocompetent individuals. Cell Biochem. Biophys. 2015, 72, 599–604. [Google Scholar] [CrossRef]
  11. Hervás-Marín, D.; Higgins, F.; Sanmartín, O.; López-Guerrero, J.A.; Bañó, M.C.; Igual, J.C.; Quilis, I.; Sandoval, J. Genome wide DNA methylation profiling identifies specific epigenetic features in high-risk cutaneous squamous cell carcinoma. PLoS ONE 2019, 14, e0223341. [Google Scholar] [CrossRef]
  12. Rodríguez-Paredes, M.; Bormann, F.; Raddatz, G.; Gutekunst, J.; Lucena-Porcel, C.; Köhler, F.; Wurzer, E.; Schmidt, K.; Gallinat, S.; Wenck, H.; et al. Methylation profiling identifies two subclasses of squamous cell carcinoma related to distinct cells of origin. Nat. Commun. 2018, 9, 577. [Google Scholar] [CrossRef]
  13. Kuss-Duerkop, S.K.; Westrich, J.A.; Pyeon, D. DNA Tumor Virus Regulation of Host DNA Methylation and Its Implications for Immune Evasion and Oncogenesis. Viruses 2018, 10, 82. [Google Scholar] [CrossRef] [Green Version]
  14. Durzynska, J.; Lesniewicz, K.; Poreba, E. Human papillomaviruses in epigenetic regulations. Mutat. Res./Rev. Mutat. Res. 2017, 772, 36–50. [Google Scholar] [CrossRef] [PubMed]
  15. Ingle, A.; Ghim, S.; Joh, J.; Chepkoech, I.; Bennett Jenson, A.; Sundberg, J.P. Novel laboratory mouse papillomavirus (MusPV) infection. Vet. Pathol. 2011, 48, 500–505. [Google Scholar] [CrossRef] [PubMed]
  16. Joh, J.; Jenson, A.B.; King, W.; Proctor, M.; Ingle, A.; Sundberg, J.P.; Ghim, S.-J. Genomic analysis of the first laboratory-mouse papillomavirus. J. Gen. Virol. 2011, 92, 692–698. [Google Scholar] [CrossRef] [PubMed]
  17. Doorbar, J.; Egawa, N.; Griffin, H.; Kranjec, C.; Murakami, I. Human papillomavirus molecular biology and disease association. Rev. Med. Virol. 2015, 25, 2–23. [Google Scholar] [CrossRef] [Green Version]
  18. Hu, J.; Cladel, N.M.; Budgeon, L.R.; Balogh, K.K.; Christensen, N.D. The Mouse Papillomavirus Infection Model. Viruses 2017, 9, 246. [Google Scholar] [CrossRef]
  19. Spurgeon, M.E.; Uberoi, A.; McGregor, S.M.; Wei, T.; Ward-Shaw, E.; Lambert, P.F. A novel in vivo infection model to study papillomavirus-mediated disease of the female reproductive tract. MBio 2019, 10, e00180-19. [Google Scholar] [CrossRef] [Green Version]
  20. Uberoi, A.; Yoshida, S.; Frazer, I.H.; Pitot, H.C.; Lambert, P.F. Role of Ultraviolet Radiation in Papillomavirus-Induced Disease. PLoS Pathog. 2016, 12, e1005664. [Google Scholar] [CrossRef]
  21. Yu, L.; Majerciak, V.; Xue, X.-Y.; Uberoi, A.; Lobanov, A.; Chen, X.; Cam, M.; Hughes, S.H.; Lambert, P.F.; Zheng, Z.-M. Mouse papillomavirus type 1 (MmuPV1) DNA is frequently integrated in benign tumors by microhomology-mediated end-joining. PLoS Pathog. 2021, 17, e1009812. [Google Scholar] [CrossRef] [PubMed]
  22. Hu, J.; Budgeon, L.R.; Cladel, N.M.; Balogh, K.; Myers, R.; Cooper, T.K.; Christensen, N.D. Tracking vaginal, anal and oral infection in a mouse papillomavirus infection model. J. Gen. Virol. 2015, 96, 3554–3565. [Google Scholar] [CrossRef] [PubMed]
  23. Handisurya, A.; Day, P.M.; Thompson, C.D.; Buck, C.B.; Pang, Y.-Y.S.; Lowy, D.R.; Schiller, J.T. Characterization of Mus musculus papillomavirus 1 infection in situ reveals an unusual pattern of late gene expression and capsid protein localization. J. Virol. 2013, 87, 13214–13225. [Google Scholar] [CrossRef] [Green Version]
  24. Lichti, U.; Anders, J.; Yuspa, S.H. Isolation and short-term culture of primary keratinocytes, hair follicle populations and dermal cells from newborn mice and keratinocytes from adult mice for in vitro analysis and for grafting to immunodeficient mice. Nat. Protoc. 2008, 3, 799–810. [Google Scholar] [CrossRef] [PubMed]
  25. Chatterjee, A.; Rodger, E.J.; Stockwell, P.A.; Weeks, R.J.; Morison, I.M. Technical considerations for reduced representation bisulfite sequencing with multiplexed libraries. J. Biomed. Biotechnol. 2012, 2012, 741542. [Google Scholar] [CrossRef]
  26. Chatterjee, A.; Rodger, E.J.; Stockwell, P.A.; Le Mée, G.; Morison, I.M. Generating Multiple Base-Resolution DNA Methylomes Using Reduced Representation Bisulfite Sequencing. In Oral Biology: Molecular Techniques and Applications; Seymour, G.J., Cullinan, M.P., Heng, N.C.K., Eds.; Springer: New York, NY, USA, 2017; pp. 279–298. [Google Scholar]
  27. Ludgate, J.L.; Wright, J.; Stockwell, P.A.; Morison, I.M.; Eccles, M.R.; Chatterjee, A. A streamlined method for analysing genome-wide DNA methylation patterns from low amounts of FFPE DNA. BMC Med Genom. 2017, 10, 54. [Google Scholar] [CrossRef] [Green Version]
  28. Stockwell, P.A.; Chatterjee, A.; Rodger, E.J.; Morison, I.M. DMAP: Differential methylation analysis package for RRBS and WGBS data. Bioinformatics 2014, 30, 1814–1822. [Google Scholar] [CrossRef]
  29. Chatterjee, A.; Stockwell, P.A.; Rodger, E.J.; Duncan, E.J.; Parry, M.F.; Weeks, R.J.; Morison, I.M. Genome-wide DNA methylation map of human neutrophils reveals widespread inter-individual epigenetic variation. Sci. Rep. 2015, 5, 17328. [Google Scholar] [CrossRef] [Green Version]
  30. Chatterjee, A.; Stockwell, P.A.; Rodger, E.J.; Morison, I.M. Genome-scale DNA methylome and transcriptome profiling of human neutrophils. Sci. Data 2016, 3, 160019. [Google Scholar] [CrossRef] [Green Version]
  31. Krueger, F.; Andrews, S.R. Bismark: A flexible aligner and methylation caller for Bisulfite-Seq applications. Bioinformatics 2011, 27, 1571–1572. [Google Scholar] [CrossRef]
  32. Thorvaldsdóttir, H.; Robinson, J.T.; Mesirov, J.P. Integrative Genomics Viewer (IGV): High-performance genomics data visualization and exploration. Brief. Bioinform. 2013, 14, 178–192. [Google Scholar] [CrossRef] [Green Version]
  33. Chen, E.Y.; Tan, C.M.; Kou, Y.; Duan, Q.; Wang, Z.; Meirelles, G.V.; Clark, N.R.; Ma’Ayan, A. Enrichr: Interactive and collaborative HTML5 gene list enrichment analysis tool. BMC Bioinform. 2013, 14, 128. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  34. Cross, S.H.; Bird, A.P. CpG islands and genes. Curr. Opin. Genet. Dev. 1995, 5, 309–314. [Google Scholar] [CrossRef]
  35. Xi, Y.; Li, W. BSMAP: Whole genome bisulfite sequence MAPping program. BMC Bioinform. 2009, 10, 232. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  36. Chujan, S.; Kitkumthorn, N.; Siriangkul, S.; Mutirangura, A. CCNA1 promoter methylation: A potential marker for grading Papanicolaou smear cervical squamous intraepithelial lesions. Asian Pac. J. Cancer Prev. 2014, 15, 7971–7975. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  37. Teschendorff, A.E.; Liu, X.; Caren, H.; Pollard, S.M.; Beck, S.; Widschwendter, M.; Chen, L. The dynamics of DNA methylation covariation patterns in carcinogenesis. PLoS Comput. Biol. 2014, 10, e1003709. [Google Scholar] [CrossRef] [Green Version]
  38. Biktasova, A.; Hajek, M.; Sewell, A.; Gary, C.; Bellinger, G.; Deshpande, H.A.; Bhatia, A.; Burtness, B.; Judson, B.; Mehra, S.; et al. Demethylation Therapy as a Targeted Treatment for Human Papillomavirus-Associated Head and Neck Cancer. Clin. Cancer Res. 2017, 23, 7276–7287. [Google Scholar] [CrossRef] [Green Version]
  39. Gozuacik, D.; Kimchi, A. DAPk protein family and cancer. Autophagy 2006, 2, 74–79. [Google Scholar] [CrossRef] [Green Version]
  40. Yanatatsaneejit, P.; Chalertpet, K.; Sukbhattee, J.; Nuchcharoen, I.; Phumcharoen, P.; Mutirangura, A. Promoter methylation of tumor suppressor genes induced by human papillomavirus in cervical cancer. Oncol. Lett. 2020, 20, 955–961. [Google Scholar] [CrossRef]
  41. Lang, F.; Li, X.; Vladimirova, O.; Hu, B.; Chen, G.; Xiao, Y.; Singh, V.; Lu, D.; Li, L.; Han, H.; et al. CTCF interacts with the lytic HSV-1 genome to promote viral transcription. Sci. Rep. 2017, 7, 39861. [Google Scholar] [CrossRef] [Green Version]
  42. Pfeifer, G.P. Mechanisms of UV-induced mutations and skin cancer. Genome Instab. Dis. 2020, 1, 99–113. [Google Scholar] [CrossRef] [Green Version]
  43. Young, L.; Listgarten, J.; Trotter, M.; Andrew, S.; Tron, V. Evidence that dysregulated DNA mismatch repair characterizes human nonmelanoma skin cancer. Br. J. Dermatol. 2008, 158, 59–69. [Google Scholar] [CrossRef]
  44. Young, L.C.; Hays, J.B.; Tron, V.A.; Andrew, S.E. DNA mismatch repair proteins: Potential guardians against genomic instability and tumorigenesis induced by ultraviolet photoproducts. J. Investig. Dermatol. 2003, 121, 435–440. [Google Scholar] [CrossRef] [Green Version]
  45. Lewis, C.J.; Stevenson, A.; Fear, M.W.; Wood, F.M. A review of epigenetic regulation in wound healing: Implications for the future of wound care. Wound Repair Regen 2020, 28, 710–718. [Google Scholar] [CrossRef]
  46. Mitra, S.; Lauss, M.; Cabrita, R.; Choi, J.; Zhang, T.; Isaksson, K.; Olsson, H.; Ingvar, C.; Carneiro, A.; Staaf, J.; et al. Analysis of DNA methylation patterns in the tumor immune microenvironment of metastatic melanoma. Mol. Oncol. 2020, 14, 933–950. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  47. Handisurya, A.; Day, P.M.; Thompson, C.D.; Bonelli, M.; Lowy, D.R.; Schiller, J.T. Strain-Specific Properties and T Cells Regulate the Susceptibility to Papilloma Induction by Mus musculus Papillomavirus 1. PLoS Pathog. 2014, 10, e1004314. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  48. Fraga, M.; Herranz, M.; Espada, J.; Ballestar, E.; Paz, M.F.; Ropero, S.; Erkek, E.; Bozdogan, O.; Peinado, H.; Niveleau, A.; et al. A mouse skin multistage carcinogenesis model reflects the aberrant DNA methylation patterns of human tumors. Cancer Res. 2004, 64, 5527–5534. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  49. Qian, H.; Xu, X. Reduction in DNA methyltransferases and alteration of DNA methylation pattern associate with mouse skin ageing. Exp. Dermatol. 2014, 23, 357–359. [Google Scholar] [CrossRef]
  50. Wei, T.; Buehler, D.; Ward-Shaw, E.; Lambert, P.F. An Infection-Based Murine Model for Papillomavirus-Associated Head and Neck Cancer. mBio 2020, 11, e00908-20. [Google Scholar] [CrossRef]
  51. Maufort, J.P.; Williams, S.M.; Pitot, H.C.; Lambert, P.F. Human papillomavirus 16 E5 oncogene contributes to two stages of skin carcinogenesis. Cancer Res. 2007, 67, 6106–6112. [Google Scholar] [CrossRef] [Green Version]
  52. Xue, X.-Y.; Majerciak, V.; Uberoi, A.; Kim, B.-H.; Gotte, D.; Chen, X.; Cam, M.; Lambert, P.F.; Zheng, Z.-M. The full transcription map of mouse papillomavirus type 1 (MmuPV1) in mouse wart tissues. PLoS Pathog. 2017, 13, e1006715. [Google Scholar] [CrossRef] [PubMed]
  53. Maglennon, G.A.; McIntosh, P.; Doorbar, J. Persistence of viral DNA in the epithelial basal layer suggests a model for papillomavirus latency following immune regression. Virology 2011, 414, 153–163. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  54. Middleton, K.; Peh, W.; Southern, S.; Griffin, H.; Sotlar, K.; Nakahara, T.; El-Sherif, A.; Morris, L.; Seth, R.; Hibma, M.; et al. Organization of human papillomavirus productive cycle during neoplastic progression provides a basis for selection of diagnostic markers. J. Virol. 2003, 77, 10186–10201. [Google Scholar] [CrossRef] [Green Version]
  55. Peh, W.L.; Middleton, K.; Christensen, N.; Nicholls, P.; Egawa, K.; Sotlar, K.; Brandsma, J.; Percival, A.; Lewis, J.; Liu, W.J.; et al. Life Cycle Heterogeneity in Animal Models of Human Papillomavirus-Associated Disease. J. Virol. 2002, 76, 10401. [Google Scholar] [CrossRef] [Green Version]
  56. Kremer, W.W.; Vink, F.J.; van Zummeren, M.; Dreyer, G.; Rozendaal, L.; Doorbar, J.; Bleeker, M.C.G.; Meijer, C.J.L.M. Characterization of cervical biopsies of women with HIV and HPV co-infection using p16ink4a, ki-67 and HPV E4 immunohistochemistry and DNA methylation. Mod. Pathol. 2020, 33, 1968–1978. [Google Scholar] [CrossRef] [PubMed]
  57. Galván, S.C.; Martínez-Salazar, M.; Galván, V.M.; Méndez, R.; Díaz-Contreras, G.T.; Alvarado-Hermida, M.; Alcántara-Silva, R.; García-Carrancá, A. Analysis of CpG methylation sites and CGI among human papillomavirus DNA genomes. BMC Genom. 2011, 12, 580. [Google Scholar] [CrossRef] [PubMed] [Green Version]
Viruses 13 02045 g001 550
Figure 1. Experimental setup for the mouse model and preliminary methylation analysis. (A) Timeline of infection model. Colour-coded mouse images represent the three groups as shown in the boxes below. Boxes show sample names of each group and their associated control samples. (B) Representative images of the experimental groups show mouse tails at sample location; (i) Active viral infection, (ii) Regressed lesion (original lesion location delineated in black), (iii) CsA treated, uninfected. (C) Hierarchical clustering dendrogram showing the difference in total methylation of all high-coverage autosomal fragments for each of the 16 samples. A and R indicate actively MmuPV1 infected and MmuPV1 infected then regressed, respectively. (C)A and (C)R indicate mock infected and mock infected then regressed, respectively. CsA and (C) indicate CsA treated and untreated, respectively. Except for one sample, the largest difference in methylation is between four mice that were mock infected and three mice that were infected with MmuPV1.
Figure 1. Experimental setup for the mouse model and preliminary methylation analysis. (A) Timeline of infection model. Colour-coded mouse images represent the three groups as shown in the boxes below. Boxes show sample names of each group and their associated control samples. (B) Representative images of the experimental groups show mouse tails at sample location; (i) Active viral infection, (ii) Regressed lesion (original lesion location delineated in black), (iii) CsA treated, uninfected. (C) Hierarchical clustering dendrogram showing the difference in total methylation of all high-coverage autosomal fragments for each of the 16 samples. A and R indicate actively MmuPV1 infected and MmuPV1 infected then regressed, respectively. (C)A and (C)R indicate mock infected and mock infected then regressed, respectively. CsA and (C) indicate CsA treated and untreated, respectively. Except for one sample, the largest difference in methylation is between four mice that were mock infected and three mice that were infected with MmuPV1.
Viruses 13 02045 g001
Viruses 13 02045 g002 550
Figure 2. Global methylation description. (A) Volcano plots of differentially methylated fragments (DMFs) in comparison of MmuPV1 infected then regressed to mock infected then regressed control tissue (left) and actively infected to mock infected control tissue (right). Horizontal dotted line represents Benjamini-Hochberg false-discovery rate of 5% to control for multiple t-tests. Vertical dotted line represents a mean difference of 0 in the proportion of methylation between groups; (B) Proportion of methylated DMFs by genomic region for regressed control tissues (blue) and regressed (orange). Data are represented in quartiles with whiskers indicating outliers outside of upper and lower quartiles. (C) Within CpG islands, differentially hypermethylated regions (n = 328) and (D) hypomethylated regions (n = 9) were located predominantly in the core.
Figure 2. Global methylation description. (A) Volcano plots of differentially methylated fragments (DMFs) in comparison of MmuPV1 infected then regressed to mock infected then regressed control tissue (left) and actively infected to mock infected control tissue (right). Horizontal dotted line represents Benjamini-Hochberg false-discovery rate of 5% to control for multiple t-tests. Vertical dotted line represents a mean difference of 0 in the proportion of methylation between groups; (B) Proportion of methylated DMFs by genomic region for regressed control tissues (blue) and regressed (orange). Data are represented in quartiles with whiskers indicating outliers outside of upper and lower quartiles. (C) Within CpG islands, differentially hypermethylated regions (n = 328) and (D) hypomethylated regions (n = 9) were located predominantly in the core.
Viruses 13 02045 g002
Viruses 13 02045 g003 550
Figure 3. Extensive hypermethylation identified in regressed skin. (A) Individual DMFs in the core promoter region of genes with >20% methylation difference between regressed and control groups by individual sample. Five genes with insufficient data to calculate methylation in at least one sample have been removed from the heatmap (B) Individual DMFs in the core promoter region of genes with >40% methylation difference between regressed and control groups. More than one DMF was identified for some genes. Error bars display standard error of the mean.
Figure 3. Extensive hypermethylation identified in regressed skin. (A) Individual DMFs in the core promoter region of genes with >20% methylation difference between regressed and control groups by individual sample. Five genes with insufficient data to calculate methylation in at least one sample have been removed from the heatmap (B) Individual DMFs in the core promoter region of genes with >40% methylation difference between regressed and control groups. More than one DMF was identified for some genes. Error bars display standard error of the mean.
Viruses 13 02045 g003
Viruses 13 02045 g004 550
Figure 4. Pathway and transcriptional function analysis of differentially methylated fragments (DMFs). Enrichr analysis with default parameters was used to identify pathways, histone modifications and transcription factor targets associated with DMFs in core gene promoter regions. Individual heatmaps show: (A) enriched KEGG and Wiki Pathways (p < 0.05); (B) enriched histone modifications from Epigenomics Roadmap and ENCODE datasets (p < 0.01); (C) Transcription factor targets from ENCODE ChIP-Seq and TRANSFAC/JASPAR positional weight matrices (p < 0.01).
Figure 4. Pathway and transcriptional function analysis of differentially methylated fragments (DMFs). Enrichr analysis with default parameters was used to identify pathways, histone modifications and transcription factor targets associated with DMFs in core gene promoter regions. Individual heatmaps show: (A) enriched KEGG and Wiki Pathways (p < 0.05); (B) enriched histone modifications from Epigenomics Roadmap and ENCODE datasets (p < 0.01); (C) Transcription factor targets from ENCODE ChIP-Seq and TRANSFAC/JASPAR positional weight matrices (p < 0.01).
Viruses 13 02045 g004
Viruses 13 02045 g005 550
Figure 5. MmuPV1 viral reads within individual samples. (A) Images of mouse tails at time of harvest. Marker pen was used to delineate the area of the lesions at the height of infection. (B,C) Coverage and distribution of viral reads visualised using IGV software. Coverage depth of reads set to group autoscale for (B) regressed (0–266,605) and (C) active groups (0–541,120).
Figure 5. MmuPV1 viral reads within individual samples. (A) Images of mouse tails at time of harvest. Marker pen was used to delineate the area of the lesions at the height of infection. (B,C) Coverage and distribution of viral reads visualised using IGV software. Coverage depth of reads set to group autoscale for (B) regressed (0–266,605) and (C) active groups (0–541,120).
Viruses 13 02045 g005
Table
Table 1. Methylation alignment and efficiency by individual sample and number of differentially methylated fragments (DMFs) using Benjamini–Hochberg correction for each group comparison.
Table 1. Methylation alignment and efficiency by individual sample and number of differentially methylated fragments (DMFs) using Benjamini–Hochberg correction for each group comparison.
Comparison
Groups		IndividualNo. of Unique AlignmentsMapping EfficiencyNo. of DMFs at 5% FDR
RegressedR1 6,942,49844.3%834
R27,861,43048.6%
R37,240,96947.0%
Regressed Control(C)R111,964,23550.6%
(C)R212,197,80749.8%
(C)R39,824,42739.5%
Actively increasing lesionsA17,449,22849.6%0
A210,712,79633.7%
A36,442,79049.5%
Active Control(C)A113,217,13838.3%
(C)A29,198,76049.8%
(C)A310,374,18351.9%
No CsA control(C)16,434,78047.8%0
(C)26,853,51949.1%
CsACsA19,428,60249.5%
CsA212,426,44752.7%
Non-scarifiedCsA1, CsA2 0
Scarified(C)A1, (C)A2, (C)A3
ActiveA1, A2, A3 0
RegressedR1, R2, R3
Table
Table 2. Total viral reads and percent viral reads within E4 region in active, regressed and control tissues.
Table 2. Total viral reads and percent viral reads within E4 region in active, regressed and control tissues.
GroupSampleTotal ReadsViral Reads% Viral Reads% in E4Infection Status
Regressed
lesions		R115,674,8491,129,4277.274Regressed lesion
R216,166,560587100Regressed lesion
R315,392,780819,1145.359.9Regressed lesion
Regressed
Control		(C)R123,666,799371300Mock-infected
(C)R224,496,535286800Mock-infected
(C)R324,877,419414300Mock-infected
Actively increasing
lesions		A115,020,374849,6775.752.1Active lesion
A231,835,71948,5210.211.4Active lesion
A313,025,548708,9765.456Active lesion
Active
Control		(C)A134,512,945624600Mock-infected
(C)A218,479,427253400Mock-infected
(C)A319,996,210247500Mock-infected
CsACsA119,044,226127100CsA/No infection
CsA223,584,967182900CsA/No infection
No CsA
Control		(C)113,470,444110600No CsA/Not inf.
(C)213,946,497128200No CsA/Not inf.
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Tschirley, A.M.; Stockwell, P.A.; Rodger, E.J.; Eltherington, O.; Morison, I.M.; Christensen, N.; Chatterjee, A.; Hibma, M. The Mouse Papillomavirus Epigenetic Signature Is Characterised by DNA Hypermethylation after Lesion Regression. Viruses 2021, 13, 2045. https://doi.org/10.3390/v13102045

AMA Style

Tschirley AM, Stockwell PA, Rodger EJ, Eltherington O, Morison IM, Christensen N, Chatterjee A, Hibma M. The Mouse Papillomavirus Epigenetic Signature Is Characterised by DNA Hypermethylation after Lesion Regression. Viruses. 2021; 13(10):2045. https://doi.org/10.3390/v13102045

Chicago/Turabian Style

Tschirley, Allison M., Peter A. Stockwell, Euan J. Rodger, Oliver Eltherington, Ian M. Morison, Neil Christensen, Aniruddha Chatterjee, and Merilyn Hibma. 2021. "The Mouse Papillomavirus Epigenetic Signature Is Characterised by DNA Hypermethylation after Lesion Regression" Viruses 13, no. 10: 2045. https://doi.org/10.3390/v13102045

APA Style

Tschirley, A. M., Stockwell, P. A., Rodger, E. J., Eltherington, O., Morison, I. M., Christensen, N., Chatterjee, A., & Hibma, M. (2021). The Mouse Papillomavirus Epigenetic Signature Is Characterised by DNA Hypermethylation after Lesion Regression. Viruses, 13(10), 2045. https://doi.org/10.3390/v13102045

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop