Skip to main content
Download PDF
  • Research article
  • Open access
  • Published:

Transcriptome analysis of HPV-induced warts and healthy skin in humans

BMC Medical Genomics volume 13, Article number: 35 (2020) Cite this article

Abstract

Background

The human papillomaviruses (HPV) are a group of viruses that, depending on the strain, can cause cancer or the formation of benign growths known as warts. Scarce information exists with regard to the genetic nature of non-genital cutaneous warts induced by the human papillomavirus (HPV).

Methods

The main purpose of this study is to investigate the differences between the gene expression profiles of common warts and healthy skin in HPV-positive individuals by RNA sequencing on the Illumina HiSeq 2500. After obtaining shave biopsies of common warts and healthy skin from twelve Arab males, we were able to analyze the transcriptomes of 24 paired cases and controls.

Results

Common warts were found to possess a highly significant and unique molecular signature. Many of the most up-regulated (KRT16, EPGN, and ABCG4) and down-regulated genes (C15orf59, CYB561A3, and FCGRT) in warts were the subject of little investigation in the published literature. Moreover, the top 500 differentially expressed genes were found to be associated with immune and autoimmune pathways, such as the neutrophil degranulation, toll-like receptor 7/8 (TLR 7/8) cascade, toll-like receptor 9 (TLR9) cascade, and toll-like receptor 10 (TLR10) pathways, among others.

Conclusions

Our findings are particularly important because they serve as the most comprehensive to date with regard to the modulation of human skin gene expression by HPV infection.

Peer Review reports

Background

The human papillomavirus (HPV) is a DNA virus that has been associated with many diseases and is the definitive cause of virtually all cervical cancer cases [1]. Nearly two hundred types of HPV have been characterized, and these types are classified as high-risk or low-risk depending on their potential to cause cancerous or benign lesions, respectively [2]. 4.5% of new cancer cases worldwide are caused by high-risk HPV infection, but low-risk types often manifest in the form of warts [3, 4]. Common warts (Verruca vulgaris) are by far the most prevalent type of wart, making up 70% of all non-genital cutaneous warts, and are thought to be benign in nature [5, 6].

Human skin is composed of two layers, the epidermis and the dermis, and HPV exclusively targets the lowermost epidermal layer, the stratum basale, in order to establish a persistent infection, after which it hijacks the keratinocyte differentiation process for the purpose of the productive viral life cycle [7, 8]. Correspondingly, HPV biology differs in basal and suprabasal keratinocytes, the latter of which are involved in different stages of epithelial differentiation [9]. Keratinocytes make up 90% of all epidermal cells, and perturbations to their differentiation make the skin more susceptible to infection and disease [10]. Currently, there is no cure for HPV itself, and its treatment focuses on the alleviation of clinical symptoms until the infection is naturally cleared by the immune system [11, 12]. However, there is no conclusive evidence as to whether the virus is completely eliminated from the body or simply lowered to undetectable levels [12].

While much research has been dedicated to the genetics of HPV-associated cancers, there is a dearth of information concerning the genetic background of non-genital cutaneous warts, which indicates a pressing need to better understand their etiology. Therefore, the aim of the present study is to provide a genome-wide comparison of the transcriptomes of common warts and normal skin using next-generation sequencing.

Methods

Subjects

The present study was approved by the Institutional Review Board (IRB) at Jordan University of Science and Technology (JUST) (Ref. 13/105/2017). All procedures performed in this study were in compliance with the IRB’s human research ethics protocol and were performed at King Abdullah University Hospital in Irbid, Jordan. Twelve Arab males presenting with common warts were recruited from the general population after giving written informed consent. Clinical diagnosis was established by a dermatologist based on typical clinical features for warts, especially keratotic surfaces and black thrombosed capillaries or pin-point bleeding upon paring. Epidermal tissue samples from warts (n = 12) and normal skin (controls) (n = 12) were obtained via a superficial shave biopsy after the application of a local anesthetic. Normal skin tissue samples were procured from a site adjacent to the area of the wart excision. Most samples came from the dorsal aspect of the hands (n = 20), while a minority were obtained from the forehead (n = 2) and the dorsal aspect of the foot (n = 2). Biopsies were immediately stored in microcentrifuge tubes containing RNA-SafeGuard reagent (GMbiolab Co., Ltd., Taiwan) and refrigerated at − 20 °C until further processing.

RNA extraction

