Epigenome-wide analysis of common warts reveals aberrant promoter methylation

Epigenetic alteration of host DNA is a common occurrence in both low- and high-risk human papillomavirus (HPV) infection. Although changes in promoter methylation have been widely studied in HPV-associated cancers, they have not been the subject of much investigation in HPV-induced warts, which are a temporary manifestation of HPV infection. The present study sought to examine the differences in promoter methylation between warts and normal skin. To achieve this, DNA was extracted from 24 paired wart and normal skin samples and inputted into the Infinium MethylationEPIC BeadChip microarray. Differential methylation analysis revealed a clear pattern of hyper- and hypomethylation in warts compared to normal skin, and the most differentially methylated promoters were found within the EIF3EP2, CYSLTR1, C10orf99, KRT6B, LAMA4, and H3F3B genes as well as the C9orf30 pseudogene. Moreover, pathway analysis showed that the H3F3A, CDKN1A, and MAPK13 genes were the most common regulators among the most differentially methylated promoters. Since the tissue samples were excised from active warts, however, this differential methylation could either be a cellular response to HPV infection or an HPV-driven process to establish the wart and/or promote disease progression. Conclusively, it is apparent that HPV infection alters the methylation status of certain genes to possibly initiate the formation of a wart and maintain its presence.


Introduction
Epigenetics is the study of heritable changes in gene expression that are not caused by changes to the DNA sequence itself, but by covalent modifications such as DNA methylation (DNA-M) [1]. Mammalian DNA-M, which primarily involves the addition of a methyl group to a cytosine base in a CpG dinucleotide, results in increased gene expression when it occurs at higher levels within the gene's body instead of its promoter region [2]. On a similar note, promoter methylation is of particular epigenetic importance because the vast majority of those located upstream of a gene contain a CpG island, the latter of which is a region with a high concentration of CpG sites [3]. In contrast to the hypermethylated CpG sites scattered throughout the human genome, CpG islands are not methylated, and the methylation of CpG islands initiates remodeling mechanisms that ultimately result in gene silencing [4,5].
The methylation status of promoters is integral to maintaining normal expression levels of the genes they regulate. In fact, promoter hypermethylation is a key part of cancer development and progression, as it results in the silencing of tumor suppressor gene expression [6]. In addition, host promoter hypermethylation has also been implicated in infections by both oncogenic and non-oncogenic Ivyspring International Publisher viruses such as the human papillomaviruses (HPV) [7]. HPV comprises a family of double-stranded DNA viruses that exclusively infect the basal epithelium of the skin and mucosa [8]. The majority of HPV infections are asymptomatic and resolve without the need for medical intervention but, depending on the individual and the HPV type, can also result in a number of malignancies and dermatological diseases [9]. One such condition is the wart, which arises due to the benign proliferation of HPV-infected epithelial keratinocytes [10]. The most prevalent type of wart is the common wart, which accounts for nearly 70% of all cutaneous warts encountered in clinical settings [11]. As a result of their benign nature, common warts are subject to a much lesser degree of scrutiny than other HPV-associated diseases.
The impermanent nature of cutaneous warts strongly suggests that epigenetic changes are involved in the mechanism of wart formation and their eventual disappearance. However, a paucity of information exists regarding the methylation status of cutaneous warts, especially in the context of the promotor regions. Therefore, the primary objective of the current study was to provide an exploratory survey of the genome-wide changes in promoter methylation patterns in cutaneous warts compared to healthy skin.

Study participants
Ethical approval to conduct this study was obtained from Jordan University of Science and Technology's (JUST) Institutional Review Board (IRB). Twelve Arab males presenting with common warts were recruited from the general population after providing written informed consent. Shave biopsies of common warts and adjacent normal skin were performed, allowing paired tissue samples (wart and normal skin) to be obtained from each participant.

Whole genome bisulfite sequencing
A QIAamp DNA Mini Kit (Qiagen, Germany) was used to perform DNA extraction, and optional RNase A digestion was incorporated. DNA purity and integrity were determined by means of the BioTek PowerWave XS2 Spectrophotometer (BioTek Instruments, Inc., USA) and agarose gel electrophoresis, respectively. Genomic DNA that fulfilled our standards for quality and quantity were shipped on dry ice to the Australian Genome Research Facility (AGRF) in Melbourne, where the quality was further ascertained by the QuantiFluor® dsDNA System (Promega, USA). The Zymo EZ DNA Methylation Kit (Zymo Research, USA) was utilized in order to perform bisulfite conversion on the 24 samples. Lastly, the samples were inputted into the Infinium MethylationEPIC BeadChip microarray (Illumina, USA) for a genome-wide interrogation of over 850,000 CpG sites.

