Genome-Wide Tiling Array Analysis of HPV-Induced Warts Reveals Aberrant Methylation of Protein-Coding and Non-Coding Regions

Abstract

The human papillomaviruses (HPV) are a group of double-stranded DNA viruses that exhibit an exclusive tropism for squamous epithelia. HPV can either be low- or high-risk depending on its ability to cause benign lesions or cancer, respectively. Unsurprisingly, the majority of epigenetic research has focused on the high-risk HPV types, neglecting the low-risk types in the process. Therefore, the main objective of this study is to better understand the epigenetics of wart formation by investigating the differences in methylation between HPV-induced cutaneous warts and normal skin. A number of clear and very significant differences in methylation patterns were found between cutaneous warts and normal skin. Around 55% of the top-ranking 100 differentially methylated genes in warts were protein coding, including the EXOC4, KCNU, RTN1, LGI1, IRF2, and NRG1 genes. Additionally, non-coding RNA genes, such as the AZIN1-AS1, LINC02008, and MGC27382 genes, constituted 11% of the top-ranking 100 differentially methylated genes. Warts exhibited a unique pattern of methylation that is a possible explanation for their transient nature. Since the genetics of cutaneous wart formation are not completely known, the findings of the present study could contribute to a better understanding of how HPV infection modulates host methylation to give rise to warts in the skin.

Keywords: HPV; cutaneous warts; methylation; skin.

Figure 1

Distributions of CpG sites (A) per genomic tiling region and (B) across genomic tiling region. In Panel B, the relative coordinates of 0 and 1 correlate to the start and end coordinates of the genomic tiling region. Those coordinates that are smaller than 0 and larger than 1 indicate flanking regions normalized by region length.

Figure 2

Comparisons of the density distributions of methylation levels (β) in warts (W) and normal skin (NS).

Figure 3

For the 1000 most differentially methylated tiling regions, scatterplot (A) and (B) volcano plot analyses were performed. In panel A, the mean of mean methylation levels (β) for warts (W) is on the y-axis, while the mean of mean methylation levels (β) for normal skin (NS) is on the x-axis. β values range from unmethylated (0) to methylated (1). In panel B, the volcano plot shows the differential methylation of genomic tiling regions as quantified by log₂ of the mean quotient in means across all sites in a region on the x-axis and the adjusted combined p-value on the y-axis between warts (W) and normal skin (NS). The color scale corresponds with the combined rank of each genomic tiling region.

Figure 4

Heatmap showing the hierarchical clustering of samples displaying only the 1000 most variable loci with the highest variance across all samples. Clustering used complete linkage and Manhattan distance. Patient identification number is shown on the bottom x-axis. The normal skin (NS) and wart (W) samples are shown on the top x-axis. Values of 0 (red color) and 1 (purple color) indicate decreased and increased methylation, respectively.

Figure 5

Scatter plot showing samples after performing Kruskal’s non-metric multidimensional scaling based on the matrix of average methylation levels and Manhattan distance.

Figure 6

Scatterplot of LOLA enrichment analysis showing the effect size (log-odds ratio) versus the significant q-value (−log10 (q-value)) of the 1000 most hypermethylated tiling regions. These regions show strong association and significant overlap with Sheffield_dnase and encode_tfbs as indicated by the large odds ratio and higher q-value. LOLA reference databases collections are shown on the right side of the plot with color coding.

Figure 7

Scatter plot of LOLA enrichment analysis showing the effect size (log-odds ratio) versus the significant q-value (−log10 (q-value)) of the 1000 most hypomethylated tiling regions. These regions show strong association and significant overlap with Sheffield_dnase, ucsc_features, codex, cistrome_epigenome, encode_tfbe, and encode_segmentation as indicated by the large odds ration values and the higher q−value. LOLA reference databases collections are shown on the right side of the plot with color coding.

Figure 8

Bar plots of LOLA enrichment analysis showing the log-odds ratios of the 1000 most hypermethylated tiling regions. Terms that exhibit statistical significance (p-value < 0.01) are shown. For encode_tfbs, c-Fos, STAT3, and c-Myc are among the top enriched terms showing large odds ratio values. Fibroblast cells, epithelial cells, muscle tissue, and skin tissue were among the most enriched cell and tissue types. Coloring of the bars reflects the putative targets of the terms.

Figure 9

Bar plots of LOLA enrichment analysis showing log-odds ratios of the 1000 most hypomethylated tiling regions. Terms that exhibit statistical significance (p-value < 0.01) are shown. The most enriched terms include Dnase weak-NHDF, Dnase-fibrobalsts, Weak Enhancer-HUVEC, NR2F2-Endometrial Stromal Cell, AFF1-leukaemia, H3K4me1-LNCaP, and androgen receptor (AR)-abl. Coloring of the bars reflects the putative targets of the terms.

Figure 10

Pathway signaling network generated from the genes located within the top 100 DM tiling regions.