Next Article in Journal
Thrombotic Events Are Unusual Toxicities of Chimeric Antigen Receptor T-Cell Therapies
Previous Article in Journal
Is Renal Ischemic Preconditioning an Alternative to Ameliorate the Short- and Long-Term Consequences of Acute Kidney Injury?
Previous Article in Special Issue
The Functional Role and Regulatory Mechanism of FTO m6A RNA Demethylase in Human Uterine Leiomyosarcoma
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Splicing Characterization and Isoform Switch Events in Human Keratinocytes Carrying Oncogenes from High-Risk HPV-16 and Low-Risk HPV-84

by
Maryam Nasiri-Aghdam
1,2,
Mariel Garcia-Chagollan
3,
Ana Laura Pereira-Suarez
4,
Adriana Aguilar-Lemarroy
1 and
Luis Felipe Jave-Suarez
1,*
1
División de Inmunología, Centro de Investigación Biomédica de Occidente, Instituto Mexicano del Seguro Social, Guadalajara 44340, Mexico
2
Departamento de Biología Molecular y Genómica, Centro Universitario de Ciencias de la Salud, Universidad de Guadalajara, Guadalajara 44340, Mexico
3
Instituto de Investigación en Ciencias Biomédicas, Centro Universitario de Ciencias de la Salud, Universidad de Guadalajara, Guadalajara 44340, Mexico
4
Department of Microbiology and Pathology, Centro Universitario de Ciencias de la Salud, Universidad de Guadalajara, Guadalajara 44340, Mexico
*
Author to whom correspondence should be addressed.
Int. J. Mol. Sci. 2023, 24(9), 8347; https://doi.org/10.3390/ijms24098347
Submission received: 28 March 2023 / Revised: 27 April 2023 / Accepted: 4 May 2023 / Published: 6 May 2023

Abstract

:
Infection of epithelial cells with high-risk HPV (HR-HPV) types, followed by expression of virus oncogenic proteins (E5, E6, and E7), leads to genomic imbalance, suppression of tumor inhibitors, and induction of oncogenes. Low-risk HPV (LR-HPV) may slow the rate at which cervical cancer spreads to an invasive stage since co-infection with LR-HPV is linked to a decreased risk of future invasive cancer than infection with HR-HPV alone. We then propose that cancer-progressing changes may be distinguished through identifying the functional differences between LR-HPV and HR-HPV. Lentiviral strategies were followed to establish HaCaT cells with constitutive expression of HPV oncogenes. RNAseq experiments were designed to analyze the transcriptome modulations caused by each of the E5, E6, and E7 oncogenes of HPV-16 and HPV-84 in HaCaT cells. We identified enhanced RNA degradation, spliceosome, and RNA polymerase pathways related to mRNA processing. ATTS (alternative transcription termination site) was discovered to be more prevalent in cells with HPV-16E5 than HPV-84E5. In HPV-16E6-infected cells, ATTS gain was significantly higher than ATTS loss. Cells with HPV-16E7 had more isoforms with intron retention (IR) than those with HPV-84E7. We identified switches in ADAM10, CLSPN, and RNPS1 that led to greater expression of the coding isoforms in HR-HPV. The results of this work highlight differences between LR-HPV and HR-HPV in mRNA processing. Moreover, crucial cervical cancer-related switch events were detected.

1. Introduction

Viral infection accounts for approximately 12–20% of all cancers [1,2]. Human papillomavirus (HPV) is one of the most common cancer-causing viruses and a significant etiological factor for many malignancies, including cervical cancer, anal cancer, vulvar and vaginal cancer, penile cancer, head and neck cancer, and skin cancer [3]. Currently, 229 different HPV types are identified and categorized as high-, low-, or unknown-risk, depending on their carcinogenic potential [4,5]. Infection of epithelial cells with high-risk HPV types (HR-HPV), followed by the expression of viral oncogenic proteins (E5, E6, and E7), leads to genomic imbalance, suppression of tumor inhibitors, and induction of oncogenes. This process eventually drives neoplastic progression, which can take years to conclude [6]. The HPV vaccine may not protect against all HR-HPV strains and may not benefit women who already have the virus [7]. Although low-risk HPV (LR-HPV) proteins can alter cell proliferation, apoptosis, and immortalization, their carcinogenic potential is lower than that of HR-HPV proteins [8]. LR-HPV infections cause benign cervical lesions, such as common genital warts, that rarely progress towards cancer [9]. LR-HPVs may slow the rate at which cervical cancer spreads to an invasive stage since co-infection with LR-HPVs decreases the risk of future invasive cancer compared to infections with HR-HPVs alone [10]. Thus, it is proposed that through identifying the functional differences between LR-HPVs and HR-HPVs, it may be possible to distinguish cancer-progressing alterations.
A crucial stage in the post-transcriptional control of gene expression is alternative splicing (AS) of mRNAs. This mechanism significantly modifies cell function and proteome diversity through regulating the expression of various isoforms at particular times and in specific cell types. In this sense, through altering the mRNA processing machinery, cancer cells can gain advantages [11]. Therefore, mRNA processing regulators are a new class of oncoproteins or tumor suppressors because they can control disease progression through altering the RNA isoforms associated with the primary cancer processes [12]. Cancer cells and each cancer type and subtype have unique splicing patterns, making AS a new hallmark of cancer [11,13]. With the advance of bioinformatic tools, it is now possible to research eight different AS mechanisms, including alternative transcription start sites (ATSS), alternative transcription termination sites (ATTS), alternative three-end acceptor sites (A3), alternative five-end donor sites (A5), exon skipping (ES), intron retention (IR), mutually exclusive exon (MEE), and multiple exon skipping (MES) [14]. Moreover, detecting isoform switch events and their predicted consequences for individual genes may help clarify splicing dysregulation during HPV infection.
Although there has been some research on the alternative splicing processes by HR-HPV types, modulation of AS through the proteins encoded with LR-HPVs receives far less attention. Splicing events generated by E5, E6, and E7 oncoproteins of HPV-84 and HPV-16 were examined for the first time to highlight the differences between HR- and LR-HPVs; these modulations were described in eight distinct splicing mechanisms and individual isoform switch events.

