Gene Expression Profiling of MicroRNAs in HPV-Induced Warts and Normal Skin

Infection with the human papillomavirus (HPV) is a common occurrence among the global population, with millions of new cases emerging on an annual basis. Dysregulated microRNA (miRNA) expression is increasingly being identified to play a role in a number of different diseases, especially in the context of high-risk HPV infection. The present study investigated the miRNA expression profiles of warts induced by low-risk HPV. In warts, miR-27b, miR-24-1, miR-3654, miR-647, and miR-1914 were downregulated while miR-612 was upregulated compared to normal skin. Using miRTargetLink Human, experimentally supported evidence was obtained showing that miR-27b targeted the vascular endothelial growth factor C (VEGFC) and CAMP-responsive element binding protein 1 (CREB1) genes. The VEGFC and CREB1 genes have been reported to be involved in tumorigenesis and wart formation, respectively. Similarly, the oxidized low-density lipoprotein receptor 1 (OLR1) gene, which plays an important role in the humoral immunity of the skin, and the plexin D1 (PLXND1) gene, which is highly expressed in tumor vasculature, were both found to be common targets of miR-27b, miR-1914, and miR-612.


Introduction
Human skin is an integumentary organ that consists of two main layers, the outermost epidermis and the underlying dermis, connected by the basement membrane [1]. The multilayered epidermis is predominantly made up of keratinocyte cells that originate in its basal layer and function to create a physical barrier against environmental and pathogenic threats [2]. Furthermore, keratinocytes faced with invasion by pathogens can induce apoptosis, trigger inflammatory responses, and produce various types of antimicrobial peptides as a means of innate immune defense [3][4][5]. One pathogen that exclusively targets keratinocytes found in the basal epidermal layer is the human papillomavirus (HPV) [6].
As the most common sexually transmitted infection, HPV infects millions of people each year and can be spread by skin-to-skin contact [7]. However, the majority of individuals infected with HPV do not exhibit any clinical symptoms due to quick resolution by the host's immune system [8,9].

RNA Sequencing (RNA-seq)
Samples that met standards for quality control were shipped on dry ice to the Australian Genome Research Facility in Melbourne, Australia, where RNA-seq was performed on the Illumina HiSeq 2500 according to the manufacturer's protocol. Real-time image analysis was carried out using the HiSeq Control Software (HCS) v2.2.68 (Illumina, Inc., San Diego, CA, USA) and Real Time Analysis (RTA) v1.18.66.3 (Illumina, Inc., San Diego, CA, USA). Afterwards, sequence data was generated through the Illumina bcl2fastq 2.20.0.422 pipeline (Illumina, Inc., San Diego, CA, USA).

Identification of Differentially Expressed miRNAs
Differentially expressed (DE) miRNAs were identified using the normalized sequence reads from all 24 samples. An miRNA was determined to be DE between warts and normal tissue samples if it had an adjusted p-value (AP) of less than or equal to 0.05.

