Consistency of mouse models with human intracerebral hemorrhage: core targets and non-coding RNA regulatory axis

Changes in gene expression profile in the mouse brain after ICH

Three days after collagenase-induced ICH, the brains of the mice were harvested. Gene expression levels were determined and compared to healthy mice. A total of 39,483 genes were identified. The results showed an increase in the expression of 13,434 genes after ICH and a significant increase in 1093 genes. Of the 18,486 genes with decreased expression, 120 were significantly downregulated (|log2FC|>1 and p.adj<0.05, Supplementary Figure 1A). The top 20 Kyoto Encyclopedia of Genes and Genomes (KEGG) enriched pathways are shown in Figure 1A. The classical pathways were the AGE-RAGE pathway, the NF-κB pathway, the TLR pathway and the TNF pathway. Table 1 lists the differentially expressed genes (DEGs) involved in more than two pathways simultaneously.

Figure 1. KEGG and GO analysis of DEGs in mice with ICH. (A) Top 20 Kyoto Encyclopedia of Genes and Genomes (KEGG) of differentially expressed genes (DEGs). (p.adj ≤ 3.86×10^-6). (B) Top 5 biological processes (BPs) of Gene ontology (GO) and DEGs involved. (p.adj ≤1.28×10^-24). (ICH 3 d vs. Con, n=3).

The DEGs were analyzed by Gene Ontology (GO). The top five biological processes (BPs) and the genes involved are shown in Figure 1B. Up to 93 genes were involved in cytokine production, which was the top-ranked BP. Several well-known genes in this group were Pf4, Tlr family genes, Hmgb2, Il6 and Cd36. We also examined the top five cellular components (CCs) and molecular functions (MFs) of all DEGs (Supplementary Figure 1B, 1C). The results showed that these genes tended to be expressed in the chromosomes and extracellular matrix. The top five MFs in which the genes were involved were immune receptor activity, integrin binding, cytokine binding, chemokine activity and toll-like receptor (TLR) binding.

Changes in miRNA and circRNA expression levels in mouse brain after ICH

A total of 1170 known miRNAs were identified in mouse brain. After data normalization, 22 miRNAs showed significant changes in the ICH 3 d group compared to the control group (|log2FC|>1 and p.adj<0.05, Supplementary Figure 2A). Among these, the expression of 21 miRNAs was upregulated and only mmu-miR-375-5p was significantly downregulated (log2FC=-1.7, p.adj = 1.25 × 10^-6). A regulatory network of differentially expressed miRNA (miR)-gene was constructed (Figure 2A). It’s worth noting that up to 311 genes were identified as the target of mmu-miR-223-3p (log2FC=2.1, p.adj=2.54×10^-16). The mmu-miR-223-3p, mmu-miR-451a and mmu-miR-155-5p were upregulated after ICH, and their important role in regulating gene expression levels has been demonstrated in the cerebrospinal fluid of patients with subarachnoid hemorrhage [12].

Figure 2. Differential regulatory networks following ICH in mice. (A) Regulatory network of differentially expressed (DE) miRNAs and DEGs in ICH mice. (B) miRNAs that may be targets of circRNA-000895.

Next, we found 100 circRNAs with significant changes in expression levels in the ICH 3 d group. Among them, 71 were upregulated and 29 were downregulated compared to the control (Supplementary Figure 2B). Since circRNAs can inhibit the function of miRNAs by binding to them, predicting their binding sites is important to provide a theoretical basis for future research. Using miRanda and psRobot to predict binding sites, we identified one notable circRNA, mmu-circRNA-000895, which can bind to 199 miRNAs (Figure 2B) and was upregulated at ICH 3 d (log2FC=1.49, p.adj=4.6×10^-4).

In addition, we summarize the DEGs and the ceRNA networks that can potentially be regulated by circRNA-000895: miR-214-3p-Ccl9, miR-761-Ccl9, miR-362-5p-Itgb1, miR-15a-5p-Nfatc1, and miR-511-5p-Tnfrsf1a. Furthermore, circRNA-000895 can regulate the expression of Lypd1 by binding to miR-486a-3p or miR-486b-3p. In the transcriptome gene screening, the expression of Lypd1 decreased after ICH. Although its p-value was 0.02 and the false discovery rate was 0.13, it was not classified as a gene with significant changes after ICH because its FC was only 1.8. It can be speculated that circRNA-000895 may regulate multiple genes after ICH in a ceRNA-dependent manner.