2. Results

2.1. Analysis of KEGG Pathways Showed High Enrichment of Pathways Involved in the mRNA Processing

To determine potential KEGG pathways and gene sets linked with HR- and LR-HPV types, the expression data of the HaCaT transduced cells with the oncogenes of HPV-16 and HPV-84 were subjected to GSEA analysis. Over-represented pathways with normalized gene set enrichment scores (NES) ≤ 1 are illustrated in Figure 1. Among them, mRNA processing pathways, including spliceosome, RNA degradation, and RNA polymerase, were repeatedly enriched between groups. Interestingly, HaCaT-16E5 cells have overexpression of genes involved with RNA degradation, whereas HaCaT-84E5 does not exhibited this pathway as an enriched one. In HaCaT cells transduced with HPV-16 and -84 oncogenes, enrichment of genes associated with spliceosome was observed when the E6 and E7 oncogenes were introduced. In addition, E6 and E7 from HPV-84 induced the expression of RNA polymerase genes, in contrast to E6/E7 from HPV-16, which induced the expression of genes related to the RNA degradation pathway. The gene set enrichment analysis results highlight the importance of mRNA processing pathways; therefore, these biological processes were further investigated as outlined in the following section.

2.2. Splicing Mechanisms Demonstrate Subtle Differences throughout the Genome

Genome-wide alternative splicing was investigated in the eight previously identified categories, and the results make it quite evident that not all alternative splicing events were equally used (Figure 2). In all comparisons, ATSS was the most prominent event, while MEE was the least frequent, occurring only in cells with the E5 and E7 oncogenes of HPV-84. In the cells carrying E5 of HPV-16, there were more significant ATTS gains than losses, with a difference of about five significant isoforms. However, the losses outweighed the gains in the A5. The most notable HPV-84 event was more significant ATSS loss than gain, which for HPV-16 was not seen. When comparing HPV-16 and HPV-84, HPV-84 had a greater loss in MES but fewer significant isoforms with ATTS, ES, and IR events than HPV-16. The most noticeable variation seen in cells with E6 was fewer isoforms with ATTS losses in HPV-16 than in HPV-84. Finally, the cells containing E7 from HPV-16 exhibited a higher number of isoforms with IR gains than those with the E7 from HPV-84, which shows a higher number of isoforms with MES losses.
In Figure 3, the uneven utilization within each comparison is shown through calculating the percentage of events that are gains for each AS type (as opposed to loss). It is evident from this plot that the overall trend in the usage of alternative splicing is the same in the two HPV types. However, there are some variances. We can, indeed, see that ATTS has a much more significant skew in HPV-16E6 compared to HPV-84E6.

2.3. Overall Predicted Switch Consequences Are Different

According to coding potential, identified domains, intron retention, sensitivity to non-sense mediated decay (NMD), and ORF length, isoform switch effects were predicted and characterized (Figure 4). The cells with HPV-16E5 oncogene exhibit increased isoforms with domain loss and sensitivity to NMD. Moreover, compared to HPV-84E5, HPV-16E5 showed a higher number of isoforms with shorter or longer ORF. HPV-16E6 induced more coding transcripts than HPV-84, where most switches are non-coding. One isoform with a domain switch was observed for HPV-16E6. Additionally, the number of isoforms with complete ORF gain and those with longer ORF are approximately higher in HPV-16E6 compared to HPV-84E6. Although the overall pattern in cells expressing E7 is identical for HPV-16 and HPV-84, HPV-84 has a more significant proportion of isoforms with domain loss, whereas HPV-16 has a higher proportion of NMD-insensitive isoforms.

2.4. Switch Overlaps and Increased Usage of Coding Isoforms Are Detected

The general overlaps in isoform switches for all three oncogenes from HPV-16 and HPV-84 are illustrated in Figure 5A–C. There was little overlap across all comparisons, indicating that each HPV type’s switch events are distinct. A minor overlap was seen between HPV-16E5 and HPV-84E5 for only one gene: ABHD14B (Figure 5D). The coding isoform ENST00000395008 of this gene had increased usage in the control group, and neither HPV-16E5 nor HPV-84E5 expressed this or any other coding isoforms of this gene. However, the overlap between HPV-16 and HPV-84 for E6 and E7 oncogenes approached three genes. Both HPV-16 and HPV-84, which are E6-transformed cells, had an elevated expression of at least one of the coding isoforms of AARS2, BZW2, and NUP205. A novel coding isoform of PCIF1 was upregulated in the control group versus HPV-16E6 and HPV-84E6. The DLST gene was the only gene where an isoform switch increased the isoform expression in both HPV-16E7 and HPV-84E7. When comparing the control group to HPV-16E7 and HPV-84E7, elevated expression of the other three genes—ABHD14B, RPUSD4, and TRAF7—was observed. Figure 5E listed all coding isoform switches when comparing HPV-16 to HPV-84 for all three oncogenes.

2.5. Isoform Switch Events in Specific Genes Are Notable

We discovered three events in the genes ADAM10, RNPS1, and CLSPN that were previously reported to play a role in the progression or prognosis of cervical cancer. The switch plots of these three genes are illustrated in Figure 6A–C. The first switch event occurred in ADAM10, where HPV-16E5 expressed the coding isoform ENST000000260408 more strongly than HPV-84E5. When comparing HPV-16E6 to HPV-84E6, two events were detected. The isoform ENST000000318121 was observed in CLSPN to be more abundant in HPV-16E6. In the third event, RNPS1, the expression of ENST000000566397 was higher in HPV-16E7.

