Integrated bioinformatics analysis for the identification of hub genes and signaling pathways related to circANRIL

Cell culture and construction of overexpression plasmid

The human cell line HEK-293T was purchased from the Cell Bank, Chinese Academy of Sciences (http://www.cellbank.org.cn) and the cells were cultured in high-glucose Dulbecco’s modified Eagle medium (DMEM) supplemented with 10% fetal bovine serum (FBS), 100 U/mL penicillin and 100 U/mL streptomycin, in a 37 °C, 5% CO₂ and saturated humidity incubator. To generate the circANRIL- OE cell line, sequcences from the circRNA hsa_circ_0008574 were obtained from http://www.circbase.org/, we constructed an expression plasmid using the following sequence: CGGAATTC TGAAATATGCTATCTTACAGTGTCCCTTTTGATGAGAAGAATAAGCCT CATTCTGATTCAACAGCAGAGATCAAAGAAAAGACTTCTGTTTTCTGGCCACCAGATAATGTTATCTGTGCTTAAAGAATTGAAAAACACACATCAAAGGAGAATTTTCTTGGAAA GAGAGGGTTCAAGCATCACTGTTAGGTGTGCTGGAATCCTTTCCCGAGTCA GCTTTCTAGAAGAAAACCGGGGAGATCTATTTGGAATGTATCTAACTCCAAAGAA ACCATCAGAGGTAACAGTAGAGACGGGGTTTCACCATGTTGGCCAGACTGGTCTT GAACTCCCGACCTCGTGATTCGCCCGCCTCGGCCTCCCAAAGTGCTGGGATTA CAGGTGTGAGACACCACGCCCGGCGGATAGAGAGAATTTTGACAGGTGAATATA TTTTTTCTTGAGGATCCCG. The restriction sites BamHI and EcoRI are underlined, and the splice sites were designed as donor sites (GT on the forward strand) and acceptor sites (AG on the reverse circularized sequence) which are shown in bold and italics, the other sites were target sequences. After digestion by the restriction enzymes EcoRI and BamHI, the PCR fragments were sequenced before cloning into the pUC57-Amp vector (Genewiz, suzhou, China) to generate the pUC57-hsa_circ_0008574 construct, which was transformed into Stbl3 competent cells (Invitrogen, Waltham, MA, USA). The status of individual colonies selected from plates was also confirmed by sequencing.

Verification and analysis of the stable cell line

Lentivirus particles were generated by cotransfection of 4.4 µg pUC57-hsa_circ_0008574, 2.4 µg psPAX2 and 2 µg pMD2.G into HEK 293T cells with Lipofectamine™ 2000 (Invitrogen).The infection efficiency was determined under a fluorescence microscope after 48 h. Stably transfected cell lines were selected using puromycin (TaKaRa, Japan) and the cells were divided into three groups: circANRIL-OE group (transfected with the constructs), blank control group (untransfected cells) and negative control group (circRNA control, transfected with the empty vector). Each group contained three samples, and each sample contained approximately 1 × 10⁶ cells.

To validate the specificity of the HEK-293T cell transduction, total RNA was extracted from HEK-293T cells by using TRIzol reagent (TaKaRa, Japan) according the manufacturer’s protocol. mircoRNA (mRNA) levels of circANRIL and lincANRIL in these groups were analyzed with a qRT-PCR kit (TianGen, Beijing, China) following the manufacturer’s instruction.

Identification of DEGs

According to the manufacturer’s protocol, a cDNA library was established by using the RNA Sequencing Library Preparation Kit (Illumina). The libraries were applied to a flowcell on a cBOT instrument (Illumina) to generate clusters, and then subjected to Illumina high-throughput RNA sequencing platform (Illumina, San Diego, CA, USA). To obtain clean data, the sequencing data were subjected to quality control, and the clean reads were mapped to the human reference genome (GRCh38) using TopHat (v2.1.0), The R function cor was utilized to obtain the Pearson correlation coefficient (Pearson’s r), the correlation coefficient ranged from 0 to 1, and higher scores indicate highly similar expression patterns among samples. Principal component analysis (PCA) was performed by R function prcomp, and the results were visualized using the ggfortify package version: 0.4.5 (https://cran.r-project.org/bin/windows/contrib/3.4/ggfortify_0.4.6.zip).

Differential expression analysis was performed using the edgeR language package version 3.4 (http://www.bioconductor.org/packages/release/bioc/html/edgeR.html) to screen for DEGs between the circANRIL-OE group and control groups. Raw counts were normalized with the TMM method, and then presented as log2 counts per million (LogCPM) values. We obtained the logFC and P value of each gene after assessing the data, and the adjusted P value (false discovery rate, FDR) was determined after multiple test correction with the Benjamini & Hochberg method, FDR < 0.05 and —logFC— > 1 were used as thresholds for DEG identification.

Gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) analyses of DEGs

All the DEGs were subjected to GO and KEGG pathway enrichment analyses. The analyses were performed using Database for Annotation, Visualization, and Integrated Discovery (DAVID, version 6.8, https://david-d.ncifcrf.gov/) functional annotation tools. GO biological process (BP) analysis was included in this study to explore the key signaling pathways. GO BP terms and KEGG pathway terms with a P value of < 0.05 and an enrichment score of >5 were defined as significantly enriched. Default DAVID database setting with medium stringency and Homo sapiens was defined as background.

Protein–protein interaction (PPI) network

To construct the PPI network, the DEGs were input for PPI analysis into the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database (version 10.0, http://www.string-db.org/). The species selected was “Homo sapiens”. The PPI score was set at 0.9 (highest confidence), Cytoscape (version 3.4.0, http://chianti.ucsd.edu/cytoscape-3.4.0/) was used to determine the PPI network pairs. The key topological properties of the PPI network were analyzed with the Cytoscape plug-in Network Analyzer (CytoNCA, version 2.1.6, http://apps.cytoscape.org/apps/cytonca) in Cytoscape (version 3.4.0). The parameter was set to ‘without weight’. All PPI pairs with a combined score >0.9 were extracted,Three parameters were chosen to evaluate every node in the network, including degree (the number of nodes directly linked to one node), betweenness (the relative importance of a vertex within the network) and closeness (the proximity of a network node to other nodes). the final combined scores were used to sequence and evaluate central proteins.

Construction of the ceRNA network and hub gene selection

To explore the underlying mechanism of circANRIL, the online database circInteractome (https://circinteractome.nia.nih.gov/) was used to identify RNA- binding proteins (RBPs) related to circANRIL and to determine whether any RBPs differentially expressed between the two groups. The TargetScan database (http://www.targetscan.org) was used to predict miRNAs within the 90th percentile. We then predicted the target genes of these miRNAs with miRWalk2.0 (http://mirwalk.umm.uni-heidelberg.de/), and the results of the following 6 databases were comprehensively considered to confirm the target genes: miRWalk (http://mirwalk.umm.uni-heidelberg.de/), MicroT4 4.0 (version 5; http://diana.imis.athena-innovation.gr), miRanda (https://anaconda.org/bioconda/miranda), miRDB (http://mirdb.org/miRDB/), PITA version 6 (https://web.archive.org/web/20200529050843/https://omictools.com/pita-tool) and RNA22 databases (https://cm.jefferson.edu/rna22/). Only those miRNA target genes included in all the gene sets were used for further analysis. Here, we reasoned that if circANRIL is upregulated, the same miRNA-regulated target genes should also be upregulated. Therefore, the potential miRNA-target gene pairs were obtained by overlapping the previously predicted target genes of the PPI and the upregulated DEGs. The circRNA-miRNA and miRNA-mRNA were combined for the construction of the circRNA-miRNA-mRNA regulatory network using Cytoscape 3.4.0. miRNAs and mRNAs with FDR <1%, absolute log2FC >2 and P < 0.01 were retained to establish the ceRNA network. Hub genes were selected according to the node degrees of the ceRNA network by calculating the number of circRNA-miRNA and miRNA-mRNA pairs. qRT-PCR was used for further validation of the selected genes. The overall procedure is summarized in the flow diagram in Fig. 1.

Statistical analysis

Data in each group were assessed using SPSS 22.0 software (SPSS Inc.) and are presented as the mean ± SD. Differences between two groups were analyzed by the chi-square test and the Mann−Whitney U test. Comparisons between multiple groups were carried out using one-way analysis of variance or the Kruskal −Wallis H-test. P < 0.05 was considered statistically significant. Graphs were generated using GraphPad Prism 5 (GraphPad Software, San Diego, CA). All experiments were performed at least three times.

Generation of a cell line stably overexpressing circANRIL

The pUC57-hsa_circ_0008574 recombinant expression vector was confirmed by sequencing analysis (Fig. 2). HEK-293T cells stably transfected with the pUC57-hsa_circ_0008574 recombinant plasmid were identified by green fluorescence (Figs. 3A, 3B). and qRT-PCR was used to determine the circANRIL overexpression efficiency in the corresponding cell line (denoted as circANRIL-OE). The circANRIL-OE group had a significantly higher circANRIL expression level than the negative and blank control group ( P < 0.05), while linear ANRIL showed a trend of decreased expression in the circANRIL-OE group compared to the control groups (Figs. 3C, 3D, Fig. S3-sheet 1).

Functional analysis of the significant DEGs

Sequencing reads from an mRNA library constructed from total RNA from the two groups of cells were mapped to the human reference genome(GRCH38) by TopHat. GENCODE annotation of the human genes was performed by htseqcount (http://htseq.readthedocs.io/en/release_0.9.1/). After filtering out low quality reads, a total of 20,313 genes were further analyzed. All the raw sequencing data have been deposited into the SRA (accession number: PRJNA756043). The Pearson correlations(r) which ranged from 0.921 to 1, and P values are shown in (Fig. 4A). The median sample similarity between groups and within groups was 0.990 and 0.954, respectively, which indicated a significant difference between the circANRIL-OE and control groups. PCA analysis revealed significant separation between the two groups, two primary components (PC1 and PC2) that explained 91.74% of the total variance (PC1 = 59.73%, PC2 = 32.01%) (Fig. 4B). According to the threshold, in total, 1745 DEGs were identified, 1052 DEGs of which were upregulated and 693 of which were downregulated. The list of DEGs between the circANRIL-OE and control groups was subjected to hierarchical clustering. A volcano plot and a heat map of the DEGs are shown in Figs. 4C, 4D (Fig. S4-sheet 2).

GO BP and KEGG analyses of the DEGs

In enrichment analysis, 26 GO BP terms were significantly enriched (P < 0.05) by the upregulated DEGs (P < 0.05) and 127 terms were enriched by the downregulated DEGs.The top 20 enriched GO BP terms are shown in Fig. 5. The results showed that the upregulated DEGs were significantly related to the following terms: RNA processing (GO: 0006396), regulation of cell proliferation (GO: 0042127), ncRNA metabolic process (GO: 0034660), ribonucleoprotein complex biogenesis (GO: 0022613), ncRNA processing (GO: 0034470), ribosome (GO: 0042254), RNA metabolic process (GO: 0016072), RNA processing (GO: 0006364), response to steroid hormone stimulus (GO: 0048545). The downregulated DEGs were predominantly associated with signal transduction (GO: 0007165), positive regulation of GTPase activity (GO: 0043547), cell adhesion (GO: 0007155), multicellular organism development(GO: 0007275), positive regulation of cell proliferation (GO: 0008284), nervous system development (GO: 0007399), negative regulation of cell proliferation (GO: 0008285), positive regulation of Apoptotic process (GO: 0043065), inflammatory response(GO: 0006954), and immune response(GO: 0006955). KEGG pathway analysis revealed that the upregulated DEGs were enriched in the MAPKsignaling pathway (hsa04010), while the downregulated genes were significantly enriched in terms 32 pathways in cancer (hsa05200), the PI3K-Akt signaling pathway (hsa04151), neuroactive ligand–receptor interaction (hsa0408), proteoglycans in cancer (hsa05205), HTLV-I infection (hsa05166), Rap1 signaling pathway (hsa04015), transcriptional misregulation in cancer (hsa05202), influenza A (hsa05164), cytokine-cytokine receptor interaction (hsa04060) and axon guidance (hsa04360) (Fig. 5; Fig. S5-sheet 3).

PPI network analysis and hub gene selection

According to previously described methods, a PPI network including 548 nodes,179 upregulated and 369 downregulated nodes was generated (Fig. 6; Fig. S6-sheet 4). The top 20 proteins ranked from the most to the least significant P values for degree, betweenness, and closeness are shown in Table 1. The results identified PSMB8, CXCL12, FN1, CXCL8, and THBS1 as potential core proteins in the PPI network because all three of their scores ranked in the top 20, However, these genes were excluded from further study because all of them were downregulated proteins.

The RBPs of 13 circRNAs were identified via circBase (http://www.circbase.org/) (Table 2). The search yielded six flanking region binding sites (corresponding to EIF4A3, AUF1, U2AF65, FUS, TIAL1, and PTB), but the corresponding genes were not differentially expressed between the groups.

Moreover, the miRNA targets of the 13 circRNAs of ANRIL were predicted using the TargetScan database (http://www.targetscan.org). A total of 197 circRNA-miRNA pairs including 13 circRNAs and 56 miRNAs were obtained, and 2456 miRNA-targets of these 56 miRNAs were predicted by miRWalk2.0 and six other databases. Based on the previously mentioned methods, the overlapping genes of these predicted mRNAs and genes encoding the upregulated proteins of the PPI network were determined for 36 miRNA-target pairs, including 22 miRNAs and 20 target genes. Furthermore, 71 circRNA-miRNA-mRNA pairs were obtained by identifying miRNAs and target genes regulated by the same miRNAs through intergration of the circRNA-miRNA and miRNA-mRNA relationships. The circRNA-miRNA-mRNA network contained 31 circRNA-miRNA pairs, 69 miRNA-mRNA coexpression relationships, and 18 miRNA-mRNA pairs. There were a total of 34 nodes, including 11 circRNAs, 10 miRNAs and 13 mRNAs. Among them, miR-1305 had the most miRNA target sites (six) while miR-1264 had two miRNA target sites, indicating that miR-1305 and miR-1264 serve as crucial miRNAs regulating the ceRNA network (Fig. 7; Fig. S7-sheet 5, sheet 6, sheet 7). There were a total of eight target genes regulated by these two miRNAs: COL5A2 (ID: 1290), DOK6 (ID: 220164), ABCA1 (ID: 19), GOLGA4 (ID: 2803), WDR3 (ID: 10885), NAA15 (ID: 80155), ROBO2 (ID:6092) and U2SURP(ID:23350). Based on their biological functions, five of these genes are possibly related to circANRIL regulation in the ceRNA network were assessed with qRT-PCR analysis:COL5A, DOK6 , ABCA1, GOLGA4 , and WDR3. ANRIL(ID: 100048912) is known to be related to altered circANRIL expression, and was thus also included in the qRT-PCR validation (Table 3). As shown in Fig. 8 (Fig. S8-sheet 8), only COL5A2 and WDR3 mRNA expression was markedly increased in the circANRIL-OE cell line, while the others showed decreased expression.