MicroRNAs in the miR-17 and miR-15 families are downregulated in chronic kidney disease with hypertension.

BACKGROUND: In older adults (aged 70-74 years), African-Americans have 4-fold higher risk of developing hypertension-attributed end-stage renal disease (ESRD) than European-Americans. A hypothesized mechanism linking hypertension and progressive chronic kidney disease (CKD) is the innate immune response and inflammation. Persons with CKD are also more susceptible to infection. Gene expression in peripheral blood can provide a view of the innate immune activation profile. We aimed to identify differentially expressed genes, microRNAs, and pathways in peripheral blood between cases with CKD and high blood pressure under hypertension treatment versus controls without CKD and with controlled blood pressure in African Americans.
METHODS: Case and control pairs (N = 15x2) were selected from those without diabetes and were matched for age, sex, body mass index, APOL1 risk allele count, and hypertension medication use. High blood pressure under hypertension treatment was defined as hypertension medication use and systolic blood pressure (SBP) ≥ 145 mmHg. CKD was defined as estimated glomerular filtration rate (eGFR) < 60 mL/min/1.73m2. Cases were selected from those with CKD and high blood pressure under hypertension treatment, and controls were selected from those without CKD (eGFR: 75-120 mL/min/1.73m2 and urine albumin-to-creatinine ratio < 30mg/g) and with blood pressure controlled by hypertension medication use (SBP < 135 mmHg and D(diastolic)BP < 90 mm Hg). We perform RNA sequencing of mRNA and microRNA and conducted differential expression and co-expression network analysis.
RESULTS: Of 347 miRNAs included in the analysis, 14 were significantly associated with case status (Benjamini-Hochberg adjusted p-value [BH p] < 0.05). Of these, ten were downregulated in cases: three of each belong to the miR-17 and miR-15 families. In co-expression network analysis of miRNA, one module, which included 13 of the 14 significant miRNAs, had significant association with case status. Of the 14,488 genes and 41,739 transcripts included in the analysis, none had significant association with case status. Gene co-expression network analyses did not yield any significant associations for mRNA.
CONCLUSION: We have identified 14 differentially expressed miRNAs in the peripheral blood of CKD cases with high blood pressure under hypertension treatment as compared to appropriate controls. Most of the significant miRNAs were downregulated and have critical roles in immune cell functions. Future studies are needed to replicate our findings and determine whether the downregulation of these miRNAs in immune cells may influence CKD progression or complications.

Introduction

About one-third of end-stage renal disease (ESRD) in the U.S. has been attributed to hypertension [1]. In older adults, African-Americans have 4-fold higher risk of developing hypertension-attributed ESRD than European-Americans, and the incidence of hypertension-attributed ESRD increases with age [2]. A hypothesized mechanism linking hypertension and kidney function decline is the innate immune response and inflammation [3]. Inflammation biomarkers in blood have been associated with kidney function decline and incident hypertension [4-7]. Additionally, older adults with chronic kidney disease (CKD) have increased risk of infection-related hospitalization or bloodstream infection [8,9]. Although immune dysfunction in persons with ESRD due to metabolic disorder and retention of uremic solute is well established [10,11], the mechanisms underlying the links between immune function, hypertension, and kidney function in earlier stages of CKD is not well understood. Gene expression analysis in peripheral blood can provide a view of the innate immune activation profile and may lead to insights into the pathophysiology of hypertension and CKD, and their complications.

Studies using immune cells have identified specific gene expression profiles of CKD patients versus hemodialyzed patients, and of patients with essential hypertension versus controls [12,13]. Recently microRNAs (miRNA) have been shown to have a potential role in influencing blood pressure and kidney function through modulating the immune response [14-16]. miRNAs are small noncoding RNAs (~22 nucleotides in length) with important post-transcriptional regulatory functions and are an attractive target of investigations.

