Small RNA deep sequencing reveals novel miRNAs in peripheral blood mononuclear cells from patients with IgA nephropathy.

Peripheral blood mononuclear cells (PBMCs) contribute to the deposition of immunoglobulin A (IgA) and progression of IgA nephropathy (IgAN). This study was performed to identify novel microRNAs (miRNAs/miRs) associated with IgAN. Small RNAs were isolated from PBMCs collected from 10 healthy participants and 10 patients with IgAN; the RNAs were then subjected to high‑throughput small RNA sequencing. The results showed that miRNAs constituted 70.33 and 69.83% of small RNAs in PBMCs from healthy participants and patients with IgAN, respectively. In total, 44 differentially expressed miRNAs were identified, of which 34 were upregulated and 10 were downregulated. Among these differentially expressed miRNAs, most showed novel associations with IgAN, except miR‑148a‑3p, miR‑184 and miR‑200a. Furthermore, Kyoto Encyclopedia of Genes and Genomes pathway analysis revealed that the target genes of the differentially expressed miRNAs were primarily enriched in cancer pathways, the PI3K‑Akt signaling pathway and MAPK pathways, all of which control cell proliferation and gene expression. Moreover, miR‑3121‑3p, miR‑203a‑3p and miR‑200a‑3p may regulate core 1 synthase, glycoprotein‑N‑acetylgalactosamine 3‑β‑galactosyltransferase 1 (C1GALT1) expression by binding to its 3' untranslated region. In conclusion, 44 differentially expressed miRNAs were discovered, 41 of which were newly found to be associated with IgAN. The differentially expressed miRNAs may regulate the progression of IgAN by controlling the behavior of PBMCs or deposition of IgA via targeting of signaling pathways or expression of C1GALT1. These findings may provide a basis for further research regarding IgAN diagnosis and therapy.

Introduction

Immunoglobulin A (IgA) nephropathy (IgAN) is the most prevalent form of primary glomerulonephritis and frequently causes end-stage renal disease, to the extent that 30–40% of patients with biopsy-proven IgAN will progress to end-stage renal disease in 20 years (1). IgAN is characterized by highly heterogeneous clinical features, ranging from asymptomatic microscopic hematuria to a rapidly progressive form of glomerulonephritis (2). Renal biopsy is the initial method of IgAN diagnosis and disease assessment (1). Because of its invasiveness, renal biopsy is difficult, and repeated procedures pose considerable inconveniences for patients. IgA1 is the main subclass of IgA; however, the sensitivity and specificity of serum IgA1 and common biomarkers (e.g. serum creatinine, glomerular filtration rate, and proteinuria) are insufficient for IgAN diagnosis. Thus, there is an urgent need to investigate the mechanism underlying the development of IgAN, and to screen for disease biomarkers and targets of IgAN for diagnosis and therapy (2,3).

Currently, multiple biomarkers for IgAN are known, such as proteins and microRNAs (miRNAs/miRs) in blood and urine (3). miRNAs constitute a class of non-coding RNAs that regulate gene expression by binding to complementary sequences in the coding portion or 3′ untranslated regions (UTRs) of mRNAs (4). miRNAs have been shown to exhibit functional dysregulation in diseases such as cancer, cardiopathy, and nephropathy (5). miRNAs are not easily degraded by RNase and are stable in blood (6); therefore, they have previously been considered to be good biomarkers for diagnosis. Previously, urinary levels of miR-146a, miR-155, miR-200a, miR-200b and miR-429 were shown to be associated with clinical and histological severity of IgAN (7,8). miR-148b and miR-let7b are known to promote the accumulation of IgA1 in IgAN, and are potential serum biomarkers for the detection of primary IgAN (9–11). However, there is not much research into the association between miRNAs and IgAN in blood, and the roles of various miRNAs remain unclear in the pathogenesis of IgAN.