3. Discussion

The maturation, processing, and alternative splicing of mRNAs are necessary for regulating gene expression and preserving proteome diversity. The development of cervical cancer is thought to be significantly influenced via HPV-mediated alternative splicing, and these alterations can be utilized as diagnostic biomarkers and targets for therapeutic strategies [15]. On the other hand, putative disease-causing splice isoforms might provide crucial details about cancer development and help develop therapeutic approaches [16]. Due to the importance of such alterations, we hypothesized that high- and low-risk HPV types might differ in splicing mechanisms and isoform switch events, enabling us to distinguish between non-cancerous and actual cancer driver events. We constructed six keratinocyte cellular models, one for each of the E5, E6, and E7 oncogenes derived from high-risk HPV-16 and low-risk HPV-84. The expression data of these models were initially used to perform gene set enrichment analysis, followed by extensive alternative splicing analysis and detection of individual and generic patterns of isoform switch events. Firstly, key pathways were identified via gene set enrichment analysis, with all being strongly related to mRNA processing and transcription regulation. Among these are RNA degradation, spliceosomes, and RNA polymerase. Except for cells containing HPV-84 E5, they were found to be enriched in all experimental groups. Several key splicing factors were found with altered expression in cervical cancer, primarily from SRSF and hnRNP protein families, which were reviewed in detail in several recently published articles [17,18,19].
RNA degradation aids ribonucleotide recycling, while also performing surveillance through eliminating aberrant RNA that might produce damaging proteins [20]. There are scarce studies on modifications of RNA degradation-related genes in cervical cancer. A 2020 study showed that deleting EDC4 (Enhancer of mRNA decapping protein 4) in cervical cancer cells enhanced sensitivity to cisplatin and inhibited cell proliferation produced via cisplatin and DNA damage [21]. This finding could imply that EDC4 is a new target to avoid chemotherapy resistance in cervical cancer. In a study on N6-methyl-adenosine (m6A) mediated via METTL3 (methyl-transferase-like 3), it was discovered that METTL3 promotes the stability of HK2 (hexokinase 2) through m6A modification, thereby promoting the Warburg effect (also known as aerobic glycolysis), which could lead to a new insight for the treatment of cervical cancer [22]. Concerning the RNA polymerase pathway, it was demonstrated that Pol III dysregulation stemming from oncogenes or tumor suppressors is present in many malignancies and contributes to carcinogenesis [23]. In summary, studies on mRNA processing pathways, specifically on RNA degradation and RNA polymerase in cervical cancer, are limited. Further research is needed to identify their function in the progression and treatment of cervical cancer.
The enrichment of mRNA processing pathways piqued our interest in further investigating mRNA processing activities in our cell models. Following the pipeline defined using IsoformSwitchAnalyzeR, we identified eight alternative splicing events in a number of significant isoforms. Notably, ATTS is more common in cells carrying HPV-16E5 than in cells harboring HPV-84E5. ATTS, on the other hand, is demonstrated to have the most notable alteration in cells with E6, with fewer isoforms with ATTS being detected in HPV-16 compared to HPV-84. Moreover, in HPV-16E6, ATTS gain is significantly higher than ATTS loss. Generally, using ATTS produces transcripts with different 3′ ends or even transcripts with different coding regions [24]. However, the stability of the transcripts generated from ATTS can vary depending on the transcript and the cellular context [25]. Alternative transcription start and termination sites were previously identified as the primary drivers of transcript isoform diversity across human tissues [26]. According to Kim et al., the distribution of ASE types differed between malignant and normal tissues, with cancer cells showing less exon skipping but more alternative start or end sites than normal cells [27]. ATTS dysregulation is a common occurrence in cancer, and it can result in increased oncogene expression and decreased tumor suppressor gene expression [28]. Cervical cancer was found to have dysregulation in two ATTS regulatory factors: cold-inducible RNA binding protein (CIRP, also known as CIRBP or A18 hnRNP), which binds and stabilizes pro-survival gene transcripts, and RBBP6, which suppresses polyadenylation [29]. Overall, ATTS seems to play an important role in cervical cancer, even distinguishing between high- and low-risk HPV types; however, further study is required to back up this finding.
A further finding in our study concerns E7-containing cells where HPV-16 has a more significant number of isoforms with IR than HPV-84. IR appears to be necessary at several phases of cell differentiation and development. Intron-retaining transcripts can be detained in the nucleus (a process known as intron detention) and destroyed via nuclear degradation mechanisms. As introns often include in-frame premature termination codons, when preserved in the mature mRNA transcript and exported to the cytoplasm, they can be identified via the cytoplasmic surveillance machinery and eliminated via nonsense-mediated decay (NMD). Furthermore, the intron retention mechanism may yield non-coding RNA(s) implicated in controlling oncogenes and tumor suppressor genes. [30]. Except for breast cancer, it has been reported that retention of alternative introns is enriched in 16 different cancer types; however, cervical cancer was not included in this study [31]. In our investigation, IR events distinguished between high- and low-risk E7; nonetheless, its role in cancer progression requires additional research.
In our study, we found a more significant number of isoforms with domain loss, sensitivity to NMD, and shorter or longer ORFs, primarily mediated by high-risk HPV16 oncogenes. It is now well understood that viral proteins influence the cellular splicing machinery, producing RNA isoforms with carcinogenic properties [17,18]. Furthermore, alternative splicing contributes to virally induced cancer progression through increasing the expression of proteins implicated in proliferation and immune response [32]. We focused on three isoform switch events in ADAM10, CLSPN, and RNPS1 that led to greater expression of the coding isoforms in high-risk HPV.
These genes have previously been linked to cervical cancer. The proteolytic cleavage of NKG2D ligands (NKG2DL) on the cancer cell surface via ADAM10 can impair the recognition of cancer cells by T or NK cells [33]. In cervical cancer, ADAM10 is the target of miR-140-5p, which is controlled by the SNHG20 lncRNA. The inhibition of SNHG20 can reduce ADAM10 protein expression, resulting in decreased cervical cancer cell proliferation [34]. The second gene—CLSPN—regulates the ATR/Chk1 signaling axis in the G2 DNA damage checkpoint [35]. In histological and cytological samples, claspin expression is substantially associated with HR-HPV infection and lesion grade, which could be therapeutically helpful in detecting HPV-related cervical lesions [36]. The third isoform switch occurred in RNPS1, which is an essential regulator of the splicing process previously found to be overexpressed in cervical cancer cells. Alternative splicing controlled by RNPS1 favors an active Rac1b/RhoA signaling axis, possibly contributing to cervical cancer cell invasion and metastasis [37]. Given the relevance of isoform switches in the control of coding and non-coding isoforms, it seems necessary to consider these switch events in future transcriptome analyses.