Given that African Americans have a strong predisposition for hypertension-attributed kidney disease, we decided to study hypertension and kidney disease jointly in older African Americans. We aimed to identify differentially expressed genes and miRNAs in the peripheral blood to gain insights into the immune profile of hypertension and kidney disease. We used a matched case-control design (n = 15x2) in the African Americans cohort of the Atherosclerosis Risk in Communities (ARIC) study [17]. We conducted RNA sequencing to quantify expression levels and performed differential expression and co-expression network analyses to identify gene expression profiles in the peripheral blood of hypertension and CKD cases.

Methods

Ethics statement

Written informed consent was obtained from all ARIC study participants, and approval was granted from the relevant institutional review boards (IRB) for the participating study centers (University of North Carolina, University of Minnesota, University of Mississippi Medical Center, and Johns Hopkins University). JHSPH IRB No. H.34.99.07.02.A1, study title “Atherosclerosis Risk in Communities (ARIC) Study—Morbidity/Mortality Follow-up Field Center.” This research was conducted in accordance with the principles described in the Declaration of Helsinki.

Study design

The ARIC Study is an ongoing prospective cohort study in four US communities [17]. A total of 15,792 participants aged 45–64 years were recruited from Forsyth County, North Carolina; Jackson, Mississippi; suburban Minneapolis, Minnesota; and Washington County, Maryland between 1987 and 1989 (Visit 1). Four follow-up examinations (visits 2–5) have been conducted. Blood samples preserved for RNA analysis using PAXgene blood tubes were available from visit 5 (2011–13).

Case and control pairs (N = 15 pairs) were selected from those individuals without diabetes but on hypertension treatment and they were matched by age, sex, body mass index, APOL1 renal risk allele count, and hypertension medications (ACE inhibitor, angiotensin receptor blocker, and calcium channel blocker) to reduce heterogeneity. The use of hypertension medication was included in the selection criteria because a high proportion (>60%) of older African Americans were on hypertension treatment [18], and some hypertension medications have been reported to influence gene expression in peripheral blood [19]. In the ARIC study, hypertension medication use was determined by inspection of medication bottles at study visit. High blood pressure under hypertension treatment was defined as systolic blood pressure (SBP) ≥ 145 mmHg with hypertension medication use. CKD was defined as estimated glomerular filtration rate (eGFR) < 60 mL/min/1.73m2. Cases were selected from those with CKD and high blood pressure under hypertension treatment. Controls were selected from those with blood pressure controlled by hypertensive medications (SBP < 135 mmHg and D(diastolic)BP < 90 mm Hg) and without CKD (eGFR: 75–120 mL/min/1.73m2 and urine albumin-to-creatinine ratio [UACR] < 30mg/g). The APOL1 renal risk variants are strongly associated with CKD progression [20], thus matching by APOL1 risk variants provided the opportunity to identify differential expressed genes that were independent of APOL1. The characteristics of the cases and controls were compared using the t-test, Kruskal-Wallis, or Fisher’s exact test, as appropriate.

mRNA and miRNA library preparation, sequencing, processing, alignment, and quantitation methods are described in S1 Text.

Surrogate variable generation

Surrogate variables (SVs) representing unknown confounders were estimated using the svaseq function from the R package sva [21], with the null model absent of variables, and case-status as the variable of interest in the full model. SVs were included as covariates in models in all subsequent expression analyses as expanded upon in the Differential Expression Analyses Methods sections.

mRNA differential expression analyses

The mRNA analyses were conducted on genes and transcripts from a list of 20,377 genes tagged as “protein-coding” in the GENCODE V19 annotation files. Results were adjusted within each analysis using the Benjamini-Hochberg [22] procedure for multiple testing.

The R package DESeq2 [23] was used to conduct gene-level analyses, with the null hypothesis of a zero log2-fold change. As DESeq2 incorporates normalization and default outlier replacement procedures, SVs were generated on DESeq2 outlier-replaced, normalized counts and included as covariates for these analyses. Gene-level analyses were carried out on14,488 genes with a normalized count ≥1 in at least 14 out of 29 samples (~50% as one sample failed library preparation, see Results), with a target FDR of alpha = 0.05 used for independent filtering. SVs were generated using 10,524 genes with a minimum normalized count of 10 in all 29 samples.

