Identifications of potential therapeutic targets and drugs in angiotensin II-induced hypertension.

This study aimed to identify the underlying therapeutic targets of angiotensin II (AngII)-induced hypertension, and screen the related drugs.The gene expression profiles of GSE93579 and GSE75815 were used to identify differentially expressed genes (DEGs) between AngII-induced hypertension and control samples based on meta-analysis. These DEGs were analyzed using Gene-Ontology (GO) function and pathway enrichment methods. Subsequently, the weighed gene coexpression network analysis (WGCNA)-based meta-analysis was applied to determine transcriptional signature with DEGs. Additionally, the functions of the modules were analyzed based on the network, and miRNAs were identified. Finally, small molecule drugs correlation with DEGs was identified.In total, 346 upregulated DEGs (e.g., Rgs7 bp) and 360 downregulated DEGs (e.g., Ebf3) were identified between AngII and control samples. In addition, a total of 150 DEGs in the brown, red, and yellow modules with higher correlation coefficient according to WGCNA, were used to construct the coexpression network, including Rgs7 bp and Ebf3, etc. in brown modules. Besides, 3 modules were obtained after the functions of the modules analysis. Moreover, 5 miRNAs were integrated in modules, including miR-124A, miR-524, miR-493, miR-323, and miR-203. Finally, anisomycin was the highest correlation with DEGs.MiR-124a might be involved in the pathogenesis of hypertension via targeting Ebf3 and Rgs7 bp, which possibly represent a novel and effective strategy for treatment of hypertension. Anisomycin might be performed to reduce blood pressure by blocking MAPK signaling pathway.

Introduction