Total RNA was extracted from the epidermal samples by means of an RNeasy Mini Kit (Qiagen, Germany), and optional on-column DNase digestion was performed. The quality (260/280), quantity (ng/ml), and RNA integrity numbers (RIN) of the purified RNA was then ascertained using the Biotek PowerWave XS2 Spectrophotometer (BioTek Instruments, Inc., USA) and the Agilent Bioanalyzer (Agilent, USA). All samples selected for sequencing had a RIN value of greater than 7.5.

RNA sequencing (RNA-Seq)

Samples were shipped on dry ice to the Melbourne node of the Australian Genome Research Facility (AGRF). RNA-Seq was carried out on the Illumina HiSeq 2500 according to the manufacturer’s protocol. Image analysis was performed in real time by the HiSeq Control Software (HCS) v2.2.68 and Real Time Analysis (RTA) v1.18.66.3 running on the instrument computer. Then, the Illumina bcl2fastq 2.20.0.422 pipeline was used to generate the sequence data. All generated data met AGRF quality standards.

Bioinformatics analysis

The sequence reads from all 24 samples were analyzed and screened according to AGRF quality control measures. edgeR version 3.22.3 was used to detect and quantify differential expression of digital gene expression data (counts of reads mapped for each gene). Unsupervised clustering of the 24 samples was carried out using a multi-dimensional scale (MDS) plot for raw data and normalized data. The trimmed mean of M-values (TMM) normalisation method was performed to overcome library size differences among samples. The TMM normalisation was done on feature counts data, the latter of which was then transformed to log cpm (counts per million) values to normalize for the different sequencing depths for each sample (Supplementary Figure S1). A generalised linear model was then used to quantify the differential expression between the groups.

