Massive analysis of cDNA Ends (MACE) and miRNA expression profiling identifies proatherogenic pathways in chronic kidney disease.

Epigenetic dysregulation contributes to the high cardiovascular disease burden in chronic kidney disease (CKD) patients. Although microRNAs (miRNAs) are central epigenetic regulators, which substantially affect the development and progression of cardiovascular disease (CVD), no data on miRNA dysregulation in CKD-associated CVD are available until now. We now performed high-throughput miRNA sequencing of peripheral blood mononuclear cells from ten clinically stable hemodialysis (HD) patients and ten healthy controls, which allowed us to identify 182 differentially expressed miRNAs (e.g., miR-21, miR-26b, miR-146b, miR-155). To test biological relevance, we aimed to connect miRNA dysregulation to differential gene expression. Genome-wide gene expression profiling by MACE (Massive Analysis of cDNA Ends) identified 80 genes to be differentially expressed between HD patients and controls, which could be linked to cardiovascular disease (e.g., KLF6, DUSP6, KLF4), to infection / immune disease (e.g., ZFP36, SOCS3, JUND), and to distinct proatherogenic pathways such as the Toll-like receptor signaling pathway (e.g., IL1B, MYD88, TICAM2), the MAPK signaling pathway (e.g., DUSP1, FOS, HSPA1A), and the chemokine signaling pathway (e.g., RHOA, PAK1, CXCL5). Formal interaction network analysis proved biological relevance of miRNA dysregulation, as 68 differentially expressed miRNAs could be connected to 47 reciprocally expressed target genes. Our study is the first comprehensive miRNA analysis in CKD that links dysregulated miRNA expression with differential expression of genes connected to inflammation and CVD. After recent animal data suggested that targeting miRNAs is beneficial in experimental CVD, our data may now spur further research in the field of CKD-associated human CVD.

Introduction

Patients with chronic kidney disease (CKD) face a tremendous burden of cardiovascular morbidity, which cannot adequately be lowered by traditional therapeutic strategies. First evidence suggests that failure in epigenetic gene regulation may contribute to this high cardiovascular disease burden in patients with CKD.-

Among mechanisms of epigenetic regulation, miRNAs are centrally involved in the maintenance of cellular homeostasis and the control of inflammatory processes., miRNAs are small (~23 nt), endogenous and non-coding RNAs that regulate gene expression at the posttranscriptional level, either by degradation of target mRNAs or by inhibition of translation.

Importantly, one single miRNA can regulate a multitude of target genes and thus control specific pathways and signaling cascades such as inflammatory responses in endothelial cells and macrophages., Consequently, dysregulation of miRNA expression centrally contributes to the onset and development of diverse pathological conditions, including cardiovascular disease (CVD)., In line, distinct miRNAs have been proposed as useful biomarkers for diagnosis and prognosis of CVD.-

Against this background, great expectations are placed in the development of synthetic anti-miRs as novel treatment strategies in cardiovascular medicine.,- Such anti-miRs, which directly inhibit distinct miRNAs, have already been introduced as therapeutic agents in other disciplines of internal medicine.

However, before such therapeutic options can be transferred to patients suffering from cardiovascular disease in the context of CKD, more knowledge on regulation of inflammatory responses and atherogenesis by distinct miRNAs in uremia is needed.

Therefore, we aimed to characterize dysregulated miRNAs and associated signaling pathways in CKD, using high-throughput miRNA sequencing in combination with the recently introduced digital gene expression profiling technique MACE (Massive Analysis of cDNA Ends) as the latest advancement of tag-based gene expression analysis methods.

Results

Generation of MACE libraries