Data processing
RnBeads, a computational R package, was adapted to process and analyze the raw intensity data (IDAT files) from the BeadChip [12]. Quality control, preprocessing, batch effects adjustment, and normalization were carried out on all probes and samples according to the RnBeads package pipeline.

Differential methylation and statistical analysis
The mean of the mean β (mean.mean β) values of all the interrogated CpG sites in each promoter were computed. The distribution of CpG sites per promoter is shown in Figure 1, while Figure 2 depicts the distribution of CpG sites across promoters. DM for each promoter was calculated using the following three measures: the mean.mean β difference between warts (W) and normal skin (NS), the log2 of the mean quotient in β means across all CpG sites in a promoter, and the adjusted combined p-value of all CpG sites in the promoter using a limma statistical test [12,13]. Furthermore, these three measures were used to create a combined ranking, in which promoters that exhibit more DM are assigned a lower combined rank [12]. Promoters were sorted from smallest to largest using the combined rank score, and the top-ranking 1000 DM promoters were selected for further analysis. In order to correct for multiple testing, the Benjamini-Hochberg procedure was utilized to set the false discovery rate (FDR) at 5%.

Gene ontology enrichment analysis
Enrichment analysis for gene ontology (GO) terms associated with the top-ranking 500 DM promoters was performed using the GO consortium [14].

Signaling pathway analysis
A signaling network of the top-ranking 1000 DM promoters was investigated using the SIGnaling Network Open Resource (SIGNOR) 2.0 [15]. Due to the large number of connections, the type of relation was selected to only include 'direct' interactions with a relaxed layout and a score of '0.0'.

Sample clustering based on methylation data
Based on all methylation values of the top-ranking 1000 DM promoters, the 24 samples showed an expected clustering pattern, as samples with similar methylation patterns or phenotypes tended to cluster together (Figure 3). In addition, the dimension reduction test was applied to the dataset using multidimensional scaling (MDS) and principal component analysis (PCA) in order to inspect for a strong signal in the methylation values of the samples (Figures 4 and 5). MDS and PCA confirmed that the difference between wart (W) and normal skin (NS) samples predominated the analysis.

Differential methylation of promoters
44,929 genomic identifiers passed the quality control and pre-processing steps, including some identifiers that did not map to gene symbols or which were not assigned (NA). Genomic identifiers without symbols were then removed, leaving 27,790 with symbols. The list of DM promoters in warts was limited to the top-ranking 1000 DM promoters using the combined rank score. Using this scoring method, a total of 576 and 424 promoters were found to be hypomethylated and hypermethylated, respectively, in warts (W) compared to normal skin (NS), with a mean β difference =>0.064 and =< -0.064 and p-value =< 0.001 (adjusted p-value =<0.007) (Figure 6). Among the 576 hypomethylated promoters, the β difference ranged from -0.064 to -0.458, while the mean β difference ranged between 0.064 and 0.367 for the 424 hypermethylated promoters. The log2 of the quotient in methylation between warts and normal skin had a maximum value of 1.633 and minimum value of -1.924 (Figure 7). The top-ranking 100 DM promoters with the lowest combined rank score are shown in Table 1.

Gene ontology enrichment analysis
Gene ontology (GO) enrichment analysis of biological process (BP) and molecular function (MF) was conducted on the top-ranking 500 DM hypermethylated promoters (Figure 8, Figure 9, Table 2, and Table 3) and the top-ranking 500 DM hypomethylated promoters (Figure 10, Figure 11, Table 4, and Table 5).

Pathway analysis
Signaling network analysis of the top-ranking 1000 DM promoters illustrated that several promoter genes were common regulators of this gene network, with a minimum of 7 direct connectivities each. These promoter genes include H3F3A, CDKN1A, MAPK13, IKBKG, CAPN2, CAMKK1 and CUL1 (Figure 12). Moreover, H3F3A was found to be the most common regulator when the signaling network analysis was carried out on the top 100 DM promoters.