Functional enrichment analysis of DE genes was performed using the Database for Annotation, Visualization and Integrated Discovery (DAVID) version 6.8 (Leidos Biomedical Research, Inc., USA). Pathways and network patterns associated with DE genes were analyzed using Reactome version 65 (https://reactome.org) and corrected for false discovery rate (adjusted p-value) using the Benjamani-Hochberg method. The signaling network of DE genes was explored using the Signaling Network Open Resource 2.0 (Signor). The ReactomePA pathway analyser was used to enrich for the list of pathways associated with the top 500 DE genes. Next, ComplexHeatmap version 1.18.1 was used to generate heatmaps for each comparison showing the clustering of each sample present in the comparison using the top 50 most differentially expressed (DE) genes for that comparison.

Moreover, cell type enrichment analysis was performed for infiltrating immune cell genes with FDR > 0.05. ClustVis (http://biit.cs.ut.ee/clustvis/) was used to create a heatmap that compared between the infiltrating immune cell genes in warts and normal skin. In addition, we performed immune cell type enrichment analysis data for genes with FDR > 0.05. The analysis was conducted using the xCell webtool (http://xCell.ucsf.edu/) and xCell gene signatures as references of immune cells [14]. The generated adjusted cell type enrichment scores were used to create a heatmap comparing infiltrating immune cell genes in warts and normal skin using the ClustVis webtool (http://biit.cs.ut.ee/clustvis/).

Lastly, HPV subtypes were identified from the RNA-seq data by utilizing the web-based Viral Genome-Targeted Assembly Pipeline (VirusTAP) [15]. VirusTAP performs quality trimming and subtracts unnecessary reads such as rRNA, bacterial genomes, and the host genome. Afterwards, de novo transcriptome assembly and megablast homology search were performed of contigs against the NCBI database. Finally, the trimmed reads were mapped back to the assembled contigs.

Results

HPV genotyping

A total of 5 HPV subtypes were detected in the RNA-seq data, including HPV 57 (Alphapapillomavirus 4), HPV 2 (Alphapapillomavirus 4), HPV 27 (Alphapapillomavirus 4), HPV 4 (Gammapapillomavirus 1), and HPV XS2 (Alphapapillomavirus 2). Tables 1 and 2 show the results of HPV genotyping for the control and wart samples. Interestingly, no contigs related to virus sequences were found in one of the wart samples, 9 W, while one of the control samples, 4C, was shown to have HPV XS2.

Table 1 HPV typing and location of wart (W) samples
Full size table
Table 2 HPV typing and location of control (C) samples
Full size table

Differentially expressed (DE) genes

Unsupervised clustering revealed clear differences between the control (normal skin) and wart samples (Fig. 1). The wart samples were grouped together and separated from control samples based on the similarity of gene expression profiles. This confirms the paired nature of the samples. The wart samples appear more heterogeneous than the normal skin samples. However, sample 4C, which is the only control sample to contain an HPV type, was further away than the rest of the control samples. In addition, sample 6 W, which contained HPV 4, was located far from the rest of the clustered wart samples. This was to be expected as HPV4 is a Gammapapillomavirus species, while the rest of the identified HPV types were Alphapapillomaviruses. Although samples 9 and 12 came from the forehead and foot, respectively, they were located within each cluster alongside the samples obtained from the hand area.

Fig. 1
figure 1

Multidimensional scaling (MDS) plot showing variation among samples based on (a) raw data and (b) normalized data. Each point represents 1 sample, and the distance between 2 points reflects the leading logFC of the corresponding RNA-Seq samples. The leading logFC (base 2 logarithm of fold change) is the average of the largest absolute logFC between each pair of samples. The plot dimension 1 (dim 1) illustrates that the control (C) and wart (W) samples form separate clusters with certain dispersion among samples

Full size image

From among 13,574 tested genes, 8007 were shown to be DE (adjusted p-value (AP) < 0.05) between warts and normal tissue. The 8007 genes were then filtered out by fold change (FC), resulting in a total of 3140 DE genes (FC > 2) (Fig. 2). Of these DE genes, a total of 1732 (55%) and 1408 (45%) genes were found to be upregulated and downregulated, respectively (Fig. 2). Supplementary Table S1 lists the top 100 DE genes in warts, while Supplementary Figure S2 shows the overall heatmap of the top 500 most DE genes. In contrast to Supplementary Figure S2, Fig. 3 shows the hierarchal clustering heatmap for the top 50 most DE genes in warts and normal skin. The similarity in clustering patterns within each of the wart and control (normal skin) groups indicate that there are distinct expression profiles for warts.

Fig. 2
figure 2

Smear plot of DE genes in warts. The logFC (base 2 logarithm of fold change) for each gene is plotted against the average abundance in logCPM (base 2 logarithm of counts per million). The two horizontal blue lines represent two fold changes. The red points are those with an AP of less than 0.05, while the black points are those with an AP of more than 0.05

Full size image
Fig. 3
figure 3

Heatmap of the top 50 DE genes between warts (W) and normal skin (C). The hierarchically clustered genes are represented by rows, and the samples are represented by columns, while the dendrograms and flat clusters are symbolized by the green and orange bars. Genes that have positively correlated logCPM values cluster together, as large positive correlations correspond to small distances. The red and blue colors indicate gene upregulation and downregulation, respectively

Full size image

Pathway and functional enrichment analysis

Enrichment analysis of the DE genes revealed 71 functional annotation clusters, and the clusters with an enrichment score > 2 and GO term p-value < 0.05 are shown in Table S2. The top enriched terms include the cornified envelope, keratinization, keratinocyte differentiation, and peptide cross-linking. Moreover, by enriching for only the top 500 DE genes and using a cut-off p-value of 0.05, 10 pathways of interest were reported (Fig. 4). Pathway analysis of the top 1000 DE genes was performed in order to determine any significant gene-protein associations. It was found that the CDKN1A (logFC = − 1.18; AP = 0.00015), MAPK14 (logFC = − 1.27; AP = 1.16 × 10− 10), MAP3K5 (logFC = − 1.43; AP = 4.97 × 10− 8), and PPP2CA (logFC = − 1.2; AP = 8.21 × 10− 8) genes had the highest number of interactions with other DE genes in wart (Fig. 5). In fact, functional enrichment analysis showed that CDKN1A is involved in the protein serine/threonine kinase activity pathway, MAPK14 is involved in the protein serine/threonine kinase and protein kinase activity pathways, and MAP3K5 is involved in the protein phosphorylation, protein serine/threonine kinase, and protein kinase activity pathways.

Fig. 4
figure 4

Top 10 pathways associated with the top 500 DE genes in terms of adjusted p-value

Full size image
Fig. 5
figure 5

Pathway network generated from the top 1000 DE genes demonstrating the various interactions between proteins. Four genes have at least five connections (MAPK14, PPP2CA, CDKN1A, and MAP3K5)

Full size image

To identify the immune gene signature, we performed cell type enrichment analysis and heatmap clustering of the genes in the present study with an FDR cut-off value of 0.05. We found 26 immune cell populations which are differentially activated in warts compared to normal skin (Fig. 6). The immune cell signature dendrogram appears to be split between natural killer T (NKT) cells, immature dendritic cells (iDCs), class-switch memory B-cells, conventional dendritic cells (cDCs), regulatory T cells (Tregs), and gamma delta T (Tgd) cells with varying expression of these in the samples.

Fig. 6
figure 6

Heatmap of immune cell infiltrates expressed in warts (W) compared to normal skin samples (C). Both rows and columns are clustered using correlation distance and average linkage. The dendrogram at the top shows the samples; the dendrogram at the right shows different immune cell signatures. Blue color indicates lower expression of a particular immune cell signature and red color indicates higher expression

Full size image

Discussion

Although warts are themselves benign, HPV infection has been identified as the definitive cause of cervical cancer and has been associated with a number of other cancers and diseases [16]. As a result, the majority of peer-reviewed studies conducted on HPV infection have been in the context of cancer research. However, HPV has been detected in a wide plethora of skin conditions other than cancer, including actinic keratoses, epidermal cysts, lichen sclerosus, psoriatic plaques, seborrheic keratoses, and skin tags, and it has also been detected on healthy skin and plucked hairs [17]. This ubiquity of HPV in human physiology and pathology mandates that its role in non-genital cutaneous warts be elucidated.

With regard to upregulated genes, KRT16 was found to be the most upregulated gene in warts. KRT16 encodes for keratin 16, a type I cytokeratin that plays an essential role in the barrier function, innate immunity, and network signaling of the skin epidermis as well as the integrated stress response shared by all eukaryotes [18, 19]. Aberrant KRT16 expression has been implicated in the as a number of inflammatory skin diseases as well as breast and pancreatic cancer [13, 20, 21]. In fact, KRT16 expression was increased by two-fold in keratinocytes infected with high-risk HPV [22]. In addition, EPGN was shown to be the second most upregulated gene in the present study and encodes for epithelial mitogen, an epidermal growth factor whose biologic function was still unknown [23]. The role of the EPGN gene in keratinocytes and warts have not yet been considered.

In terms of downregulated genes, this study found that the chromosome 15 open reading frame 59 (C15orf59) is the most downregulated gene in warts. Very little is known about the C15orf59 gene either in terms of its function or its role in disease. The second most down-regulated gene is the cytochrome b561 family member A3 (CYB561A3) gene, which is normally expressed intercellularly as well as at the cell membrane and was found to be down-regulated in sepsis [24]. Compared to the C15orf59 and CYB561A3 genes, much more is known about the third most down-regulated gene, Fc fragment of IgG receptor and transporter (FCGRT), which encodes for the neonatal Fc receptor and is normally expressed in human epidermal keratinocytes [25]. In cell lines infected with high-risk HPV, FCGRT expression on epithelial keratinocytes was found to modulate the uptake of IgG antibodies [26].

In this study, the top wart-associated pathways were involved in the immune system, pointing towards the effects of HPV infection on immune function. The most significantly associated pathway was neutrophil degranulation, a component of the innate immune system which was found to be a typical feature of inflammatory disorders such as asthma and acute lung injury [27]. The second and third most significantly associated pathways, formation of the cornified envelope and keratinization, are involved in the formation of wart hyperkeratosis, which in itself is a form of programmed cell death or apoptosis [28]. The toll-like receptor 7/8 (TLR7/8) cascade, the fourth most significantly associated pathway, plays an integral role in innate immunity, but its aberrant expression was significantly associated with persistent HPV16 infection [29]. In addition, the TLR9, TLR10, and TLR5 cascades were associated with the DE genes in our study, and they were similarly dysregulated in HPV infection and cervical cancer progression [29, 30]. Two of the top ten pathways involved the myeloid differentiation factor 88 (MyD88), a protein involved in innate immune function that was upregulated in HPV-transformed cervical cancer cell lines [31].

Moreover, the immune gene expression signature showed a list of immune cells which are differentially activated in warts compared to normal skin. Natural killer T (NKT) cells, a key component of the innate immune system, are important in killing certain tumour cells as well as virally infected cells, and they have been associated with HPV infection [32]. Similarly, regulatory T (Treg) cells are considered to aid in HPV growth by playing a role in immunosuppression [33]. Additionally, three dendritic cell (DC) subsets, namely immature dendritic cells (iDCs), conventional dendritic cells (cDCs) and plasmacytoid dendritic cells (pDCs), showed low activity in wart samples compared to normal skin. DCs have an essential role in initiating and regulating innate and adaptive immune responses to viral infection [34].

A pathway analysis of all 3140 DE genes found that four genes were central to the pathway: cyclin dependent kinase inhibitor 1A (CDKN1A), mitogen-activated protein kinase 14 (MAPK14), mitogen-activated protein kinase kinase kinase 5 (MAP3K5), and protein phosphatase 2 catalytic subunit alpha (PPP2CA). Our study showed that CDKN1A had the highest number of interactions with the other DE genes in this study and was itself downregulated in warts. The CDKN1A gene, which encodes for the cyclin-dependent kinase inhibitor p21, is a transcriptional target of the tumor suppressor protein p53 and functions to induce cell cycle arrest in the event of DNA damage [35, 36]. In the context of HPV type, it appears to be that the alphapapillomaviruses, including the low-risk HPV2 and HPV57 and the high-risk HPV15 and HPV18, target the p53-p21 pathway regardless of oncogenic potential [37, 38].

Additionally, MAPK14 has been documented as having an important role in cell survival and apoptosis as well as differentiation and proliferation [39]. In fact, MAPK14 plays a dual role in the skin, where it is responsible both for evoking an inflammatory reaction and for containing the aforementioned inflammation [40]. Moreover, MAP3K5, which is also known as apoptosis signal-regulating kinase 1 (ASK1), responds to inflammatory and physical stressors by activating the JNK and p38 pathways, and it plays a key role in hair regrowth in wounded skin [41]. Lastly, PPP2CA has been determined to play a critical role in the development of both the epidermis and the hair follicles [42].

Several strengths can be attributed to the present study that consolidate the significance and accuracy of its findings. All of our patients were Jordanian-Arab males ranging from 18 to 27 years old, and both the case (wart) and control (normal skin) biopsies were obtained from the same patient. None of the warts biopsied for this study have been modified by any type of treatment or medicine. Moreover, case and control biopsies for ten out of twelve patients in this study were carried out on the same area of the body, i.e. the hands, further reducing the amount of variability between subjects. Finally, all the warts included in this study were of the same type, i.e. Verruca vulgaris.

However, some limitations exist regarding the current study. Although only common warts (Verruca vulgaris) were investigated, HPV-induced warts are heterogeneous in nature, and several different HPV types have been associated with common warts, including HPV 2, 4, 7, 27, and 57 [43]. In addition, the HPV genome is itself highly heterogenous due to the fact that each HPV type has evolved to maximize viral fitness to its biological ecosystem within the host, i.e. niche adaptation [44, 45]. In the present study, HPV typing revealed that the investigated warts were caused by HPV types belonging to the Alphapapillomavirus 2 (HPV XS2), Alphapapillomavirus 4 (HPV 2/27/57), and Gammapapillomavirus 1 (HPV4) species. Another limitation is the fact that histological analysis of the wart samples was not possible, as the sample biopsies were small and entirely used up during the RNA extraction process. Since HPV biology differs between epidermal layers, this means that the present RNA-seq data contains a mixture of basal and suprabasal expression patterns. A final limitation is the relatively small sample size, but this is in line with the sample sizes of other RNA-seq studies, which included as little as four patients compared to our twelve [46].

Conclusion

The present study served to clarify the genetic background of HPV infection in the context of non-genital cutaneous warts. The findings of the current study highlighted dysregulation of multiple genes which might play a critical role in wart formation, including KRT16, EPGN, C15orf59, CYB561A3, and FCGRT. Furthermore, several immune system related pathways were found to be associated with identified DE genes including neutrophil degranulation, myeloid differentiation factor 88, and toll-like receptors.

Our data also suggests that several genes involved in the immune response may play a critical role in facilitating the HPV infection process. Indeed, immune cell signature analysis showed different activity levels of multiple cells involved in the immune system including NKT, DC, and Treg cells. Many of the genes found to be implicated in common warts were not the subject of previous peer-reviewed investigation, suggesting future directions and lines of research targeting those individual genes as well as the expression patterns of individual HPV types.

Availability of data and materials

The data generated over the course of the present study has been deposited to NCBI’s Gene Expression Omnibus (GEO) under the GEO accession number of GSE136347.

Abbreviations

AGRF:

Australian Genome Research Facility

BCC:

Basal cell carcinoma

DE:

Differentially expressed

DNA:

Deoxyribonucleic acid

FC:

Fold change

HPV:

Human papilloma virus

IL:

Interleukin

IRB:

Institutional Review Board

PCR:

Polymerase chain reaction

RIN:

RNA Integrity Number

RNA:

Ribonucleic acid

RNA-Seq:

RNA-Sequencing

SNP:

Single nucleotide polymorphism

References

  1. Ganguly P, Ganguly N. Transcriptomic analyses of genes differentially expressed by high-risk and low-risk human papilloma virus E6 oncoproteins. VirusDisease. 2015;26:105–16.

    Article  Google Scholar 

  2. Egawa N, Doorbar J. The low-risk papillomaviruses. Virus Res. 2017;231:119–27.

    Article  CAS  Google Scholar 

  3. Serrano B, Brotons M, Bosch FX, Bruni L. Epidemiology and burden of HPV-related disease. Best Pract Res Clin Obstet Gynaecol. 2018;47:14–26.

    Article  Google Scholar 

  4. Doorbar J, Egawa N, Griffin H, Kranjec C, Murakami I. Human papillomavirus molecular biology and disease association. Rev Med Virol. 2015;25:2–23.

    Article  CAS  Google Scholar 

  5. Mammas IN, Spandidos DA, Sourvinos G. Genomic diversity of human papillomaviruses (HPV) and clinical implications: an overview in adulthood and childhood. Infect Genet Evol. 2014;21:220–6.

    Article  Google Scholar 

  6. Mulhem E, Pinelis S. Treatment of nongenital cutaneous warts. Am Fam Physician. 2011;84:288–93.

    PubMed  Google Scholar 

  7. Yousef YA, Tbakhi A, Al-Hussaini M, AlNawaiseh I, Saab A, Afifi A, et al. Mutational analysis of the RB1 gene and the inheritance patterns of retinoblastoma in Jordan. Familial Cancer. 2018;17:261–8.

    Article  CAS  Google Scholar 

  8. McBride AA. Mechanisms and strategies of papillomavirus replication. Biol Chem. 2017;398:919–27.

    Article  CAS  Google Scholar 

  9. White EA. Manipulation of Epithelial Differentiation by HPV Oncoproteins. Viruses. 2019;11:369.

  10. Di Meglio P, Perera GK, Nestle FO. The multitasking organ: recent insights into skin immune function. Immunity. 2011;35:857–69.

    Article  Google Scholar 

  11. Perfetto L, Briganti L, Calderone A, Cerquone Perpetuini A, Iannuccelli M, Langone F, et al. SIGNOR: a database of causal relationships between biological entities. Nucleic Acids Res. 2016;44:D548–54. https://doi.org/10.1093/nar/gkv1048.

    Article  CAS  PubMed  Google Scholar 

  12. Miranda PM, Silva NNT, Pitol BCV, Silva IDCG, Lima-Filho JL, Carvalho RF, et al. Persistence or clearance of human papillomavirus infections in women in Ouro Preto. Brazil Biomed Res Int. 2013;2013:578276.

    CAS  PubMed  Google Scholar 

  13. D’Erme AM, Wilsmann-Theis D, Wagenpfeil J, Hölzel M, Ferring-Schmitt S, Sternberg S, et al. IL-36γ (IL-1F9) is a biomarker for psoriasis skin lesions. J Invest Dermatol. 2015;135:1025–32.

    Article  Google Scholar 

  14. Aran D, Hu Z, Butte AJ. xCell: digitally portraying the tissue cellular heterogeneity landscape. Genome Biol. 2017;18:220.

    Article  Google Scholar 

  15. Yamashita A, Sekizuka T, Kuroda M. VirusTAP: Viral Genome-Targeted Assembly Pipeline. Front Microbiol. 2016;7:32.

  16. Chow LT, Broker TR. Human papillomavirus infections: warts or cancer?. Cold Spring Harb Perspect Biol. 2013;5:a012997.

  17. Ljubojevic S, Skerlev M. HPV-associated diseases. Clin Dermatol. 2014;32:227–34.

    Article  Google Scholar 

  18. Lessard JC, Piña-Paz S, Rotty JD, Hickerson RP, Kaspar RL, Balmain A, et al. Keratin 16 regulates innate immunity in response to epidermal barrier breach. Proc Natl Acad Sci U S A. 2013;110:19537–42.

    Article  CAS  Google Scholar 

  19. Evstafieva AG, Kovaleva IE, Shoshinova MS, Budanov AV, Chumakov PM. Implication of KRT16, FAM129A and HKDC1 genes as ATF4 regulated components of the integrated stress response. PLoS One. 2018;13:e0191107.

    Article  Google Scholar 

  20. Mao Y, Shen J, Lu Y, Lin K, Wang H, Li Y, et al. RNA sequencing analyses reveal novel differentially expressed genes and pathways in pancreatic cancer. Oncotarget. 2017;8:42537–47.

    PubMed  PubMed Central  Google Scholar 

  21. Joosse SA, Hannemann J, Spotter J, Bauche A, Andreas A, Muller V, et al. Changes in keratin expression during metastatic progression of breast Cancer: impact on the detection of circulating tumor cells. Clin Cancer Res. 2012;18:993–1003.

    Article  CAS  Google Scholar 

  22. Vliet-Gregg PA, Hamilton JR, Katzenellenbogen RA. NFX1-123 and human papillomavirus 16E6 increase notch expression in keratinocytes. J Virol. 2013;87:13741–50.

    Article  CAS  Google Scholar 

  23. Schneider MR, Yarden Y. Structure and function of epigen, the last EGFR ligand. Semin Cell Dev Biol. 2014;28:57–61.

    Article  CAS  Google Scholar 

  24. Zhang J, Cheng Y, Duan M, Qi N, Liu J. Unveiling differentially expressed genes upon regulation of transcription factors in sepsis. 3 Biotech. 2017;7:46.

    Article  Google Scholar 

  25. Kuo TT, Baker K, Yoshida M, Qiao S-W, Aveson VG, Lencer WI, et al. Neonatal fc receptor: from immunity to therapeutics. J Clin Immunol. 2010;30:777–89.

    Article  CAS  Google Scholar 

  26. Wang H-K, Wei Q, Moldoveanu Z, Huh WK, Vu HL, Broker TR, et al. Characterization of serum antibodies from women immunized with Gardasil: a study of HPV-18 infection of primary human keratinocytes. Vaccine. 2016;34:3171–7.

    Article  CAS  Google Scholar 

  27. Lacy P. Mechanisms of degranulation in neutrophils. Allergy Asthma Clin Immunol. 2006;2:98–108.

    Article  Google Scholar 

  28. Eckhart L, Lippens S, Tschachler E, Declercq W. Cell death by cornification. Biochim Biophys Acta - Mol Cell Res. 1833;2013:3471–80.

    Google Scholar 

  29. Daud II, Scott ME, Ma Y, Shiboski S, Farhat S, Moscicki A-B. Association between toll-like receptor expression and human papillomavirus type 16 persistence. Int J Cancer. 2011;128:879–86.

    Article  CAS  Google Scholar 

  30. DeCarlo CA, Rosa B, Jackson R, Niccoli S, Escott NG, Zehbe I. Toll-like receptor transcriptome in the HPV-positive cervical cancer microenvironment. Clin Dev Immunol. 2012;2012:785825.

    Article  Google Scholar 

  31. Morale MG, da Silva AW, Silva AM, Villa LL, Boccardo E. HPV-transformed cells exhibit altered HMGB1-TLR4/MyD88-SARM1 signaling axis. Sci Rep. 2018;8:3476.

    Article  Google Scholar 

  32. Orange JS, Ballas ZK. Natural killer cells in human health and disease. Clin Immunol. 2006;118:1–10.

    Article  CAS  Google Scholar 

  33. Molling JW, de Gruijl TD, Glim J, Moreno M, Rozendaal L, Meijer CJLM, et al. CD4+CD25hi regulatory T-cell frequency correlates with persistence of human papillomavirus type 16 and T helper cell responses in patients with cervical intraepithelial neoplasia. Int J Cancer. 2007;121:1749–55.

    Article  CAS  Google Scholar 

  34. Belz G, Mount A, Masson F. Dendritic Cells in Viral Infections. In: Dendritic Cells. Berlin, Heidelberg: Springer; 2009. p. 51–77.

    Chapter  Google Scholar 

  35. Benson EK, Mungamuri SK, Attie O, Kracikova M, Sachidanandam R, Manfredi JJ, et al. p53-dependent gene repression through p21 is mediated by recruitment of E2F4 repression complexes. Oncogene. 2014;33:3959–69.

    Article  CAS  Google Scholar 

  36. Bunz F, Dutriaux A, Lengauer C, Waldman T, Zhou S, Brown JP, et al. Requirement for p53 and p21 to sustain G2 arrest after DNA damage. Science. 1998;282:1497–501.

    Article  CAS  Google Scholar 

  37. Fischer M, Uxa S, Stanko C, Magin TM, Engeland K. Human papilloma virus E7 oncoprotein abrogates the p53-p21-DREAM pathway. Sci Rep. 2017;7:2603.

  38. Gheit T. Mucosal and cutaneous human papillomavirus infections and cancer biology. Front Oncol. 2019;9:355.

  39. Young PR. Perspective on the discovery and scientific impact of p38 MAP kinase. J Biomol Screen. 2013;18:1156–63.

    Article  CAS  Google Scholar 

  40. Kim C, Sano Y, Todorova K, Carlson BA, Arpa L, Celada A, et al. The kinase p38α serves cell type–specific inflammatory functions in skin injury and coordinates pro- and anti-inflammatory gene expression. Nat Immunol. 2008;9:1019–27.

    Article  CAS  Google Scholar 

  41. Osaka N, Takahashi T, Murakami S, Matsuzawa A, Noguchi T, Fujiwara T, et al. ASK1-dependent recruitment and activation of macrophages induce hair growth in skin wounds. J Cell Biol. 2007;176:903–9.

    Article  CAS  Google Scholar 

  42. Fang C, Li L, Li J. Conditional knockout in mice reveals the critical roles of Ppp2ca in epidermis development. Int J Mol Sci. 2016;17:756.

    Article  Google Scholar 

  43. Bacaj P, Burch D. Human papillomavirus infection of the skin. Arch Pathol Lab Med. 2018;142:700–5.

    Article  CAS  Google Scholar 

  44. Harari A, Chen Z, Burk RD. Human papillomavirus genomics: past, present and future. Current Problems in Dermatology (Switzerland). 2014;45:1–18.

    Article  Google Scholar 

  45. Chen Z, DeSalle R, Schiffman M, et al. Niche adaptation and viral transmission of human papillomaviruses from archaic hominins to modern humans. PLoS Pathog. 2018;14:e1007352.

  46. Swindell WR, Sarkar MK, Liang Y, Xing X, Baliwag J, Elder JT, et al. RNA-seq identifies a diminished differentiation gene signature in primary monolayer keratinocytes grown from lesional and uninvolved psoriatic skin. Sci Rep. 2017;7:18045.

    Article  Google Scholar 

Download references

Acknowledgements

The authors are grateful to all the participants of this study for their invaluable contribution. This study was funded by Jordan University of Science and Technology. The authors would like also to express their gratitude to King Khalid University, Saudi Arabia for providing administrative and technical support.

Funding

This work was supported by the Deanship of Research at Jordan University of Science and Technology under grant number (Ref # 177/2017).

Author information

Authors and Affiliations

  1. Department of Applied Biological Sciences, Jordan University of Science and Technology, Irbid, 22110, Jordan

    Laith N. AL-Eitan & Amneh H. Tarkhan

  2. Department of Biotechnology and Genetic Engineering, Jordan University of Science and Technology, Irbid, 22110, Jordan

    Laith N. AL-Eitan

  3. Department of Anatomy, College of Medicine, King Khalid University, Abha, 61421, Saudi Arabia

    Mansour A. Alghamdi

  4. Genomics and Personalized Medicine Unit, College of Medicine, King Khalid University, Abha, 61421, Saudi Arabia

    Mansour A. Alghamdi

  5. Department of Internal Medicine, Jordan University of Science and Technology, Irbid, 22110, Jordan

    Firas A. Al-Qarqaz & Hadeel S. Al-Kofahi

  6. Division of Dermatology, Department of Internal Medicine, King Abdullah University Hospital, Jordan University of Science and Technology, Irbid, 22110, Jordan

    Firas A. Al-Qarqaz & Hadeel S. Al-Kofahi

Authors
  1. Laith N. AL-Eitan
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Amneh H. Tarkhan
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Mansour A. Alghamdi
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Firas A. Al-Qarqaz
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Hadeel S. Al-Kofahi
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

LNA-E designed and supervised the study. LNA-E, AHT, FAA-Q and HSA collected the biological samples and clinical data. LNA-E, AHT and MAA lead the implementation of the method and performed the data analysis. LNA-E, AHT and MAA helped with the interpretation, description of the results and drafted the manuscript. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Laith N. AL-Eitan.

Ethics declarations

Ethics approval and consent to participate

Ethical approval was obtained from the IRB committee at Jordan University of Science and Technology (Ref. # 19/105/2017). All participants gave written informed consent before taking part in this study.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

Additional file 1 Figure S1.

Distribution of gene expression values. Box plots show the distribution of log counts (A) before and (B) after normalization. Figure S2. Overall heatmap of the top 500 most DE genes. The hierarchically clustered genes are represented by rows, and the samples are represented by columns, while the dendrograms and flat clusters are symbolized by the green and orange bars. Genes that have positively correlated logCPM values cluster together, as large positive correlations correspond to small distances. The red and blue colors indicate gene upregulation and downregulation, respectively.

Additional file 2 Table S1.

Top 100 DE genes in warts (FC > 2 and AP < 0.05) as sorted by the quasi-likelihood F test (QLF).

Additional file 3 Table S2.

Functional annotation clustering of the DE genes only showing DE genes with a p-value < 0.05. BP stands for biological processes, MF for molecular function, CC for cellular component, and GO for gene ontology.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

AL-Eitan, L.N., Tarkhan, A.H., Alghamdi, M.A. et al. Transcriptome analysis of HPV-induced warts and healthy skin in humans. BMC Med Genomics 13, 35 (2020). https://doi.org/10.1186/s12920-020-0700-7

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12920-020-0700-7

Keywords

BMC Medical Genomics

ISSN: 1755-8794

Contact us