DEGs in human peripheral blood 72 h and 24 h after ICH

We searched the GEO database for a 0-24 h and 72-96 h human peripheral blood transcriptome and re-analyzed the raw data using bioinformatics (GEO accession number: GSE 125512). Data were obtained from eleven patients (four African Americans and seven white Americans) with a mean age of 57.5 years. Blood samples were taken 24 and 72 hours after the onset of ICH and transcriptome sequencing was performed. A total of 16,640 genes were identified in human peripheral blood. As all data are from post-onset blood, we have relaxed the fold change accordingly [13]. For genes with p.adj<0.05, log2FC≤-0.5 indicated that the genes were downregulated in ICH 72 h, and log2FC≥0.5 indicated that the genes were upregulated. Using these criteria, 456 genes were significantly upregulated and 183 were significantly downregulated in ICH 72 h (Supplementary Figure 3A).

GO analysis was also performed on these DEGs. These DEGs were involved in 171 biological processes (Figure 3A). The first four BPs were related to neutrophil activity and these BPs involved the same 35 DEGs; These genes are expressed in the bone marrow under physiological conditions and were significantly up-regulated 72 h after ICH, indicating that ICH injury triggered a systemic inflammatory response on the third day after the onset of the disease. The results of the analysis of the top 20 CCs and MFs are shown in Supplementary Figure 3B, 3C.

Figure 3. KEGG and GO analysis of DEGs in human patients with ICH. (A) Top 20 BPs enrichment of GO of DEGs. (p.adj ≤ 7.0×10^-6). (B) Top 5 of KEGG enrichment of DEGs (p.adj ≤0.04). (C) The DEGs involved in top 5 KEGG. (ICH 72h vs. ICH 24h, n=11).

KEGG analysis was then performed. Among the top five KEGG enrichment pathways in human peripheral blood DEGs (Figure 3B, 3C), four of them also appeared in the KEGG enrichment pathways of mouse DEGs: ‘PI3K-Akt signaling pathway’, ‘Malaria’, ‘Viral protein interaction with cytokine and cytokine receptor’, and ‘Hematopoietic cell lineage’. To further understand the regulation mechanism of these genes, networks of the top five KEGG DEGs were constructed (Figure 4).

Figure 4. Network of TF/miRNAs and DEGs (top five KEGG pathways) in ICH patients.

Seven genes significantly differentially expressed after human/mouse ICH

After analyzing the human and mouse ICH transcriptome data, the DEGs of the two species were integrated and seven identical genes were identified. These genes are: Tlr4, Il1r2, Csf1, Cx3cr1, Cxcl5, Pf4 and Ppbp were all upregulated in mouse brain in ICH 3 d. However, in human peripheral blood, the expression levels of TLR4 and IL1R2 were higher at 24 h than at 72 h after ICH. The regulatory networks of these genes were then constructed from three aspects: TFs and miRNAs (Figure 5A). According to the predictions, 56 miRNAs target CXCL5, of which hsa-miR-765 and hsa-miR-148b-3p simultaneously target CSF1, and hsa-miR-7-5p and hsa-miR-942-5p target TLR4. In addition, hsa-miR-26b-5p targets TLR4 and IL1R2, and hsa-miR-335-5p targets CSF1, PPBP and TLR4. In addition, NFKB1 and RUNX1 are both TFs of CSF1 and PF4.

Figure 5. Regulate network and enrichment analysis for the seven DEGs of human/mouse ICH. (A) Genes differentially expressed in both human peripheral blood and mouse brain after ICH and their regulated network, including miRNAs and TFs. (B) The top 20 KEGG enrichment pathways involved with these genes and their TFs (p.adj ≤ 0.001). (C) The functional networks of these genes and TFs.

