Lack of Conserved miRNA Deregulation in HPV-Induced Squamous Cell Carcinomas

Squamous cell carcinomas (SCCs) in the anogenital and head and neck regions are associated with high-risk types of human papillomaviruses (HR-HPV). Deregulation of miRNA expression is an important contributor to carcinogenesis. This study aimed to pinpoint commonly and uniquely deregulated miRNAs in cervical, anal, vulvar, and tonsillar tumors of viral or non-viral etiology, searching for a common set of deregulated miRNAs linked to HPV-induced carcinogenesis. RNA was extracted from tumors and nonmalignant tissues from the same locations. The miRNA expression level was determined by next-generation sequencing. Differential expression of miRNAs was calculated, and the patterns of miRNA deregulation were compared between tumors. The total of deregulated miRNAs varied between tumors of different locations by two orders of magnitude, ranging from 1 to 282. The deregulated miRNA pool was largely tumor-specific. In tumors of the same location, a low proportion of miRNAs were exclusively deregulated and no deregulated miRNA was shared by all four types of HPV-positive tumors. The most significant overlap of deregulated miRNAs was found between tumors which differed in location and HPV status (HPV-positive cervical tumors vs. HPV-negative vulvar tumors). Our results imply that HPV infection does not elicit a conserved miRNA deregulation in SCCs.


Introduction
Human papillomaviruses (HPV), species-specific viruses infecting mucosal and cutaneous epithelium of the host, can cause benign diseases such as papillomata or warts; however, some HPV types have oncogenic potential and are reported to be responsible for up to 5% of all tumors worldwide and up to 30% of those whose development is attributed to infectious agents [1]. The most common carcinoma associated with HPV infection is cervical cancer, which is caused by HPV in almost 100% of cases. HPV is also associated with other tumors of the anogenital region such as vaginal tumors, where HPV participates in 78% of cases, vulvar tumors with 25% of HPV-positive cases, anal tumors with almost 90% of HPV-positive cases, and penile tumors with 50% of HPV-positive cases [2]. The second region where HPVs play the etiological role in cancer development is head and neck, especially the oropharyngeal location [3]. The proportion of HPV-associated oropharyngeal cancers is increasing around the world, reaching 40-80% of cases in the USA and 20-90% of cases throughout Europe [4][5][6].
The proper molecular profiling of HPV-positive and HPV-negative tumors and comparisons of differences between the two etiologies is necessary for expanding our understanding of the mechanisms of carcinogenesis. There are only few publications addressing the differences in molecular profiles between HPV-associated and HPV-independent tumors. Koncar et al. determined the expression of 24 proteins by immunohistochemistry, mutations of 48 genes by sequencing, and copy number alterations for six genes by in situ hybridization for cohorts of vulvar, anal, cervical, and oropharyngeal tumors, reporting the HPV-positive tumors to have similar molecular profiles [7]. However, the HPV positivity in their samples was determined based on the p53 status, which is not an established surrogate marker for the detection of HPV. Nevertheless, Tuna and Amos compared The Cancer Genome Atlas (TCGA) database data from studies analyzing HPV-driven tumors; they found similar genomic alterations in HPV-negative and HPV-positive cases but also detected distinct epigenomic and transcriptomic profiles between these groups [8]. Genomic comparison of HPV-positive and HPV-negative oral squamous cell carcinomas (SCCs) was performed by Gillison et al. and they specified a number of unique genetic features for virus-associated tumors [9]. Similarly, significant genetic differences between HPV-positive and HPV-negative tumors were revealed in several other studies [10,11]. So far, only two studies using whole-exome sequencing of vulvar tumors with regard to HPV positivity have been published [12,13]. They have shown that genetic alterations in vulvar carcinomas encompass mutations and copy number alterations that differ between HPV-positive and HPV-negative cases, in addition to common alterations observed irrespective of the HPV status.
Promising specific, sensitive, and clinically significant biomarkers are miRNAs, small single-stranded non-coding RNAs that regulate gene expression and play an important role in cell development, growth, and differentiation. Their expression has been reported to be deregulated in tumors. MicroRNAs play a role in stimulating tumor growth by negative regulation of tumor suppressor genes or by positive regulation of oncogenes. Specific miRNA profiles have been reported in many studies of various tumors, for example, of leukemia, breast cancer, lung cancer, or prostate cancer [14][15][16][17][18]. Numerous studies focused on the analysis of miRNA profiles in cervical cancer [19][20][21][22]; however, miRNA expression in tumors of the rest of anogenital regions is less researched, with only two studies analyzing miRNA profiles in vulvar cancer [23,24]. MicroRNA expression in head and neck cancers has been increasingly studied, and several studies on cell lines or tumors have been published [25][26][27][28][29]; however, they did not examine their association with HPV status. HPV status was addressed in two studies of Lajer et al. [19,30] who defined the group of core HPV miRNAs specific for two types of HPV-positive tumors, cervical and head and neck. Furthermore, Miller et al. analyzed miRNA expression in HPV-positive and HPV-negative oropharyngeal carcinomas and validated the data against The Cancer Genome Atlas, identifying HPV-associated oncogenic miRNAs [31]. In our previous study, we focused on miRNA profiles in HPV-negative and HPV-positive tonsillar tumors and in cervical tumors and on the comparison with a model system of keratinocyte clones [32] and defined the core HPV miRNAs, but they did not overlap with those from the study by Lajer et al. The lack of comparability between the outcomes of these studies can be most likely attributed to the anatomical heterogeneity of the analyzed head and neck tumors and to the different methodological approaches.
In our study, we analyzed the miRNA expression profiles in a set of anogenital tumors associated with HPV, namely, cervical, vulvar, and anal tumors, and in a set of tonsillar tumors by next-generation sequencing (NGS) in an attempt to identify commonly deregulated miRNAs exclusive for HPV-dependent SCCs. To our knowledge, this is the first study to analyze miRNA profiles in four HPV-driven types of tumors which employs a unified methodological approach. Furthermore, this study provides miRNA expression analysis of insufficiently researched anal and vulvar carcinomas.