Discussion
To the best of the authors' knowledge, this is the first study to investigate the genome-wide changes in promoter methylation patterns associated with HPV-induced cutaneous warts. The present findings provide an exploratory analysis that creates clear lines of future research on this topic, especially with regard to validation studies involving larger sample sizes.
In the present study, the most differentially methylated (DM) promoter in warts compared to normal skin was found within the eukaryotic translation initiation factor 3 subunit E pseudogene 2 (EIF3EP2) gene, a pseudogene with no function or association previously reported in the literature. Likewise, little is known about the second most DM gene in warts, the chromosome 9 open reading frame 30 (C9orf30) pseudogene. In contrast, the third most DM gene is the protein-coding cysteinyl leukotriene receptor 1 (CYSLTR1) gene, which is normally involved in allergic and hypersensitive reactions [16].
Variation in the CYSLTR1 gene modulates asthma risk as well as adenoid hypertrophy progression, and it has been implicated in the disease outcome of colorectal, prostate, and squamous cell carcinoma [17][18][19][20][21]. Moreover, CYSLTR1 is highly expressed in the normal human skin epidermis, but its expression was found to be even higher in atopic dermatitis [22]. Table 2 depicts all the protein-coding genes containing DM promoters from among the top-ranking 100 listed in Table 1.             Figure 13. Pathway signaling network generated from the top-ranking 100 DM promoters.
Among the protein-coding genes, C10orf99 and KRT6B promoters exhibited high levels of differential methylation in warts. The chromosome 10 open reading frame 99 (C10orf99) gene encodes for an antimicrobial peptide that is widely expressed in the skin and digestive tract [23]. In a pathologic context, C10orf99 was determined to contribute to psoriasis development by promoting keratinocyte proliferation [24,25]. Likewise, the keratin 6B (KRT6B) gene encodes for a type II keratin that is normally present in mammalian epithelial cells and is rapidly induced in human keratinocytes after skin wounding [26]. KRT6B has been identified as a potential biomarker for differentiating between lung adenocarcinoma and lung squamous cell carcinoma, and its increased expression is associated with lower disease-free survival rates in young breast cancer patients [27,28]. Mutations in the KRT6B gene result in an autosomal dominant skin disorder known as pachyonychia congenita, which involves plantar keratoderma and pain alongside thickened toenails [29]. In contrast, two of the most differentially methylated protein-coding promoters, namely the kallikrein related peptidase 2 (KLK2) and Izumo sperm-egg fusion 1 (IZUMO1) genes, are integral for sperm function. KLK2 over-expression has been associated with the promotion of prostate cancer cell growth [30].
As previously mentioned, the ephemeral nature of warts hints towards the involvement of an epigenetic component. Correspondingly, some of the most DM promoters were found within the laminin subunit alpha 4 (LAMA4) and H3 histone family member 3B (H3F3B) genes, which are responsible for cell differentiation and nucleosomal displacement, respectively [31,32]. In certain breast cancer subtypes, increased LAMA4 expression was noted to contribute to the chromatin remodeling mechanisms that are a part of cancer progression [33]. Moreover, atypical HF3B expression was reported to be associated with colorectal cancer and chondroblastoma [34,35]. On a similar note, epigenetic modifications have been linked to changes in metabolism in a number of different non-communicable diseases, including cancer and diabetes [36]. In the present study, promoters were differentially methylated within the 17β-hydroxysteroid dehydrogenase type 14 (HSD17B14), leukotriene C4 synthase (LTC4S), folate receptor 3 (FOLR3), alcohol dehydrogenase 7 (ADH7), and adiponectin receptor 2 (ADIPOR2) genes that are involved in steroid, eicosanoid, folate, retinol, and glucose and lipid metabolism, respectively [37][38][39][40][41]. Like the CYSLTR1 gene, LTC4S polymorphisms were associated with asthma risk and drug responsiveness in different ethnic populations [42][43][44][45].
Pathway analysis demonstrated that the most common regulator among the top-ranking 1000 DM promoters was the H3 histone family member 3A (H3F3A) gene. Like the H3F3B gene, H3F3A encodes for a histone variant and is involved in transcriptional regulation [46]. Aberrant H3F3A expression has been associated with the promotion of pediatric and adolescent cancers as well as lung cancer cell migration [46,47]. The second most common regulator was the cyclin dependent kinase inhibitor 1A (CDKN1A) gene, which is mostly involved in CDK2 inhibition and is a primary target of p53 activity [48]. The CDKN1A gene was associated with better patient survival in HPV-related oropharyngeal squamous cell carcinoma [49]. The third most common regulator in HPV-induced warts is the mitogen-activated protein kinase 13 (MAPK13) gene. MAPK13 is a member of the MAP kinase family and functions to regulate cellular responses to a range of different stimuli, especially in the context of keratinocyte apoptosis and skin homeostasis [50]. Analysis of genome-wide promoter methylation revealed that MAPK13 was hypermethylated in the majority of primary and metastatic melanomas [51]. MAPK13 was also found to be hypermethylated in esophageal squamous cell carcinoma [52].
In summary, it is apparent that HPV-induced warts possess certain patterns of promoter methylation that could be responsible for their formation and maintenance. One limitation of the current study is that it is not possible at this stage to determine whether the differential methylation occurred as a result of the host cells' response to infection or due to HPV-driven processes responsible for wart formation and progression. Future research is required in order to assess the functional and clinical importance of the hypo-and hypermethylated promoter sites as well as to determine the exact nature of this differential methylation.