We performed genome-wide gene expression profiling of all HD patients (n = 10) and all control subjects (n = 10) by MACE in order to identify differentially expressed genes in HD patients. Therefore, RNA from PBMCs was used to generate 20 independent MACE libraries (10 “control libraries” and 10 “HD libraries”). In total 147 204 502 reads (= MACE tags) were generated across all libraries. After removal of duplicates, trimming of polyA-sequences and elimination of low-quality reads, of reads that were mapped to introns or intergenic regions, and of reads that could not be mapped at all, the total number of reads in coding regions was 51 167 601, comprising 21 308 090 reads from the control libraries and 29 859 511 reads from the HD libraries (Tables S1A and S1B). Across all libraries, these 51 167 601 reads accounted for 17 723 different genes (Table S2).

Differentially expressed genes in HD patients

Statistical data analysis (R) comprised hierarchical clustering of differentially expressed transcripts (Fig. 1A) as well as principal component analysis (Fig. 1B) and Pearson correlations (Fig. 1C) of all transcripts. We found HD patients and control subjects to cluster differentially; moreover, in line with their broad spectrum of comorbidity, HD patients displayed much higher heterogeneity in gene expression than control subjects.

Figure 1. Analysis of MACE results. (A) Unsupervised hierarchical cluster analysis with Euclidean distance measure of the 80 differentially expressed genes derived from MACE analysis. (B) Principle component analysis (PCA) of all expressed genes from the 20 examined samples. The first (x-axis) and second principal component (y-axis) accounted for 16% and 9%, respectively, of the total variation in the data. (C) Pearson product-moment correlation coefficient (PCC) for all samples compared within the control and patient group as well as between both groups.

The DESeq R/Bioconductor package was used to estimate the biological variance between replicates and to test for differential expression. We identified 80 differentially expressed genes (P < 0.05) in HD patients compared with healthy controls, of which 60 were upregulated and 20 were downregulated. Figure 2 visualizes differences in gene expression between HD patients and control subjects, whereas Table 1A and Table 1B list the top 15 upregulated as well as the top 15 downregulated genes in HD patients (based on their FDR [False Discovery Rate] adjusted P value, according to Benjamini and Hochberg).

Figure 2. Differences in gene expression between control subjects and hemodialysis patients. An MA plot was created to visualize the relation between log2 fold-change between controls and hemodialysis patients and log2 fold-average gene expression. Genes with a FDR (False Discovery Rate) adjusted P value < 0.05 are indicated in red.

Table 1A. Top 15 upregulated genes in hemodialysis patients

Gene symbol	Gene name	Control NEV	Dialysis NEV	P	FC	Biologic function
IER2	Immediate early response 2	262.5	525.8	1.8E-05	1.0	-
ZFP36	ZFP36 ring finger protein	689.3	1505.1	2.7E-05	1.1	Posttranscriptional regulation of tumor necrosis factor (TNF)
JUNB	Jun B proto-oncogene	707.3	1524.6	4.3E-05	1.1	Transcription factor, regulating gene activity following primary growth factor response
SOCS3	Suppressor of cytokine signaling 3	43.3	140.8	5.2E-05	1.7	Negative regulation of cytokines that signal through the JAK/STAT pathway
KLF6	Kruppel-like factor 6	1052.9	1546.8	0.000 12	0.6	Transcription factor
DUSP6	Dual specificity phosphatase 6	736.5	1461.1	0.000 17	1.0	Negatively regulate MAP kinases; inactivates ERK2
VIM	Vimentin	4174.7	6242.8	0.000 26	0.6	Class-III intermediate filament
PPP1R15A	Protein phosphatase 1, regulatory subunit 15A	170.1	315.4	0.000 38	0.9	Regulates growth arrest after stress; induces apoptosis by TP53 phosphorylation
NPC2	Niemann-Pick disease, type C2	136.0	210.4	0.000 55	0.6	Regulates transport of cholesterol
LMO2	LIM domain only 2 (rhombotin-like 1)	213.6	381.8	0.000 67	0.8	Hematopoietic development
FOS	FBJ murine osteosarcoma viral oncogene homolog	693.7	2118.4	0.001 41	1.6	Transcription factor; regulation of cell proliferation, differentiation, transformation
PLAGL2	Pleiomorphic adenoma gene-like 2	46.6	86.5	0.001 41	0.9	Transcriptional regulator
JUND	Jun D proto-oncogene	583.7	815.6	0.001 41	0.5	Transcription factor
FOSB	FBJ murine osteosarcoma viral oncogene homolog B	47.8	319.4	0.001 72	2.7	Transcription factor
EGR1	Early growth response 1	262.7	734.4	0.001 97	1.5	Transcriptional regulator; involved in differentitation and mitogenesis

