Profiling and initial validation of urinary microRNAs as biomarkers in IgA nephropathy

Statement on ethics

The study was carried out in accordance with the Declaration of Helsinki for Human Research and approved by the Ethics Committee of the Chinese PLA General Hospital (approval number S2014-004-02). Written informed consent for inclusion was obtained from each participant.

Study design

A screening and a validation phase were designed for the present study. Besides healthy controls, we recruited patients with membranous nephropathy (MN) or minimal change disease (MCD) as other glomerulopathy controls. In the screening phase, 32 urine samples collected from 18 patients with IgAN and 14 controls were subjected to Affymetrix GeneChip miRNA 4.0 Array to identify miRNAs that were significantly differentially expressed. In the validation phase, the expression of different miRNAs was validated by real-time quantitative polymerase chain reaction (RT-qPCR) in samples from 102 patients with IgAN, 41 patients with MN, 27 patients with MCD and 34 healthy controls.

Enrollment of subjects

Patients who received renal biopsy for the first time from December 2013 to November 2014 at the Nephrology Department of the Chinese PLA General Hospital were consecutively recruited for the present study. The diagnosis was confirmed based on clinical findings and renal biopsies. A total of 322 consecutive patients were diagnosed as IgAN, MN or MCD. Among them, we excluded 41 patients because they were <20 or >50 years old at renal biopsy, 26 patients because they had been treated with steroids within 6 months prior to renal biopsy, and 7 patients because there were less than eight glomeruli in light microscopy sections. We also excluded patients with a concomitant diagnosis of chronic hepatic disease (n = 14), diabetes (n = 11), urinary, respiratory or gastrointestinal tract infection (n = 7), rheumatoid arthritis (n = 3), Henoch-Schonlein purpura (n = 7) and systemic lupus erythematosus (n = 10).

Healthy volunteers who were between 20 and 50 years old were recruited at the Physical Examination Center of the same hospital. All of them had normal renal function, normal urinalysis results, and no personal or family history of nephropathy or other serious illness.

Finally, 120 patients with IgAN, 45 patients with MN, 31 patients with MCD and 40 healthy volunteers were enrolled in this study. Among the IgAN patients, the first 6 patients with grade I–II according to Lee’s grading system who met the inclusion criteria and a sex ratio of 1–1 were selected as a screening cohort, as well as 6 patients with grade III and 6 patients with grade IV–V. In addition, 14 age- and gender-matched controls which included 4 patients with MN, 4 patients with MCD and 6 healthy volunteers were selected in the screening cohort.

The demographic and clinical data, such as age, gender, mean aortic pressure (MAP), serum creatinine (Scr) and 24-hour urinary protein excretion (UPE), of all included participants were recorded at the time of kidney biopsy. The estimated glomerular filtration rate (eGFR) was calculated using the Chronic Kidney Disease Epidemiology Collaboration (CKD–EPI) equations (Levey et al., 2009).

Evaluation of renal pathological lesion

Renal tissues obtained by biopsy were stained with haematoxylin and eosin (H &E), periodic acid Schiff (PAS), periodic acid silver methenamine (PASM) and Masson’s trichrome for light microscopy and with IgG, IgA, IgM, C3 and C1q for immunofluorescence. Two pathologists who were blinded to patients’ data evaluated the renal biopsy slides separately.

For patients with IgAN, the severity of renal pathological lesions were evaluated according to Lee’s grading system (Lee et al., 1982) and Oxford classification (Roberts et al., 2009). First, all cases of IgAN were categorized into grade I–V according to Lee’s grading system. Then, four variables of Oxford classification were further evaluated; that is, mesangial proliferation (M1—more than half the glomeruli have more than three cells in a mesangial area, M0—not meet the criteria of M1), endocapillary hypercellularity (E1—present, E0—absent), segmental glomerulosclerosis (S1—present, S0—absent ), and tubular atrophy/interstitial fibrosis (T0—0%–25%, T1—26%–50%, T2—>50%).

Urine sample preparation

A whole-stream early morning urine specimen was collected at the day of renal biopsy. The urine sample was centrifuged at 3,000 g for 30 min and at 13,000 g for 5 min at 4 °C; supernatant was discarded (Wang et al., 2012; Wang & Szeto, 2013). Then the urinary cell pellet was lysed by RNA lysis buffer MZ (catalog number DP501; Tiangen Biotech, Beijing, China) and stored at −80 °C until used.

