Integrated Analysis of Virus and Host Transcriptomes in Cervical Cancer in Asian and Western Populations

doi:10.21203/rs.3.rs-42875/v1

Download PDF

Research

Integrated Analysis of Virus and Host Transcriptomes in Cervical Cancer in Asian and Western Populations

https://doi.org/10.21203/rs.3.rs-42875/v1

This work is licensed under a CC BY 4.0 License

Journal Publication

published 01 May, 2021

Read the published version in Genomics →

Version 1

posted

You are reading this latest preprint version

Background: Race may influence vulnerability to HPV and have an effect on cervical cancer survival independent of socioeconomic status. Integration of the virus and host transcriptomes from different populations provides an unprecedented opportunity to understand these racial disparities in the prevalence of HPV and its associated cervical cancers in addition to fundamentally advancing the knowledge of HPV-induced cervical carcinogenesis and progression.

Methods: We performed RNA-Seq analysis of 90 tumors and 39 adjacent normal tissues from cervical cancer patients at Zhejiang University (ZJU) in China, and conducted a comparative analysis with RNA-Seq data of 286 cervical cancers from TCGA.

Results: We found a significantly higher rate of HPV positives and HPV integrations in TCGA than in ZJU. In addition to LINC00393 and HSPB3 as new common integration hotspots in both populations, we found new hotspots such as SH2D3C and CASC8 in TCGA, and SCGB1A1 and ABCA1 in ZJU. We described the first, to our knowledge, virus-transcriptome-based classification of cervical cancer that is highly predictive of clinical outcome. Particularly, patients with expressed E5 performed better than those without E5 expression. However, the constituents of these virus-transcriptome-based tumor subtypes differ dramatically between the two populations. We further characterized the immune infiltration landscapes between different HPV statuses and revealed significantly elevated levels of regulatory T cells and M0 macrophages in HPV positive tumors, which were associated with poor prognosis.

Conclusions: These findings increase our understanding of the racial disparities in the prevalence of HPV and its associated cervical cancers between Asian and Western populations, and also have important implications in the classification of tumor subtypes, prognosis, and anti-canc

Cancer Biology

Asian

Cervical cancer

HPV

Immune infiltration

RNA-Seq

Survival

Virus integration

Virus genes expression profiles

Western populations

Cervical cancer is the second most common cancer affecting women worldwide and accounts for around 530,000 cases per year, with the prevalence varying between different races [1, 2, 3, 4]. Persistent infection with high-risk human papilloma viruses (HPVs) causes most cases [3]. HPVs contain a 7.9-kb circular double-stranded DNA genome that consists of four parts: an early region (E1, 2, 4–7 genes), a late region (L1, 2 genes), a long control region (LCR), and a small, highly variable, non-coding region (NCR) between E5 and L2. HPVs infect basal epithelial cells located in the cervical transformation zone and participate in tumor development and cervical malignant transformation partly through expressed viral oncogenes, in which E6 and E7 are most focused on for their roles in inactivating p53 and pRb [3].

However, systematic expression profiles of viral genes and their association with molecular and clinical features of cervical cancer are poorly understood. In some cases, HPVs integrate into the host genome and trigger substantial host genome alterations [5]. The viral integration patterns have been well studied and several integration hotspots have been identified in several independent studies [5, 6], while it is unknown whether the above features are consistent across patients from different races and ethnicities. There may be nonnegligible differences in genetic susceptibility to HPV infection among races and ethnicities. Asian/Pacific Islanders have the lowest HPV prevalence of 42.2% compared to blacks (66.2%), and whites (71.5%) [2]. More black and Hispanic women get HPV-associated cervical cancer than women of other races and ethnicities [4]. Race may influence vulnerability to HPV [7] and have an effect on cervical cancer survival independent of socioeconomic status [1]. Therefore, comparative studies are needed to understand these racial disparities in HPV prevalence and the rate of death from cervical cancer.

Though effective HPV vaccines are available for primary prevention, there are still people who are already infected with HPV or diagnosed with cervical cancer, who cannot benefit from vaccines, and for whom curative therapeutic approaches are needed. Patients with tumor recurrence or advanced stage are limited to chemotherapy and radiation. Emerging combination immunotherapy is an attractive option for patients with advance-stage cancer or metastatic disease. Characterizing tumor immune infiltration and its differences between HPV-positive and -negative cervical cancers will help in designing an effective treatment of combination immunotherapy for these patients [8].

Further studies among Asian populations are required to understand racial disparities in HPV prevalence and its associated cervical cancer. Here, we performed RNA sequencing (RNA-Seq) of 90 tumors and 39 adjacent normal tissues from patients with cervical cancer recruited at Zhejiang University (ZJU), China (Table 1), and conducted a comparative analysis with RNA-Seq data of 286 cervical cancers from The Cancer Genome Atlas (TCGA) project (Table S1). The primary aim of this study was to compare HPV infection, integration and virus gene expression profiles from patients with cervical cancer between ZJU and TCGA populations. Our secondary aim was to characterize virus-induced host response and immune infiltration landscapes in the two populations. Integration of virus and host transcriptomes from the two populations provides an unprecedented opportunity to understand these racial disparities in the prevalence of HPV and its associated cervical cancer in addition to fundamentally advancing the knowledge of HPV-induced cervical carcinogenesis and progression.

Table 1