Peripheral blood mononuclear cells (PBMCs) are known to participate in IgA deposition and IgAN development, and previous studies have shown that miRNAs in PBMCs could promote the progression of IgAN (9,10,12,13). The miRNA expression profile of PBMCs in patients with IgAN has been studied by using a miRNA microarray, and miR-148b was found to be a potential promoter of IgAN (9). In the present study, small RNAs from PBMCs were sequenced to identify novel miRNAs associated with IgAN. Furthermore, the potential roles of differentially expressed miRNAs were investigated through bioinformatics analysis. In total, 44 differentially expressed miRNAs were identified, and this study is the first to report an association for most of these miRNAs with IgAN, indicating that these miRNAs may serve as biomarkers or targets for IgAN.

Materials and methods

Patients

In this study, 10 healthy participants and 10 patients with IgAN were enrolled during the period from January 2018 to December 2018 in The Second Affiliated Hospital of Guangzhou University of Chinese Medicine. The admission criteria of IgAN nephropathy patients were as follows: i) IgAN was confirmed by renal biopsy; ii) normal renal function (blood creatinine, blood urea nitrogen); iii) ≥18 years old; iv) no previous hormone, immunosuppressant or kidney transplantation treatments; v) the patients did not experience secondary IgAN such as purpura nephritis, lupus nephritis and hepatitis B-related nephritis. The exclusion criteria for healthy participants were as follows: i) Those with chronic diseases, such as coronary heart disease, hypertension, acute and chronic cerebrovascular disease or diabetes; ii) those with infectious diseases and fever; iii) those with mental illnesses; iv) those with any metabolic syndromes. The research protocol was approved by the Medical Ethical Committee of Guangdong Provincial Hospital of Chinese Medicine, and all patients provided written informed consent. All patients were diagnosed with IgAN by an experienced pathologist, and the clinical features are shown in Table I.

Table I.

Demographic and baseline clinical data of the healthy participants and patients with IgAN.

Characteristic	IgAN	Healthy control
Number of cases	10	10
Sex, male:female	5:5	5:5
Age, years	29.30±4.24	27.2±9.31
IgA, g/l	3.53±1.11	2.36±0.52
Serum creatinine, µmol/l	74.40±20.04	80.8±14.48
GFR, ml/min/1.73 m2	108.26±19.52	102.01±10.56

Data are presented as the mean ± SD. IgA, immunoglobin A; IgAN, IgA nephropathy; GFR, glomerular filtration rate.

PBMC small RNA isolation

Human PBMC isolation buffer (TBD sciences; Tianjin Haoyang Biological Products Technology Co. Ltd.) was used to isolate PBMCs from 2 ml peripheral blood samples of patients, in accordance with the manufacturer's protocol. TRIzol® (Invitrogen; Thermo Fisher Scientific, Inc.) was used to isolate total RNA from PBMCs, in accordance with the manufacturer's protocol. RNA molecules of 18–40 nucleotides in length were collected for the analysis of small RNAs.

Small RNA sequencing and differentially expressed miRNA screening

Small RNA library preparation was performed using the NEB Next® Multiplex Small RNA Library Prep Set for Illumina® (New England BioLabs, Inc.). Briefly, RNA molecules were ligated with a 3′ RNA adapter, followed by ligation with a 5′ adapter. Subsequently, adapter-ligated RNAs were subjected to reverse transcription into cDNA at 50°C for 60 min and amplified. PCR was performed using the following thermocycling conditions: Initial denaturation at 94°C for 30 sec; followed by 15 cycles of denaturation at 94°C for 15 sec, annellation at 62°C for 30 sec, extension at 70°C for 15 sec; and a final extension at 70°C for 5 min. The PCR products were then size-selected by polyacrylamide gel electrophoresis, in accordance with the manufacturer's instructions. Purified library products were evaluated using the Agilent 2200 TapeStation (Agilent Technologies, Inc.) and diluted to 2 pM for cluster generation in situ on the HiSeq 2500 single-end flow cell (Illumina, Inc.), followed by sequencing on the HiSeq 2500 at 1×50 base pairs.