Next, KEGG enrichment analysis was performed for the seven genes and their TFs, and the data were visualized for the top 20 pathways (Figure 5B). A functional network of the seven genes and their TFs was also constructed (Figure 5C). NFKB1 and RELA, the TFs of CSF1 and CXCL5, were involved in all pathways shown in Figure 5C, except cytokine receptor interaction. In addition, miR-30b/c-5p also have binding sites on CSF1. CSF1 and TLR4 share the same TF named SPL1. CXCL5 and CSF1 can be transcribed by the TFs NFKB1 and RELA and both have targeted binding sites with hsa-miR-148b-3p; this miRNA is also homologous in human and mouse. We also found that hsa-miR-335-5p targets PPBP, CSF1 and TLR4. miR-335-5p is highly homologous in humans and mice. Notably, mmu-miR-335-5p has potential binding sites with mmu-circRNA-000895. It should be noted that the presence of hsa-SEMA3F_0001 has been detected in human urine exosomes (ID no. exo_circ_52661 in exoRBase v2). Therefore, the circRNA-000895-miR-335-5p-PPBP/CSF1/TLR4 signaling pathways are highly relevant for translational medicine in the context of ICH.

RT-qPCR experiments were performed to investigate the expression levels of these genes after ICH in mice. The relative expression levels of Pf4, Ppbp, Cx3cr1, Cxcl5, Csf1 and Tlr4 were consistent with the transcriptome data. Compared to the control group and ICH 1d, all these genes were upregulated in ICH 3 d (Figure 6). The expression level of Il1r2 peaked at ICH 1d and decreased in ICH 3 d, but was still higher than in control mice. This expression trend was most consistent with human ICH patients.

Figure 6. Relative expression levels of the seven genes after ICH in mice. (A) TLR4: ICH 3 d vs. Con, ICH 3 d vs. ICH1 d, ICH 1 d vs. Con, ***P< 0.0001. (B) IL1R2: ICH 3 d vs. Con, **P < 0.01; ICH 1 d vs. Con, ***P< 0.0001. (C) CSF1: ICH 3 d vs. Con, ***P< 0.001. (D) CX3CR1: ICH 3 d vs. Con, *P < 0.05. (E) CXCL5: ICH 3 d vs. Con, ***P< 0.0001; ICH 3 d vs. ICH1 d, **P < 0.01. (F) PF4: ICH 3 d vs. Con, ICH 1 d vs. Con, ***P< 0.0001. (G) PPBP: ICH 3 d vs. Con, ICH 3 d vs. ICH1 d, ***P< 0.001. (n = 6/group/time point).

Limitation

This article detects differentially expressed mRNAs, circRNAs and miRNAs after ICH. Although the sample sizes of human and mouse were small, we were able to identify the DEGs common to both species after the disease. Further research is needed to develop a more detailed regulatory network of DEGs with translational significance after ICH. The expression levels of these genes in healthy individuals and those at high risk of ICH should be included to confirm whether they can be used as early biological warning signals; this will promote research into their regulation mechanism after ICH and determine their potential as drug targets for human non-infectious diseases.

Animals

Eight-week-old male C57BL/6 mice weighing 20-25 g were purchased from Liaoning Changsheng Biotechnology Company (Liaoning, China). Mice were housed in the laboratory for at least 1 week. Mice food was purchased from the Animal Experimental Center of the First Affiliated Hospital of Harbin Medical University.

Grouping and ICH models

Mice were randomly divided into two groups: Control group (Con, n=3) and ICH group (ICH, n=3). Collagenase was used to induce a murine model of ICH. The experimental method is described in a previous study by Wang et al. [41]. Briefly, anesthetized mice were fixed to a stereotaxic instrument and, after skin preparation, a longitudinal incision was made between the two ears. After positioning the skull with the stereotaxic instrument, a hole was drilled in the skull and 0.4μL collagenase was injected into the striatum with a micro-syringe. After waiting for absorption, the gap in the skull was filled with wax, the skull wiped with saline and the skin sutured. Successful establishment of the model was confirmed by the behavior of the mice; neurobehavioral injury was the inclusion criterion for the experiment.

Sample preparation

Samples for gene expression profiling and mRNA, miRNA and circRNA analysis were prepared using the same treatment. C57BL/6 mice were deeply anesthetized with inhalation anesthetic machine using isoflurane at ICH 3 d. After the heart was exposed, the systemic circulation was perfused with saline. Brain tissue was removed and stored at -80° C.

Messenger RNA isolation and microarray scanning

Brain tissue was pulverized with liquid nitrogen. TRIzol was used for mRNA extraction. RNA quantity and quality were measured, and RNA integrity was assessed by standard denaturing agarose gel electrophoresis. Total RNA from each sample was linearly amplified and labeled with Cyanine 3 (Cy3)-UTP. Subsequently, the labeled cRNA was then fragmented and hybridized, and 100 μL of hybridization solution was added to the gene expression microarray slide. The hybridized arrays were washed, fixed and scanned using a DNA microarray scanner (Agilent G2505C Scanner, Agilent Technologies, USA).