4. Materials and Methods

4.1. Obtaining Keratinocytes Carrying the E5, E6, and E7 Oncogenes of HPV-16 and HPV-84

Briefly, E5, E6, and E7 were amplified using specific primers from genomic DNA extracted from cervical biopsies of women infected with HPV16 or 84. The amplified ORFs were subcloned into the lentiviral expression vector pLVX-Puro and verified via sequencing. Lentiviral particles were produced using the acquired plasmids pLVX-16E5, pLVX-16E6, pLVX-16E7, pLVX-84E5, pLVX-84E6, and pLVX-84E7. In the following step, HaCaT cells were individually infected with the lentiviral particles containing each viral gene. Transduced cells were selected through puromycin addition for three weeks. Transduced HaCaT cells were named HAC16E5, HAC16E6, HAC16E7, HAC84E5, HAC84E6, HAC84E7, and HACPLVX [38]. HaCaT cells were initially kept frozen in liquid nitrogen until use. When the cultures were prepared, the cells were thawed at 37 °C, centrifuged at 1000× g rpm for 10min to remove excess DMSO (dimethyl sulfoxide), suspended in fresh medium, and cultured in 75 cm2 culture flasks. Cell lines were grown in Invitrogen’s GIBCO™ DMEM (Dulbecco’s Modified Eagle Medium) medium with L-glutamine (584 mg/L), sodium pyruvate (110 mg/L), and D-glucose (4.5 g/L). The medium was supplemented with penicillin (100 U/mL), streptomycin (100 μg/mL), and 10% fetal bovine serum. The cultures were maintained at a confluence of 70–80%; upon reaching this level, the cells were detached with 0.25% trypsin (EDTA 380 mg/L) and counted to make the necessary dilutions for the experiments. Culture passages were created according to the cell duplication rate. Cell cultures were maintained in a humidified incubator at 37 °C in a 5% CO2 atmosphere.

4.2. RNA Extraction and Sequencing

Total RNA was obtained from 5 × 106 cells of each cell line in three replicates using the PureLinkTM RNA Mini Kit extraction system (Ambion, Naugatuck, CT, USA). Cells were washed two times with PBS, and cell lysis was performed. The lysate obtained was column purified, and the RNA was eluted in water. The RNA concentration and integrity were evaluated with the Agilent 2100 bioanalyzer using RNA nano chips. RT-qPCR confirmed the expression of E5, E6, and E7 through specific primers. Sequencing of total mRNA was performed using the Illumina Nova-Seq 6000 platform with a 150 bp paired-end method, where most reads have a Phred quality score greater than 35.

4.3. Differential Expression and Gene Set Enrichment

Data quality was verified using the FastQC v.0.11.9 tool [39] on the Galaxy platform available at https://usegalaxy.org/, accessed on 5 April 2021 [40]. Reads were mapped to the hg38v38 human genome using STAR v.2.7.9a [18] tool. Mapped reads were filtered based on the mapping quality utilizing the Qualimap v.2.2.1 tool [41]. The overall mapping rate was about 92%. The STAR-aligned reads were used separately as DESeq2 v.1.36.0 [42] input files for differential expression calculation between our experimental groups. After data pre-processing, the normalized counts from DESeq2 were uploaded to the GSEA v. 4.3.2 desktop application (https://www.gsea-msigdb.org/, accessed on 21 December 2021 [43,44], and enrichment analyses were performed for the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways [45]. The number of permutations was set to 1000, and the permutation type was selected as “gene_set”. Pathways with a normalized enrichment score >1 or <−1 at FDR < 10% were considered significant.

4.4. Transcript Assembly

For accurate reconstruction of all expressed isoforms of each gene, and to estimate the relative abundance of those isoforms, the StringTie v2.2.0 transcript assembler was used [46]. For this aim, the STAR alignment file of all reads from each sample in BAM format was taken as the input for StringTie to assemble transcripts. All the genes found in any samples were merged to create a single assembly. After merging, a second run of StringTie was performed to recalculate the abundance of the merged transcripts.

4.5. Isoform Switch Events and Alternative Splicing Analysis

