Abstract

Background. Hundreds of millions of people worldwide suffer from chronic hepatitis B (CHB), making it one of the most serious public health problems. CHB can lead to serious acute-on-chronic liver failure (ACLF). Since ACLF has a high mortality rate, predicting the risk of ACLF is a critical issue in clinical practice. Methods. To investigate the occurrence of ACLF in CHB, we used high-throughput RNA sequencing to detect the expression of a novel endogenous noncoding RNA (circRNA) in healthy controls (HC) as well as CHB and ACLF patients. Differentially expressed circRNAs were selected and analyzed by Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) biological pathway analyses, circRNA-miRNA coexpression networks, and statistical analyses. Results. A total of 21,101 circRNAs were identified in HC, CHB, and ACLF subjects. Compared with the HCs, 4 circRNAs were upregulated and 14 circRNAs were downregulated in ACLF subjects, with a consistent trend in all three groups. GO analysis revealed that regulation of lymphocyte activation and macrophage tolerance induction were the most important biological processes in ACLF progression. KEGG pathway analysis revealed that primary immunodeficiency and NOD-like receptor signaling pathways were the most enriched terms. The circRNA-miRNA coexpression network and statistical analyses showed that circRNA_07734, circRNA_08533, and circRNA_16083 were the most important circRNAs in the process of ACLF. Conclusions. Immune dysfunction and deficiency may play a key role in the development of ACLF, especially circRNA_07734, circRNA_08533, and circRNA_16083.

1. Introduction

Chronic hepatitis B (CHB) is one of the most serious public health problems worldwide. More than 2 billion people were infected with hepatitis B virus (HBV), 360 million of whom were chronically infected [1]. Although CHB is initially asymptomatic, it can lead to serious acute-on-chronic liver failure (ACLF, HBV-ACLF) after several years, which has a poor prognosis and a high mortality rate (>70%) [2]. ACLF is serious in CHB patients, who constitute approximately 70% of all ACLF cases in most Eastern countries [2, 3]. In China, due to the high prevalence of chronic HBV infection, HBV-associated ACLF patients account for approximately 80% of all ACLF cases [4]. Since ACLF has a very high mortality rate, prediction of ACLF risk is a major issue in the medical field.

At present, the model for end-stage liver disease (MELD) is recognized as a reliable method to measure the mortality risk of patients with liver failure. Shen et al. assessed six widely used short- and long-term prognostic models for patients with ACLF. They concluded that the integrated MELD model may be optimal [5]. Liu et al. suggested that incorporating the γ-glutamyl transpeptidase to platelet ratio into the MELD model may allow more accurate prediction of survival in ACLF patients [6]. Ma et al. established a scoring system offering superior predictive performance, in terms of both specificity and sensitivity, for ACLF patients compared with the MELD [7]. Wei et al. reported that serum Golgi protein 73 (GP73) may be useful in the diagnosis of ACLF in populations with chronic HBV infections [8]. Cui et al. evaluated the prognostic utility of serum Hcy levels, which may serve as a biomarker for 3-month mortality rate in ACLF patients [9]. Although many studies predicting the risk of ACLF have been performed, no method had been proven effective and the pathogenesis of ACLF remains unclear.

Circular RNA (circRNA) is a novel endogenous noncoding RNA characterized as a widespread, abundant, stable, conserved, and tissue-specific molecule in mammalian cells [10]. Different from linear RNA, circRNA is not degraded by RNase R; it is composed of a class of RNA developing covalently closed loop structures without 5-3 polarities [10, 11]. circRNAs are important competing endogenous RNAs. They act as sponges for microRNA (miRNA) and regulate the expression of miRNA-targeted transcripts. Moreover, circRNAs have the potential to serve as biomarkers for the diagnosis and prognosis of several cancers [12, 13]. For example, circRNAs associated with the occurrence and development of hepatocellular carcinoma (HCC) [14], such as circRNA_100338 and circRNA_0005075, serve as biomarkers for HCC and as targets for HCC therapeutics associated with HBV [1517].

The numbers of circular RNAs are different between CHB and ACLF and increase with symptom severity [18]. So far, there is no report on circRNAs in CHB and ACLF, and whether aberrant expression of circRNAs plays a role in ACLF is still unknown. Our study is aimed at investigating the numbers of circRNAs in CHB and ACLF and at identifying the important circRNAs in ACLF development.

2. Material and Method

2.1. Instrument and Database

Total RNA was extracted using QIAamp® RNA Blood Mini (Qiagen, No. 52304). The concentration and purity of RNA were detected on a NanoDrop ND-1000 spectrophotometer (Thermo Fisher Scientific, Inc.). The quality of extracted RNA was detected by using Agilent 4200 TapeStation system (G2991AA) according to the Agilent RNA ScreenTape Assay Quick Guide for 4200 TapeStation System. RNA integrity was evaluated using the Agilent 4200 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA).

The DESeq package (https://ngs.csr.uky.edu/DESeq) , which uses a model based on the negative binomial distribution method, was applied to detect differentially expressed circRNAs and mRNA. DAVID (Database for Annotation, Visualization and Integrated Discovery; http://www.david.abcc.ncifcrf.gov/) was used to analyze the potential functions of the linear transcripts. The homemade miRNA target prediction software, Arraystar (Rockville, MD, USA), is based on miRanda (http://microrna.org/microrna/home.do) and Target Scan (http://www.targetscan.org/vert_71) .

2.2. Ethics Statement and Sample Preparation

The current study was approved by the Medical Ethics Committee of The First Affiliated Hospital, School of Medicine, Zhejiang University (Registry number 201743), China. All clinical examinations and data collection procedures were conducted in accordance with the Declaration of Helsinki.

Three CHB and three ACLF patients admitted to the Liver Center were included. Three healthy controls (HCs) were recruited from the Physical Examination Center. The clinical information is listed in Supplemental table 1.

Blood samples (5 mL) were collected into heparinized tubes and diluted 2 : 1 with phosphate-buffered saline. Mononuclear cells (MNCs) were obtained with equal amounts (1.077 g/cm3) of Ficoll-Paque (Ficoll-Paque™ PLUS; GE Healthcare, Chicago, IL, USA) and density gradient centrifugation. Samples with an RNA 7 were subsequently subjected to analysis.

2.3. circRNA Library Construction and Identification of Differentially Expressed circRNAs

For the circRNA library construction, we used the circRNA analysis software CIRI. First, raw reads were filtered to remove connector sequences, sequences containing more N, and low-quality reads, and clean reads were obtained. Then, by comparing with the ribosomal database, the ribosomal RNA sequences were removed and effective reads were obtained. On the one hand, effective reads were compared with the reference genome by HISAT2 (mapping) and the mapping reads were assembled into mRNA by string to generate the mRNA expression matrix. Effective reads, on the other hand, predicts the back-splicing junction of circRNA through CIRI, The predicted circRNA was obtained, and then, the predicted circRNA was rematched back to the reads of the reverse shearing site to obtain the expression matrix of circRNA. Meanwhile, the circRNA shear ratio could be obtained by the combined analysis of the two aspects.

We use circSeq data to compare and analyze the differential expression of the same circRNA in two samples, and two criteria were selected: First is fold change, that is, the change multiple of the expression level of the same circRNA in two samples. The second is value or FDR (adjusted value). The value of each circRNA is calculated by FDR error control method, and then, the value is corrected by multiple hypothesis testing. The screening condition of differentially expressed circRNAs is and .

2.4. Identification of circRNAs and mRNA

The concentration and purity of total RNA extracted from MNCs were measured on a NanoDrop ND 1000 spectrophotometer (Thermo Fisher Scientific, Inc.) at 260 and 280 nm wavelengths. Before the synthesis of first- and second-strand cDNA, Ribo-Zero (ribosomal RNA) was depleted and RNA-fragmented. cDNA was purified followed by processing of adenylate 3 ends, ligation of adapters, and enrichment of DNA fragments. Finally, 1 μL was loaded onto a 2100 Bioanalyzer (Agilent Technologies) for validation. High-throughput RNA sequencing was performed according to the method of Memczak et al. [19].

Raw data were extracted using Agilent Feature Extraction software and predicted by circBase. circRNAs and mRNA were quantitatively analyzed by Shanghai OE Biotech Ltd., Co. (Shanghai, China). The circRNAs and mRNA of the HC, CHB, and ACLF subjects were verified with CIRI software (Shanghai OE Biotech), and the identified sequence was selected for further analysis. The expression levels of circRNAs among the three groups were measured by mapped back-splice junction reads per million mapped reads (RPM), and the mRNA_expression was measured by FPKM (Fragments Per Kilobase per Million). The location of the chromosome where the circRNA sequence overlapped was annotated. The criteria were chosen according to the fold change and value of the circRNA and mRNA in the same sample between two groups.

2.5. Pathway Analyses and miRNA Target Prediction

Gene functions were classified into three subgroups: biological process (BP), cellular component (CC), and molecular function (MF). The most enriched GO terms among the CHB and ACLF groups (ranked by enrichment score) were selected, and the enriched GO terms between HC and ACLF were ranked by the number of differentially expressed linear transcripts. The different biological pathways of linear transcripts were analyzed by the KEGG pathway.

circRNA can be used as a target molecule of miRNA and regulated by miRNA. Analysis of the circRNA-miRNA interaction could help elucidate the function and mechanism of circRNAs in contact with miRNAs. The circRNA-miRNA interactions were predicted by Arraystar. The target interaction network diagram of the circRNA–miRNA/gene interaction network was built by Cytoscape (available on http://www.cytoscape.org/).

Spearman’s correlation was performed to evaluate the relation between circRNA with mRNA and liver function indexes. All statistical analyses were conducted with SPSS software (ver. 22.0; IBM Corp., Armonk, NY, USA), and was considered statistically significant.

3. Results

3.1. Overview of circRNA Profiles

To investigate the expression profile of circRNAs involved in HC, CHB, and ACLF, tissues of these three groups were analyzed by circular RNA high-throughput sequencing. Results showed that approximately 86.61% of circRNAs had a predicted spliced length of less than 2,000 nt; 48.44% had a length less than 500 nt, and 26.53% had a length between 500 nt and 1,000 nt (Figure 1(a)). Moreover, Venn diagram of circRNA expression showed that a total of 21,101 circRNAs in the HC, CHB, and ACLF were detected. Among the 21,101 circRNAs, 3581 circRNAs were differentially expressed in CHB tissues compared with HC and ACLF; 3214 circRNAs were differentially expressed in ACLF tissues compared with HC and ACLF (Figure 1(b)). There are no outliers of circRNA expression in the HC, CHB, or ACLF group due to instrument or operational anomalies. The distribution of circRNAs on the genome is shown in Supplemental Figure 1.

3.2. Identification of Differentially Expressed circRNAs and mRNA

According to the circRNA criteria ( and ), 550 differentially expressed circRNAs in the CHB and ACLF groups were identified compared with HC. Among them, 176 differentially expressed circRNAs were detected in the CHB group (89 upregulated and 87 downregulated) and 374 differentially expressed circRNAs were detected in the ACLF group (181 upregulated and 193 downregulated) compared with HC. The differentially expressed circRNAs between ACLF and HC groups are shown in a volcano plot (Figure 2(a)). The unsupervised hierarchical cluster analysis results are shown in Figure 2(b).

The differentially expressed candidate circRNAs among the three groups showing progressive up- or downregulation were further analyzed. Compared with HC, there were 4 upregulated and 14 downregulated circRNAs (Table 1, Figure 3) in ACLF. mRNAs were analyzed by the same method. The results showed that there were 7 upregulated and 9 downregulated mRNA in ACLF compared with HC. The expression information and statistic differences of circRNA and mRNA are shown in the Supplemental Table 2, 3.

3.3. GO and KEGG Analyses

The functions of the target genes of differentially expressed circRNAs in the ACLF group were further analyzed by GO enrichment and KEGG pathway analysis. Significant terms in the GO and KEGG pathway analyses () were selected and ranked by value.

After obtaining the differentially expressed circRNA, we conducted GO enrichment analysis for the differentially expressed circRNA using the information of circRNA source genes. The number of differential transcripts included in each GO entry was counted, and the significance of differential circRNA enrichment in each GO entry was calculated using hypergeometric distribution test. The results of calculation will return a value of enrichment significance, small value indicates enrichment of circRNA in this GO item, and its calculation formula is as follows:

Enrichment score is calculated as follows:

is the number of circRNAs with GO annotation in all circRNAs; is the number of circRNAs with GO annotation in source genes of differentially expressed circRNAs in ; is the number of circRNAs whose source genes are annotated as a specific GO term among all circRNAs; is the number of differentially expressed circRNAs annotated by the source gene for a specific GO term. The biological functions of circRNA-derived genes can be understood according to the results of GO analysis combined with the biological significance.

The top-10 downregulated GO terms, classified by BP, CC, and MF, are shown in Supplemental Figure 2A. The most enriched BP, CC, and MF terms were regulation of lymphocyte activation, integrin alphaL-beta2 complex, and RNA polymerase I core binding, respectively. Among the upregulated GO terms, the most enriched BP, CC, and MF terms were positive regulation of macrophage tolerance induction, Golgi lumen, and extracellular matrix structural constituent, respectively (Supplemental Figure 2B).

Identified circRNAs were annotated using the KEGG public database. Pathway analysis was carried out on differentially circRNA-derived genes (combined with KEGG annotation results) using the KEGG database, and the significance of differentially circRNA enrichment in each pathway item was calculated using the hypergeometric distribution test. The results of the calculation returned a value of enrichment significance, and a small value indicated that the differential transcript was enriched in this pathway. Refer to GO enrichment analysis for the corresponding calculation formula.

KEGG pathway analysis (Table 2, Supplemental Figure 3) revealed that the most enriched term in downregulated circRNAs was primary immunodeficiency, while in upregulated circRNAs, NOD-like receptor signaling pathway was the most enriched term. As shown in Table 2, the KEGG pathway hsa05340 was the most important.

3.4. circRNA-miRNA Coexpression Network

A network map of the most significant relationships between circRNAs and miRNAs is shown in Figure 4. Among the differentially expressed circRNAs selected from the ACLF group, circRNA_03646, circRNA_06400, circRNA_09331, circRNA_17523, circRNA_05781, circRNA_07734, circRNA_08533, circRNA_17526, circRNA_15699, and circRNA_16083 underpinned the circRNA-miRNA coexpression network, as determined by Arraystar. Among them, thee expression of circRNA_03646, circRNA_06400, circRNA_07734, circRNA_08533, and circRNA_16083 had statistical significance between the ACLF and HC group. Moreover, circRNA_07734, circRNA_08533, and circRNA_16083 also had statistical significance in KEGG analysis. circRNA_07734 only contacted two miRNAs in the network, which is less than circRNA_08533 and circRNA_16083.

3.5. Correlation Analysis

Spearman’s correlation analysis indicated that the RPM values ofcircRNA_07734, circRNA_08533, and circRNA_16083 were related to ALT, AST, TBIL, and IBIL (Table 3). circRNA_07734 and circRNA_08533 showed strong negative correlation with ALT, AST, TBIL, and IBIL, which indicated if ALT, AST, TBIL, and IBIL increased higher, the expression of circRNA_07734 and circRNA_08533 would be lower. circRNA_16083 is different with circRNA_07734 and circRNA_08533, which had positive correlation with ALT, AST, TBIL, and IBIL. In clinical practice, ALT, AST, TBIL, and IBIL are used to evaluate the liver function. Commonly, the higher the ALT, AST, TBIL, and IBIL levels were, the more severe the liver injury will be. Therefore, circRNA_07734, circRNA_08533, and circRNA_16083 may be useful as the sensitive indicators of liver injury.

4. Discussion

circRNAs have received increasing attention due to their highly stable and specific spatiotemporal expression patterns. circRNAs are a transcriptional product in various tissue and cell types in humans, mice, drosophila, and nematodes, among other species. circRNAs can be classified into three types: exonic (circRNA arising from the exons of the linear transcript), intergenic (circRNA located outside of a known gene locus), and sense-overlapping (circRNA whose gene locus overlaps with linear RNA but is transcribed from the opposite strand). The differentially expressed hosting genes of the circRNAs and the disordered splicing machinery will lead to circRNAs differentially expressed. So far, numerous circRNAs have been reported in different animals. The majority of circRNAs showing specific expression patterns are related to tissue development and disease [19]. Therefore, circRNA has potential as a biomarker for some diseases [11]. Thus far, many studies have reported abnormal expression levels of circRNAs in a number of cancers, such as gastric [20, 21], colorectal [22, 23], breast [24, 25], and lung cancers [26, 27]. However, few studies have focused on the development of ACLF.

In this study, we detected several circRNAs in HC, CHB, and ACLF subjects. The circRNAs were initially identified by reference to a database or website. circBase, (http://circbase.org/) which includes merged genomic data on circRNAs from humans, mice, nematode worms, pike, and coelacanth, provides evidence of their expression [28]. Therefore, we identified the circRNAs detected herein using circBase. circRNA-seq showed that there were 4 upregulated and 14 downregulated circRNAs in ACLF compared with HC. Dysregulation of circRNAs in patients with chronic hepatitis B prompted us to investigate the functional role of these circRNAs.

CHB is mainly caused by aberrant cellular immunity, due, for example, to natural killer (NK) cells, cytotoxic T cells, and antibody-dependent lymphocytes. When infected with hepatitis B, three types of lymphocytes target antigens of the hepatocyte membrane: HBsAg, HBcAg, and other immunoregulatory cells. Auxiliary and inhibitory T cells take part in the immune response, which can lead to immune disorder hypofunction. Many researchers have reported dysfunction and inhibition of T cells in CHB. For example, Raziorrouh et al. [29] reported that tetramer+CD4+ T cell numbers were reduced during CHB. Zhang et al. reported that HBV-specific T cell responses to recombinant HBV core protein were reduced in CHB patients and positively correlated with HBV viral load in coinfected, chronic HBV patients [30]. Therefore, the progression of ACLF from CHB may be related to the BP of immunodeficiency or dysfunction, which can in turn be caused by HTLV-I infection, Staphylococcus aureus infection, or other diseases. The process of immunity dysfunction may be related to the NOD-like receptor signaling pathway, the NF-kappa B signaling pathway, and/or the interleukin-1-mediated signaling pathway.

Since circRNAs contain multiple miRNA-binding sites, the function of individual circRNAs can be predicted by their miRNA target gene and annotated according to the function of their miRNA target gene. The circRNA-miRNA coexpression network identified two upregulated circRNAs and nine downregulated circRNAs in this study. After combining the results from the GO and KEGG analyses, circRNA_00389, circRNA_07734, circRNA_08533, circRNA_ 15699, and circRNA_16083 were identified as the most important circRNAs during the progression of ACLF. Among them, circRNA_00389 had the highest enrichment score in the KEGG analysis, consistent with the GO analysis. The top term in BP analysis was regulation of lymphocyte activation, which corresponded to downregulated circRNA_00389. However, the expression of circRNA_00389 had no statistical significance between the ACLF and HC group. Therefore, combined with the results of Spearman’s correlation analysis, circRNA_07734, circRNA_08533, and circRNA_16083, which are related to immune dysfunction and correlate with ALT, AST, TBIL, and IBIL, are important circRNAs in ACLF.

This study was conducted to initially explore the role of circRNAs in CHB and acute liver failure due to CHB. The etiology of ALCF is CHB. Compared with healthy subjects, many circRNAs changed in CHB and ACLF. We speculated that circRNAs was an intermediate regulatory molecule that alters the activity of downstream miRNA and thus regulates mRNA levels of some genes. Therefore, we did not find significant changes in circRNAs between the CHB and ACLF groups.

The limitation of this study is that there were only three samples in each group, but the selected circRNAs all had the same expression tendency in CHB and ACLF (Supplemental Table 1). The results of test analysis between groups were consistent with the expression tendency of selected circRNAs. Besides, the Spearman correlation analysis of circRNA_07734, circRNA_08533, and circRNA_16083 showed that there was statistical significance with differentially expressed mRNA and ALT, AST, TBIL, and IBIL (Table 3). Therefore, our results were credible and can be used as the basis for further studies on the regulatory role of circRNA downstream target genes in CHB and ACLF. In further work, we should perform qRT-PCR to determine the expression levels of these circRNAs and carry out research on mechanisms.

5. Conclusion

In this work, 21,101 circRNAs were detected in CHB, HC, and ACLF groups. In total, 4 upregulated and 14 downregulated circRNAs were identified in the ACLF group compared with the HC group, which showed a consistent trend of change in the three groups. GO and KEGG biological pathway analyses and circRNA-miRNA coexpression network and statistical analyses revealed that circRNA_07734, circRNA_08533, and circRNA_16083 are the most important circRNAs in ACLF. Their RPM values correlated with ALT, AST, TBIL, and IBIL; thus, immune dysfunction or immunodeficiency may be critical for ACLF development.

Data Availability

The datasets generated and/or analyzed during the current study are not publicly available but are available from the corresponding author on reasonable request.

Disclosure

An earlier version, named as “Circular RNA expression profiles in the development of acute-on-chronic liver failure in hepatitis B virus-infected patients,” has been presented as preprint in http://www.researchsquare.com, in the following link: https://www.researchsquare.com/article/rs-38521/v3. The funder had no role in the study design, data collection, analysis, decision to publish, or manuscript preparation.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Authors’ Contributions

Jiangshan Lian and Feiyan Lin contributed equally to this work.

Acknowledgments

This work was supported by the opening foundation of the State Key Laboratory for Diagnosis and Treatment of Infectious Diseases, The First Affiliated Hospital, College of Medicine, Zhejiang University, grant No. SKLID2020KF01, and the Scientific Research Projects of Wenzhou Municipal Science and Technology Bureau (Y20180106).

Supplementary Materials

The data in the supplementary tables provide the clinical parameters in the HC, CHB, and ACLF groups and the expression of selected circRNA (RPM) and mRNA (FPKM) in three groups. The supplementary figures show the distribution of circRNAs on each chromosome. Gene Oncology (GO) analysis of ACLF risk status according to the down- and upregulated circRNAs. KEGG enrichment analysis of downregulated circRNAs and upregulated circRNAs (Supplementary Materials)