Differentially Expressed miRNAs
Among the 1830 miRNAs included in the analysis, a total of 6 differentially expressed (DE) miRNAs were found in warts. Of these, 5 miRNAs were found to be upregulated and 1 miRNA was downregulated in warts compared to normal skin (Table 1). The fold change for each of the 12 control and 12 wart samples was individually plotted for the 6 DE miRNAs, and an overall consistent DE of the latter was shown in the wart samples ( Figure S1).
Interaction networks for each candidate miRNA and overlapping target genes were created using miRTargetLink Human. Interactions with strong experimental evidence are shown in Figure 1. The term "strong evidence" indicates that the interaction of miRNAs with the target genes is supported by strong experimental methods from the literature, such as reporter gene assays. By contrast, "weak evidence" indicates that the interactions are supported by weaker experimental methods like microarrays. Information about the target genes was extracted from miRTarBase as well as from in-house generated data from miRTargetLink Human databases [30].
The target gene overlap between the 6 candidate miRNAs is shown in Figure 2 and Table 2. Three miRNAs, namely hsa-miR-27b, hsa-miR-1914-3p and hsa-miR-612, are at the center of the interaction network, having more than 12 interactions each. These three miRNAs have two common target genes (OLR1 and PLXND1). Moreover, the results showed that two genes (VEGFC and CREB1) are likely to be potential targets of hsa-miR-27b, a finding that is supported by previous reports. The miRNet webtool was used to predict diseases associated with the candidate miRNAs, and it revealed a total of 33 associated diseases for hsa-miR-24-1, 29 diseases for hsa-mir-27b, and 2 diseases for hsa-mir-612. By contrast, no associated diseases were predicted for hsa-mir-647, hsa-mir-3654, and hsa-mir-1914 (Table S1). 58 genes for hsa-mir-647, 34 genes for hsa-mir-3654, 88 genes for hsa-mir-1914-5p, 141 genes for hsamir-1914-3p, 37 genes for hsa-mir-24-1, and 422 genes for hsa-mir-27b. Further, potential target genes were predicted using TargetScan, it revealed a total of 31 genes for hsa-mir-612, 4469 genes for hsamir-647, 3000 genes for hsa-mir-3654, 3946 genes for hsa-mir-1914, 1392 genes for hsa-mir-24-1, and 2737 genes for hsa-mir-27b. Interaction networks for each candidate miRNA and overlapping target genes were created using miRTargetLink Human. Interactions with strong experimental evidence are shown in Figure 1. The term "strong evidence" indicates that the interaction of miRNAs with the target genes is supported by strong experimental methods from the literature, such as reporter gene assays. By contrast, "weak evidence" indicates that the interactions are supported by weaker experimental methods like microarrays. Information about the target genes was extracted from miRTarBase as well as from in-house generated data from miRTargetLink Human databases [30].
The target gene overlap between the 6 candidate miRNAs is shown in Figure 2 and Table 2. Three miRNAs, namely hsa-miR-27b, hsa-miR-1914-3p and hsa-miR-612, are at the center of the interaction network, having more than 12 interactions each. These three miRNAs have two common target genes (OLR1 and PLXND1). Moreover, the results showed that two genes (VEGFC and CREB1) are likely to be potential targets of hsa-miR-27b, a finding that is supported by previous reports. The miRNet webtool was used to predict diseases associated with the candidate miRNAs, and it revealed a total of 33 associated diseases for hsa-miR-24-1, 29 diseases for hsa-mir-27b, and 2 diseases for hsamir-612. By contrast, no associated diseases were predicted for hsa-mir-647, hsa-mir-3654, and hsamir-1914 (Table S1).

Functional Enrichment Analysis
GO enrichment analyses for miRNA target genes were performed using the DAVID online database webtool. The top 10 most significant GO terms of each criteria are presented in Figure 3, which shows that the target genes of the candidate miRNAs were mainly enriched in protein binding and poly(A) RNA binding on the MF level, enriched in nucleoplasm and nucleus on the CC level, and enriched in regulation of transcription, DNA-templated and transcription, DNA-templated on the BP level. The top 20 most significant KEGG pathway terms are presented. The miRNA target genes were mainly enriched in focal adhesion, ErbB signaling pathway, and adherens junction.

Functional Enrichment Analysis
GO enrichment analyses for miRNA target genes were performed using the DAVID online database webtool. The top 10 most significant GO terms of each criteria are presented in Figure 3, which shows that the target genes of the candidate miRNAs were mainly enriched in protein binding and poly(A) RNA binding on the MF level, enriched in nucleoplasm and nucleus on the CC level, and enriched in regulation of transcription, DNA-templated and transcription, DNA-templated on the BP level. The top 20 most significant KEGG pathway terms are presented. The miRNA target genes were mainly enriched in focal adhesion, ErbB signaling pathway, and adherens junction.