MicroRNA isolation and microarray scanning

Sequencing was supported by Annoroad Gene Technology Co, Ltd. (China). The mirVana RNA Isolation Kit was used to isolate low-molecular-weight RNA from mouse brain tissue in ICH 3 d and control mice. The Illumina HiSeq 4000 platform was used to sequence miRNA expression profiles according to the manufacturer’s instructions.

CircRNA isolation and scanning

After total RNA was extracted from mouse brain tissue, the purity and concentration of the total RNA samples were determined using the NanoDrop ND-1000 to confirm that the samples could be used in the next steps. Total RNA was digested with ribonuclease R (RNase R) (Epicentre Technologies, USA) to remove linear RNA and enrich circRNAs. The circRNAs were then amplified and transcribed into fluorescent cRNAs using a random priming method (Super RNA Labeling Kit, Arraystar, USA). The labeled cRNAs were purified using the RNeasy Mini Kit (Qiagen, USA). The concentrations and specific activities of the labeled cRNAs were measured using the NanoDrop ND-1000. Labeled cRNAs were fragmented and hybridized, and 50μL of hybridization solution was added to the circRNA expression microarray slide. The hybridized arrays were washed, fixed and scanned using the Agilent G2505C scanner, and the Agilent Feature Extraction software (version 11.0.1.1) was used to analyze the array images.

RT-qPCR

After obtaining brain tissue from mice at different ICH time points (n=6), total RNA was extracted using a total RNA extraction kit (DP431, TIANGEN), and the concentration and purity were determined using a NanoVue Plus Micro-Volume UV-Vis spectrophotometer. The total RNA was then reverse transcribed into cDNA using a reverse transcription kit (FSQ-201, TOYOBO), and the concentration and purity of the cDNA were determined. Finally, 20μL of amplification reaction mix (QPS-201, TOYOBO) was prepared and the samples were introduced into an ABI StepOne Real-Time PCR system. The primer sequences are listed in Table 2. The ΔΔCT method was used to calculate the results. The data obtained were analyzed by one-way ANOVA. Tukey multiple comparison was used to test for statistical significance between groups.

Statistical analysis

Functional annotation of differentially expressed genes

The R package edgeR was used for differential expression analysis. Differentially expressed genes (DEGs) with |log2FC|>1 and adjusted p-values(p.adj)<0.05 were identified. Gene ontology (GO) annotation (BP: biological process, CC: cellular component, and MF: molecular function) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis were used to select the major functional and biological pathways of DEGs with Benjamini-Hochberg (BH) p.adj <0.05. These genes were enriched using the R package clusterProfiler, and the R package GOplot was used to identify the top five BPs, CCs and MFs.

Construction of regulatory network for differentially expressed miRNAs and genes

The edgeR package was used to identify differentially expressed miRNAs, and miRTarBase was used to find miRNA-gene interactions. A regulatory network of differentially expressed miRNAs and genes was then constructed. This network was visualized using Cytoscape software.

Construction of regulatory network for differentially expressed circRNAs and miRNAs

The edgeR package was used to identify differentially expressed circRNAs, and miRanda and psRobot were used to estimate the miRNA binding sites of sheared circRNAs. The circRNA-miRNA regulatory network was then constructed using Cytoscape software.

Bioinformatic analysis of a human disease model

Transcriptomes from human disease models were obtained from the Gene Expression Omnibus (GEO) database (GSE125512) and a volcano map was used to display DEGs with |log2FC|>0.5 and p.adj <0.05. Similarly, GO annotation and KEGG enrichment analysis were used to identify the primary functions and biological pathways of enriched DEGs with BH p.adj<0.05. The R package GOplot was used to highlight the five most significant KEGG pathways. In addition, miRTarBase was used to identify miRNA-gene pairs and the TRRUST (v2.0) database to identify transcription factor (TF)-gene pairs. A multifactorial regulatory network was then constructed.

Integration of mouse and human regulatory networks

Homology analysis was used to identify a primary regulatory network containing seven genes, miRNAs and TFs in a combination of critical disease-associated functional molecules in mice and humans. KEGG enrichment analysis was then performed on the combined network. Finally, Cytoscape software was used to generate a functional network map based on homologous genes and key TFs.