Clinical Samples
Samples of SCCs were obtained from the Department of Pathology and Molecular Medicine, 2nd Faculty of Medicine, Charles University and Motol University Hospital, Prague, and from the Institute of Pathology, 3rd Faculty of Medicine, Charles University and University Hospital Kralovske Vinohrady, Prague. All tissue samples originated from unique patients. Formalin-fixed paraffin-embedded (FFPE) samples prepared from tumor tissues were stained and macrodissected [32]. The HPV status of FFPE samples was determined by P16 immunochemistry and HPV DNA detection [33]. All HPV-positive tumors contained high-risk HPV16. For each anatomical location (tonsillar, cervical, vulvar, anal), the sample set consisted of three HPV-positive tumors, three HPV-negative tumors, and three samples of normal tissue. Cervical HPV-negative tumors were not included due to their unavailability.

Sample Processing for NGS
Total RNA was isolated from macrodissected sections of FFPE samples using AllPrep ® DNA/RNA FFPE (Qiagen, Germantown, MD, USA). RNA concentration and quality were measured by a Nanodrop ™ spectrophotometer (Thermo Scientific, Waltham, MA, USA) and Experion chip electrophoresis (Bio-Rad, Hercules, CA, USA). From purified RNA, sequencing libraries were prepared using the TruSeq Small RNA Library prep kit (Illumina, San Diego, CA, USA). The libraries were sequenced on NextSeq 500 (Illumina, San Diego, CA, USA) at the Division of Molecular Medicine, Ruder Bošković Institute (Zagreb, Croatia), or the EMBL's Genomics Core Facility (GeneCore) (Heidelberg, Germany) using the NextSeq 500/550 High Output Kit (Illumina, San Diego, CA, USA). Sequencing reads were deposited in the SRA, GenBank (BioProject ID: PRJNA718204).

Analysis of NGS Data
Sequencing reads were trimmed off adapter sequences using the FastQ toolkit Basespace App (Illumina). The trimmed reads were mapped onto human reference genome assembly GRCh38.p13 using the Geneious mapper [34] (max. mismatches: 10%, max. gaps: 10%). The complete analysis of DE between replicate tumor samples and normal tissues (expression normalization, DE calculation, calculation of statistical significance) was carried out using the DESeq2 package [35]. The miRNAs exhibiting fold change of expression (FC) ≥ 2.0 and statistical support of padj (p-value adjusted for multiple testing) ≤ 0.1 were extracted and further analyzed. The Adonis-Bray test [36] with 10,000 permutations was used to test the association of global miRNA expression with the sample metadata categories (tissue type, anatomic location, HPV status).