NEV indicates normalized expression value; FC: fold change (log2[dialysis/control ratio]).

Table 1B. Top 15 downregulated genes in hemodialysis patients

Gene symbol	Gene name	Control NEV	Dialysis NEV	P	FC	Biologic function
CXCL5	Chemokine (C-X-C motif) ligand 5	33.6	15.7	0.001 37	-1.1	Inflammatory chemokine; neutrophil activation
DERL3	derlin 3	44.3	20.9	0.001 37	-1.1	Degradation of misfolded glycoproteins in the ER
ZAP70	Zeta-chain (TCR) associated protein kinase 70kDa	764.4	531.1	0.001 54	-0.5	T-cell development and lymphocyte activation
ATHL1	ATH1, acid trehalase-like 1	759.5	413.3	0.001 69	-0.9	-
ZBTB16	Zinc finger and BTB domain containing 16	61.8	27.3	0.007 64	-1.2	Zinc finger transcription factor; cell cycle progression
CARD11	Caspase recruitment domain family, member 11	263.9	181.2	0.013 72	-0.5	T-cell receptor mediated T-cell activation; activator of NF-kappa-B
ZC3H6	zinc finger CCCH-type containing 6	262.9	154.4	0.017 31	-0.8	-
PVRIG	Poliovirus receptor related immunoglobulin domain containing	237.5	148.8	0.024 90	-0.7	-
VPRBP	Vpr (HIV-1) binding protein	40.3	21.8	0.025 53	-0.9	Component of E3 ubiquitin-protein ligase complex; role in HIV infection
TBC1D10C	TBC1 domain family, member 10C	1389.0	1032.7	0.032 46	-0.4	Inhibitor of Ras signaling pathway and calcineurin
ZNF75A	Zinc finger protein 75a	49.4	31.6	0.034 53	-0.6	Transcriptional regulation
MZF1	Myeloid zinc finger 1	111.4	75.7	0.036 05	-0.6	Transcription regulator; hemopoietic development
RPL35	Ribosomal protein L35	5097.5	3846.3	0.036 43	-0.4	Ribosomal protein
ZBTB32	Zinc finger and BTB domain containing 32	10.3	3.8	0.03647	-1.5	Transcriptional regulator
MAN1A2	Mannosidase, α, class 1A, member 2	187.5	134.5	0.036 47	-0.5	Involved in the maturation of Asn-linked oligosaccharides

NEV indicates normalized expression value; FC: fold change (log2[dialysis/control ratio])

We next separately compared gene expression between HD patients with prevalent cardiovascular disease (CVD) and HD patients without prevalent CVD (Table S3). These subgroups differed in the expression of several genes. However, when applying the same strict statistical criteria as in the analysis of the total cohort, these differences lost statistical significance, given the small patient numbers in theses subgroups (n = 5).

Finally, we compared gene expression in healthy controls separately either to HD patients with CVD (Table S4), or to HD patients without prevalent CVD (Table S5). Differences in gene expression were more pronounced between healthy control subjects and HD patients with prevalent CVD (33 differentially expressed genes, Table S4), than between control subjects and HD patients without prevalent CVD (13 differentially expressed genes, Table S5).

Interaction network analyses