For miRNA analysis, clean reads were mapped to hg19. miRNAs and Piwi-interacting (pi)RNAs were identified using miRBase version 21 (www.mirbase.org) (14) and pirnabank (pirnabank.ibab.ac.in) (15), respectively. Transfer (t)RNA, ribosomal (r)RNA, small nuclear (sn)RNA and small nuclear non-coding (sno)RNA data were from Rfam12.1 (rfam.xfam.org) (16). The software edgeR (17) was used to compare the expression of miRNAs between PBMCs from healthy participants and those from patients with IgAN and to calculate the P-value. miRNA reads were normalized to miRNA reads per million reads, and the average values were examined. miRNAs with a fold change ≥2 and P<0.05 were considered to be indicate a statistically significant difference in the expression profile (Table SI).

Kyoto Encyclopedia of Genes and Genomes (KEGG), Gene Ontology (GO) and miRNA-mRNA regulatory networks analyses

KEGG and GO functional annotations were used to highlight the most overrepresented GO terms and KEGG pathways that were closely associated with top 10 differentially expressed miRNAs. KEGG pathways analysis was based on the KEGG database (http://www.genome.jp/) (18) to obtain the signal transduction and disease pathway annotation information for the candidate target genes. The P-value was calculated by a Fisher's exact test, and the signal transduction and disease pathways with statistical significance were compared with the background were obtained with P<0.05 as the significant threshold. GO functional analysis was based on the GO resource (http://geneontology.org) (19) to provide molecular function, biological process and cellular component annotation for candidate target genes. P<0.05 was considered to indicate a statistically significant difference of the enrichment degree of GO annotation. miRNA-target gene network analysis was performed using the mirWalk program (http://mirwalk.umm.uni-heidelberg.de/) (20).

miRNA binding sites prediction on the C1GALT1 3′ UTR

The potential miRNAs that may regulate C1GALT1 were predicted using TargetScan (http://www.targetscan.org) (21) and miRDB (http://mirdb.org) (22).

Results

Small RNA sequencing and analysis in healthy participants and patients with IgAN

To study the aberrant expression of miRNAs in the PBMCs of patients with IgAN, small RNAs were isolated and cDNA libraries were prepared. The cDNA libraries were sequenced, and the data were analyzed following the workflow shown in Fig. 1A. Small RNAs isolated from the PBMCs of healthy participants and patients with IgAN were primarily composed of 18–24 bp sequences (Fig. 1B). The numbers of 22 and 21 bp nucleotides were highest in these cells, which indicated that mature miRNAs were the most abundant class of small RNAs in PBMCs (Fig. 1B). To annotate the RNA species in PBMCs, the experimental sequences were mapped to databases containing human miRNAs, rRNAs, tRNAs, snRNAs and snoRNAs. Fig. 1C shows that miRNAs, tRNAs, rRNAs, snRNAs, snoRNAs, piRNAs and Y RNAs were detected in PBMCs in this study. miRNAs were the most common RNA species in both healthy participants (70.33% of small RNAs) and patients with IgAN (69.83% of small RNAs); these levels were not significantly different between groups.

Figure 1.

Small RNA sequencing analysis in healthy participants and patients with IgAN. (A) Workflow of the small RNA sequencing experiment. (B) Small RNA read length distributions of PBMCs from healthy participants and patients with IgAN. (C) Percentages of each non-coding RNA category in PBMCs from healthy participants and patients with IgAN by comparing reads with miRNA/rRNA/tRNA/snRNA/snoRNA database. IgAN, immunoglobin A nephropathy; PBMC, peripheral blood mononuclear cells; miRNA, microRNA; rRNA, ribosomal RNA; tRNA, transfer RNA; snRNA, small nuclear RNA; snoRNA, small nuclear non-coding RNA; piRNA, piwi-interacting RNA.

Significantly differentially expressed miRNAs in PBMCs of patients with IgAN

To identify miRNAs associated with the development of IgAN, miRNAs that were significantly differentially expressed were assessed (P<0.05, fold change ≥2 or ≤0.5). As shown in Fig. 2A and B, there were 44 differentially expressed miRNAs identified: A total of 34 miRNAs were upregulated and 10 miRNAs were downregulated in PBMCs from patients with IgAN, compared with PBMCs from healthy participants. The chromosomal distribution of these differentially expressed miRNAs showed that most originated from chromosomes 1, 7, 8 and 11. In particular, chromosome 1 had 8 miRNAs, which constituted the highest number (Fig. 2C).

Figure 2.

Differentially expressed miRNA profiles in healthy participants and patients with IgAN. (A) Heatmap of 44 differentially expressed miRNAs in HCs and patients with IgAN. High expression level is showed in red and lower levels in blue (fold change ≥2 or ≤0.5 and P<0.05). (B) Numbers of miRNAs upregulated and downregulated in patients with IgAN compared with HCs. (C) Distributions of differentially expressed miRNAs on each chromosome. miRNA/miR, microRNA; IgAN, immunoglobin A nephropathy; HC, healthy control; hsa, Homo sapiens; chr, chromosome.

Functional prediction for upregulated genes

To elucidate the potential roles of differentially expressed miRNAs in IgAN, the top 10 upregulated miRNAs were initially identified. Among them, miR-6843-3p showed the most significant upregulation in patients with IgAN; miR-671-3p and hsa-miR-5100 both showed high abundance in healthy participants and patients with IgAN (Table II). The 3′ UTRs of mRNAs are targeted by miRNAs, and the mRNA stabilities are reduced or translation efficiency is inhibited. In this study, from the top 10 upregulated miRNAs, >800 potential target genes were identified; subsequently, an interaction network of miRNAs and mRNAs was constructed by mirWalk program (Fig. 3). miR-1287-5p had the largest subnetwork, which included 397 potential target genes. KEGG pathway analyses of enriched target genes were conducted, which revealed that ‘pathways in cancer’, the ‘PI3K-Akt signaling pathway’, the ‘MAPK signaling pathway’, the ‘cGMP-PKG signaling pathway’, the ‘oxytocin signaling pathway’ and the ‘adrenergic signaling in cardiomyocytes’ were likely targets of the identified miRNAs (Fig. 4A). GO term analyses showed that most genes were enriched in ‘binding’ and ‘protein binding’ in the molecular function component. In the cellular component, genes were mainly enriched in the ‘organelle’, ‘cytoplasm’ and ‘membrane-bounded organelle’ (Fig. 4B). In the biological component, genes were mainly enriched in ‘cellular process’, ‘biological regulation’, ‘regulation of biological process’ and ‘single-organism process’ (Fig. 4B).

Table II.

Top 10 upregulated miRNAs in small RNA sequencing analysis.

Order	miRNA	IgAN RPM	Healthy control RPM	log2 (fold change)	P-value
1	hsa-miR-6843-3p	0.054	0.399	2.8854	0.000744547
2	hsa-miR-3675-3p	0.079	0.349	2.1433	0.00775704
3	hsa-miR-3120-5p	0.673	2.901	2.1079	0.000941252
4	hsa-miR-873-5p	0.786	3.282	2.062	0.005967363
5	hsa-miR-5100	28.728	113.291	1.9795	0.000493526
6	hsa-miR-4661-5p	1.323	4.252	1.6843	0.006082614
7	hsa-miR-204-5p	0.717	2.122	1.5654	0.005360923
8	hsa-miR-184	0.57	1.687	1.5654	0.007100943
9	hsa-miR-671-3p	51.661	140.942	1.448	0.006420685
10	hsa-miR-1287-5p	7.066	15.25	1.1098	0.009113073

RPM, reads per kilobase of transcript, per million mapped reads; IgAN, immunoglobin A nephropathy; miR/miRNA, microRNA; hsa, Homo sapiens.

Figure 3.

miRNA-mRNA regulatory network of top 10 upregulated miRNAs (blue) and mRNAs (green) containing predicted miRNA targets. TargetScan and miRDB were used for predicting miRNA targets. miRNA, microRNA.

Figure 4.

Characteristics of potential target genes for the top 10 upregulated miRNAs in patients with IgAN. (A) Top 30 most enriched pathways of potential target genes for the upregulated miRNAs were analyzed using the Kyoto Encyclopedia of Genes and Genomes database. The size of the midpoint in the graph indicates the number of differential genes enriched, and the color indicates the P-value. (B) Gene Ontology enrichment analysis was used to annotate the function of potential target genes for the upregulated miRNAs. IgAN, immunoglobin A nephropathy; miRNA, microRNA.

Functional prediction for downregulated genes

The top 10 downregulated miRNAs are listed in Table III; miR-6782-5p showed the most significant downregulation in patients with IgAN, and miR-122-5p showed the highest abundance. The target genes of the 10 downregulated miRNAs were predicted, and an interaction network of miRNAs and mRNAs was constructed (Fig. 5). This demonstrated that miR-7705 had the largest subnetwork, which included 242 potential target genes. To identify the potential roles of the downregulated miRNAs, GO and KEGG enrichment analyses were performed. The ‘PI3K-Akt signaling pathway’, ‘microRNAs in cancer’, ‘pathways in cancer’, ‘Ras signaling pathway’ and ‘MAPK signaling pathway’ were identified as likely targets of the miRNAs (Fig. 6A). GO analysis showed that most genes in the molecular function group were involved in ‘organic cyclic compound binding’ and ‘protein binding’. Within the cellular component, ‘organelle’, ‘cytoplasm’, and ‘nucleus’ had the largest groups of enriched genes (Fig. 6B). In the biological component, ‘cellular process’, ‘single-organism process’, ‘metabolic process’, ‘regulation of biological process’ and ‘primary metabolic process’ were the top 5 enriched groups of target genes (Fig. 6B). Because aberrant O-glycosylation, which is catalyzed by core 1 synthase, glycoprotein-N-acetylgalactosamine 3-β-galactosyltransferase 1 (C1GALT1), is critical for IgA accumulation, potential miRNAs that can regulate it were predicted using TargetScan and miRDB. As shown in Fig. 7, the C1GALT1 3′ UTR possesses binding regions for miR-3121-3p (4 regions), miR-203a-3p (4 regions), and miR-200a-3p (1 region).

Table III.

Top 10 downregulated miRNAs in small RNA sequencing analysis.

Order	miRNA	IgAN RPM	Healthy control RPM	log2 (fold change)	P-value
1	hsa-miR-6782-5p	0.149	0	−3.8972	0.000324358
2	hsa-miR-3121-3p	0.152	0.029	−2.3899	0.005884351
3	hsa-miR-203a-3p	1.039	0.306	−1.7636	0.015739813
4	hsa-miR-548j-3p	0.54	0.182	−1.569	0.024266267
5	hsa-miR-122-5p	19.627	6.864	−1.5157	0.009925759
6	hsa-miR-490-3p	0.324	0.117	−1.4695	0.013063963
7	hsa-miR-6754-3p	0.295	0.115	−1.3591	0.020103511
8	hsa-miR-200a-3p	2.488	1.091	−1.1893	0.026159717
9	hsa-miR-4742-5p	1.616	0.734	−1.1386	0.00109596
10	hsa-miR-7705	0.32	0.154	−1.0551	0.042978753

RPM, reads per kilobase of transcript, per million mapped reads; IgAN, immunoglobin A nephropathy; miR/miRNA, microRNA; hsa, Homo sapiens.

Figure 5.

miRNA-mRNA regulatory network of top 10 downregulated miRNAs (blue) and mRNAs (green) containing predicted miRNA targets. TargetScan and miRDB were used for predicting miRNA targets. miRNA, microRNA.

Figure 6.

Characteristics of potential target genes for the top 10 downregulated miRNAs in patients with IgAN. (A) Top 30 most enriched pathways of potential target genes for the downregulated miRNAs were analyzed using the Kyoto Encyclopedia of Genes and Genomes database. The size of the midpoint in the graph indicates the number of differential genes enriched, and the color indicates the P-value. (B) Gene Ontology enrichment analysis was used to annotate the function of potential target genes for the downregulated miRNAs. IgAN, immunoglobin A nephropathy; miRNA, microRNA.

Figure 7.

miRNA sequences and potential target sites on C1GALT1 3′ UTR. TargetScan and miRDB were used for predicting miRNA targets. Sequences in red and italic indicate the predicted complementary binding sequences. miRNA/miR, microRNA; C1GALT1, core 1 synthase, glycoprotein-N-acetylgalactosamine 3-β-galactosyltransferase 1.

Discussion

The exact pathogenesis of IgAN is not well-defined, although it is characterized by the aberrant glycosylation and deposition of IgA in mesangial areas. Previous studies found that abnormal immune regulation is closely associated with the development of IgAN (3,9,12,13). Proliferating cell nuclear antigen expression increased in PBMCs, and the percentage of circulating B lymphocytes (CD19) increased in patients with IgAN (12,23), indicating abnormal proliferation of PBMC in the development of IgAN. PMBCs from patients with IgAN produced more IgA than PBMCs from healthy controls (13). Interleukin 4, mainly produced by T cells, can induce IgA synthesis in PBMCs (13). Degradation of C1GALT1 by miR-148b mimics can attenuate IgA glycosylation in PBMCs (9). The results of these studies indicate that abnormally regulated PBMCs proliferation along with IgAN production and glycosylation have important roles in IgA accumulation and IgAN development.

miRNAs play important roles in cell behaviors such as cell proliferation, migration and apoptosis by regulating gene expression (5). miRNAs are ubiquitous and abundant in body fluids and cells; importantly, they are also stable (24). Therefore, miRNAs can be good biomarkers and drug targets (5,25). Previously, several studies have demonstrated that miRNAs could be involved in IgAN development, including miR-148b, let7b, miR-146a, miR-155, miR-200a, miR-200b and miR-429 (7–10). However, the role of miRNA expression in the pathogenesis of IgAN has not been well explored. To identify additional functional miRNAs, small RNA sequencing was performed; this analysis showed that miRNAs constituted nearly 70% of small RNAs in PBMCs from both healthy participants and patients with IgAN. Furthermore, 44 differentially expressed miRNAs were identified, of which 34 were upregulated and 10 were downregulated. Among them, miR-148a-3p, miR-184, and miR-200a were previously reported to be associated with IgAN (26,27), the other 41 were newly found to be associated with IgAN. In addition, miR-135a is involved in podocyte injury, and miR-203 is known to promote diabetic nephropathy (28,29). To the best of the authors' knowledge, this is the first report that the remaining miRNAs are potentially associated with kidney disease.

Among the top 10 upregulated miRNAs, miR-6843-3p exhibited the greatest fold change, followed by miR-3675-3p and miR-3120-5p. miR-671-3p showed the most abundant expression in both healthy participants and patients with IgAN, followed by miR-5100 and miR-1287-5p. KEGG pathway analysis showed enriched target genes in the PI3K-Akt signaling pathway and MAPK pathways, which were reportedly involved in IgAN (30,31), indicating a potential role in IgAN for the top 10 upregulated miRNAs in PBMCs in this study.

Among the top downregulated miRNAs, miR-6872-5p exhibited the greatest fold change, followed by miR-3121-3p and miR-203a-3p. miR-122-5p showed the most abundant expression in both healthy participants and patients with IgAN, followed by miR-200a-3p and miR-4742-5p. Previously, miRNAs were shown to promote cell proliferation and the production of aberrant glycosylated IgA1 in IgAN (9,12). C1GALT1, which is a critical enzyme involved in IgA glycosylation, may be regulated by miRNAs, such as miR-148b (9). Among the members of the miR-200 family, the intrarenal expression of miR-200c was downregulated in renal tissue; however, miR-200a and miR-200b were upregulated (7). In the present study, miR-200a-3p was identified in PBMCs, and 4 potential binding sites were predicted. Moreover, miR-3121-3p and miR-203a-3p, the top 2 downregulated miRNAs, were newly identified and predicted to regulate C1GALT1 expression. miR-200a-3p and miR-203a-3p were previously demonstrated to inhibit cell proliferation (28,32,33). KEGG pathway analysis of downregulated miRNA target genes suggested that the PI3K-Akt signaling pathway, RAS signaling pathway and MAPK signaling pathways may also be involved in the role of miRNAs in IgAN by regulating cell proliferation or gene expression involved in IgAN development. These findings suggested that differentially expressed miRNAs may be involved in the progress of IgAN by controlling PBMC biological behaviors and IgA deposition through direct targeting of genes or indirect regulation of various signaling pathways.

In conclusion, the present study identified 44 differentially expressed miRNAs associated with IgAN, of which numerous associations were novel. In addition, these miRNAs may provide biomarkers and drug targets for IgAN diagnosis and therapy.