Global Characterization of miRNA Expression
Since deregulation of the cellular miRNA network is commonly observed in carcinogenesis, we aimed to characterize the patterns of miRNA expression in a set of tissue samples encompassing tumors from different anatomical locations and differing in HPV status. Samples of HPV-associated tumors (anogenital and tonsillar tumor tissues), of the corresponding HPV-negative tumors (if available), and of healthy tissues were collected from three patients per sample type and processed for NGS. Sequencing yields ranged from 338,068 to 4,273,955 total reads per sample; the percentage of miRNA-derived reads varied from 6.3% to 82.9% (Table S1). The expression level values were calculated for individual miRNAs as their proportion among the total miRNA reads (Table S2).
The similarity of tissue samples based on their miRNA expression profiles is visualized in Figure 1 by principal component analysis (PCA) and a heatmap. In both plots, the clustering and thus the similarity of samples correlated with their anatomical locations. Tonsillar samples formed a separated cluster. Among the anogenital samples, anal samples clustered separately from genital samples. A notable separation was detected between the genital (i.e., cervical and vulvar) samples. All cervical tumors, all vulvar HPV-negative tumors, and one vulvar HPV-positive tumor were distinct from normal cervical and vulvar tissues and two vulvar HPV-positive tumors; the extent of dissimilarity of their miRNA expression profiles was readily discernable in both PCA and the heatmap (Figure 1). The Adonis-Bray test revealed that the miRNA expression was most significantly associated with anatomical location (anal/cervical/vulvar/tonsillar; p = 0.0001) and tissue type (tumor/nonmalignant; p = 0.001) of samples. This is in accordance with miRNA expression being long recognized to exhibit tissue-specific and tumor-specific differences [37][38][39][40]. The association with HPV status (positive/negative) was weaker and not significant (p = 0.074).

Differential Expression of miRNA in SCCs
We compared miRNA expression between tumor and normal tissues of the same locations, identifying miRNAs which exhibited differential expression (DE). This type of comparison is particularly suitable for pointing out the miRNAs whose deregulation is linked with carcinogenesis. Figure 2a shows the total numbers of deregulated miRNAs in the seven tumor types differing in anatomical location and HPV status. The total numbers of deregulated miRNAs varied greatly among tumor types. Cervical HPV-positive tumors exhibited the highest number of deregulated miRNAs, followed by vulvar HPV-negative tumors. This is in accordance with the marked separation of the miRNA expression profiles of these tumor samples from all other tissues as seen in Figure 1. On the opposite end of the spectrum of miRNA deregulation, HPV-positive vulvar tumors yielded only a single downregulated miRNA (see below). Again, this was mirrored in the clustering of miRNA expression profiles; 2/3 of vulvar HPV-positive tumor samples were intermingled between normal vulvar tissues ( Figure 1).
The cross-tumor comparison revealed that the largest proportion of deregulated miRNAs was specific to a single tumor type (149/239 and 94/210 upregulated and downregulated miRNAs, respectively; Figure 2b). The numbers of miRNAs which exhibited shared deregulation decreased with the increasing extent of common deregulation among the seven tumor types. Both upregulated and downregulated miRNAs followed this pattern ( Figure 2b). Downregulated miRNAs were more abundant than upregulated miRNAs among the most widely deregulated miRNAs. Specifically, 58 vs. 68, 28 vs. 28, 4 vs. 14, and 0 vs. 6 miRNAs (upregulated vs. downregulated) exhibited shared deregulation in two, three, four, and five tumor types, respectively (Figure 2b). No miRNA was deregulated in more than five tumor types. The miRNAs whose deregulation was detected in most tumor types of the set are likely to play a conserved role in tumorigenesis. A literature search confirmed that all six miRNAs downregulated in 5/7 tumor types were determined to function as tumor suppressors in SCCs (Table 1). Notably, the sole miRNA downregulated in HPV-positive vulvar cancers (MIR451A) was present among these most widely deregulated miRNAs.   Table S3. Venn diagrams were created using InteractiVenn [79].