Next we generated interaction networks of gene products of differentially expressed genes in HD patients by using the String 9.0 software (http://string-db.org/; Fig. 3). Thereby, differentially expressed genes could be annotated to distinct pathways, comprising the T cell receptor signaling pathway (PAK1, RHOA, FOS, CARD11, ZAP70), the MAPK signaling pathway (DUSP1, IL1B, FOS, JUND, HSPA1A, DUSP6, PAK1), the Toll-like receptor signaling pathway (FOS, IL1B, MYD88, TICAM2), the chemokine signaling pathway (RHOA, PAK1, CXCL5, GNG5), or the Jak-STAT signaling pathway (SOCS3, PIM1, IFNGR1, CSF2RB). Notably, several differentially expressed genes are coding for transcription regulators such as bzip transcription factors (FOSB, JUNB, FOS), zinc finger proteins (KLF4, PLAGL2, EGR1) or proteins belonging to the small GTPase superfamily, Rho type (RHOA, RHOB).

Figure 3. Interaction networks of differentially expressed genes. Protein interaction networks were generated with the String 9.0 database using the 80 differentially expressed genes in hemodialysis patients.

Gene ontology analysis for biological and functional differences

To confirm these associations, we performed gene ontology (GO) enrichment analyses in order to annotate differential gene expression in HD patients to biological and functional alterations of immune cells: within the GO term “biological process,” the most significant differences were found for the immune system process (P = 0.003) including the subcategories immune system development (e.g., JUNB, EGR1, KLF4), regulation of immune system process (e.g., ZFP36, FOS, LMO2), and immune response (e.g., DUSP6, SOCS3, EGR1). Within the GO term “molecular function,” the most significant differences were found for nucleic acid binding transcription factor activity (P = 0.005) with the subcategory sequence-specific DNA binding transcription factor activity (e.g., JUNB, FOSB, KLF4).

Dysregulated atherosclerosis-related genes in HD patients

We next analyzed whether those 80 genes differentially expressed between HD patients and controls could be linked to cardiovascular disease and—given the central role of inflammation in CKD-associated atherogenesis—to infection / immune disease. The Genetic Association Database (accessible from the National Institutes of Health; http://geneticassociationdb.nih.gov) linked 22 out of these 80 genes to cardiovascular disease (e.g., KLF6, DUSP6, VIM, CXCL5, KLF4), and 34 genes to infection / immune disease (e.g., ZFP36, SOCS3, FOS, JUND, ZAP70) (Table S6A and B).

Generation of miRNA expression libraries

In order to analyze whether these changes in gene expression among HD patients can be linked to differential miRNA expression, we performed ultra-deep miRNA sequencing with pooled small RNA from HD patients and controls. Thus, we generated two independent miRNA expression libraries (“control library” and “HD library”). After eliminating low-quality reads and tags which were detected only once, the total number of reads was 10 967 102, comprising 5 964 297 reads (109 834 unique sequences, “UniTags”) from the control library and 5 002 805 reads (89 794 UniTags) from the HD library. Applying the bioinformatics pipeline omiRas, UniTags from the control library were mapped to 287 342 genomic loci, whereas UniTags from the HD library were mapped to 244 072 loci. The distribution of mapping regions and the relative frequency of non-coding RNA (ncRNA) classes is depicted in Figure S1.

Differentially expressed miRNAs in HD patients

In total, 650 different miRNAs were identified in the combined libraries, 182 of which were differentially expressed (P < 0.001) between HD patients and control subjects (Table S7). Of these 182 differentially expressed miRNAs, 75 were upregulated in HD patients, whereas 107 were downregulated. The top 15 upregulated and the top 15 downregulated (both based on p-value) miRNAs in HD patients are presented in Table 2A and Table 2B. Notably, several of these highly differentially expressed miRNAs, have previously been linked with cardiovascular disease such as miR-21, miR-26b, miR-146b, or miR-155. Table S8 provides an overview of current study findings in the context of earlier data on miRNA expression in human cardiovascular and kidney disease.

Table 2A. Top 15 upregulated miRNAs in hemodialysis patients

miRNA	Control NEV	Dialysis NEV	FC	P
hsa-miR-451a	5529.9	11 650.0	1.1	0
hsa-miR-21-5p	73 186.8	90 972.8	0.3	0
hsa-miR-223-3p	16 541.0	22 304.6	0.4	0
hsa-miR-144-3p	99.8	643.5	2.7	3.0E-306
hsa-miR-142-3p	22 816.2	27 078.9	0.2	2.0E-260
hsa-miR-144-5p	132.2	483.9	1.9	2.9E-149
hsa-miR-16-5p	27 536.5	30 696.7	0.2	3.6E-124
hsa-miR-146b-5p	15 657.3	17 844.5	0.2	1.9E-102
hsa-miR-15a-5p	4689.4	5875.9	0.3	8.9E-95
hsa-miR-424-5p	559.2	964.1	0.8	2.5E-77
hsa-miR-26b-5p	24 301.3	26 623.2	0.1	4.3E-77
hsa-miR-29b-3p	5484.4	6560.7	0.3	6.1E-69
hsa-miR-27a-5p	388.1	693.2	0.8	2.6E-62
hsa-miR-19b-3p	13 219.3	14 615.1	0.1	5.1E-51
hsa-miR-486-5p	702.1	1037.9	0.6	4.9E-47

NEV indicates normalized expression value; FC: fold change (log2[dialysis/control ratio]).

Table 2B. Top 15 downregulated miRNAs in hemodialysis patients

miRNA	Control NEV	Dialysis NEV	FC	P
hsa-miR-150-5p	71 417.0	56 863.9	-0.3	0
hsa-miR-140-3p	15 405.0	12 346.2	-0.3	9.3E-239
hsa-miR-342-3p	13 503.9	10 770.4	-0.3	8.2E-218
hsa-miR-423–3p	9818.1	7508.6	-0.4	3.0E-217
hsa-miR-320a	1801.9	944.5	-0.9	2.4E-190
hsa-miR-181a-5p	24 848.5	21 607.1	-0.2	3.5E-162
hsa-miR-126-3p	11 834.5	9670.2	-0.3	1.9E-154
hsa-miR-92a-3p	19 326.6	16 657.9	-0.2	1.9E-141
hsa-miR-30d-5p	19 128.0	16 619.5	-0.2	5.3E-126
hsa-miR-199a-3p	2899.5	2041.7	-0.5	1.1E-105
hsa-miR-23b-3p	4815.6	3730.1	-0.4	5.8E-98
hsa-miR-146a-5p	2922.6	2191.1	-0.4	9.8E-75
hsa-miR-155-5p	6880.2	5766.9	-0.3	3.1E-70
hsa-miR-708-5p	208.9	60.6	-1.8	1.7E-61
hsa-miR-151a-3p	2135.4	1590.7	-0.4	2.8E-57

NEV indicates normalized expression value; FC: fold change (log2[dialysis/control ratio]).

Regulation of differentially expressed genes by miRNAs

Finally, we aimed to analyze whether CKD-specific miRNA dysregulation may explain differences in gene expression between HD patients and controls. Using the module for combined analysis of miRNA and mRNA data of omiRas, which integrates information of 8 miRNA-mRNA interaction databases, we found 155 interactions between 68 differentially expressed miRNAs and 47 differentially expressed target genes (Table S9). The interaction networks are presented in Figure 4A . Importantly, genes that were upregulated in HD patients could be linked to miRNAs of which the expression was downregulated (Fig. 4A); in line, genes that were downregulated in HD patients could be linked to miRNAs with upregulated expression (Fig. 4B). Moreover, among those genes differentially expressed between HD patients and controls, 13 out of 22 genes linked to cardiovascular disease, and 20 out of 34 genes linked to infection / immune disease by Genetic Association Database analysis, could be connected to dysregulated miRNA expression (Tables S7 and S9).

Figure 4. Interaction networks of differentially expressed genes and miRNAs. (A) Interaction networks between upregulated genes and downregulated miRNAs in hemodialysis patients. (B) Interaction networks between downregulated genes and upregulated miRNAs in hemodialysis patients.

Discussion

Chronic kidney disease (CKD) patients suffer from a dramatically elevated cardiovascular event rate, which is mainly driven by accelerated arteriosclerosis. Traditional cardiovascular risk factors cannot fully explain this disease burden, and the implication of CKD specific risk factors is widely acknowledged.

We, and others, recently claimed that failure in epigenetic gene regulation centrally contributes to this high cardiovascular disease burden in CKD patients. However, earlier studies in this field were constrained to DNA methylation analysis, and the impact of miRNA dysregulation in CKD-associated atherosclerosis has not been analyzed until now, even though miRNAs are central regulators of virtually all processes in human (patho)physiology, including cardiovascular health and disease.-,

Against this background, we now aimed to identify dysregulated miRNAs in CKD by high-throughput sequencing and quantification of miRNAs. This allowed us to characterize 182 miRNAs differentially expressed between CKD patients and age- and gender-matched controls. Of note, this is the first report that applied next-generation technology for miRNA analysis in CKD. Compared with more conventional techniques such as qPCR and microarray platforms, next-generation sequencing has the advantage of not being biased by thermodynamics and thus allowing exact and digital quantification of miRNA expression.

To analyze the implications of miRNA dysregulation in the context of cardiovascular disease, we performed genome-wide gene expression analysis of all study participants by Massive Analysis of cDNA Ends (MACE). This allowed us to link miRNA dysregulation to differential gene expression, given that one single miRNA can simultaneously regulate a multitude of target genes and as one gene can be regulated by many different miRNAs.

MACE has recently been introduced as a digital gene expression method, in which each cDNA molecule is represented by one cDNA fragment (tag) of 94 bps, originating from a region 100–500 bps from the 3′ (poly-A) end of the transcript. High throughput sequencing of tags provides high resolution gene expression values and reveals differential expression of low-abundant transcripts beyond the scope of microarrays and RNaseq, if sequenced at similar depth. Moreover, the applied TrueQuant method for elimination of PCR copies ensured correct quantification of the low-complexity sequences.

Thereby, analysis of 147 204 502 reads identified 80 genes differentially expressed in CKD patients. Interaction and gene ontology analyses allowed us to specifically identify dysregulation within pathways linked to immune system processes as well as to inflammation (e.g., the Toll-like receptor signaling pathway, the MAPK signaling pathway, and the Chemokine signaling pathway). Additionally, according to the Genetic Association Database, many of these differentially expressed genes could be linked to cardiovascular disease and also to infection/immune disease. Finally, interaction network analysis of dysregulated miRNAs and target transcripts connected 68 differential expressed miRNAs with 47 reciprocal expressed target genes thus pointing toward a central role of miRNA dysregulation in CKD-associated inflammation and/or arteriosclerosis.

Beyond the field of nephrology, a substantial contribution of miRNA dysregulation to cardiovascular disease has widely been acknowledged in recent years. Here, most studies focused on arteriosclerosis, as the single most important pathophysiological cause of cardiovascular disease in humans. Arteriosclerosis has been characterized as an inflammatory disease, which is mainly driven by monocytes and monocyte-derived macrophages. Of note, miRNA networks centrally regulate the response of these cells during inflammation, and distinct miRNAs are differentially expressed in macrophages in response to certain inflammatory stimuli such as pro-inflammatory cytokines, TLR ligands, or oxLDL.-,- Among these, miR-21, which is related to key processes in atherogenesis and which is highly overexpressed in human atherosclerotic plaques compared with non-atherosclerotic arteries,, has been identified in the present study as one of the most upregulated miRNAs in CKD. Similar, in our CKD patients, we found upregulation of miR-34a, miR-146b, miR-26b, and miR-30e, all of which have been found to be overexpressed in atherosclerotic tissues from CAD patients in earlier studies.,,

Two earlier studies compared miRNA profiles of blood cells from patients with prevalent cardiovascular disease to healthy controls subjects (neither of which included patients with overt kidney disease). By real-time PCR profiling, Hoekstra and colleagues found an increased expression of miR-135a and a reduced expression of miR-147 in coronary artery disease (CAD) patients. A miRNA microarray approach by Taurino and colleagues identified upregulation of miR-140–3p and miR-182 among CAD patients. Notably, in the current analysis, we found upregulation of miR-182, but downregulation of miR-140–3p, among CKD patients.

Given their high stability, miRNAs can be detected and analyzed in any body fluid such as plasma or urine. Thus, circulating miRNA levels entered the field of emerging biomarkers in cardiovascular medicine, and a first epidemiologic study recently identified three miRNAs (miR-126, miR-223, and miR-197) as predictors of myocardial infarction.

In murine studies, Donners and colleagues recently demonstrated that hematopoietic deficiency of miR-155 enhances the development of atherosclerotic plaques in hyperlipidemic mice by limiting leukocyte recruitment to the atherosclerotic lesions. In this context, they found a skewed monocyte subset distribution with an increase in the Ly6C(hi) monocyte subset, which are considered the murine counterpart of human classical and intermediate monocytes. Interestingly, patients with CKD are characterized by profound shifts of monocyte subsets with high counts of intermediate monocytes (CD14++CD16+ monocytes),, which are discussed as major cellular contributors to atherosclerotic disease.- Notably, in our present study, miR-155 expression was strongly reduced in CKD patients, and further clinical studies will have to discern in how far distinct miRNA expression may contribute to monocyte subset biology in man.

Very few data on miRNA expression in CKD patients have been available before, which again mainly focused on circulating miRNAs. Neal and colleagues found levels of circulating miRNAs—among them the miR-155—to be strongly reduced in CKD patients. In patients with acute kidney injury, elevated plasma levels of miR-210 were characterized as independent predictor of mortality. Of note, due to different underlying mechanisms of biogenesis, data on circulating miRNA should not uncritically be compared with cellular miRNA expression. Moreover, although circulating miRNAs may undoubtedly serve as valuable biomarkers for distinct pathological conditions, the characterization of cellular miRNA profiles may allow a more comprehensive understanding of epigenetic pathways in pathophysiological conditions such as atherosclerosis. Such understanding is further amplified by interaction network analysis, as done in the present manuscript, which provides first evidence that miRNA networks coordinate the expression of inflammatory related target genes in CKD. These network analyses underscore the notion that single miRNAs regulate multiple genes within specific signaling pathways, leading to a profound dysregulation of downstream gene networks when miRNA expression is disturbed.

In summary, we provide a first comprehensive miRNA analysis in CKD which allows to link dysregulated miRNA expression with differential expression of genes connected to inflammation and atherosclerosis. Since targeting of distinct miRNAs has been successfully applied in other fields of internal medicine and since animal studies suggested similar therapy options against atherosclerotic vascular disease,,, our data may hopefully spur further research on epigenetic dysregulation in CKD-associated cardiovascular disease.

Patients and Methods

Study participants

Ten clinically stable male CKD patients (56.1 ± 3.9 y) undergoing standard hemodialysis (HD) therapy thrice a week were recruited from the Department of Internal Medicine IV, Nephrology and Hypertension of the Saarland University Medical Center. From each patient, 20 ml EDTA (EDTA) anticoagulated blood was drawn before a dialysis session after the long-interdialytic interval. Five CKD patients had prevalent cardiovascular disease, defined as the presence of coronary artery disease (prior myocardial infarction or coronary revascularization), cerebrovascular disease (prior stroke or carotid revascularization), and / or peripheral artery disease (prior revascularization of lower-limb arteries). The remaining five CKD patients had no prevalent cardiovascular disease, which allowed post-hoc comparisons of gene expression between subgroups of CKD patients.

Ten healthy sex- and age-matched employees from Saarland University Medical Center served as controls. Mean age of controls was 53.5 ± 2.4 y, and their mean estimated glomerular filtration rate was 86.0 ± 14.3 ml/min/1.73 m2. No control subject had prevalent cardiovascular disease.

All participants gave written informed consent in accordance with the Declaration of Helsinki. The study protocol was approved by the local Ethics Committee.

RNA isolation and preparation of MACE (massive analysis of cDNA ends) and miRNA libraries

Peripheral blood mononuclear cells (PBMCs) were immediately isolated from anticoagulated blood by Ficoll-Paque (Lymphocyte Separation Medium; PAA) gradient density centrifugation. Total RNA was isolated from peripheral blood mononuclear cells using the NucleoSpin® miRNA kit (Macherey Nagel), which allows to separate the large and small fraction of the total RNA. The large fraction of the total RNA was used for the preparation of the MACE libraries, the small RNA fraction for miRNA library preparation.

For preparation of MACE libraries, poly-adenylated mRNA was isolated from 1 µg of the large fraction of total RNA using Dynabeads® mRNA Purification Kit (Life Technologies GmbH) and cDNA was produced by first and second strand synthesis using the SuperScript® III First-Strand Synthesis System (Life Technologies GmbH), but with modified, barcoded poly-T adapters that are suitable to bind on the Illumina Hiseq2000 flow cell and that are biotinylated at the very 5′ end. After cDNA production, the biotinylated cDNA was random-fragmented to reach an average size of 250 bps, the 3′ends were captured by streptavidin beads and GenXPro’s TrueQuant adapters were ligated. Ten barcoded samples were sequenced simultaneously in one lane of an Illumina Hiseq2000 with 1 × 100 bps.

For preparation of miRNA libraries, small RNA (17–30 bps) was size-selected from total RNA using the flashPAGE™ Small Nucleic Acid Purification System (Life Technologies). TrueQuant Adapters were ligated directly to the small RNA, basically as described by Hafner et al. The miRNA population containing flanking p5 and p7 adapters was sequenced on a Hiseq2000 machine (Illumina). The MACE and miRNA library preparation was performed by GenXPro GmbH.

Data Analysis

Quantification of mRNA expression

MACE reads in each patient library were PolyA-trimmed and the marked quality region was removed. Mapping of reads to the human genome (hg19) was performed with novoalign (http://novocraft.com/). The mapping coordinates of each read were overlapped with the refseq annotation track from the UCSC table browser (http://genome.ucsc.edu/cgi-bin/hgTables?command=start) to quantify mature mRNA expression. Normalization and test for differential expression was performed in the statistical programming language R (www.r-project.org/) with the DEseq package.

GO-enrichment analysis was performed using the GO-enrichment toolkit from http://genxpro.ath.cx, which is based on the Fisher’s exact test among transcripts that are differentially expressed at a P value of 0.05.

Quantification of miRNA expression

Small RNA-Seq libraries were analyzed with omiRas (detailed workflow is presented at http://tools.genxpro.net/omiras). For each small RNA-Seq library, data processing started with 3′ adaptor clipping by a local alignment of the adaptor sequence to each read. Subsequently, Illumina's marked quality region was removed and the reads were summarized to UniTags. Singletons were removed from the data set and the remaining tags were mapped to the human genome (hg19) with bowtie. The mapped tags were annotated with the help of various models of coding and non-coding RNAs retrieved from the UCSC table browser. Tags mapping to exonic regions of coding genes were excluded from further analysis. Non-coding RNAs were quantified in each library. For tags mapping to multiple genomic loci the number of reads corresponding to the tag was divided by the number of mapping loci.

To account for differences in sequencing depth, tag-counts were normalized (NEV, normalized expression value) and differential expression was detected with the DEGseq bioconductor package.