Hypertension (systolic blood pressure (BP) ≥140 mm Hg and diastolic BP ≥ 90 mm Hg) is an important public-health challenge and affects over a third of the people worldwide.[ In China, hypertension is a leading cause of premature death for cardiovascular disease, and its prevalence and mortality are more than 30% and 43%, and costing is about 6.6% of total health expenditure.[ Although the treatment of hypertension has been improved, the treatment rate in China was only about 20%.[ Thus, more drugs are usually required. However, the pathophysiology processes of hypertension are still complex.

In general, hypertension is caused by the interplay of genetic and environmental factor, such as inheritance, sodium chloride intake, smoking.[ The rennin–angiotensin system (RAS) has an important role in regulating blood pressure and the homeostasis of sodium and potassium.[ Angiotensin II (AngII) in RAS is known to have effects beyond the elevation of blood pressure by mitogen-activated protein kinase (MAPK) signaling pathway.[ Therefore, angiotensin-converting enzyme inhibitor and angiotensin receptor blocker are used to reduce the blood press. Therefore, the study of target biomarker is essential for hypertension.

MicroRNAs (miRNAs or miRs) as small non-coding RNAs, are well-known post-transcriptional regulation of mRNA expression. Some studies find that miRNAs are related to the etiology of hypertension, such as miRs20a-5p, 103a-2-5p, 30c-5p, 4763-5p, 4709-3p, and 4717-3p.[ In hypertensive patients, miR-126 and miR-9 expression levels display a significant positive correlation.[ Additionally, miRNA-130/301 members regulated endothelin-1 for pulmonary vascular cross-talk in pulmonary hypertension.[ However, the therapeutic targets in hypertension remain undefined.

In the present study, we aimed to identify the underlying therapeutic targets of AngII-induced hypertension, and screen the related drugs. The gene expression profiles of GSE93579 and GSE75815 were used to identify differentially expressed genes (DEGs) between AngII-induced hypertension and control samples based on meta-analysis. These DEGs were analyzed using Gene-Ontology (GO) function and pathway enrichment methods. Subsequently, the weighed gene coexpression network analysis (WGCNA)-based meta-analysis was applied to determine transcriptional signature with DEGs in brown, red, and yellow modules. In addition, three modules were obtained after the functions of the modules analysis based on the network, and miRNAs were integrated in the modules. Finally, small molecule drugs correlation with DEGs was identified.

Methods

Microarray data and preprocessing

The gene expression profiles of GSE93579 and GSE75815 were downloaded from the Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/) database through retrieving “hypertension and angiotensin.” Overall, the microarray data of 15 aorta samples from hypertensive mice were deposited at GEO server with GSE93579 by Rippe et al,[ including 8 AngII samples treated with AngII (1 μg/kg/min) and 7 control samples treated with saline. In their study, the protocols were approved by the animal care and use committee in Malmö/Lund (permit numbers M46-13 and M57-14) according to national guidelines and the European Communities’ Council Directive 86/609/EEC. The raw data was background corrected and normalized using oligo package in R software (version 3.3.2, https://cran.r-project.org/bin/windows/base/).[

In addition, the microarray data of 23 vascular tissue samples from hypertensive mice were deposited at GEO server with GSE75815 by Siedlinski,[ including 12 AngII samples from AngII-infused hypertensive mice at speed of 490 ng/kg/min and 11 control samples. This study was performed using relevant guidelines and Local Ethics Committee no. 1 in Kraków (Poland) approved the protocols employed (permissions nos. 151/2012, 100/2013, 254/2015, and 157/2016). The probes were corresponded to gene symbol according to the latest annotations file from NCBI gene database. When more than 1 probe was corresponded to the same gene symbol, the median value would be the expression level of gene. Subsequently, these data fitted log-normal distribution using log2 function, and normalized using median function of limma package in R software (Linear Models for Microarray Data, http://www.bioconductor.org/packages/release/bioc/html/limma.html). Finally, heatmap.sig.genes function in MetaDE package (https://cran.r-project.org/web/packages/MetaDE/) was used to show the correlation matrices of DEGs in GSE93597 and GSE75815.

Identification of DEGs using meta-analysis

DEGs between were identified using MetaDE package, which contains many advanced genomic meta-analysis methods to identify.[ First, heterogeneity test was performed. The threshold value of homogeneity test was set as tau2 = 0 and Qpval > 0.05, which means that these genes expression levels are homogeneous without bias in each gene expression profile. Then the threshold value of DEGs was set as |log fold change (FC)| > 1 and P value adjustment < 0.05.

Gene ontology (GO) enrichment function and pathway analysis

GOstat is helpful to analyze biological processes (BP) of genes from high-throughput experiments, such as microarrays.[ GO as the framework for the model of biology is used to describe gene function, and relationships between these concepts including BP, cellular components (CC), and molecular functions (MF). To comprehend the DEGs involved in BP, CC, and MF, the GO enrichment analyses were performed with GOstat (http://gostat.wehi.edu.au/) using the hypergeometric distribution method. The formula is , where N is the total number of genes enriched in all GO terms; n is the number of feature genes enriched in all GO terms; M is the number of genes in a specific GO term; and m is the number of feature genes in a specific GO term. P value < 0.05 was significantly different.

WGCNA of DEGs based on meta-analysis

WGCNA is a systems method for identifying target gene of disease in bioinformatics applications, which can be used for describing the correlation among genes by finding significant modules in the high-throughput sequencing data.[ First, meta-analysis was used to integrate the above DEGs expression levels and compare gene–gene correlation in the 2 microarray data. In the present study, GSE93597 as master dataset and GSE75815 as slave dataset were used to identify the modules related to AngII-induced hypertension. As follows: the correlation of the expression values of DEGs in the 2 datasets was analyzed. The higher the correlation, the higher the consistency of gene expression in each dataset. This is a prerequisite for constructing a WGCNA network using meta-analysis. The correlation matrix of gene coexpression was constructed according to Smn = |cor(m,n), where Smn indicates correlation coefficient of coexpression patterns between genes m and n. The adjacency is defined as amn = power(S, which measures pairwise correlations between gene expressions. Adjacency functions for both weighted and unweighted networks require the user to choose threshold parameters, correlation coefficient between log2 k (node conut) and log2 p(k) (frequency of node) ≥0.9 was set as threshold. The correlation matrix Smn was transformed to adjacency matrix amn. Then adjacency matrix amn was transformed to topological matrix using the following equation:

where lmn indicates the sum of adjacency coefficient of the common edge between genes m and n; km indicates sum of connection strengths of m with the other network genes. Gene significance (GS) measures were used to incorporate external information into the coexpression network. Module significance is determined through the average absolute GS for all genes in a module. Subsequently, coexpression network was visualized using the Cytoscape 3.3 software.

Identification of modules in the coexpression network

MCODE (degree cutoff = 2, node score cutoff = 0.2, and K-core = 2) and BINGO (adjusted P < .05) plugin in Cytoscape software (version 3.3, http://cytoscape.org/) were used to analyze the functions of the modules based on the network.[ Modules indicate the topological structures of genes interrelation in the network.[

Identification of miRNA related to DEGs in modules

WebGestalt was performed to search miRNA related to the modules in the coexpression network. The cutoff of P value was less than 0.05. Finally, the miRNA regulatory network was constructed using miRNAs and target genes based on the modules.

Small molecule drugs related to DEGs in coexpression network

These DEGs were uploaded to connectivity map, and then the small molecule drugs connected disease and DEGs were identified. The threshold value of small molecule drugs was set as |score| > 0.8.[

Results

DEGs in hypertension mice exposure to AngII

In total, 706 DEGs were identified between AngII and control samples based on meta-analysis, including 346 upregulated DEGs (e.g., Ebf3) and 360 downregulated DEGs. These DEGs could be well distinguished in AngII and control samples (Fig. 1). The top 10 upregulated and downregulated DEGs were shown in Table 1, such as Fos, Gprasp2, Bmper, and Lrg1.

Figure 1

The heatmap of differentially expressed genes (DEGs) between Angiotensin II (AngII)-induced hypertension and control samples. Red colors indicate overexpressed genes and green colors indicate downexpressed genes in AngII-induced hypertension mice compared with the corresponding control samples.

Table 1

The top 10 upregulated and downregulated DEGs between AngII and control samples based on meta-analysis.

The result of GO enrichment function and pathway analysis

In total 22 GO terms including 8 BP terms (regulation of transcription and regulation of RNA metabolic process, etc.), 7 CC terms (extracellular region and membrane-enclosed lumen, etc.), and 8 MF terms (ion binding and DNA binding, etc.) were enriched (Fig. 2A). In addition, 13 pathways, such as MAPK signaling pathway, focal adhesion, and regulation of actin cytoskeleton, were also identified (Fig. 2B).

Figure 2

The Gene-Ontology terms and pathway enrichment results of DEGs under the cutoff of P value < .05. A, The biological processes (BP), cellular components (CC), and molecular functions (MF) terms of DEGs were involved in hypertension mice exposure to AngII. The horizontal axis indicates the count of genes. The colors of column indicate the different P value. B, The pathway enrichment results of DEGs.

The construction of coexpression network

The consistency analysis showed that the correlation between GSE93597 and GSE75815 data was high (correlation coefficient = 0.71 and P = 5.2E−48) (Fig. 3). Gene cluster dendrogram is obtained by average linkage hierarchical clustering (Fig. 4A and B). In addition, the brown, red, and yellow modules had a higher correlation coefficient among these modules (P < .05, correlation coefficient > 0.6) (Fig. 4C and D). Subsequently, a total of 150 DEGs in the 3 modules were used to construct the weighed gene coexpression network, including 67 DEGs in brown modules (42 downregulated DEGs and 25 upregulated DEGs), 40 DEGs in red modules (such as 16 downregulated DEGs and 24 upregulated DEGs), and 43 DEGs in yellow modules (such as 21 downregulated DEGs and 22 upregulated DEGs) (Fig. 5).

Figure 3

The correlation of DEGs expression value between GSE93597 and GSE75815 data.

Figure 4

Weighed gene coexpression network analysis of DEGs data in GSE93597 and GSE75815 data. The color row underneath the dendrogram indicates the module in gene cluster dendrogram based on GSE93597 (A) and GSE75815 (B) data. The correlation between modules and hypertension in GSE93597 (C) and GSE75815 (D) data was calculated. The horizontal axis indicates modules with different colors.

Figure 5

The weighed gene coexpression network. The different colors indicate different modules including 150 DEGs.

The miRNA regulatory network

Three modules were obtained after the functions of the modules analysis, such as modules 1, 2, and 3 (Fig. 6A–C). In addition, the network on miRNAs and target genes based on the modules was constructed (Fig. 6A), including miR-124A, miR-524, miR-493, miR-323, and miR-203. Besides, the module 1 was involved in cellular process, macromolecule metabolic process, and so on; the module 2 was involved in intracellular signaling pathway, regulation of Rho protein signal transduction, and so on; the module 3 was involved in cobalamin transport, regulation of transcription termination, and so on.

Figure 6

Modules based on coexpression network. The miRNA regulatory network was constructed with miRNAs and target genes in the module 1 (A). The modules 2 (B) and 3 (C) were identified.

Small molecule drugs competing against AngII

Overall, 3 drugs among total 6 drugs connected DEGs in the coexpression network were negatively correlated with disease (Table 2). In other words, the 3 drugs could compete against the effect of AngII. Anisomycin has the highest correlation in the 3 drugs.

Table 2

Total 6 drugs connected DEGs in the coexpression network were negatively correlated with hypertension.

Discussion

In the present study, a total of 706 DEGs were identified between AngII and control samples based on meta-analysis, including 346 upregulated DEGs (e.g., Rgs7 bp) and 360 downregulated DEGs (e.g., Ebf3). In addition, a total of 150 DEGs in the brown, red, and yellow modules with higher correlation coefficient according to WGCNA were used to construct the weighed gene coexpression network, including 67 DEGs in brown modules (including Rgs7 bp and Ebf3, etc.), 40 DEGs in red modules, and 43 DEGs in yellow modules. Besides, 3 modules (modules 1, 2, and 3) were obtained after the functions of the modules analysis based on the coexpression network. Moreover, 5 miRNAs were integrated in the module 1, including miR-124A, miR-524, miR-493, miR-323, and miR-203. Finally, anisomycin was the highest correlation with DEGs in the coexpression network.

Overexpression of EBF3 has previously been served as tumor suppressors in several cancer types.[ Since dopamine plays an important role in regulating renal function and blood pressure in humans with essential hypertension, lack of dopamine receptors results in hypertension.[ In addition, overexpression of Ebf3 induces a significant increase in the terminal differentiation of dopamine neurons.[ Furthermore, dopamine receptor subtypes, such as D1 receptor (D1R) and D5R, closely related to AngII type 1 receptor (AT1R).[ For example, the D1R inhibits AT1R function by completing abolition of AT1R signaling in the short term (min) or decreasing the AT1R receptor in the long term (24 h).[ Therefore, overexpression of Ebf3 might be served as a candidate biomarker in AngII-induced hypertension. In addition, Rgs7 bp (or R7 bp) named regulator of G-protein signaling 7 binding protein is expressed in the mammalian nervous system.[R7BP strongly inhibits the effect on homo-oligomerization of regulator of G-protein signaling 7 by Tayou.[ Moreover, in a previous study, G-protein-coupled receptor 7 may be a promising target in hypertension via mediating the PLC/PKC pathway to regulate arterial myogenic tone.[ Therefore, Rgs7 bp might be served as a candidate biomarker in hypertension.

As reported in the previous studies, miRNAs are involved in the pathogenesis of various cardiovascular diseases, including pulmonary arterial hypertension characterized by profound elevation of blood pressure in the peripheral arteries of the lung.[ In addition, endothelial dysfunction is closely related to pathogenesis of hypertension. Meanwhile, hypertension can lead to endothelial dysfunction.[ In a review, miRNAs are defined as the potential therapeutic marker in many cardiovascular disorders through the complex interactions linking miRNAs, genes, and pathways leading to endothelial dysfunction.[ Therefore, miRNAs might be potential therapeutic marker in hypertension. Remarkably, miR-124a expression is regulated by miR-150,[ and both of them are usually tumor suppressors.[ Cox regression analysis confirms downregulated miR-150 levels are a significant predictor of poor survival in pulmonary arterial hypertension by decreased levels of KLF2.[ Moreover, miR-124a is reduced and its target genes expression levels (such as PTPB1 and PKM2) are increased in a rat model of pulmonary arterial hypertension.[ In this study, overexpression of Ebf3 and downexpression of Rgs7 bp were the target genes of miR-124a. Therefore, miR-124a was involved in the pathogenesis of hypertension via targeting Ebf3 and Rgs7 bp.

Anisomycin was the highest correlation with DEGs in the coexpression network. It was an inhibitor of protein synthesis through blocking synthesis of peptide band and was a specific antagonist of JNK (MAPK subfamilies) with a concentration of 25 ng/mL.[ In this study, MAPK signaling pathway was enriched in AngII-induced hypertension mice, including 20 DEGs. Similarly, AngII-induced hypertension production is possibly mediated by p38 MAPK.[ Thus, MAPK inhibitor might achieve clinical outcome for the treatment of hypertension. Therefore, anisomycin might offer an effective strategy to treat hypertension and prevent the elevation of blood pressure.

However, the study had also several limitations. As follows, sample sizes were not affluent; experimental verification was insufficient. In addition, the human samples should be used instead of mouse samples. All of these were the sources of potential bias or imprecision. However, the algorithms in the study have been identified in many studies and uniform. Therefore, the miR-124a expression levels would be measured, and the PCR of Ebf3 and Rgs7 bp will be used to verify the predicted results, and anisomycin blocked the effect of AngII by MAPK pathway in hypertension will be identified in the next study.

In conclusion, Ebf3 and Rgs7 bp might be served as candidate biomarkers in AngII-induced hypertension. Furthermore, miR-124a might be involved in the pathogenesis of hypertension via targeting Ebf3 and Rgs7 bp. MiR-124a, or its targets, possibly represent a novel and effective strategy for the treatment of hypertension. In addition, anisomycin potentially prevented the elevation of blood pressure by blocking MAPK signaling pathway.