IsoformSwitchAnalyzeR was used to detect and visualize the mechanisms of alternative splicing and individual isoform switch events [14]. In the pipeline, the transcript assembly files for each sample, the merged-GTF file from StringTie, and the reference GTF file (Homo_sapiens.GRCh38.102.gtf) were used to construct the isoformSwitchAnalysisPart1 object [47]. For the filtering step, a gene expression cutoff of 3 and an isoform expression cutoff of 1 were applied. The Isoform switch test was implemented using DEXSeq via IsoformSwitchAnalyzeR [48]. To predict the consequence of each isoform switch event in coding, non-coding, and nonsense-mediated decay (NMD) sensitive categories, data from external annotation sources, including Pfam for protein domains [49], SignalP for signal peptides [50], CPAT for potential coding sequences [51], and IUPred2A for prediction of intrinsically disordered regions [52], were then imported and integrated to the IsoformSwitchAnalyzeR object. To measure the effect size in isoform usage, IsoformSwitchAnalyzeR uses isoform fraction (IF) values, which quantify the fraction of the parent gene expression originating from a specific isoform (calculated as isoform_exp/gene_exp). Finally, the isoform(s) used frequently (positive dIF) is compared to the isoform(s) used infrequently (negative dIF) to identify potential functional consequences of the isoform switch. Only switch events in the coding isoform with a significant IF were considered when filtering out the isoform switch events. Additionally, all switch events were verified twice in the ENSEMBL database, and those in coding isoforms that translated to the same number of amino acids were eliminated. The events that did not affect the protein domains, signal peptides, or intrinsically disordered areas were also removed from the final table of coding isoforms. The following steps involved performing an alternative splicing analysis for the following eight categories: alternative 3′ acceptor sites (A3), alternative 5′ splice sites (A5), alternative transcription start sites (ATSS), alternative transcription termination sites (ATTS), exon skipping (ES), intron retention (IR), mutually exclusive exons (MEE), and multi-exon skipping (MES), both at genome-wide and isoform levels.

5. Conclusions

Our findings point to key differences produced by the oncogenes E5, E6, and E7 of HPV-16 and HPV-84. We focused primarily on mRNA processing pathways, where we discovered specific enhanced splicing processes and isoform switch events in HPV-16 that may favor cancer progression. These changes, in addition to providing information on biological processes implicated in cervical cancer pathways, may be beneficial for identifying prospective biomarkers and therapeutic targets. Nonetheless, only large prospective studies may determine the actual clinical and diagnostic utility of these alterations in cervical cancer.

Author Contributions

Conceptualization, M.N.-A., A.A.-L. and L.F.J.-S.; methodology, M.N.-A., M.G.-C., A.L.P.-S., A.A.-L. and L.F.J.-S.; investigation, M.N.-A., M.G.-C., A.L.P.-S., A.A.-L. and L.F.J.-S.; software, M.N.-A.; formal analysis, M.N.-A. and L.F.J.-S.; resources, M.G.-C., A.L.P.-S. and L.F.J.-S.; data curation, M.N.-A.; writing—review and editing, M.N.-A. and L.F.J.-S.; visualization, M.N.-A.; supervision, L.F.J.-S.; funding acquisition, A.A.-L. and L.F.J.-S. All authors have read and agreed to the published version of the manuscript.

Funding

This research was granted by IMSS.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

All raw and processed data for this project was deposited in the Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo, accessed on 30 August 2022) under accession number GSE228187. Scripts for differential expression and splicing analysis were released in the following GitHub repository: https://github.com/MaryamNasiriAghdam/MNasiri_HPV, accessed on 14 February 2022.

Acknowledgments

M.N.A. expresses her gratitude to the Mexican National Council for Technology (CONACYT) for financing her postgraduate study.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analysis, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