RNA extraction

Total RNA was extracted using a miRcute miRNA Isolation Kit (catalog number DP501; Tiangen Biotech, Beijing, China) according to the manufacturer’s protocol. The quantity (ng/ml) and purity (ratio of absorbances of the RNA isolates at 260 nm and 280 nm [A260/280]) of the RNAs obtained was evaluated by NanoDrop 2000 spectrophotometer (Thermo Scientific, Waltham, Massachusetts, USA).

MiRNA microarray analysis

Affymetrix GeneChip miRNA 4.0 Array (Affymetrix, Santa Clara, California, US, catalog number 902411), which covers all of the 2,578 mature human miRNAs available in miRBase version 20 (June 24, 2013, www.mirbase.org/) (Kozomara & Griffiths-Jones, 2011), was used to profile miRNA expressions. Briefly, 1 µg of RNA was polyA tailed and labelled with a FlashTag Biotin HSR RNA Labeling Kit (Affymetrix, Santa Clara, California, US, catalog number 901911). The labelled RNA was hybridized at 48 °C for 16 h on the miRNA arrays, which were washed and stained with Affymetrix Fluidics Station 450 (Affymetrix, Santa Clara, California, US) and scanned with an Affymetrix GeneChip Scanner 3000 (Affymetrix, Santa Clara, California, US) using the Command Console software (Affymetrix, Santa Clara, California, US). The data were analyzed with miRNA QCTool (Affymetrix, Santa Clara, California, US) using the Affymetrix default analysis settings and quantile as the normalization method. The microarray data have been deposited in Gene Expression Omnibus (GEO) at NCBI (http://www.ncbi.nlm.nih.gov/geo/) with accession number GSE64306.

RT-qPCR analysis

RT-qPCR was carried out in compliance with the MIQE guidelines (Bustin et al., 2009) to verify the candidate miRNAs revealed by microarray. Briefly, reverse transcription was performed using 500 ng of total RNA and a miRcute miRNA First-Strand cDNA Synthesis Kit (catalog number KR201; Tiangen Biotech, Beijing, China) according to the manufacturer’s protocol. RT-qPCR was performed on an Applied Biosystems 7500 Fast Dx Real-Time PCR Instrument (Applied Biosystems, Carlsbad, Calfornia, USA) using a miRcute miRNA qPCR Detection Kit (catalog number FP401; Tiangen Biotech, Beijing, China), according to the manufacturer’s instructions. All of the primers were purchased from Tiangen Biotech Company (Beijing, China). All PCR reactions were performed in triplicate, followed by melt curve analysis to verify their specificity and identity. U6 was selected as the endogenous reference control (Mestdagh et al., 2009). Relative miRNA expression levels were calculated using the ΔΔCt method as previously described (Livak & Schmittgen, 2001).

MiRNA target prediction and function analysis

The candidate miRNAs were imported into the miRWalk algorithm (Dweep et al., 2011) (http://www.umm.uni-heidelberg.de/apps/zmf/mirwalk/), and the prediction of their target genes was performed using nine additional algorithms, including TargetScan, DIANAmT, miRanda, miRDB, RNAhybrid, PICTAR4, PICTAR5, PITA and RNA22.

To explore the functional annotation and pathway enrichment of those predicted genes, the Gene Ontology (GO) and KEGG (Kyoto Encyclopedia of Genes and Genomes) database analyses were conducted using a DAVID (the Database for Annotation, Visualization and Integrated Discovery) online analysis tool (Huang da, Sherman & Lempicki, 2009) (http://david.abcc.ncifcrf.gov/tools.jsp), in which we focused on the GO biology process (BP) feature.

Statistical analysis

The microarray data were analyzed by the algorithm of SAM (significance analysis of microarrays, www-stat.stanford.edu/~tibs/SAM/) (Olson, 2006). The false discovery rate (FDR) was set to <0.05 and the minimum fold change (FC) was set to >2.0 or <0.5. Hierarchical clustering was carried out using the MeV 4.9 software (Multi Experiment Viewer, http://www.tm4.org/mev.htm) to generate both miRNA and sample trees based on Pearson correlation. For the RT-qPCR data, statistical analysis was performed by SPSS version 20.0 (SPSS, USA). Data were compared by Mann–Whitney U test, Kruskal–Wallis test and Spearman’s rank order correlations as appropriate. In bioinformatics analysis, the Fisher’s exact test and χ² test were used to select the significant GO category or KEGG pathway, and the FDR was calculated to correct the P value. A P-value <0.05 was considered to be statistically significant.

Patients characteristics

Demographic and clinical characteristics of the screening and validation cohorts at the time of renal biopsy are provided in Tables 1 and 2 respectively. There were no significant differences in age and MAP among different groups in both cohorts.

In the screening cohort, the levels of UPE in IgAN grade I–II and III subgroups were lower than that in MCD group (all P < 0.01) and higher than that in healthy control (HC) group (all P < 0.01); the levels of Scr in IgAN grade IV–V subgroup were higher than that in both MN and HC group (all P < 0.05), whereas the levels of eGFR showed an opposite trend.

In the validation cohort, the levels of UPE in IgAN group were lower than that in MN and MCD group (all P < 0.01), but higher than that in HC group (P < 0.001); the levels of Scr were higher than that in all control groups (all P < 0.01), whereas the levels of eGFR showed an opposite trend. The histopathological parameters are shown in Table 3.

The yield of RNA from urinary sediments

RNA yield from urinary sediments (ng of RNA isolated from per ml of urine) was determined and showed no significant differences among patients with IgAN, MN, MCD and healthy subjects. The overall median was 22.73 ng/ml (interquartile range [IQR] = 10.75–40.35). A260/280 was also recorded and showed no significant differences among different groups in the present study. The overall mean of A260/280 was 1.95 (range = 1.82–2.07, SD = 0.06).

Differentially expressed miRNAs in the screening phase

Global miRNA profiling from microarray analysis

To identify miRNAs differentially expressed in IgAN, miRNA microarrays were used to analyze their global expression profile of the four groups (IgAN, MN, MCD and HC). The 32 urine samples of the screening cohort were profiled on Affymetrix GeneChip miRNA 4.0 Arrays.

Among 2,578 mature human miRNAs represented on the microarray, 780 (30.3%) miRNAs were identified as expressed in at least one urine sample. Unsupervised hierarchical clustering analysis demonstrated that miRNA profiling clearly differentiated IgAN samples from healthy controls; there were, however, one MN sample and two MCD samples mixed with IgAN, indicating that the differences between IgAN patients and disease controls were not as clear (Fig. 1).

Differentially expressed miRNAs in IgAN subgroups

In the screening phase, we grouped IgAN grade I–II as the early pathological change group; grade III as the mild pathological change group; and grade IV–V as the severe group according to Lee’s grading system. Then we compared the different subgroups of IgAN (grade I–II, III and IV–V) with both healthy and disease controls. The differentially expressed miRNAs are shown in the form of volcano plot (Fig. S1) and Venn diagram (Fig. 2). All these miRNAs presented a FC >2 or <0.5 and P value <0.05.

In the IgAN grade I–II subgroup, two miRNAs (hsa-miR-223-3p, hsa-miR-629-5p) were in high levels and two (hsa-miR-3613-3p, hsa-miR-4668-5p) were in low levels compared to all of the control groups (Fig. 2A). In the IgAN grade III subgroup, three miRNAs (hsa-miR-150-5p, hsa-miR-572 and hsa-miR-371b-5p) were in high levels and three (hsa-miR-3613-3p, hsa-miR-4668-5p and hsa-miR-6750-5p) were in low levels compared to all control groups (Fig. 2B). IgAN grade I–II and III subgroups shared dysregulation of only two miRNAs (miR-3613-3p and miR-4668-5p). Interestingly, there were no overlaps in changes to miRNAs for IgAN grade IV–V subgroup as shown in the Venn diagram (Fig. 2C). It turned out that miRNAs were in different levels only in specific IgAN subgroups.

Analysis of candidate miRNAs in the validation phase

Validation of candidate miRNAs using RT-qRCR

For confirmation purposes, the two miRNAs (miR-3613-3p and miR-4668-5p) that were differentially expressed in both IgAN grade I–II and III subgroups in the screening phase were selected as candidates and analyzed in the validation cohort, which included 102 patients with IgAN, 41 patients with MN, 27 patients with MCD and 34 healthy controls.

RT-qPCR results showed that the levels of miR-3613-3p were significantly lower in IgAN as compared with that in HC, MN and MCD, respectively [0.27 (0.19–0.38) vs. 1.10 (0.9–1.26), 0.58 (0.32–1.00) and 0.61 (0.39–0.77), respectively, all P < 0.01] (Fig. 3A); the levels of miR-4668-5p were significantly lower in IgAN than that in HC [0.54 (0.38–0.60) vs. 1.11 (0.73–1.22), P = 0.001]; there were no significant differences for miR-4668-5p, however, between IgAN and the disease controls (Fig. 3B).

Correlations between candidate miRNAs and clinical parameters

Correlations between urinary levels of candidate miRNAs and clinical parameters, such as UPE, Scr, eGFR and Lee’s grades, were further analyzed within the IgAN group (Table 5). Urinary levels of miR-3613-3p were positively correlated with both UPE (r = 0.50, P = 0.034) and Lee’s grades (r = 0.57, P = 0.014), while negatively correlated with eGFR (r = − 0.48, P = 0.043). Similarly, significant correlations were found between urinary levels of miR-4668-5p with eGFR (r = − 0.50, P = 0.034) and Lee’s grades (r = 0.57, P = 0.013).

Table 5:

Correlations between candidate miRNAs and clinical parameters in patients with IgA nephropathy.

	miR-3613-3p			miR-4668-5p
	r	95% CI	P value	r	95% CI	P value
Age	0.12	−0.37∼0.54	0.631	0.26	−0.28∼0.67	0.29
MAP	0.42	−0.06∼0.72	0.084	0.34	−0.19∼0.73	0.173
UPE	0.501	0.03∼0.83	0.034	0.44	−0.05∼0.76	0.066
Scr	0.451	−0.06∼0.79	0.06	0.42	−0.06∼0.73	0.083
eGFR	−0.48	−0.84∼0.02	0.043	−0.5	−0.84∼0.08	0.034
Lee’s grades	0.57	0.07∼0.84	0.014	0.57	0.04∼0.85	0.013

DOI: 10.7717/peerj.990/table-5

Notes:

r

Spearman correlation coefficient

95% CI

95% confidence interval

MAP

mean aortic pressure

UPE

24-hour urinary protein excretion

Scr

serum creatinine

eGFR

estimated glomerular filtration rate

In addition, relationships between expressions of candidate miRNAs and four components of Oxford classification are shown in Fig. 4. For segmental glomerulosclerosis, patients scored as S0 had significantly lower levels of urinary miR-3613-3p and miR-4668-5p than those scored as S1 [0.23 (0.12–0.33) vs. 0.56 (0.29–2.50), 0.23 (0.18–0.25) vs. 0.53 (0.35–1.90), respectively, all P < 0.05].

Bioinformatics analysis of candidate miRNAs

Functions and pathway analysis of miRNAs targets

To explore in-depth biological information, the GO and KEGG database analyses were conducted using a DAVID online analysis tool. All of the significant GO Biological Processes for miR-3613-3p and miR-4668-5p are listed in Table S4 and Table S5 respectively; the top 20 GO Biological Processes according to P values and all significant KEGG pathways for miR-3613-3p and miR-4668-5p are shown in Fig. 5 and Fig. 6 respectively.

There were 254 and 49 significant GO terms for miR-3613-3p and miR-4668-5p respectively (P < 0.05). The most significantly enriched GO terms for miR-3613-3p were the regulation of transcription (GO: 0045449, GO: 0006350) and macromolecule metabolic process (GO: 0010604, GO: 0010557) as shown in Fig. 5A; it also had significant enrichment in B cell differentiation (GO: 0030183, P = 0.006) and activation (GO: 0042113, P = 0.016) as listed in Table S4. Similarly, the most significantly enriched GO terms for miR-4668-5p (as shown in Fig. 6A) also included transcription (GO: 0006350, GO: 0045449) and macromolecule metabolic process (GO: 0051252, GO: 0010558), as well as nuclear transport (GO: 0006607, GO: 0051169).

In addition, as shown in Figs. 5B and 6B, there were 18 and 2 significant KEGG pathways for miR-3613-3p and miR-4668-5p respectively (P < 0.05); both included Wnt signaling pathway (hsa04310, P = 0.012 and P = 0.013 respectively).