Limited Role of HPV in miRNA Deregulation
Since infection with high-risk HPV16 is a key factor which delineates tumors between and within different anatomical locations, we investigated whether infection with HPV is functionally translated into a common set of deregulated miRNAs. Core HPV miRNAs were focused on in two previous studies [19,32]; however, both analyses were limited to two groups of HPV-positive tumors (cervical and head and neck SCCs). Figure 2c shows that the miRNAs which exhibited deregulation only in HPV-positive tumors were present in minute numbers. Specifically, a single miRNA (miR-146a-5p) was upregulated exclusively in HPV-positive anal and cervical tumors, and a single miRNA (miR-210-3p) was upregulated exclusively in HPV-positive anal and tonsillar tumors. Both miRNAs have been previously shown to be involved in tumorigenesis. The deregulation of miR-146a-5p was observed in many types of cancers (reviewed in [80]), including cervical tumors where it functions as oncomiR, stimulating proliferation of cells [81]. MiR-210-3p was characterized as a hypoxia-regulated miRNA involved in many biological processes and was shown to be overexpressed in many tumors including head and neck cancer [82]. Five miRNAs (miR-663a, miR-769-5p, miR-1307-5p, miR-3196, miR-4800-3p) were upregulated and three miRNAs (miR-139-5p, miR-142-5p, miR-574-3p) were downregulated exclusively in HPV-positive tonsillar and cervical tumors. Downregulation of miR-139-5p and miR-574-3p in HPV-positive tonsillar and cervical tumors is in agreement with the results reported by Lajer et al. [19]. The deregulation of both miRNAs was observed in many other cancer types; the role of miR-139-5p as a tumor suppressor was confirmed in cervical tumors and head and neck tumors [83,84]; however, the function of miR-574-3p was not studied in these types of tumors. No miRNA was deregulated in more than two types of exclusively HPV-positive tumors (Figure 2c, Table S3). Therefore, our results do not support the existence of universal "Core HPV" miRNAs whose deregulation is dependent on HPV infection.
In our set of seven SCC types, the highest similarity in miRNA deregulation was detected between cervical and HPV-negative vulvar tumors. These tumors exhibited the largest total numbers of deregulated miRNAs (Figure 2a). Of all the miRNAs deregulated in HPV-negative vulvar tumors, 69% and 76% (upregulated and downregulated miRNAs, respectively) were also deregulated in cervical tumors. Twenty-nine upregulated miRNAs and 42 downregulated miRNAs exhibited exclusive deregulation in both cervical and vulvar HPV-negative tumors (Figure 2c). For these miRNAs, we mined available literature for experimental evidence of their role in carcinogenesis of SCCs, i.e., of the effect on cell proliferation and invasiveness. While information was lacking for most upregulated miRNAs, 81% of downregulated miRNAs were experimentally characterized ( Table 2). Of these 34 miRNAs, 30 were determined to function as tumor suppressors (17-in cervical carcinomas, 13-in other types of SCCs); their downregulation observed here is thus in line with their carcinogenic role in SCCs. Since the widespread downregulation of tumor suppressor miRNAs in vulvar HPV-negative tumors is independent of HPV infection, it should be considered convergent with cervical tumors.
Unless diagnosed early, cervical cancers are notorious for their poor prognosis; the otherwise high mortality is kept in check only by strenuous screening and HPV vaccination. In our set of SCCs, cervical carcinomas exhibited the most extensive miRNA deregulation (Figure 2a). Among vulvar carcinomas, the prognosis (survival rate) of HPV-negative tumors is far worse than that of HPV-positive tumors [85,86]; due to additional differences other than HPV status, both types of vulvar carcinomas are regarded as separate rather than related cancer types [87]. With respect to the total extent of miRNA deregulation, vulvar SCCs with differing HPV status are fundamentally different (extensive vs. negligible deregulation; Figure 2a); our findings thus tentatively suggest that the pattern of miRNA deregulation might be linked to the severity of gynecological malignancies. The current understanding of molecular mechanisms of vulvar carcinogenesis is extremely limited [88]. The downregulation of multiple tumor suppressor miRNAs detected in this study suggests a possible mechanism contributing to a severe course of vulvar HPV-negative cancers. Interestingly, differences in the prognosis [33] and miRNA deregulation pattern [32] were also previously found between HPV-positive and HPV-negative tonsillar tumors, although their potential connection has not been investigated yet.  In conclusion, this study demonstrates that miRNA deregulation in various SCC types is largely tumor-specific and that the HPV status of SCCs is not a determinant of either the extent or the composition of the deregulated miRNA pool. Our findings may provide useful directions for further studies on the pathogenesis of individual types of SCCs with respect to their differences in anatomical location and HPV status.