References

  1. Akram, N.; Imran, M.; Noreen, M.; Ahmed, F.; Atif, M.; Fatima, Z.; Waqar, A.B. Oncogenic Role of Tumor Viruses in Humans. Viral Immunol. 2017, 30, 20–27. [Google Scholar] [CrossRef] [PubMed]
  2. Mui, U.N.; Haley, C.T.; Tyring, S.K. Viral Oncology: Molecular Biology and Pathogenesis. J. Clin. Med. 2017, 6, 111. [Google Scholar] [CrossRef] [PubMed]
  3. Francies, F.Z.; Dlamini, Z. Aberrant Splicing Events and Epigenetics in Viral Oncogenomics: Current Therapeutic Strategies. Cells 2021, 10, 239. [Google Scholar] [CrossRef]
  4. Karolinska Institute. International Human Papillomavirus Reference Center. Available online: https://www.hpvcenter.se/human_reference_clones/ (accessed on 1 September 2020).
  5. de Villiers, E.-M.; Fauquet, C.; Broker, T.R.; Bernard, H.-U.; zur Hausen, H. Classification of papillomaviruses. Virology 2004, 324, 17–27. [Google Scholar] [CrossRef]
  6. Olusola, P.; Banerjee, H.N.; Philley, J.V.; Dasgupta, S. Human Papilloma Virus-Associated Cervical Cancer and Health Disparities. Cells 2019, 8, 622. [Google Scholar] [CrossRef] [PubMed]
  7. Liao, J.B. Viruses and human cancer. Yale J. Biol. Med. 2006, 79, 115–122. [Google Scholar]
  8. Aranda-Rivera, A.K.; Cruz-Gregorio, A.; Briones-Herrera, A.; Pedraza-Chaverri, J. Regulation of autophagy by high- and low-risk human papillomaviruses. Rev. Med. Virol. 2021, 31, e2169. [Google Scholar] [CrossRef]
  9. Brentjens, M.H.; Yeung-Yue, K.A.; Lee, P.C.; Tyring, S.K. Human papillomavirus: A review. Dermatol. Clin. 2002, 20, 315–331. [Google Scholar] [CrossRef]
  10. Sundström, K.; Ploner, A.; Arnheim-Dahlström, L.; Eloranta, S.; Palmgren, J.; Adami, H.-O.; Helm, N.Y.; Sparén, P.; Dillner, J. Interactions Between High- and Low-Risk HPV Types Reduce the Risk of Squamous Cervical Cancer. J. Natl. Cancer Inst. 2015, 107, djv185. [Google Scholar] [CrossRef]
  11. Biamonti, G.; Catillo, M.; Pignataro, D.; Montecucco, A.; Ghigna, C. The alternative splicing side of cancer. Semin. Cell Dev. Biol. 2014, 32, 30–36. [Google Scholar] [CrossRef]
  12. Urbanski, L.M.; Leclair, N.; Anczuków, O. Alternative-splicing defects in cancer: Splicing regulators and their downstream targets, guiding the way to novel cancer therapeutics. Wiley Interdiscip. Rev. RNA 2018, 9, e1476. [Google Scholar] [CrossRef] [PubMed]
  13. Ladomery, M. Aberrant alternative splicing is another hallmark of cancer. Int. J. Cell Biol. 2013, 2013, 463786. [Google Scholar] [CrossRef] [PubMed]
  14. Vitting-Seerup, K.; Sandelin, A. IsoformSwitchAnalyzeR: Analysis of changes in genome-wide patterns of alternative splicing and its functional consequences. Bioinformatics 2019, 35, 4469–4471. [Google Scholar] [CrossRef] [PubMed]
  15. Francies, F.Z.; Bassa, S.; Chatziioannou, A.; Kaufmann, A.M.; Dlamini, Z. Splicing Genomics Events in Cervical Cancer: Insights for Phenotypic Stratification and Biomarker Potency. Genes 2021, 12, 130. [Google Scholar] [CrossRef] [PubMed]
  16. Bergsma, A.J.; van der Wal, E.; Broeders, M.; van der Ploeg, A.T.; Pijnappel, W.P. Alternative Splicing in Genetic Diseases: Improved Diagnosis and Novel Treatment Options. Int. Rev. Cell Mol. Biol. 2018, 335, 85–141. [Google Scholar] [CrossRef] [PubMed]
  17. Cerasuolo, A.; Buonaguro, L.; Buonaguro, F.M.; Tornesello, M.L. The Role of RNA Splicing Factors in Cancer: Regulation of Viral and Human Gene Expression in Human Papillomavirus-Related Cervical Cancer. Front. Cell Dev. Biol. 2020, 8, 474. [Google Scholar] [CrossRef]
  18. Boudreault, S.; Roy, P.; Lemay, G.; Bisaillon, M. Viral modulation of cellular RNA alternative splicing: A new key player in virus-host interactions? Wiley Interdiscip. Rev. RNA 2019, 10, e1543. [Google Scholar] [CrossRef]
  19. Basera, A.; Hull, R.; Demetriou, D.; Bates, D.O.; Kaufmann, A.M.; Dlamini, Z.; Marima, R. Competing Endogenous RNA (ceRNA) Networks and Splicing Switches in Cervical Cancer: HPV Oncogenesis, Clinical Significance and Therapeutic Opportunities. Microorganisms 2022, 10, 1852. [Google Scholar] [CrossRef]
  20. Arraiano, C.M.; Andrade, J.M.; Domingues, S.; Guinote, I.B.; Malecki, M.; Matos, R.G.; Moreira, R.N.; Pobre, V.; Reis, F.P.; Saramago, M.; et al. The critical role of RNA processing and degradation in the control of gene expression. FEMS Microbiol. Rev. 2010, 34, 883–923. [Google Scholar] [CrossRef]
  21. Wu, X.; Zhong, Y.; Chen, Q.; Zhang, X.; Zhang, H. Enhancer of mRNA Decapping protein 4 (EDC4) interacts with replication protein a (RPA) and contributes to Cisplatin resistance in cervical Cancer by alleviating DNA damage. Hereditas 2020, 157, 41. [Google Scholar] [CrossRef]
  22. Wang, Q.; Guo, X.; Li, L.; Gao, Z.; Su, X.; Ji, M.; Liu, J. N6-methyladenosine METTL3 promotes cervical cancer tumorigenesis and Warburg effect through YTHDF1/HK2 modification. Cell Death Dis. 2020, 11, 911. [Google Scholar] [CrossRef]
  23. Yeganeh, M.; Hernandez, N. RNA polymerase III transcription as a disease factor. Genes Dev. 2020, 34, 865–882. [Google Scholar] [CrossRef] [PubMed]
  24. de Klerk, E.; ’t Hoen, P.A.C. Alternative mRNA transcription, processing, and translation: Insights from RNA sequencing. Trends Genet. 2015, 31, 128–139. [Google Scholar] [CrossRef] [PubMed]
  25. Yuan, F.; Hankey, W.; Wagner, E.J.; Li, W.; Wang, Q. Alternative polyadenylation of mRNA and its role in cancer. Genes Dis. 2021, 8, 61–72. [Google Scholar] [CrossRef] [PubMed]
  26. Reyes, A.; Huber, W. Alternative start and termination sites of transcription drive most transcript isoform differences across human tissues. Nucleic Acids Res. 2018, 46, 582–592. [Google Scholar] [CrossRef]
  27. Kim, E.; Goren, A.; Ast, G. Insights into the connection between cancer and alternative splicing. Trends Genet. 2008, 24, 7–10. [Google Scholar] [CrossRef]
  28. Zhang, Y.; Liu, L.; Qiu, Q.; Zhou, Q.; Ding, J.; Lu, Y.; Liu, P. Alternative polyadenylation: Methods, mechanism, function, and role in cancer. J. Exp. Clin. Cancer Res. 2021, 40, 51. [Google Scholar] [CrossRef]
  29. Monteuuis, G.; Schmitz, U.; Petrova, V.; Kearney, P.S.; Rasko, J.E.J. Holding on to Junk Bonds: Intron Retention in Cancer and Therapy. Cancer Res. 2021, 81, 779–789. [Google Scholar] [CrossRef]
  30. Dvinge, H.; Bradley, R.K. Widespread intron retention diversifies most cancer transcriptomes. Genome Med. 2015, 7, 45. [Google Scholar] [CrossRef]
  31. Qi, F.; Li, Y.; Yang, X.; Wu, Y.-P.; Lin, L.-J.; Liu, X.-M. Significance of alternative splicing in cancer cells. Chin. Med. J. 2020, 133, 221–228. [Google Scholar] [CrossRef]
  32. Rossello, A.; Steinle, A.; Poggi, A.; Zocchi, M.R. Editorial: ADAM10 in Cancer Immunology and Autoimmunity: More Than a Simple Biochemical Scissor. Front. Immunol. 2020, 11, 1483. [Google Scholar] [CrossRef]
  33. Guo, H.; Yang, S.; Li, S.; Yan, M.; Li, L.; Zhang, H. LncRNA SNHG20 promotes cell proliferation and invasion via miR-140-5p-ADAM10 axis in cervical cancer. Biomed. Pharmacother. 2018, 102, 749–757. [Google Scholar] [CrossRef] [PubMed]
  34. Spardy, N.; Covella, K.; Cha, E.; Hoskins, E.E.; Wells, S.I.; Duensing, A.; Duensing, S. Human papillomavirus 16 E7 oncoprotein attenuates DNA damage checkpoint control by increasing the proteolytic turnover of claspin. Cancer Res. 2009, 69, 7022–7029. [Google Scholar] [CrossRef] [PubMed]
  35. Benevolo, M.; Musio, A.; Vocaturo, A.; Donà, M.G.; Rollo, F.; Terrenato, I.; Carosi, M.; Pescarmona, E.; Vocaturo, G.; Mottolese, M. Claspin as a biomarker of human papillomavirus-related high-grade lesions of uterine cervix. J. Transl. Med. 2012, 10, 132. [Google Scholar] [CrossRef] [PubMed]
  36. Deka, B.; Chandra, P.; Yadav, P.; Rehman, A.; Kumari, S.; Kunnumakkara, A.B.; Singh, K.K. RNPS1 functions as an oncogenic splicing factor in cervical cancer cells. IUBMB Life 2023, 75, 514–529. [Google Scholar] [CrossRef] [PubMed]
  37. Artaza-Irigaray, C.; Molina-Pineda, A.; Aguilar-Lemarroy, A.; Ortiz-Lazareno, P.; Limón-Toledo, L.P.; Pereira-Suárez, A.L.; Rojo-Contreras, W.; Jave-Suárez, L.F. E6/E7 and E6* From HPV16 and HPV18 Upregulate IL-6 Expression Independently of p53 in Keratinocytes. Front. Immunol. 2019, 10, 1676. [Google Scholar] [CrossRef]
  38. Andrews, S. FastQC: A Quality Control Tool for High Throughput Sequence Data. 2010. Available online: http://www.bioinformatics.babraham.ac.uk/projects/fastqc/ (accessed on 5 April 2021).
  39. Afgan, E.; Baker, D.; Batut, B.; van den Beek, M.; Bouvier, D.; Čech, M.; Chilton, J.; Clements, D.; Coraor, N.; Grüning, B.A.; et al. The Galaxy platform for accessible, reproducible, and collaborative biomedical analyses: 2018 update. Nucleic Acids Res. 2018, 46, W537–W544. [Google Scholar] [CrossRef]
  40. Dobin, A.; Davis, C.A.; Schlesinger, F.; Drenkow, J.; Zaleski, C.; Jha, S.; Batut, P.; Chaisson, M.; Gingeras, T.R. STAR: Ultrafast universal RNA-seq aligner. Bioinformatics 2013, 29, 15–21. [Google Scholar] [CrossRef]
  41. García-Alcalde, F.; Okonechnikov, K.; Carbonell, J.; Cruz, L.M.; Götz, S.; Tarazona, S.; Dopazo, J.; Meyer, T.F.; Conesa, A. Qualimap: Evaluating next-generation sequencing alignment data. Bioinformatics 2012, 28, 2678–2679. [Google Scholar] [CrossRef]
  42. Love, M.I.; Huber, W.; Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014, 15, 550. [Google Scholar] [CrossRef]
  43. Mootha, V.K.; Lindgren, C.M.; Eriksson, K.-F.; Subramanian, A.; Sihag, S.; Lehar, J.; Puigserver, P.; Carlsson, E.; Ridderstråle, M.; Laurila, E.; et al. PGC-1alpha-responsive genes involved in oxidative phosphorylation are coordinately downregulated in human diabetes. Nat. Genet. 2003, 34, 267–273. [Google Scholar] [CrossRef] [PubMed]
  44. Subramanian, A.; Tamayo, P.; Mootha, V.K.; Mukherjee, S.; Ebert, B.L.; Gillette, M.A.; Paulovich, A.; Pomeroy, S.L.; Golub, T.R.; Lander, E.S.; et al. Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles. Proc. Natl. Acad. Sci. USA 2005, 102, 15545–15550. [Google Scholar] [CrossRef] [PubMed]
  45. Kanehisa, M.; Goto, S. KEGG: Kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000, 28, 27–30. [Google Scholar] [CrossRef] [PubMed]
  46. Pertea, M.; Pertea, G.M.; Antonescu, C.M.; Chang, T.-C.; Mendell, J.T.; Salzberg, S.L. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat. Biotechnol. 2015, 33, 290–295. [Google Scholar] [CrossRef] [PubMed]
  47. Soneson, C.; Love, M.I.; Robinson, M.D. Differential analyses for RNA-seq: Transcript-level estimates improve gene-level inferences. F1000Research 2015, 4, 1521. [Google Scholar] [CrossRef]
  48. Ritchie, M.E.; Phipson, B.; Wu, D.; Hu, Y.; Law, C.W.; Shi, W.; Smyth, G.K. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015, 43, e47. [Google Scholar] [CrossRef]
  49. Finn, R.D.; Bateman, A.; Clements, J.; Coggill, P.; Eberhardt, R.Y.; Eddy, S.R.; Heger, A.; Hetherington, K.; Holm, L.; Mistry, J.; et al. Pfam: The protein families database. Nucleic Acids Res. 2014, 42, D222–D230. [Google Scholar] [CrossRef]
  50. Petersen, T.N.; Brunak, S.; von Heijne, G.; Nielsen, H. SignalP 4.0: Discriminating signal peptides from transmembrane regions. Nat. Methods 2011, 8, 785–786. [Google Scholar] [CrossRef]
  51. Wang, L.; Park, H.J.; Dasari, S.; Wang, S.; Kocher, J.-P.; Li, W. CPAT: Coding-Potential Assessment Tool using an alignment-free logistic regression model. Nucleic Acids Res. 2013, 41, e74. [Google Scholar] [CrossRef]
  52. Mészáros, B.; Erdos, G.; Dosztányi, Z. IUPred2A: Context-dependent prediction of protein disorder as a function of redox state and protein binding. Nucleic Acids Res. 2018, 46, W329–W337. [Google Scholar] [CrossRef]