Discussion
In the present study, miR-27b was found to be upregulated in warts compared to normal skin, and this miRNA displayed the highest number of interactions (n = 27) with target genes. miR-27b, one of two miR-27 homologs, modulates adipocyte differentiation as well as adipogenesis regulation by targeting peroxisome proliferator-activated receptor gamma (PPARγ) [31]. High-risk HPV infection is reported to upregulate miR-27b levels in HPV-positive cervical cancer tissues which, in turn, increases cell proliferation and decreases apoptosis by inhibiting PPARγ expression [32][33][34]. By contrast, the inhibition of PPARγ expression was found to enhance the response of cervical cancer cells to radiation treatment, while ligand activation of the PPARγ nuclear receptor resulted in the induction of differentiation and apoptosis in non-small cell lung cancer cells [35,36].
In addition, strong evidence illustrated that miR-27b-3p targets the vascular endothelial growth factor C (VEGFC) and CAMP-responsive element binding protein 1 (CREB1) genes ( Table 2). The VEGFC protein functions mainly in the process of lymphangiogenesis but has also been implicated in the promotion of tumor metastasis and growth [37]. VEGFC gene transcription was significantly increased in allergen-stimulated keratinocytes, but the VEGFC protein itself was not reported to directly affect the in vitro proliferation of epidermal keratinocytes [38,39]. In the context of high-risk HPV infection, VEGFC expression was stimulated by cigarette smoke and associated with the grade of cervical intraepithelial neoplasia [40,41]. On the other hand, the inhibition of CREB family members led to the reduction of papilloma formation in the murine epidermis via induction of apoptosis [42]. Correspondingly, CREB1 knockdown was found to promote apoptosis in murine follicular cells [43]. Moreover, in bladder cancer cells, CREB1 has been implicated in epithelial to mesenchymal transition, a carcinogenic process which is non-significantly associated with the HPV status of squamous cell oropharyngeal carcinoma [44,45]. miR-24-1 was found to be upregulated in warts compared to normal skin. miR-24 was reported to induce apoptosis and inhibit growth in laryngeal squamous cell carcinoma cells [46]. Moreover, keratinocyte differentiation normally promotes miR-24 expression, but this induction does not occur in keratinocytes infected with high-risk HPV [47]. Moreover, aberrant miR-24 expression was reported in HPV-positive cervical cancer cell lines [48]. Elevated miR-24 levels were associated with secreted frizzled-related protein 4 (SFRP4) levels in obese and diabetic subjects [49]. miR-24-1 was strongly upregulated during white adipocyte differentiation in a murine model, and its overexpression helped enhance adipocyte differentiation in murine mesenchymal stem cells [50,51].
Compared to normal skin, miR-3654, miR-647, and miR-1914 were found to be upregulated in warts. Both miR-647 and miR-1914 were found to be elevated in colorectal cancer tissues and cell lines as miR-647/1914 was found to promote the proliferation and migration of cancer cells by the downregulation of nuclear factor I/X [52]. By contrast, another study found that miR-1914 helped suppress the chemoresistant abilities of colorectal cancer cells also by nuclear factor I/X (NFIX) downregulation [53]. Moreover, miR-3654 has been identified for its involvement in prostate cancer progression [54].
Contrastingly, miR-612 was found to be downregulated in warts compared to normal skin. miR-612 plays a critical role in the development of colorectal cancer cells, and it was also found to have a suppressive role in reducing the stemness and tumor metastasis of hepatocellular carcinoma cells [55,56]. In hepatocellular carcinoma patients, miR-612 inhibits the metastasis and epithelial-mesenchymal transition of cancer cells by targeting the AKT2 gene [57]. In addition, this miRNA has been reported as a potential regulator of TP53 and CD40 expression in adults with metabolic syndrome [58]. As can be seen, it appears that miR-612 possess antitumor functions in malignant lesions, which might explain why it was found to be downregulated in benign lesions such as warts.
With regard to interactions with other genes, miR-27b, miR-1914-3p, and miR-612 were found to have the highest number of interactions. The oxidized low-density lipoprotein receptor 1 (OLR1) and plexin D1 (PLXND1) genes were identified as common targets of the three miRNAs. The ORL1 gene encodes for the lectin-type oxidized LDL receptor 1 (LOX-1), a protein that acts as the main oxidized LDL receptor on several different types of cells [59,60]. In the skin, the LOX-1 proteins on the surface of keratinocytes can be acted upon by nociception to produce leukotriene B(4), the latter of which induces an itch-associated response in a murine model [61]. LOX-1 also plays a role in the humoral immunity of the skin, as it is expressed on the surfaces of dermal dendritic cells [62,63]. On the other hand, the PLXND1 gene is involved in body fat distribution and is abundantly expressed in tumor vasculature [64][65][66].
Despite the ubiquity of low-risk HPV infection, there is a dearth of information available regarding the expression profiles of non-genital warts. While miRNA profiles have been investigated in HPV-associated cancers, little is known about the changes in miRNA expression that occur in warts. As a result, the present study aimed to shed light on this issue and provide a much-needed comparison between the miRNA expression patterns of high-and low-risk HPV infection. One limitation of the present study was that all participants were male. However, as this is an exploratory study, only male subjects were included in order to reduce any gender-specific variation that might arise from sex-biased expression patterns.

Conflicts of Interest:
The authors declare no competing interests.