Gene-level kidney-focused analyses were also conducted, in which 397 kidney-expressed proteins from the Human Protein Atlas [24] (v14) (http://www.proteinatlas.org/humanproteome/kidney) were examined. Of these, 392 were available in our data; four were not available in the ENSEMBL GRCh37 assembly, and one was not in our final GTF of protein-coding genes. Genes meeting the normalized count threshold defined as above for expression were subset to produce the final kidney-specific results.

The Ballgown [25] R package was used to analyze 41,739 transcripts with fragments per kilobase of transcript per million mapped reads (FPKM) ≥ 0.3 in at least 14 samples, and SVs were generated on transcripts with FPKM ≥ 0.5 in at least 26 samples (90%) (17,085 transcripts).

miRNA differential expression analyses

The miRNA analyses were conducted on mature miRNAs from miRBase [26] v20 and novel miRNAs, detected as described in S1 Text. Similar to the mRNA analyses, the R package DESeq2 was used for miRNA analyses with the null hypothesis of a zero log2-fold change, and SVs were generated on the outlier-replaced normalized counts and included as covariates in the model. These analyses were carried out on 347 miRNAs with a normalized count ≥ 1 in at least 15 out of 30 samples (50%) with target FDR of alpha = 0.05 used for independent filtering, and SVs were generated using 270 miRNAs with a minimum normalized count of 3 in at least 27 samples (90%). Counts were summed across all precursors for each mature miRNA. Results were adjusted using the Benjamini-Hochberg procedure for multiple testing. We performed hierarchical clustering of the normalized counts of the miRNAs with nominal p-value < 0.01 to visualize the expression levels of these miRNAs between cases and controls.

Co-expression network analysis

We used the Weighted Gene Co-Expression Network Analysis (WGCNA) R package to construct correlated signed network and tested for associations between the eigengene (the first principal component) of each module and case status [27]. The genes, transcripts, and miRNAs used for co-expression network analysis and the generation of surrogate variables was the same as those for differential expression analyses. The gene and miRNA counts were first normalized by the size factor with outlier replacement using DESeq2 [23]. The transcript counts were normalized to FPKM using Ballgown [25]. Next, the gene, transcript, or miRNA expression levels were transformed using natural logs after adding 1. Residuals of the log transformed expression levels were generated with adjustment for surrogate variables. These residuals were used as input for co-expression analysis. In network construction, we used bi-weight mid-correlation as a measure of co-expression to minimize the influence from outliers. For soft thresholding power, we used the first power with adjusted R square for linear fit > 0.8. The minimum module size was set as 25, 30, and 8 for genes, transcripts, and miRNA, respectively. The significant threshold for module association was defined as 0.05 divided by the number of correlated modules in each specific analysis.

miRNA target prediction and gene ontology analyses

The procedure for miRNA target prediction and gene ontology analyses of differentially expressed miRNAs is similar to the procedure described by the authors of the empiricalGO [28] software in Bleazard et al. 2015. Briefly, miRanda [29] v3.3a was used for miRNA target prediction with parameters free energy < −20 kcal/mol and score > 155 [28]. The 3’UTR sequences for target prediction were obtained from GRCh37 ENSEMBL BioMart [30]. Subsequently, gene ontology (GO) annotations were obtained from GRCh37 ENSEMBL BioMart, and the empiricalGO python script was used in “basic” mode to produce a list of GO terms with a one-tailed permutated p-value for each term. The final results were subset to include only terms with a minimum size of 10, and the Benjamini-Hochberg multiple testing correction was applied to adjust for the number of GO terms in each list.

Results

Study population characteristics

The mean age of the cases and controls was 77 years, and 67% were female. No significant differences were detected between cases and controls for all of the matching characteristics (Table 1). The cases had significantly lower eGFR (mean of 46 min/mL/1.73m2 vs 88 min/mL/1.73m2, p < 0.0001) and higher SBP (mean of 156 mm Hg vs. 115 mm Hg, p < 0.0001). UACR was higher in cases (median of 13.7mg/g vs 7.8mg/g), but the difference was not significant (p = 0.13). On hypertension medications, none of the participants were on both angiotensin-converting-enzyme (ACE) inhibitor and angiotensin receptor blocker while two participants were on ACE inhibitor and beta blocker, and four participants were on ACE inhibitor and calcium channel blocker. These participants were split evenly between the case and control groups. Since gene expression levels were measured in whole blood, we also compared cases and controls for white blood cell count and percentage of lymphocyte, monocyte, and granulocyte, and did not observe significant differences.

Table 1

Study population characteristics.

Variable	Case	Control	p
Selection variables
eGFR, mL/min/m2, mean (SD)	46 (12)	88 (10)	N/A
SBP, mm Hg, mean (SD)	156 (12)	115 (33)	N/A
Matching variables
Age, year, mean (SD)	77 (4.5)	77 (5.4)	0.97
Female, % (n)	67 (10)	67 (10)	1.00
BMI, kg/m2, mean (SD)	28 (3.5)	30 (3.3)	0.09
APOL1 high-risk, % (n)	33 (5)	33 (5)	1.00
Hypertension medications
ACE inhibitor, % (n)	53.3 (8)	60 (9)	1.00
Angiotensin receptor blocker, % (n)	27 (4)	27 (4)	1.00
Beta blocker, % (n)	33.3 (5)	33.3 (5)	1.00
Calcium channel blocker, % (n)	47 (7)	60 (9)	0.71
Blood cell type variables
White blood cell count, 1000 per mm3, median (1st, 3rd quartile)	5.4 (4.1, 6.2)	4.8 (4.4, 5.9)	0.80
Lymphocyte, %, median (1st, 3rd quartile)	33.6 (29.1, 39.0)	30.4 (27.8, 40)	0.72
Monocyte, %, median (1st, 3rd quartile)	13.8 (12.8, 17.2)	13.4 (12.5, 16.0)	0.48
Granulocyte, %, median (1st, 3rd quartile)	50.9 (45.1, 57.4)	54.8 (46.3, 59.8)	0.48
Other variables
DBP, mm Hg, mean (SD)	76 (12)	70 (19)	0.31
Albuminuria, median (1st, 3rd quartile)	13.7 (5.9, 66.7)	7.8 (5.7, 13.3)	0.13
High sensitive C-reactive Protein,	1.4 (1.01, 6.87)	5.3 (2.5, 6.3)	0.31
Serum creatinine, mg/dL	1.6 (0.9)	0.8 (0.1)	0.003
Diuretic use, % (n)	73.3 (11)	64.3 (9)	0.70

p, p-value, N/A, not applicable

Sequencing and processing of mRNAs and miRNAs

Sequencing of mRNA was successful for 29 out of 30 samples; one sample (case) failed library preparation. All 30 samples were successfully sequenced for miRNA. Quality control (QC) and mapping statistics for the 29 samples with mRNA data (depth: 18.7M-45.1M paired-end reads) and the 30 samples with miRNA data (depth: 6.2M-11.5M single-end reads) are listed in S1 and S2 Tables, respectively. As S1 Table indicates, External RNA Controls Consortium (ERCC) transcripts were poorly detected in one sample, but this individual was retained for analysis because the distribution of counts in other genomic features aligned with those of the other samples, indicating a possible issue with the ERCC spike-in itself.

Differential expression analysis of genes, transcripts, and miRNAs

The gene-level differential expression analysis of 14,448 protein-coding genes considered “expressed” (normalized count ≥ 1 in ≥50% of samples) produced no significant genome-wide results after Benjamini-Hochberg (BH) multiple testing correction (S3 Table). The relevance of the kidney to both hypertension and CKD led us to examine a set of 154 genes meeting our expression threshold and expressed in the kidney, based on data available from the Human Protein Atlas [24], which contains protein expression data from four kidney samples, with proteins in the glomeruli, proximal tubules, distal tubules and the collecting ducts. In this subset of genes, none were significant after BH correction. The gene with the lowest p-value was SMIM24 (unadjusted p = 5.67x10-4, BH p = 0.087).

Transcript-level analyses were conducted on 41,739 transcripts with FPKM ≥ 0.3 in ≥ 50% of samples. After BH correction, there were no differentially expressed transcripts identified; results with p < 10−3 are presented in S4 Table.

We tested 347 miRNAs, including 12 novel miRNAs, for differential expression. We detected 14 significant miRNAs. Four were upregulated, and 10 were downregulated in the cases as compared to the controls (Table 2). The most significant miRNA was miR-17-5p (log2 fold change = -0.77, BH p = 6.7E-4). Two other miRNAs in the miR-17 family (miR-106a-5p, miR-106b-3p) were also significantly downregulated. The miR-15 family has three members that were significantly downregulated in the cases (miR-15a-5p, miR-15b-5p, and miR-16-5p). Twenty members of the miR-17 and miR-15 families were found in our data and all, but two, were downregulated in cases, although the associations of some were not significant (S5 Table). S1 Fig presents a heat map of the normalized count of the 32 miRNAs with nominal p-value < 0.01.

Table 2

Differentially expressed miRNAs (BH p<0.05).

miRNA	Mean of normalized count	Log2 Fold Change	SE	p	BH p
hsa-miR-17-5p	71.97	-0.77	0.16	1.94x10-6	6.76x10-4
hsa-miR-130a-3p	492.19	-0.60	0.15	4.31x10-5	7.47x10-3
hsa-miR-15b-5p	215.32	-0.72	0.19	1.24x10-4	1.17x10-2
hsa-miR-106b-3p	1425.98	-0.65	0.17	1.35x10-4	1.17x10-2
hsa-miR-106a-5p	5.46	-1.03	0.28	2.12x10-4	1.34x10-2
hsa-miR-16-5p	6983.54	-0.64	0.17	2.33x10-4	1.34x10-2
hsa-miR-181a-5p	15301.70	-0.59	0.16	2.80x10-4	1.39x10-2
hsa-miR-1285-3p	561.01	-0.44	0.13	4.92x10-4	2.14x10-2
hsa-miR-15a-5p	430.52	-0.76	0.23	8.36x10-4	3.18x10-2
hsa-miR-29c-5p	73.85	0.55	0.17	9.18x10-4	3.18x10-2
hsa-miR-345-5p	715.07	0.60	0.18	1.04x10-3	3.27x10-2
hsa-miR-142-3p	89.46	0.65	0.20	1.23x10-3	3.52x10-2
hsa-miR-339-3p	507.26	0.36	0.11	1.32x10-3	3.52x10-2
hsa-miR-210-3p	440.49	-0.39	0.12	1.50x10-3	3.72x10-2

SE, standard error; p, p-value; BH p, Benjamini-Hochberg adjusted p-value

Experimentally validated targets in humans were available for 11 of the 14 significant miRNAs in miRTarBase [31] 6.0 (S6 Table). The number of validated target genes for each significant miRNA ranged from 4 (miR-1285-3p) to 134 (miR-17-5p), with a median of 56. Altogether 272 unique genes have experimental evidence of being regulated by the 11 significant miRNAs. Of these, 239 were detected in our data, and 8 were associated with case status at p < 0.05 (S7 Table).

Co-expression network analysis of genes and transcripts

To investigate whether patterns in gene co-expression may be related to case-status, we conducted gene co-expression analysis. Of the 14,488 genes included in the analysis, 178 were assigned into four co-expression modules (S8 Table). The rest were pruned by WGCNA due to low correlation (< 0.3) with the eigengene in each module. The association between the eigengene of the brown module (consisting of 35 genes) and case status was nominally significant (p = 0.03, S9 Table). Of the 41,739 transcripts included in the analysis, no correlated modules were detected.

Co-expression network analysis of miRNAs

Of the 347 miRNA included in co-expression analysis, 182 were assigned into five co-expression modules with correlated miRNAs (S10 Table). The rest (n = 165) were pruned due to low correlation with the eigengene in each module. The eigengene of the turquoise module (108 miRNAs) was significantly associated with case status (p = 0.005, Table 3). The eigengene of the blue module (29 miRNAs) was nominally associated with case status (p = 0.03). For both modules, case status was associated with lower eigengene values suggesting an association with lower miRNA expression in these modules. The turquoise module included 13 of the 14 miRNAs that were significantly downregulated in the cases (Table 2).

Table 3

Association between the eigengene of each miRNA co-expression module and case status.

Module	Betab	SE	pa	Number of miRNAs
Turquoisec	-0.97	0.32	5.38x10-3	108
Blue	-0.80	0.34	2.63x10-2	29
Green	-0.06	0.37	8.63x10-1	11
Brown	0.03	0.37	9.32x10-1	23
Yellow	-0.01	0.37	9.84x10-1	11

SE, standard error; p, p-value

a p-values represent association of the module eigengenes, the first principal component of expression levels of a module, with case status.

b Beta values are in log expression level units.

c 13 out of 14 differentially expressed miRNAs (excluding hsa-miR-339-3p), belong to the turquoise module.

Gene ontology analysis of significant miRNAs

We then proceeded to predict target genes of the significant miRNAs using miRanda. The predicted targets were analyzed for enrichment of gene ontology annotations using empiricalGO with all predicted targets of the 347 expressed miRNAs as the “universe.” empiricalGO counted 5,576 targets for the 14 differentially expressed miRNAs from the target genes, of which 4,431 had associated GO terms. No GO terms with minimum size of 10 genes in our data were significant after Benjamini-Hochberg multiple test correction, while 25 terms had a one-tailed permutated p < 0.05 (S11 Table). All results for mRNA and miRNA analyses are summarized in S2 Fig.

Discussion

Main findings

We have identified 14 miRNAs that were significantly associated with CKD and high blood pressure under hypertension treatment. Ten of these miRNAs were downregulated in the cases, and 13 were grouped in one module in co-expression network analysis. Three of the 14 miRNAs belong to the miR-17 family, and three belong to the miR-15 family.

The miRNAs that were significantly associated with case status have critical roles in immune cell function. First, miR-17-5p, the most significantly downregulated miRNA, belongs to a cluster of miRNAs located in intron 3 of C13orf25 at chromosome 13 [32]. This cluster of miRNAs (miR-17/92) has a wide range of functions in immune cell development and differentiation [32]. Specifically, the miR-17 cluster has been found to promote T cell survival, regulate Th1 response, and interleukin 10 (IL10) production in regulatory T cells [33,34]. Members of two other clusters (miR-106a/363 and miR-106b/25) of the miR-17 family were also significantly downregulated in cases. miR-106a on the X chromosome has been shown to downregulate IL10 expression [35], and all three miR-17 family members differentially expressed in our study (106a, 106b, 17-5p) are known to be upregulated in activated T lymphocytes [36]. Second, three members of the miR-15 family (miR-15a-5p, miR-15b-5p, miR-16-5p) were also significantly downregulated in cases. Expression of miR-15 has been found to enhance the induction of regulatory T cells from naïve CD4+ T cells [37]. Additionally, miR-210 suppresses proinflammatory cytokines in murine macrophages [38]. Finally, in considering experimentally verified targets of the differentially expressed miRNAs in this study, the BCL2, CCND1, and VEGFA genes are each targeted by five miRNAs from at least two miRNA families. These genes are known factors in apoptosis and cell survival, having previously been studied in the context of cancer pathways [39,40]. Taken together, the downregulation of the above miRNAs in cases suggests a lower immune activation state. Whether this lower activation state has implications for CKD progression or complications requires further investigation.

Population-based studies have reported increased risk of infection among persons with early stages of CKD. In older adults (age ≥ 65 years), reduced kidney function (eGFR < 60 mL/min/1.73m2) is associated with increased risk of infection-related hospitalization or bloodstream infection [8,9]. In persons with ESRD, immune dysfunction due to metabolic disorder and retention of uremic solute is well established [10,11]. The decline in kidney function occurs in a continuum. Thus, immune dysfunction may play a role in earlier stages of CKD. The increased risk of infection in persons with earlier stages of CKD is consistent with our findings of the downregulation of some miRNAs that are critical for immune cell function.

On kidney disease, the plasma levels of miR-15b were found to be 2-fold lower in hemodialysis patients versus non-CKD controls [41]. This result is consistent with our study. The associations between levels of miR-15b and miR-17 in kidney tissues and acute kidney injury in human and in animal models have been reported [42]. Since gene expression levels are highly tissue specific, it is uncertain whether studies of gene expression in one tissue may be generalizable to a different tissue [43-47].

Strengths and limitations

One of the strengths of this study is its design. The cases and controls were carefully selected and matched for important potential confounders of gene expression, including hypertension medication use. In addition, we performed deep sequencing of the miRNA pool that allowed the discovery of novel miRNAs. However, our study has some limitations as well. First, the differentially expressed miRNAs need to be replicated in an independent study and validated by alternate laboratory methods, such as real-time polymerase chain reaction. Second, gene expression levels were measured in whole blood and not in specific types of immune cells. Thus, although cases and controls did not differ in major white blood cell types, and we used surrogate variables to control for unmeasured confounders, we cannot exclude the possibility that the differentially expressed miRNAs arise from differences in the distributions of white blood cell subpopulations. Third, our study is a cross-sectional study. We cannot distinguish whether the differential expression of miRNAs might have influenced disease development or was a consequence of the disease condition. Fourth, in contrast to miRNA co-expression network analysis, mRNA co-expression network analysis did not detect any significant modules although one module was nominally significant. This suggests larger sample size may be required for mRNA differential expression analysis. Finally, our cases and controls were under hypertension treatment, therefore our results cannot be generalized to persons without hypertension. Since the majority (>60%) of older adults in the U.S. are hypertensive [18], our results are, however, relevant to a large proportion of older adults.

Supplementary methods for mRNA and miRNA library preparation, sequencing, processing, alignment, and quantitation.

(DOCX)

Click here for additional data file.

A heat map of the normalized count of miRNAs with nominal p-value < 0.01 in differential expression analysis.

Hierarchical clustering was conducted using Euclidean distance and the Ward’s criterion [48].

(PDF)

Click here for additional data file.

Summary of mRNA and miRNA analyses and results in this study.

(PDF)

Click here for additional data file.

mRNA quality control and mapping statistics.

Mapping statistics, Quality Control statistics from PicardTools modules, and ERCC statistics.

(XLSX)

Click here for additional data file.

miRNA quality control and mapping statistics.

Mapping and Quality Control filtering statistics.

(XLSX)

Click here for additional data file.

Top genes from differential expression analysis using DESeq2 (p<1x10-4).

(XLSX)

Click here for additional data file.

Top transcripts from differential expression analysis using Ballgown (p<1x10-3).

(XLSX)

Click here for additional data file.

Association of miRNAs in the miR-17 and the miR-15 families with case status.

(XLSX)

Click here for additional data file.

Experimentally validated target genes for 11 differentially expressed miRNAs from miRTarBase.

Three out of 14 differentially expressed miRNAs were not available in the database.

(XLSX)

Click here for additional data file.

Association between experimentally validated target genes of the significant miRNA and case status.

(XLSX)

Click here for additional data file.

mRNA in each module from weighted gene co-expression network analysis.

(XLSX)

Click here for additional data file.

Association between the eigengene of each mRNA coexpression module and case status.

The eigengene is the first principal component of expression levels of a module.

(XLSX)

Click here for additional data file.

miRNAs in each correlated module from weighted gene co-expression network analysis.

(XLSX)

Click here for additional data file.

Top GO Terms for targets of differentially expressed miRNAs (one-tailed p<0.05).

These terms surpass nominal significance, but do not reach significance after BH correction of p-values for number of terms.

(XLSX)

Click here for additional data file.