Figure 1. Gene set enrichment analysis (GSEA). Plot shows normalized enrichment scores of top KEGG pathways. Pathways related to mRNA processing are highlighted in dark grey.
Figure 1. Gene set enrichment analysis (GSEA). Plot shows normalized enrichment scores of top KEGG pathways. Pathways related to mRNA processing are highlighted in dark grey.
Ijms 24 08347 g001
Figure 2. Alternative transcription events in isoforms. Plot shows eight AS mechanisms, including alternative transcription start sites (ATSS), alternative transcription termination sites (ATTS), alternative 3-end acceptor site (A3), alternative 5-end donor site (A5), exon skipping (ES), intron retention (IR), mutually exclusive exon (MEE), and multiple exon skipping (MES). HAC: HaCaT cells.
Figure 2. Alternative transcription events in isoforms. Plot shows eight AS mechanisms, including alternative transcription start sites (ATSS), alternative transcription termination sites (ATTS), alternative 3-end acceptor site (A3), alternative 5-end donor site (A5), exon skipping (ES), intron retention (IR), mutually exclusive exon (MEE), and multiple exon skipping (MES). HAC: HaCaT cells.
Ijms 24 08347 g002
Figure 3. Gains or losses in alternative splicing events.
Figure 3. Gains or losses in alternative splicing events.
Ijms 24 08347 g003
Figure 4. Overall isoform switch consequences. (A) for E5 from HPV-16 and HPV-84 versus control. (B) for E6 from HPV-16 and HPV-84 versus control. (C) for E7 from HPV-16 and HPV-84 versus control.
Figure 4. Overall isoform switch consequences. (A) for E5 from HPV-16 and HPV-84 versus control. (B) for E6 from HPV-16 and HPV-84 versus control. (C) for E7 from HPV-16 and HPV-84 versus control.
Ijms 24 08347 g004
Figure 5. Isoform switch events reveal a little overlap between isoforms induced by each HPV oncogene. Overlaps between HPV16 and HPV84 are separated based on oncogene type (A) for E5, (B) for E6, and (C) for E7. List of genes with increased coding isoform usage (D) control versus HPV16 and HPV84. (E) HPV16 versus HPV84.
Figure 5. Isoform switch events reveal a little overlap between isoforms induced by each HPV oncogene. Overlaps between HPV16 and HPV84 are separated based on oncogene type (A) for E5, (B) for E6, and (C) for E7. List of genes with increased coding isoform usage (D) control versus HPV16 and HPV84. (E) HPV16 versus HPV84.
Ijms 24 08347 g005
Figure 6. Isoform switch plots for genes previously related to cervical cancer when comparing HPV-16 to HPV-84: (A) ADAM10; (B) CLSPN; (C) RNPS1; *—p-value < 0.05, ***—p-value <0.001, ns—no significant difference.
Figure 6. Isoform switch plots for genes previously related to cervical cancer when comparing HPV-16 to HPV-84: (A) ADAM10; (B) CLSPN; (C) RNPS1; *—p-value < 0.05, ***—p-value <0.001, ns—no significant difference.
Ijms 24 08347 g006
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Nasiri-Aghdam, M.; Garcia-Chagollan, M.; Pereira-Suarez, A.L.; Aguilar-Lemarroy, A.; Jave-Suarez, L.F. Splicing Characterization and Isoform Switch Events in Human Keratinocytes Carrying Oncogenes from High-Risk HPV-16 and Low-Risk HPV-84. Int. J. Mol. Sci. 2023, 24, 8347. https://doi.org/10.3390/ijms24098347