Demographic and clinicopathologic characteristics of ZJU cervical cancer patients.
Variable Name	No. of patients
Age
< 40	14
40–50	45
> 50	31
Stage
Stage I	57
Stage II	30
Stage III	3
Pathology
Squamous cell carcinoma	73
Adenocarcinoma	10
Adenosquamous	5
Others	2
Tumor Size
Small	43
Large	18
Unknown	29

Sample acquisition and processing

A total of 90 fresh cervical cancer tissues and 39 adjacent noncancerous cervix tissues were collected from September 2014 to December 2016. All samples were obtained at the time of diagnosis before any treatment was administered. Subjects in our RNA sequencing study included 73 squamous cell carcinomas (SCCs), 10 adenocarcinomas, 5 adenosqumous carcinomas and 2 mixtures of the above cancer types (Table 1). Fresh tissues were collected, snap-frozen in liquid nitrogen, and stored at -80 °C. Clinical information (e.g., tumor stage, tumor size, and smoking status) was extracted from pathology reports generated at the time of tissue collection. Review of hematoxylin and eosin (H&E) slides was performed by a gynecologic pathologist to confirm the diagnosis. Tumor cell content of sections from frozen cervical tissue samples embedded in optimal cutting temperature (OCT) medium was confirmed to be > 70% by H&E staining. As to tumor size, clinically visible lesion 4.0 cm or less in greatest dimension by ultrasonography was defined as small, otherwise it was defined as large.

Library preparation for RNA sequencing

Total RNA was extracted from 10 to 20 mg of tissue with the RNeasy Mini Kit (Qiagen, Hilden, Germany). The RNA integrity of each sample was assessed using agarose gel electrophoresis. The concentration of total RNA was determined by NanoDrop 2000 (Thermo Fisher Scientfic, Wilmington, DE). The full transcriptome sequencing libraries were prepared using the NEBNext Ultra RNA Directional Library Prep Kit (NEB#420) according to the manufacturer’s instructions. The 150 bp paired-end sequencing was performed on a HiSeq 3000 sequencer (Illumina, San Diego, CA).

RNA-Seq analysis

Before the analysis, adapter sequences were first removed from the output RNA-Seq sequence reads using cutadapt (http://code.google.com/p/cutadapt/). Reads with low base quality (< 13) were further trimmed and, if less than 25 base pairs, removed by FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). The trimmed sequence reads were then aligned to the manually curated reference genome and gene annotation database using TopHat2 (v2.0.13) with mammalian default parameters [9]. The reference genome was combined from human reference (hg19) and reference sequences of 143 types of HPVs downloaded from a web resource of Papillomavirus Episteme (PaVE, http://pave.niaid.nih.gov). The gene annotation database was combined from Ensembl genes v75 and gene annotation in PaVE. Transcript construction, quantification and normalization of transcript abundance were performed using Cufflinks (v2.2.1) [10]. Transcript or gene expression levels were presented in fragments per kilobase of transcript per million fragments mapped (FPKM).

RNA-Seq data of 286 cervical cancers from TCGA (Table S1) were downloaded from the GDC Data Portal (https://portal.gdc.cancer.gov/). The same analyses as those for ZJU were performed on RNA-Seq data in TCGA.

HPV typing and detection of integration

The threshold of HPV expression for defining HPV positives and negatives was generated from the density plot of HPV expression in each dataset separately. Both density plots had two peaks, one big and one small, with the shared threshold being taken in the middle of the two peaks (Figure S1). For HPV positive samples, their HPV genotypes were determined according to the type of expressed HPV. For those samples infected with multiple HPV types, the genotype with the highest expression was taken.

A search for HPV integration events was performed using the ChimeraScan (v0.4.5) [11], deFuse (v0.6.2) [12], and SOAPfuse (v1.28) [13] programs. The trimmed sequence reads were then aligned to the reference genome and gene annotation database using the above gene fusion discover tools [11, 12, 13]. All HPV integration events in tumors that were detected in any normal samples were filtered out. We further filtered HPV integration sites with less than 2 reads spanning the breakpoint between human and virus genomes (chimeric reads). Then, the predicted integration events from the three programs were integrated and combined with genomic annotation to generate a list of candidates for experimental validations.

Since there may be deviation of positions of detected breakpoints in HPV integration, the breakpoints in the ± 5 bp of splice sites were considered as the splice donor or acceptor. Virus-host transcripts whose breakpoint pair consisted of a splice donor and a splice acceptor were considered as spliced. The identity of virus and human sequences around integration sites were calculated from the sequence provided by deFuse [12], and sequences considered to have been spliced as mentioned above were filtered.

Unsupervised hierarchical cluster of tumors based on HPV gene expression

The FPKM values of HPV genes were transformed by log2(FPKM + 1). Then unsupervised hierarchical clustering was performed using a Euclidean distance metric and ward.D2 clustering without reorder of row dendrograms. The Kruskal–Wallis test was used to identify differentially expressed genes among five subgroups in the two datasets separately. Gene ontology (GO) enrichment analysis based on Fisher’s Exact test was performed on genes of subclass A and B separately in the two datasets. GO terms was filtered using fold enrichment > 2 as well as Benjamini–Hochberg corrected p value < 0.05 [14, 15].

Host gene expression analysis

The ranks of expression levels of genes at or near the HPV integration sites in tumors were generated by comparing their expression levels with those of the same genes across the remaining tumors. Then these ranks were normalized as 0 to 1 by dividing total tumor sample size, in which 1 represented the highest and 0 represented the lowest expression level among samples. Genes near the integration sites were those genes that were overlapped with the 500 kb upstream or downstream of the integration sites. The chromosomal locations of 89 common fragile sites were used[16]. The genomic coordinates of regulatory regions were obtained from Reg2Map (https://personal.broadinstitute.org/meuleman/reg2map/HoneyBadger2-intersect_release/).

Inference of immune cell ratios

RNA-Seq data in TCGA and ZJU, and two microarray data sets (GSE6791 and GSE63514) were used for inferring tumor immune cell infiltrations. The raw CEL files were downloaded from the Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/) and were processed by the MAS5 algorithm. RNA-Seq and microarray data were quantile normalized in each dataset separately before the analysis. A linear support vector regression based method, CIBERSORT (v1.03) [17], was applied to estimate relative ratios of 22 leukocytes for each sample. Samples with empirical p value greater than 0.05 were excluded for further analysis. A permutation test was then applied to evaluate differential distribution of these inferred immune cell subtypes between groups using the R package “permute” (https://www.r-project.org).

Survival analyses

Survival associated immune cells or transcripts were assessed by univariate Cox proportional hazards model. Survival distributions in different HPV subsets were visualized by Kaplan–Meier curves, and the significance was assessed by a log-rank test. The events of overall survival were defined as death, while the recurrence-free interval was ended by any disease recurrence, but death was not included. Survival analyses were performed using the R package “survival”.

PCR verification experiment of HPV integration

cDNA synthesis was performed with 1 ug of total RNA with a reverse transcription kit (Takara, Tokyo, Japan). PCR and Sanger sequencing were used to verify the HPV integration. PCR primers were designed on the basis of the paired-end assembled fragments, in which one primer was located in the human genome and the other in the HPV genome. PCR was performed using a ProFlex™ PCR thermal cycler (Life technology, Carlsbad, CA). PCR products were purified and sequenced using the Sanger sequencing method (Table S2).

Characterization of HPV integration in cervical cancer

The HPV positive rates significantly differed between the two populations with 90.0% (81 of 90) in the ZJU group and 94.1% (269 of 286) in the TCGA group (Table S1 and Figure S2a, binomial test; p = 4.7e-2). Overall, the proportion of HPV subtypes also differed between the ZJU and TCGA samples, such as the rate of HPV subtypes 16 and 58 which were much higher in the ZJU while the HPV18 samples accounted for more in the TCGA (Figure S2b). These results imply that the HPV-based tumor types vary in the two populations.

HPV integration events were observed in 57 of 81 (70%) and 239 of 269 (88%) HPV-positive tumor tissues in ZJU and TCGA, respectively (Fig. 1a). The binomial test proved that the integration rate was significantly different between the two populations (p = 1.3e-5). This difference is not negligible as longer reads and larger amounts of sequence data were used in the ZJU group, which favor detection of HPV integration (Figure S2c). The overall different rate of un-cleared HPV infection and HPV integration in the two populations may be due to the local virus sublineage diversity by geographical disparities through evolution or different genetic susceptibility among different races or ethnicities. In both populations, HPV integration occurred in all viral genes and affected all human chromosomes. In addition, MIPOL1, PTPN13, VMP1 and TP63 which have been previously reported were identified in both the ZJU and TCGA samples [5, 18], and we also found LINC00393 and HSPB3 as new common integration hotspots in both populations. Several integration breakpoints such as PVT1, RAD51B, MYC and ERBB2 that were previously reported in TCGA samples were also identified in our analysis. Furthermore, we found new hotspot sites such as SH2D3C and CASC8 in the TCGA samples, and SCGB1A1 and ABCA1 in the ZJU samples. Validation experiments proved the reliability of these HPV integration events (Table S2). In contrast to a previous study that revealed that common fragile sites (CFS) are preferential targets for HPV integration [19], we did not observe such a preference for HPV integration in CFS (Figure S3).

Genomic loci affected by HPV integration were significantly associated with higher gene expression (Mann–Whitney two-sided test; p = 1.1e-4 in ZJU and p = 2.3e-21 in TCGA), which is in good agreement with a previous study [5]. To explore how HPV integration affects the expression of genes around integration sites, we divided the integration sites into five categories according to the relative expression levels of integrated genes and investigated the expression of adjacent genes in each category. Density plots highlighted that the higher the expression of an integrated gene, the higher the expression of adjacent genes (Fig. 1b). This phenomenon was not likely ascribed to the direct effect of integration on transcription factor binding sites or the overall expression levels of samples with HPV integration (Figure S4a). A similar trend was observed in TCGA samples (Figure S4b). These findings revealed several novel characteristics of host genome alteration resulting from HPV integration.

Microhomologous sequences between human and HPV at integration sites

Then we evaluated the identity of human and virus sequences of ± 15 bp around all of integration sites identified in ZJU and TCGA samples (Fig. 2, and Figure S5.). We first randomly selected the human sequence and virus sequence to draw scatter plots for the number of identical nucleotides and the maximum number of consecutive identical nucleotides. As for the results from the human and virus sequences around the integration sites, the points drifted more to the right top of the coordinate system than any two randomly selected human and virus sequences. This supported the idea that HPV may integrate into the host genome through a microhomology-mediated DNA repair pathway during host DNA damage response [6].

The relevance of HPV gene expression to cervical cancer patient outcome

Unsupervised hierarchical clustering of HPV genes detected five distinct clusters consistently in the ZJU and TCGA tumor samples: group 0 with an overall low expression level consistent with defined HPV negative; group 1 expressing mainly E6 and E7; group 2 additionally expressing E6SI compared with group 1; group 3 with higher expression of E1 and E5, partly with expressed E4; and group 4 expressing mainly E5 and E6 but not E4 (Fig. 3a and 3b, and Table S1). E6SI is a spliced transcript encoding active E6 oncoprotein. As to histological types, group 0 contained a mixture of adenocarcinoma, adenosquamous and squamous cell carcinoma samples, whereas squamous cell carcinoma accounted for the majority of samples in the other groups. Except for the above features, we did not find any significant difference in stage, tumor size, and smoking status among these five groups.

Then we investigated whether there were any molecular features and clinical outcomes associated with these groups. Consistently in these two populations, we found enrichment of HPV18 samples in group 2 (hypergeometric test, p = 1.1e-2 in ZJU and p = 1.9e-16 in TCGA) and HPV58 samples in group 3 (hypergeometric test, p = 3.4e-6 in ZJU and p = 1.5e-2 in TCGA). HPV58 is relatively prevalent in China and other Asian countries. Notably, our data showed that E6 and E6SI were not detectable in most HPV58 samples, whereas E1, E4, E5 and E7 were expressed in the HPV58 samples. Meanwhile, HPV18 subtypes were correlated with high expression levels of E6SI (t-test, p = 1.6e-2 in ZJU and p = 1.7e-6 in TCGA) as previously reported [18]. These characteristics suggested that the expression of HPV genes not only varies within an HPV subtype such as HPV16, but also depends in part on a specific tumor subtype, giving us new insights into the potential role of HPV infection via its gene expression pattern beyond HPV subtypes. Furthermore, HPV integration frequency was low in groups 3 and 4 compared with nearly 100% in groups 1 and 2, implying that lack of E5 expression is perhaps propitious to HPV integration.

The E2/E6 ratio using quantitative PCR is often used to detect HPV integrations. The rationale behind the detection is that HPV integration can lead to the disruption of the viral E2 gene, whereas E6 and E7 are retained [20]. However, our data showed it is not always true that samples with HPV integration lack E2 expression. E2 expression was observed in some samples with HPV integration in our analyses. On the other hand, E6 expression was disrupted in some cases. Therefore, the E2/E6 ratio may not be effective for the detection of HPV integrations.

Kaplan–Meier curves revealed that the overall survival rate and recurrence-free internal survival rate were significantly different among the above five HPV-based groups (Fig. 3c; p = p = 9.14e-3 and 1.52e-3, respectively) in the TCGA samples. Group 4 exhibited the most favorable outcome, yet poorer prognosis could not be distinguished among the other groups in terms of overall survival. The five groups showed distinct outcomes in disease recurrence, with group 4 performing the best and group 2 displaying the worst outcome. Overall, the groups with expressed E5 performed better than the groups without E5 expression. These results highlighted the relevance of E5 to cancer recurrence, which not only offered a powerful prognosis biomarker, but also revealed the potential role of E5 in the development of cervical cancer. A specific HPV gene expression pattern is a better prognostic predictor than a single HPV gene; E5 seems to have a better prognostic value than well-focused E6 and E7.

Although there were similar HPV-gene based tumor subgroups and their associated host gene expression profiles in the two populations, the constituents of subgroups varied substantially. The proportions of group 2 in the two populations were nearly identical, while group 4 with the best survival outcome accounted for the most part of cases (35.5%) in ZJU, compared to 15.0% in TCGA, and group 1 with the worse outcome accounted for the most part of the cases (34.6%) in TCGA, compared to 20.0% in ZJU (Figure S2d).

To assess whether the above HPV subgroups are recapitulated at the gene level on the host genome, we detected genes associated with the HPV subgroups by Kruskal–Wallis test. Among 6,095 and 10,360 genes associated with HPV subgroups in ZJU and TCGA, respectively, there were 3,569 overlapped genes between the two populations. Both populations showed the most distinct clusters between group 0 and other groups, indicating there was a mostly different gene expression pattern between HPV negative and positive samples (Fig. 3a and 3b). Further comparison of enriched GO terms revealed the functional similarity of genes in subclass A between the two populations, as well as genes in subclass B between the two populations (Fig. 3d). Overall, genes in subclass B were mainly related to cell cycle and DNA repair. In total, these results represented the robust host response to viral gene expression within specific HPV subgroups.

RNA splicing affects HPV gene expression

We also explored potential factors affecting HPV gene expression. We first examined the relationship between integration and HPV gene expression among HPV positive samples. E7 was significantly up-regulated and E5 was significantly down-regulated in tumor samples with integration (Fig. 4a and Figure S6a). In most integration events, the virus gene is expressed as a virus-host fusion transcript and some of them undergo RNA splicing. Such a virus-host transcript can be spliced from a splice donor in the viral genome to a splice acceptor in the host DNA [21]. In our HPV16 positive tumors, nucleotides (nt) 226 and 880 were the most common splice donors during the splicing of virus-host transcripts (Fig. 4b and Figure S6b). Although both breakpoints of nt 226 and 880 were detected among the samples of the four groups (68% and 74.5% in ZJU, 52.2% and 59.4% in TCGA), these integration events were more apt to occur in group 1 and group 2 (hypergeometric test, p = 5.3e-4 and 2.4e-6 in ZJU, p = 7.5e-7 and 1.5e-7 in TCGA). Interestingly, we also found the complementary situation that the transcript was spliced from a splice donor in the host DNA to a splice acceptor in the virus genome, though its incidence was less common, of which nt 3,358 was the most common splice acceptor during splicing in both populations (17.0% in ZJU and 25.7% in TCGA). Notably, breakpoints of nt 3,358 were significantly enriched in group 3 and group 4 (hypergeometric test, p = 0 in ZJU and 1.2e-11 in TCGA). These results suggest that both breakpoints on the virus genome during integration and subsequent splice events of virus-host transcripts have an important impact on HPV gene expression. HPV integration not only causes dysregulation of host genes, but also shapes distinct HPV gene expression patterns influencing patient survival outcome.

Landscapes of tumor-infiltrating immune cells in cervical cancer

We investigated the immune cell infiltrative levels from RNA sequencing data in the ZJU and TCGA tumor samples and two additional microarray datasets of cervical cancer including GSE6791 [22] and GSE63514 [23]. We first compared compositional differences in immune cell populations between tumors and adjacent normal tissues in ZJU and the other two microarray datasets. Generally, tumor tissues detected a higher fraction of T lymphocytes and yet a lower fraction of B lymphocytes and mast cells than their adjacent normal tissues. For the same dendritic cells, the activated form has a higher fraction in tumor tissues than adjacent normal tissues, whereas the resting form is just the opposite. Immune cells that have higher fractions in tumors than in their adjacent tissues are more likely associated with poorer outcomes (Fig. 5). CD8 T cells and neutrophils were the most favorable and adverse prognostic populations, respectively, which are consistent with previous reports in anal squamous cell carcinoma and gastric cancer [24, 25].

Except for the significantly higher proportion of plasma cells in ZJU tumors, there was no other significant disparity in the immune cell compositions between ZJU and TCGA tumors (Fig. 6a). This suggests that overall levels of immune infiltration are consistent and reproducible across tumors of the two populations. To determine whether immune cell infiltrative levels are associated with HPV infection or integration, we compared immune cell compositions in HPV negative samples, and HPV positive samples with integration and without integration (Fig. 6b). There were more kinds of immune cells showing higher fractions in HPV positive samples compared with HPV negative samples, indicating stronger immune infiltration in tumors with un-cleared HPV infection. For instance, dendritic cells and regulatory T cells were significantly and consistently elevated in HPV positive samples compared with HPV negative samples in both populations. Except for the naïve B cell in TCGA, there was no significant compositional differences in immune cells between HPV positive samples with and without integration, especially in ZJU. These results suggest that HPV integration has little impact on immune infiltration in tumor microenvironments.

Race may influence vulnerability to HPV and have an effect on cervical cancer survival independent of socioeconomic status [7]. Integration analysis of virus and host transcriptome from different populations provides deeper insights into these racial disparities in the prevalence of HPV and its associated cervical cancer. Our analyses found that the rate of HPV positives and HPV integration was significantly higher in TCGA than in ZJU samples. Besides several common integration hotspots, each population has its own unique virus integration sites. Molecular classification of tumor samples using HPV gene expression profiles revealed distinct cancer subtypes that have significant clinical relevance to patients’ outcome; whereas the relative proportions of tumor subtypes were dramatically disparate between the two populations. These results demonstrated that not only the prevalence of HPV, but also the constituents of the cervical cancer patients differ among races and ethnicities.

Our integration analysis also revealed several novel characteristics of the host genome affected by HPV integration. HPVs integrate into the host genome and trigger substantial host genome alterations, leading to increased expression in target genes around the virus integration sites [5, 18]. We observed an increased expression of genes around the virus integration sites compared with the overall expression levels of genes without integration nearby. Furthermore, there is a strong corresponding relationship between the expression of target genes and adjacent ones. The higher the expression of target genes on integration sites, the higher the expression of adjacent genes. Since the range of affected genes expanded much beyond the integration sites, this correspondence is not ascribed to the direct effect of integration on transcription factor binding sites in the host genome. Our observation supports an increased expression of genes in the host genome that is likely driven by integrated viral promoters [5].

In addition to virus integration causing dysregulation of the host gene, it also shapes HPV gene expression patterns influencing the clinical outcome of patients. In most integration events, virus genes are expressed as virus-host fusion transcripts and some of them undergo RNA splicing. Virus-host fusion transcripts are often more stable than their viral counterparts and thus lead to increased expression of HPV genes [21, 26]. Our analysis revealed that both breakpoints on the virus genome during integration and subsequent splice events of virus-host transcripts have an important impact on HPV gene expression. In HPV16 positive tumors, nt 226 and 880 were the most common splice donor sites during the splicing of virus-host transcripts. We also identified a novel, complementary mode of splicing in tumor samples expressing E5, in which the transcript is spliced from a splice donor in the host DNA to a splice acceptor in the virus genome. In the complementary mode, nt 3,358 was the most common splice acceptor (Fig. 4b and Figure S6b). Future studies are required to understand the complex regulatory mechanisms underlying HPV RNA splicing. Understanding the regulation of HPV RNA splicing might provide valuable indications of its potential implications in the long-term persistence of HPV and its involvement in the development of cervical cancer.

For the first time, our analysis demonstrated that molecular classification using an HPV gene expression profile is highly predictive of clinical outcome for patients with cervical cancer. A virus genome-wide expression pattern is a better prognostic predictor than a single HPV gene. In particular, E5 has a better prognostic value than well-described E6 and E7; patients with expressed E5 performed better than those without E5 expression. While the oncogenic activities of E6 and E7 are well characterized, the role of E5 is still rather nebulous. The E5 gene is frequently deleted when the HPV genome is integrated during malignant progression [27]. Lack of E5 expression is perhaps propitious for HPV integration. Future studies are required to investigate if there is a negative feedback control of oncogene expression such as E6 and E7 by depleting E5 expression.

The immune infiltration landscape from the RNA-seq data revealed that the overall levels of tumor immune cell populations were largely consistent between ZJU and TCGA. Comparison of immune infiltrative levels between samples of different HPV statuses revealed the significantly elevated levels of dendritic cells and regulatory T cells in HPV positive tumors compared with HPV negative tumors, as well as significantly higher ratios of activated dendritic cells, monocytes and macrophages in tumors compared to normal samples. However, HPV integration has little impact on immune infiltration in tumor microenvironments.

It is well known that dendritic cells normally present in cervical tissues. In an inflammatory milieu, dendritic cells can be activated and subsequently up-regulate a series of molecules like major histocompatibility complex class I and II to facilitate the priming of naïve CD4 + and CD8 + T cells in lymph nodes [28]. Thus, the dendritic cells can be potent as therapeutic vaccines, which aim to activate the antigen-specific immune response and eradicate HPV-infected cells by simulating cytotoxic T cells [29]. Several clinical trials based on E7-pulsed dendritic cells have been carried out and received encouraging outcomes [29]. The increased levels of dendritic cells were found in tumors, especially in HPV positive tumors, hence the strategy to promote the migration and maturation of dendritic cells are of great potential for most cervical cancer patients [8].

Regulatory T cells (Tregs) play an important role in the maintenance of immunological self-tolerance to self-antigens and prevention of immune pathologies, which are related to the suppression of effective tumor immunity. A large number of tumor-infiltrating Tregs is substantiated in various cancers and is associated with poor clinical prognosis [30]. We observed higher levels of Tregs infiltration in HPV positive tumors than that in HPV negative tumors. Patients with higher Tregs infiltration level tend to have worse prognosis. A similar result has been illustrated in head and neck squamous cell carcinomas [31]. These findings suggested HPV positive tumors may benefit more from the immunotherapy targeting Tregs.

Overall elevated levels of monocytes and macrophages were found in cervical cancers compared to normal samples, in which the fractions of M0 and M1 macrophages were significantly increased. Accumulating evidence has suggested that the heterogeneous phenotypes of tumor-associated macrophages (TAM) are informative for patient prognosis. The M1 macrophages are considered involved in inflammatory response and antitumor immunity. In contrast, M2 macrophages exert anti-inflammatory and pro-tumorigenic activities[32]. Consistent with previous studies on other cancers [33], the correlation between increasing levels of M1 macrophages and superior survival was observed in cervical cancer, whereas M0 macrophages were found to be associated with poorer survival. The M0 and M2 macrophages increase cancer invasion abilities in human-original and mouse-original lung cancer models [34]. Future studies are required for further investigation of M0 macrophages as a promising immunotherapeutic target.

In summary, we performed integrated analysis of virus and host transcriptomes to characterize HPV infection, integration, RNA splicing, virus genes expression profiles, virus-induced host response, and tumor immune infiltration landscapes from patients with cervical cancer in the two populations. Our analysis found that the rate of HPV positives and HPV integration was significantly higher in patients from Western populations than an Asian population. We revealed the significant correspondence of increased expression levels between an integrated gene and its adjacent genes. We described the first, to our knowledge, virus transcriptome-based classification of cervical cancer that is highly predictive of clinical outcome. In particular, patients with expressed E5 performed better than those without E5 expression. Tregs and M0 macrophages are significantly elevated in HPV positive tumor samples and are associated with poor prognosis. Integration of virus and host transcriptomes from the two populations provides an unprecedented opportunity to understand these racial disparities in the prevalence of HPV and its associated cervical cancer in addition to fundamentally advancing the knowledge of HPV-induced cervical carcinogenesis and progression.

CFS: common fragile sites

FPKM: fragments per kilobase of transcript per million fragments mapped

GEO: Gene Expression Omnibus

GO: Gene ontology

H&E: hematoxylin and eosin

HPV: human papilloma virus

LCR: long control region

NCR: non-coding region

nt: nucleotides

OCT: optimal cutting temperature

PaVE: Papillomavirus Episteme

RNA-Seq: RNA sequencing

SCCs: squamous cell carcinomas

TAM: tumor-associated macrophages

TCGA: The Cancer Genome Atlas

Tregs: Regulatory T cells

ZJU: Zhejiang University

Ethics approval and consent to participate

This study was reviewed and approved by the Ethnics Committees of Women’s Hospital of Zhejiang University School of Medicine (Hangzhou, China). The study approval number is 20170124. The study was conducted in accordance with the International Ethical Guidelines for Biomedical Research Involving Human Subjects. All subjects provided informed consent to participate in the study.

Consent for publication

Not applicable.

Availability of data and material

The data that support the findings of this study are available from the corresponding author upon reasonable request. RNA-Seq raw fastq files for these samples were deposited in the National Center for Biotechnology Information (NCBI) database (https://www.ncbi.nlm.nih.gov/) under the BioProject accession number PRJNA644250 (https://dataview.ncbi.nlm.nih.gov/object/PRJNA644250?reviewer=bn1objg9igle1qsicqjg3f2849).

Competing interests

The authors have declared that no conflict of interests exist.

Funding

This work has been support by National Key Research and Development Program of China (2016YFC1302900), Medical Health Science and Technology Key Project of Zhejiang Provincial Health Commission (WKJ-ZJ-2007), the Key Project of Zhejiang Provincial Natural Science Foundation of. China (LZ20H160001), and National Natural Science Foundation of China (81772766, and 81871864).

Authors' contributions

YL, PL and WGL considered and designed the study. QQ performed data analysis and interpreted the analysis. XFL, QZ and MA provided assistance in the data analysis. QZ and ARL performed RNA-Seq library preparation and PCR validation. WL, BL, KL, LH and WFL collected tissues and performed pathological analysis. QQ and YL wrote the manuscript. All of the authors discussed and commented the study.

Acknowledgements

We thank TCGA for providing access to RNA-Seq data of cervical cancer. We thank the Bioinformatics Core Facility at Zhejiang University School of Medicine for providing computing capacity and anonymous reviewers for reading and commenting on the manuscript.

Adams SA, Fleming A, Brandt HM, Hurley D, Bolick-Aldrich S, Bond SM, Hebert JR. Racial disparities in cervical cancer mortality in an African American and European American cohort in South Carolina. J S C Med Assoc. 2009;105: 237-44.
Akogbe GO, Ajidahun A, Sirak B, Anic GM, Papenfuss MR, Fulp WJ, Lin HY, Abrahamsen M, Villa LL, Lazcano-Ponce E, et al. Race and prevalence of human papillomavirus infection among men residing in Brazil, Mexico and the United States. Int J Cancer. 2012;131: E282-91.
Emma J Crosbie MHE, Silvia Franceschi, Prof Henry C Kitchener. Human papillomavirus and cervical cancer. THE LANCET. 2013;382: 889–99.
Viens LJ, Henley SJ, Watson M, Markowitz LE, Thomas CC, Thompson TD, Razzaghi H, Saraiya M. Human Papillomavirus-Associated Cancers - United States, 2008-2012. MMWR Morb Mortal Wkly Rep. 2016;65: 661-6.
Ojesina AI, Lichtenstein L, Freeman SS, Pedamallu CS, Imaz-Rosshandler I, Pugh TJ, Cherniack AD, Ambrogio L, Cibulskis K, Bertelsen B, et al. Landscape of genomic alterations in cervical carcinomas. Nature. 2014;506: 371-5.
Hu Z, Zhu D, Wang W, Li W, Jia W, Zeng X, Ding W, Yu L, Wang X, Wang L, et al. Genome-wide profiling of HPV integration in cervical cancer identifies clustered genomic hot spots and a potential microhomology-mediated integration mechanism. Nat Genet. 2015;47: 158-63.
Burk RD, DeSalle R. The tango and tangle of human papillomavirus and the human genome. J Natl Cancer Inst. 2006;98: 1026-7.
Lin C-T, Wang C-N, Lai C-H. Immunotherapy for advanced or relapsed cervical cancer. Gynecology and Minimally Invasive Therapy. 2013;2: 3-7.
Kim D, Pertea G, Trapnell C, Pimentel H, Kelley R, Salzberg SL. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol. 2013;14: R36.
Trapnell C, Hendrickson DG, Sauvageau M, Goff L, Rinn JL, Pachter L. Differential analysis of gene regulation at transcript resolution with RNA-seq. Nat Biotechnol. 2013;31: 46-53.
Maher CA, Kumar-Sinha C, Cao X, Kalyana-Sundaram S, Han B, Jing X, Sam L, Barrette T, Palanisamy N, Chinnaiyan AM. Transcriptome sequencing to detect gene fusions in cancer. Nature. 2009;458: 97-101.
McPherson A, Hormozdiari F, Zayed A, Giuliany R, Ha G, Sun MG, Griffith M, Heravi Moussavi A, Senz J, Melnyk N, et al. deFuse: an algorithm for gene fusion discovery in tumor RNA-Seq data. PLoS Comput Biol. 2011;7: e1001138.
Jia W, Qiu K, He M, Song P, Zhou Q, Zhou F, Yu Y, Zhu D, Nickerson ML, Wan S, et al. SOAPfuse: an algorithm for identifying fusion transcripts from paired-end RNA-Seq data. Genome Biol. 2013;14: R12.
Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, Davis AP, Dolinski K, Dwight SS, Eppig JT, et al. Gene ontology: tool for the unification of biology. The Gene Ontology Consortium. Nat Genet. 2000;25: 25-9.
The Gene Ontology C. Expansion of the Gene Ontology knowledgebase and resources. Nucleic Acids Res. 2017;45: D331-D38.
Schwartz M, Zlotorynski E, Kerem B. The molecular basis of common and rare fragile sites. Cancer letters. 2006;232: 13-26.
Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, Hoang CD, Diehn M, Alizadeh AA. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. 2015;12: 453-7.
Cancer Genome Atlas Research Network, Albert Einstein College of Medicine, Analytical Biological Services, Barretos Cancer Hospital, Baylor College of Medicine, Beckman Research Institute of City of Hope, Buck Institute for Research on Aging, Canada's Michael Smith Genome Sciences Centre, Harvard Medical School, Helen F. Graham Cancer Center, et al. Integrated genomic and molecular characterization of cervical cancer. Nature. 2017;543: 378-84.
Thorland EC, Myers SL, Gostout BS, Smith DI. Common fragile sites are preferential targets for HPV16 integrations in cervical tumors. Oncogene. 2003;22: 1225-37.
Peitsaro P, Johansson B, Syrjanen S. Integrated human papillomavirus type 16 is frequently found in cervical cancer precursors as demonstrated by a novel quantitative real-time PCR technique. J Clin Microbiol. 2002;40: 886-91.
Wentzensen N, Ridder R, Klaes R, Vinokurova S, Schaefer U, Doeberitz M. Characterization of viral-cellular fusion transcripts in a large series of HPV16 and 18 positive anogenital lesions. Oncogene. 2002;21: 419-26.
Pyeon D, Newton MA, Lambert PF, den Boon JA, Sengupta S, Marsit CJ, Woodworth CD, Connor JP, Haugen TH, Smith EM, et al. Fundamental differences in cell cycle deregulation in human papillomavirus-positive and human papillomavirus-negative head/neck and cervical cancers. Cancer Res. 2007;67: 4605-19.
den Boon JA, Pyeon D, Wang SS, Horswill M, Schiffman M, Sherman M, Zuna RE, Wang Z, Hewitt SM, Pearson R, et al. Molecular transitions from papillomavirus infection to cervical precancer and cancer: Role of stromal estrogen receptor signaling. Proc Natl Acad Sci U S A. 2015;112: E3255-64.
Hu WH, Miyai K, Cajas-Monson LC, Luo L, Liu L, Ramamoorthy SL. Tumor-infiltrating CD8(+) T lymphocytes associated with clinical outcome in anal squamous cell carcinoma. J Surg Oncol. 2015;112: 421-6.
Li TJ, Jiang YM, Hu YF, Huang L, Yu J, Zhao LY, Deng HJ, Mou TY, Liu H, Yang Y, et al. Interleukin-17-Producing Neutrophils Link Inflammatory Stimuli to Disease Progression by Promoting Angiogenesis in Gastric Cancer. Clin Cancer Res. 2017;23: 1575-85.
Jeon S, Lambert PF. Integration of human papillomavirus type 16 DNA into the human genome leads to increased stability of E6 and E7 mRNAs: implications for cervical carcinogenesis. Proc Natl Acad Sci U S A. 1995;92: 1654-8.
Venuti A, Paolini F, Nasir L, Corteggio A, Roperto S, Campo MS, Borzacchiello G. Papillomavirus E5: the smallest oncoprotein with many functions. Mol Cancer. 2011;10: 140.
Patel S, Chiplunkar S. Host immune responses to cervical cancer. Curr Opin Obstet Gynecol. 2009;21: 54-9.
Menderes G, Black J, Schwab CL, Santin AD. Immunotherapy and targeted therapy for cervical cancer: an update. Expert Review of Anticancer Therapy. 2015;16: 83-98.
Sakaguchi AT, Shimon. Regulatory T cells in cancer immunotherapy. Cell Research. 2016;27: 109-18.
Mandal R, Senbabaoglu Y, Desrichard A, Havel JJ, Dalin MG, Riaz N, Lee KW, Ganly I, Hakimi AA, Chan TA, et al. The head and neck cancer immune landscape and its immunotherapeutic implications. JCI Insight. 2016;1: e89829.
Chanmee T, Ontong P, Konno K, Itano N. Tumor-associated macrophages as major players in the tumor microenvironment. Cancers. 2014;6: 1670-90.
Richards DM, Hettinger J, Feuerer M. Monocytes and Macrophages in Cancer: Development and Functions. Cancer Microenviron. 2013;6: 179-91.
Yuan A, Hsiao YJ, Chen HY, Chen HW, Ho CC, Chen YY, Liu YC, Hong TH, Yu SL, Chen JJ, et al. Opposite Effects of M1 and M2 Macrophage Subtypes on Lung Cancer Progression. Scientific reports. 2015;5: 14273.

Download PDF

Journal Publication

published 01 May, 2021

Read the published version in Genomics →

Version 1

posted

You are reading this latest preprint version

Integrated Analysis of Virus and Host Transcriptomes in Cervical Cancer in Asian and Western Populations

Status:

Journal Publication

Version 1

Abstract

Figures

Background

Materials And Methods

Sample acquisition and processing

Library preparation for RNA sequencing

RNA-Seq analysis

HPV typing and detection of integration

Unsupervised hierarchical cluster of tumors based on HPV gene expression

Host gene expression analysis

Inference of immune cell ratios

Survival analyses

PCR verification experiment of HPV integration

Results

Characterization of HPV integration in cervical cancer

Microhomologous sequences between human and HPV at integration sites

The relevance of HPV gene expression to cervical cancer patient outcome

RNA splicing affects HPV gene expression

Landscapes of tumor-infiltrating immune cells in cervical cancer

Discussion

Conclusion

Abbreviations

Declarations

References

Supplementary Files

Status:

Journal Publication

Version 1