Quantitative analysis of the splice variants expressed by the major hepatitis B virus genotypes

Abstract

Hepatitis B virus (HBV) is a major human pathogen that causes liver diseases. The main HBV RNAs are unspliced transcripts that encode the key viral proteins. Recent studies have shown that some of the HBV spliced transcript isoforms are predictive of liver cancer, yet the roles of these spliced transcripts remain elusive. Furthermore, there are nine major HBV genotypes common in different regions of the world, these genotypes may express different spliced transcript isoforms. To systematically study the HBV splice variants, we transfected human hepatoma cells, Huh7, with four HBV genotypes (A2, B2, C2 and D3), followed by deep RNA-sequencing. We found that 13-28 % of HBV RNAs were splice variants, which were reproducibly detected across independent biological replicates. These comprised 6 novel and 10 previously identified splice variants. In particular, a novel, singly spliced transcript was detected in genotypes A2 and D3 at high levels. The biological relevance of these splice variants was supported by their identification in HBV-positive liver biopsy and serum samples, and in HBV-infected primary human hepatocytes. Interestingly the levels of HBV splice variants varied across the genotypes, but the spliced pregenomic RNA SP1 and SP9 were the two most abundant splice variants. Counterintuitively, these singly spliced SP1 and SP9 variants had a suboptimal 5' splice site, supporting the idea that splicing of HBV RNAs is tightly controlled by the viral post-transcriptional regulatory RNA element.

Keywords: HBV; pgRNA; shotgun sequencing; transcriptome assembly.

Fig. 1.

RNA-seq analysis of the HBV and host transcriptomes. QC checking of the paired-end RNA-seq libraries was carried out using fastqc. Adapter sequences were trimmed from the RNA-seq reads using skewer. Trimmed reads were aligned to the human genome and HBV pgRNAs using star in 2-pass mode. Duplicated and multi-mapped reads were discarded from the binary alignment map (BAM) files using samtools. A post-alignment QC check was performed using picard tools. PacBio CCS reads were aligned to the HBV pgRNA using minimap2. HBV splice junctions were extracted and corrected using 2passtools and flair, respectively. Reference-based transcriptome assembly and quantification were carried out using stringtie, with a post-processing step focusing on the HBV spliced transcript isoforms. Splice site sequence contexts were scored using maxentscan. Completeness of RNA splicing was evaluated using ipsa. Reads mapped to human genes were quantified using mmquant, followed by differential gene expression analysis using deseq2. A list of differentially expressed genes was submitted to the david webserver for functional annotation analysis.

Fig. 2.

HBV genotypes expressed a wide variety of spliced transcript isoforms. (a) Proportions of the spliced transcripts in HBV RNAs. Only the spliced transcripts present in both biological replicates are shown. (b) Relative abundance of the HBV splice variants in genotypes A to D. See also Tables S4 and S5.

Fig. 3.

Distinct splicing profiles were observed across the HBV genotypes. The lollipop plot indicates the positions of splice sites relative to the EcoRI site of genotype C2. Blue and red colours indicate 5′ and 3′ splice sites, respectively. SP and pSP denote the known and putative spliced pgRNA transcripts, respectively (splice variants panel). These splice variants were reproducibly detected across the independent biological replicates of HBV-transfected Huh7. Grey dotted lines denote the positions of initiation codons of C, P, preS1 and X reading frames (ORFs panel). Read coverage is shown in grey (coverage panel). Arcs represent RNA-seq reads mapped across the splice junctions (supporting read counts in red colour). Only the splice junctions supported by ≥100 reads are shown for readability purposes. Blue and red vertical lines indicate the MaxEntScan scores of the 5′ and 3′ splice sites, respectively (coverage panel). A positive MaxEntScan score predicts a good splice site sequence context, whereas a negative score predicts a poor splice site sequence context. Three main scenarios were observed. ① The presence and absence of spliced reads at position 2087 were predicted by MaxEntScan scores, in which reads were found to map across the 5′ splice sites with strong positive scores (B2 and C2), but not those with strong negative scores (A2 and D3). ② Varying spliced read counts could not be explained by similar scores. ③ Most spliced reads were mapped across a weak splice donor site. See also Fig. S4, Table S8.

Fig. 4.

Expression profiles of HBV splice variants in HBV-transfected Huh7 cells, HBV-infected PHHs and biopsy samples. The heatmaps show the mean percentages of HBV RNAs that were spliced. The RNA-seq libraries that had ≥5 splice variants are shown. The known (SP) and putative (pSP) splice variants were reproducibly detected across the independent biological replicates of HBV-transfected Huh7. Other splice variants are represented with PacBio CCS read names (see Methods). See also Figs S4–S6, Table S1.

Fig. 5.

HBV 5′ splice sites are more likely to be spliced than 3′ splice sites. (a) More spliced reads were mapped across the 5′ splice sites of HBV than 3′ splice sites. Similar results were obtained from Welch two-sample t-test (one-sided) and permutation test (e.g. P values of 0.04 and 0.06 were obtained for genotype B2, respectively). Solid black lines indicate median values. (b) Completeness of splicing at the 5′ and 3′ splice sites. Only the splice sites supported by ≥10 reads were included for comparison.

Fig. 6.

Most frequently used splice sites differed between HBV and the host. The nucleotide frequencies surrounding the splice sites are represented by the spliced reads. Exon boundaries are shaded in grey. Only the splice sites supported by ≥10 reads were included.

Fig. 7.

MA plots [log ratio versus mean expression (log scale)] show differential gene expression between the HBV-treated and control samples. Normalized counts indicate the counts divided by the normalization factors (as computed using the deseq2 default function). Red points denote the FDR-adjusted P value of <0.05. Unfilled triangles denote the genes that have undergone twofold changes in expression. See also Table S9.

Fig. 8.

Significantly dysregulated genes in the HBV-treated cells. A total of 12 genes were differentially expressed in the B2-treated sample. Circled numbers denote the ranking based on FDR-adjusted P values. See also Tables S9 and S10.