AMA Style

Nasiri-Aghdam M, Garcia-Chagollan M, Pereira-Suarez AL, Aguilar-Lemarroy A, Jave-Suarez LF. Splicing Characterization and Isoform Switch Events in Human Keratinocytes Carrying Oncogenes from High-Risk HPV-16 and Low-Risk HPV-84. International Journal of Molecular Sciences. 2023; 24(9):8347. https://doi.org/10.3390/ijms24098347

Chicago/Turabian Style

Nasiri-Aghdam, Maryam, Mariel Garcia-Chagollan, Ana Laura Pereira-Suarez, Adriana Aguilar-Lemarroy, and Luis Felipe Jave-Suarez. 2023. "Splicing Characterization and Isoform Switch Events in Human Keratinocytes Carrying Oncogenes from High-Risk HPV-16 and Low-Risk HPV-84" International Journal of Molecular Sciences 24, no. 9: 8347. https://doi.org/10.3390/ijms24098347

APA Style

Nasiri-Aghdam, M., Garcia-Chagollan, M., Pereira-Suarez, A. L., Aguilar-Lemarroy, A., & Jave-Suarez, L. F. (2023). Splicing Characterization and Isoform Switch Events in Human Keratinocytes Carrying Oncogenes from High-Risk HPV-16 and Low-Risk HPV-84. International Journal of Molecular Sciences, 24(9), 8347. https://doi.org/10.3390/ijms24098347

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop