Identification and characterization of circular RNAs expression profiles in obstructive sleep apnea-induced liver injury

Animals

Six-week-old male balb/c mice (17–21 g) were purchased from Beijing Weitong Lihua Experimental Animal Technology Co., Ltd. All mice were housed with standard mouse diet and tap water. The animal protocol of our study was approved by the Experimental Animal Ethics Committee of the Second Affiliated Hospital of Fujian Medical University (QZ5312408).

Liver CIH injury protocol and circRNAs sequencing

An CIH model was established by chronic intermittent hypoxia (IH) as previously described by us [3, 10]. Briefly, we used the intermittent hypoxia system to create a CIH environment. Specifically, the oxygen and nitrogen were automatically flowed into the chamber by a gas control system and calculated with an oxygen analyzer, creating an intermittent hypoxia condition. After the chamber reached 6% O₂ for 1 min with a 2 cycle (1 min/cycle), compressed air was pumped into the chamber to achieve 21% O₂ for another 1 min. Six mice were randomly allocated into two groups (n = 3, each group): (1) an CIH group, in which mice were placed daily in the chamber for 30 cycles/h, 8 h/day for 7 days/week, for eight consecutive weeks; (2) a control group, in which mice were housed in the chamber with 21% O₂ during the entire experiment. Following pentobarbital euthanasia, the liver tissues of mice were collected at the end of CIH exposure. All mice experiments were performed in duplicate and repeated at least three times.

CircRNAs microarray hybridization and circRNAs microarray analysis

The general flowchart of data processing and detailed methods were described in Figure 1. Briefly, the process detail of circRNA microarray hybridization can be found in our previous work [3, 7]. Agilent Microarray Scanner was used to scan the arrays and acquired array images analyzed using Agilent Feature Extraction Software 9.5.1.1. Quantile normalization and subsequent data processing were performed using R project and the Bioconductor package of Limma [11, 12].

Figure 1. The flowchart in this study.

Changes in circRNAs expression with p-values < 0.05 and fold-changes ≥1.5 were considered statistically significant. Fold change filtering was calculated for screening DECs, which were visualized with Volcano Plot. Hierarchical clustering was performed to display the distinguishable DECs patterns between samples.

Bioinformatic analyses

Gene Ontology (GO) enrichment analysis was performed to determine functional annotations of DECs-target genes, including biological processes, cellular components and molecular functions, using the clusterProfiler R package (v4.0.0) [13, 14]. Furthermore, Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis was also further performed to identify the enriched pathways among these DECs-target genes (Kanehisa and Goto, 2000). The threshold of significance was defined as P ≤ 0.05 for both GO and KEGG analyses. Moreover, we ranked the top 10 enriched GO categories and the top 10 enriched KEGG pathways of the DECs-target genes.

Validation by qRT-PCR

We used qRT-PCR to validate the differential expression level of six selected circRNAs. The total RNA from liver tissues was extracted by using the TRIzol Reagent (Takara, Dalian, China). All primers, spanning the distal ends of circRNAs, were designed using Primer 5 software (Table 1). Briefly, cDNA was synthesized using the PrimeScript™ RT Reagent Kit (Takara, China). Subsequently, the qRT-PCRs were run and analyzed on an ABI Q2Real-time PCR system (Applied Biosystems, USA) using a TB Green™ Premix Ex Taq™ II (Takara, China). Finally, relative circRNAs expression were calculated using the 2^−ΔΔCt method, with β-actin serving as an internal control.

Construction of circRNA-miRNA-mRNA networks

Based on the database of miRanda and TargetScan, we randomly selected four significant DECs to construct a competitive network among circRNA, miRNA and mRNA by a software (Arraystar’s home-made miRNA target prediction software). In addition, the binding capacity of circRNA and microRNA, as well as the capacity and number of microRNA-mRNA binding sites should be considered when we used the software. Finally, we estimated a circRNA-miRNA-mRNA network to analyze the ceRNA manner of the circRNAs and the network was visualized with Cytoscape software (Version 3.7.2) [15].

Statistical analysis

Each measurement was independently repeated three times. Continuous data were expressed as mean ± standard deviation (SD). All qPCR statistical analyses using unpaired student’s t-test. P < 0.05 were considered marginally significant.

CIH treatment affected liver histology

Normal hepatic architecture was observed in the liver tissues of NC group (Figure 2A, 2B). As shown in Figure 2C, 2D, CIH caused liver injury evidenced by some liver cells slightly edematous with unclear borders and a few hepatocyte nuclei slightly enlarged infiltrated (Figure 2C, 2D).

Figure 2. Histopathological alterations involving liver tissues. (A, B) Representative images of HE-stained liver tissues of the NC group; (C, D) Representative images of HE-stained liver tissues of the CIH group. NC group had normal hepatic architecture; CIH group showed some liver cells were slightly edematous with unclear borders. A few hepatocyte nuclei were slightly enlarged. Some cells were in sheets, and the demarcation of the hepatic cord was unclear. All images H&E with original magnification of 200X. NC, normal control; CIH, chronic intermittent hypoxia.

Overview of circRNAs expression

CircRNA microarray was conducted to screen the dysregulated circRNAs from 3 CIH-liver tissues and 3 normal liver tissues from the control group. With the above p-value and fold-change threshold, we detected a total of 80 DECs. Of these DECs, 1 were upregulated, while 79 were downregulated. The results of the quantile regression in Figure 3A show that the circRNA median of the six samples were normalized by quantile normalization (Figure 3A). Importantly, hierarchical clustering, scatter plots and a volcano plot efficiently show that clearly DECs expression profiles in two groups (Figure 3B–3D). To classify the dysregulated circRNAs in CIH-liver tissue samples, we listed the circRNAs into five categories, including 5567 exonic, 1026 sense-overlapping, 430 intronic, 240 intergenic regions and 140 antisense of the 7403 upregulated circRNAs, while the 6588 downregulated circRNAs comprised 5456 exonic, 733 sense-overlapping, 258 intronic, 74 intergenic and 66 antisense (Figure 4A). In addition, chromosomal distribution analysis revealed that most circRNAs were located at chromosome 1,2,4,5,9,11, while few located at the Y chromosome (Figure 4B).

Figure 3. Overview of the microarray signatures. (A) Box plots show the normalized intensities from the two groups. (B) Hierarchical clustering of all DECs with a fold-change ≥1.5 and p-value < 0.05. (C) CircRNAs in the Scatter-Plot located above the top green line and below the bottom green line indicated more than a 1.5-fold change of circRNAs. (D) Volcano plot filtering revealed that 80 DECs, including 1 up-regulated and 79 down-regulated circRNAs. Abbreviation: DECs: differentially expressed circRNAs.

Figure 4. Classification and distribution of the DECs. (A) Classification of the DECs was listed. (B) The number of DECs was determined in each chromosome. DECs, differentially expressed circRNAs.

Validation of circRNA expression

To verify the veracity of microarray results, we conclusively validated the selected circRNAs by qRT-PCR. There are 5 downregulated circRNAs (mmu_circRNA_000469, 37851, 38959, 38983, 31665) and 1 upregulated circRNAs (mmu_circRNA_27565) between the CIH group and the control group. These six circRNAs are selected according to a comprehensive evaluation of the p-values, fold-changes, length and original signal values of the circRNA. The qRT-PCR data were in agreement with the above microarray results (Figure 5).

Figure 5. Relative fold changes of six randomly selected circRNAs from microarray results and validated by qRT-PCR. The downward and upward histogram represent down-expression and up-expression, respectively.

GO and KEGG analyses

To explore the functional roles of circRNAs target genes, GO term enrichment and KEGG pathway analyses of these validated circRNAs related genes were performed. Specifically, the most significant GO functions were related to carboxylic acid metabolic process, negative regulation of mitochondrial membrane potential, small molecule metabolic process, cellular anatomical entity, organelle, intracellular, catalytic activity, acting on a tRNA, histone deacetylase and binding magnesium ion binding (Figure 6A–6C). Additionally, KEGG analysis identified numerous signaling pathways, associated with the dysregulated circRNAs-related genes, that were significantly altered in the OSA-induced liver injury. The KEGG pathways were significantly associated with Natural killer cell mediated cytotoxicity, one carbon pool by folate, Glycerolipid metabolism, HIF-1 signaling pathway, Renal cell carcinoma, cGMP-PKG signaling pathway, Central carbon metabolism in cancer, Thyroid hormone signaling pathway, Tuberculosis, β-Alanine metabolism (Figure 6D).

Figure 6. GO and KEGG Pathway analysis of validated circRNAs related genes. (A–C) Top 10 enriched GO terms in BP, CC and MF. (D) Top 3 enriched KEGG pathways. Abbreviations: GO: Gene Ontology; CC: cellular components; BP: biological processes; MF: molecular functions.

Construction of a circRNA-miRNA-mRNA interaction network

For predictions of potential targets of the selected circRNAs, miRWalk, miRanda and Targetscan were used. To further explore the biological functions of circRNAs in OSA-induced liver injury, we constructed a lncRNA-miRNA-mRNA network, including four selected circRNAs (mmu_circRNA_000469, 38959, 38983, 27565), and their targeted miRNAs and downstream mRNAs. And the network was visualized by Cytoscape software (version 3.6.1) (Figure 7). There are 227 target miRNAs and corresponding 321 target mRNAs of the circRNAs by bioinformatics predictions. The interaction network will help us improve the understanding of underlying mechanism of mmu_circRNA_000469, 38959, 38983, 27565 by revealing its potential connections between circRNAs and miRNAs. Our finding suggested that these circRNAs may play a critical role in the biological process of OSA-induced liver injury.

Figure 7. A network diagram shows the circRNA-miRNA-mRNA network for four candidate circRNAs. Yellow rounds, red rounds and green rounds represent circRNA, miRNA and mRNA, respectively.