Integrated profiling of microRNA expression in membranous nephropathy using high-throughput sequencing technology

The present study analyzed microRNA (miRNA) expression profiles in peripheral blood lymphocyte cells (PBLCs) from patients with membranous nephropathy (MN) and normal controls (NC), in an effort to improve the understanding of the pathogenesis of MN. High-throughput sequencing was performed on 30 MN patients and 30 healthy individuals (NC group). Known and novel miRNAs were analyzed and the results were confirmed by quantitative reverse transcription PCR (qRT-PCR). In total, 326 miRNAs showed a significant difference in expression between the MN and NC groups. This included 286 downregulated miRNAs and 40 upregulated miRNAs. In addition, there were 6 novel miRNAs that presented differential levels of expression between the MN and NC groups. The miRNAs were mapped to the genome, using a short oligonucleotide alignment program (SOAP), to analyze their expression and distribution. Twenty-five percent of the unique miRNAs in the MN group and 52.1% in the NC group were mapped to the genome. One hundred and eight mismatches were identified. Seventy-seven mismatches were detected in a higher proportion of the MN samples, compared with the NC samples. Twenty-five mismatches were detected in a higher proportion of the NC samples than the MN samples. Differential miRNA expression was also detected between 10 randomly selected pair groups, as depicted in a cluster analysis diagram. These data indicate that differential miRNA expression may be involved in the pathogenesis of MN. In addition, the discrepancies between the MN and NC groups, in the mismatched miRNAs that were mapped to the genome, strongly suggest that miRNAs play an important role in the pathogenesis of human disorders. miRNAs may provide a potential breakthrough in the research of MN and may provide a novel biomarker for the diagnosis and treatment of the disease.


Introduction
Membranous nephropathy (MN) is a renal disease that is histologically characterized by the uniform thickening of the glomerular capillary wall. This is caused by subepithelial deposits of immune complexes, which appear as granular deposits of immunoglobulin G when visualized using immunofluorescence, and as an electron-dense deposit when observed under an electron microscope. MN is a frequent cause of nephrotic syndrome in adults and can eventually progress to end-stage renal failure in many patients (1)(2)(3). The MN disease site is the glomerular visceral epithelial cell, or podocyte, which is a highly specialized and terminally differentiated cell that is found on the outside of the glomerular basement membrane (4). The disease, similar to most immune glomerular diseases, manifests itself as an immune response to the self-antigens expressed on the podocyte cell membrane (5). There is no specific treatment for MN. Diuretics and angiotensin-converting enzyme inhibitors have limited effects, whilst immunosuppressants cause side-effects and their use in treatment is controversial (6).
Currently, the diagnosis of kidney disease depends mainly on renal biopsy and immunofluorescence (7). There is an urgent need for a validated biomarker that is easy to measure and able to accurately predict long-term outcomes, to aid both in the diagnosis and treatment of the disease (8). The clinical management of patients with MN may improve with the better understanding of the underlying mechanisms of the disease and the availability of biomarkers. microRNAs (miRNAs) regulate gene expression and have been found to modulate crucial biological processes, including differentiation, proliferation and apoptosis (9). They function through various mechanisms, such as targeted miRNA degradation and translational inhibition (10,11). Hundreds of miRNAs have been identified, using molecular cloning and bioinformatics prediction strategies, in worms, flies, fish, frogs, mammals and flowering plants (12). Several studies have indicated a possible link between miRNAs and kidney disease (13)(14)(15). Although the expression results of the various array analyses are inconsistent, Integrated profiling of microRNA expression in membranous nephropathy using high-throughput sequencing technology the data indicate that the dysregulation of miRNAs may play a pivotal role in the pathogenesis of kidney disease. However, there have been few studies on the association between miRNAs and MN. In this study, we analyzed the pathogenesis of MN at the miRNA level. Next generation high-throughput sequencing is an effective technique that enables miRNA profiling at unprecedented quantitative and qualitative levels (16). The data presented in this study may enhance our knowledge of miRNA expression profiles in patients with MN, and may provided insight into the pathogenesis, diagnosis and treatment of MN.

Patients and methods
Patient samples. Peripheral blood samples were collected from 30 patients at the Second Clinical Medical College, Jinan University, Shenzhen People's Hospital (Shenzhen, China) between 2011 and 2012. The diagnosis of MN was confirmed by renal biopsy and immunofluorescence. Thirty specimens for the control group were obtained from individuals who underwent annual body check-ups and were confirmed as healthy in 2012, at Shenzhen People's Hospital.
Prior written informed consent was obtained from all patients. The project was approved by the Shenzhen People's Hospital Ethics Committee. This study was performed in accordance with the guidelines of Jinan University, which abides by the Helsinki Declaration on ethical principles for medical research involving human subjects.
Sample processing. Samples were collected in EDTA tubes and then transferred to depletion filters to separate the lymphocytes, in accordance with the LeukoLock Total RNA isolation system protocol (Ambion, Austin, TX, USA). Total RNA was extracted from the cells using the miRNeasy Mini kit (Qiagen, Hilden, Germany), in accordance with the manufacturer's instructions. The integrity of the RNA and the presence of miRNAs were assessed by micro-capillary electrophoresis, using an RNA 6000 kit and small RNA kit (Agilent Technologies Inc., Santa Clara, CA, USA), respectively. The concentration and quality of the RNA were assessed by absorbance spectrometry on a NanoDrop 2000 spectrophotometer (Thermo, Waltham, MA, USA). miRNA microarray. miRNA microarrays, composed of 455 human, 236 rat and 344 mouse miRNAs, were used to detect all human, rat, mouse and other miRNAs in the Sanger miRNA database (release 8.1). This study used miRNA power labeling (Exiqon, Vedbaek, Denmark) to 3'-or 5'-end label 0.5 µl of a sample or a human universal reference total RNA, with the cy3-like HY3 or cy5-like HY5 dye, respectively. Locked nucleic acids (LNAs) were used, as their superior sensitivity over conventional DNA-based miRNA arrays can discriminate between closely related miRNA family members. miRNA microarray analysis. To create an MN library and normal control (NC) library, total miRNAs from each subject underwent miRNA library construction and sequencing. The process is illustrated in Fig. 1. The 50 nucleotide sequence tags used for the high-throughput sequencing were obtained through data cleansing to remove low quality tags and several types of contaminants. The length distribution of the clean tags was summarized and used in a standard bioinformatics analysis. The clean tags were annotated and placed into different categories to predict novel miRNAs and edit known miRNAs. The tags that could not be annotated were placed in a separate category.
Statistical analysis. Statistical analyses were performed after the libraries were established. The miRNAs were mapped to the genome and a summary was produced for the known miRNA alignments. Cluster analysis, differential expression and base edit analysis were performed for the miRNAs.
Quantitative reverse transcription PCR (qRT-PCR) verification of miRNA results. qRT-PCR was performed to verify the results of the deep sequencing analysis for 5 differentially expressed miRNAs. Total RNA (2 µg) was reverse transcribed into cDNA using a reverse transcription kit, in accordance with the manufacturer's instructions (Promega, Madison, WI, USA). The cycle parameters for the PCR reaction were as follows: 95˚C for 5 min followed by 40 cycles of a denaturing step at 95˚C for 10 sec and an extension step at 60˚C for 60 sec. All reactions were run in triplicate. The relative quantification 2 -∆∆Ct method was used to determine the changes in expression levels between the MN and NC groups. The miRNA expression levels were normalized to the reference RNA, RUN6. The -∆∆Ct values were calculated using the following formula: -∆∆Ct = -Ct (MN -NC) , where Ct is the cycle threshold provided by the Rotor-Gene 6000 Series Software 1.7 (Qiagen).

Results
Small RNA expression and distribution in each genome. Highthroughput sequencing produced 11,436,003 and 18,696,751 Figure 1. MicroRNA (miRNA) microarray analysis process. The data analysis process used 50 nucleotide tags that were produced by high-throughput sequencing. The low quality reads were filtered out in accordance with the base quality value. The adaptor sequence was trimmed at the 3'-terminus. The 5' adaptor contaminants that were formed by ligation were removed. The clean, high-quality reads were used for standard analysis. Cluster analysis contained the data on small RNA expression and distribution in each genome. Differential expression included miRNA and novel miRNA expression in the membranous nephropathy (MN) and normal control (NC) groups. Target prediction was analyzed using the summary of known miRNA alignments and miRNA base edits. rRNA, ribosomal RNA; tRNA, transfer RNA; snRNA, small nuclear RNA; snoRNA, small nucleolar RNA; sRNA, small RNA. high quality sequence reads from the MN and NC groups, respectively. Following the removal of the contaminated reads, the total number of clean, small RNA reads was 11,109,127 and 16,001,191, in the MN and NC groups, respectively. The number of unique small RNAs was 154,580 in the MN group and 1,034,806 in the NC group. The small RNA tags were mapped to the genome using the short oligonucleotide alignment program (SOAP) program, as previously described (17), to analyze their expression and distribution. The results are presented in Table Ⅰ. The percentage of unique small RNAs mapped to the genome in the NC group was 52.10%, which was greater than the 25% observed for the MN group. By contrast, the MN percentage for the total small RNAs was 85.00%, whilst for the NC group the value was 83.10%. Fig. 2 shows the number of small RNA tags that were located on   Summary of known miRNA alignments. We aligned the small RNAs to miRNA precursors in the corresponding species to obtain the miRNA count. This also allowed the identification of bases at each position within the miRNAs. The results are shown in Table ⅠⅠ. The percentage occurence of each base (A/U/C/G; A, adenine, U, uracil, C, cytosine, D, guanine) at each position in the small RNAs was calculated. A large variation was observed between the MN and NC groups. For example, 83.99% of the bases in the second position were A in the MN group. Only 22.99% of the bases at this position were A in the NC group. Significant differences were also observed at the 9th base (10.10% U in the MN group, 63.04% U in the NC group) and 14th base (83.09% C in the MN group, 16.62% C in the NC group). However, some similarities were observed between the two groups at base positions 1, 6, 10, 12, 13, 15, 17, 18, 22, 23 and 24. At these positions, the most common base was the same in each group. The outcomes are depicted in Fig. 3.

Differential expression of miRNAs in the MN and NC groups.
Four procedures were completed to compare the miRNA expression between the MN and NC libraries: ⅰ) relative expression analysis was used to normalize the data against the number of miRNAs and the total number of small RNA reads. This was used to define the expression preferences of individual miRNAs, between the two libraries. ⅱ) The outcome of the relative expression analysis was multiplied by a constant, set at 1x10 6 : Normalized expression = (actual miRNA count/total count of clean reads) x10 6 . ⅲ) The fold change and P-values were calculated from the normalized expression data: fold change = log 2 (MN/NC), P(y/x) = (N 2 /N 1 ) y × (x+y)!/x!y!(1+N 2 /N 1 ) (x+y+1) (x and y indicate the number of reads of a miRNA in the NC and MN libraries, respectively. N 1 and N 2 represent the total number of clean reads in the NC and MN libraries, respectively). ⅳ) The fold change and Audic-Claverie method (18) were used to define the differential expression of miRNAs between the two groups. Fold changes (log 2 MN/NC) with P≤0.01 were considered to indicate a statistically significant result. Differential expression between the MN and NC libraries was found in 326 miRNAs. These consisted of 286 miRNAs that were downregulated and 40 miRNAs that were upregulated   Table Ⅲ.
Differential expression of novel miRNAs in the MN and NC groups. There were 15 novel miRNAs in the MN group and 22 novel miRNAs in the NC group. Only 6 of these showed a significantly different level of expression between the MN and NC groups. The novel-miR-82, novel-miR-98, novel-miR-89 and novel-miR-84 miRNAs were downregulated, whilst the novel-miR-152 and novel-miR-15 miRNAs were upregulated ( Table Ⅳ) (Fig. 5). The samples were clustered in accordance with their similarities in expression patterns, i.e., fold change and P-value. Red indicated that the miRNA had a higher expression level in the MN specimens, whilst green indicated that the miRNA had a higher expression level in the 10 specimens from the NC group. Gray indicated that the miRNA was not expressed in at least one sample.
miRNA base edits. Nucleotide base positions 2-8 of a mature miRNA are highly conserved and are known as the seed region.
The target of the miRNA may be dependent on this region. In our analysis, seed region base changes were detected by aligning unannotated small RNAs with mature miRNAs from the miRBase 18 database (http://www.mirbase.org/). Mismatches in the alignment were assumed to be associated with the mechanism of the disease. In this study, 108 miRNAs, which were common between the MN and NC groups were analyzed for base edits. The proportion of miRNAs with base edits was used to obtain the ratio of base edits between the two groups (Table Ⅴ). There were 77 miRNAs in which there were more base edits in the MN group than in the NC group (MN/NC >1), 6 miRNAs that were equivalent between the two groups (MN/ Validation of miRNA expression by qRT-PCR. The expression levels of 5 randomly selected miRNAs: hsa-miR-7-5p, hsa-miR-615-3p, hsa-miR-577, hsa-miR-98 and hsa-miR-375, were compared. The hsa-miR-98 and hsa-miR-375 miRNAs were upregulated, whilst hsa-miR-7-5p, hsa-miR-615-3p and hsa-miR-577 were downregulated. The qRT-PCR data were obtained using the 2 -∆∆Ct method and normalized using RNA RUN6 as a reference. The log 2 (MN/NC) value was compared with the fold change value (Table Ⅵ). The log 2 (MN/ NC) values of hsa-miR-7-5p, hsa-miR-615-3p and hsa-miR-577 were -5.06, -2.40 and -1.12, respectively. These data confirmed that these miRNAs were downregulated. The log 2 (MN/NC) values of hsa-miR-98 and hsa-miR-375 were 4.28 and 2.63, respectively, again confirming the previous data, indicating that they were upregulated.

Discussion
The high through-put sequencing used in this study is suitable for the analysis of small RNA molecules as it is able to decrease the loss of nucleotides in the reads, caused by secondary structure. The technology is also ideal as it does not require a large sample quantity (19). Such an analysis can obtain millions of small RNA sequence tags in one run and can identify the differential expression of small RNAs between two samples (20). We performed high through-put sequencing on a large number of peripheral blood lymphocytes from individuals separated into the MN and NC groups. The aim was to identify dysregulated miRNAs that may serve as reliable diagnostic markers and potential therapeutic targets. The data confirmed that dysregulated miRNAs may play an important role in the pathogenesis of nephropathy, which is consistent with previous studies (21)(22)(23). We determined the unique and total number of small RNAs in the MN and NC groups, and positioned the small RNAs within the genome. Expression analysis and distribution of the small RNAs was also performed. The number of total and  unique small RNAs was greater in the NC group than in the MN group. The miRBase 18 database (http://www.mirbase.org/) provides a range of data to facilitate studies of miRNA genomics. This has been used previously to map all miRNAs to their genomic coordinates (24), allowing a network of genomewide miRNA expression to be produced (25,26). Analysis of the properties of miRNA targets is a promising approach to the prediction of miRNA function. If the targets of specific miRNAs are enriched with genes associated with a particular biological process, it is reasonable to infer that the miRNA is also involved in the same process (27). The function of miRNAs was not predicted in this study; however, a statistical analysis was performed, revealing a discrepancy between a diseased and a normal group of specimens, regardless of the quantity and distribution of miRNAs. This provided strong evidence for the function of miRNAs in the pathogenesis of disease.
Our data demonstrated the bias in miRNA nucleotides at each base position. We found that U was the dominant nucleotide in miRNAs. This was particularly noticeable at positions 1, 6, 8, 18, 19, 20, 22 and 24 in the MN group, and at positions 1, 6, 9, 14, 22 and 24 in the NC group. These results are consistent with those of Zhang et al (28), who demonstrated that the distribution of nucleotides indicated an important role for U at the boundaries of the seed region and termini. In addition, there was a large discrepancy in the proportion of the 4 bases at each position. For example, U accounted for 96.61 and 84.27% of bases at the 1st position in the MN and NC groups, respectively. The A base accounted for 92.45 and 51.50% of bases at the 17th position in the MN and NC groups, respectively. Genes with a higher level of expression showed stronger signals, which indicated that these nucleotides were responsible for the regulation of translation initiation. The diversity of nucleotide sequences surrounding the initiation codon has been explained by differences in relative contributions from two distinct patterns (29), and preferred nucleotide sequences varied between different eukaryotic species (30). We speculate that the reasons for the different nucleotide bias, between the MN and NC groups, resulted from the individual differences in evolution or the role of miRNAs. This requires further research.
In the present study, the results of the miRNA differential expression analysis were unexpected as there were more downregulated (n=286) miRNAs than upregulated (n=40) miRNAs. By contrast, there are several reports on miRNA pathogenesis in nephropathy that have reported more upregulated miRNAs, compared with downregulated miRNAs (13,31). This is also the case for diseases, such as congenital disorders (32), cancer (33) and immunological diseases (34). However, certain studies have reported an increase in downregulated miRNAs, compared with upregulated miRNAs; however, the difference in numbers is small. Chen et al reported more downregulated (n=41) miRNAs than upregulated (n=33) ones in urothelial cell carcinoma (35). Osanto et al reported more downregulated (n=41) miRNAs than upregulated (n=29) ones in clear cell renal cell carcinoma (36). miRNAs that are more abundant in the kidneys, compared with other organs, include miR-192, miR-194, miR-204, miR-215 and miR-216 (37). The miRNA-30 family (hsa-miR-30e-5p, hsa-miR-30e-3p, hsa-miR-30d-5p, hsa-miR-30c-5p, hsa-miR-30c-2-3p, hsa-miR-3b-5p, hsa-miR-30b-3p, hsa-miR-3a-5p and hsa-miR-30a-3p) and the miR-133 family (hsa-miR-133b and hsa-miR-133a) have been linked to the connective tissue growth factor (CTGF), which is a key molecule in the process of fibrosis (38). It is also a key molecule in the process of nephropathy (39). Our study demonstrated that the miRNA-30 family was downregulated and that the miR-133 family was upregulated. The loss of miR-23b, miR-24 and miR-26a resulted in the rapid progression of marked glomerular and tubular injury. Their existence has been shown to be critical in maintaining glomerular filtration (40). These miRNAs were also downregulated in our study. We deduced that the key involvement of miRNAs in the pathogenesis of MN was an outcome of downregulation. This could explain why the downregulation of miRNAs was more common than the upregulation. However, our inference requires a more in-depth study.
We were also able to predict novel miRNAs, with 6 that showed a significant difference in expression between the MN and NC groups. Four of the 6 were downregulated and 2 were upregulated. The outcome was in agreement with the higher numbers of miRNAs that were downregulated, compared with the number that were upregulated, as described earlier in this study. Certain studies have reported that the read number for most novel miRNAs is much lower than that for the conserved miRNAs, which indicates that non-conserved miRNAs are usually expressed at a lower level (35,41). Despite the limited number of novel miRNAs in this study, we found that the fold change was relatively large, with values >6 or <-6. This indicated that the novel miRNAs may have functional relevance in the pathogenesis of MN. The identification of novel miRNA genes is important as it may reveal putative genes that exert a regulatory effect on different types of cancer (42). Dhahbi et al (43) found 20 novel miRNAs that were differentially expressed between young and senescent fibroblasts. Three novel miRNAs have been shown to exhibit relative sequence counts of >10 and are likely to be involved in the development of prostate cancer (41). These studies demonstrated that novel miRNAs have are closely associated with the occurrence of disease. The targets and functions of novel miRNAs, which had not previously been investigated in MN, have yet to be determined.
The nucleotides at positions 2-8 of mature miRNAs are known to be highly conserved. The targets of miRNAs may be altered by a change in the nucleotides in this region (44). We calculated the percentage of base edits and compared the percentages between the MN and NC groups. A discrepancy was observed, which indicated that this region may be implicated in the pathogenesis of the disease. Blow et al found that 6 of 99 surveyed pre-miRNAs were edited in at least 1 of 10 human tissues (45). The hsa-miR-1269b and hsa-miR-1034-3p miRNAs had the largest percentage (100%) of edited bases in the NC group. There were 9 miRNAs in the MN group where the percentage of edited bases reached 100%. There was a similar base edit percentage in both the MN and NC groups. For example, hsa-let-7a-5p, 85.52% in the MN group and 64.99% in the NC group; hsa-miR-1304-4p, 100% in the MN group and 100% in the NC group; hsa-miR-130-3p, 5.54% in the MN group and 2.10% in the NC group. We identified that the number of miRNAs with a ratio (MN)/ratio(NC) >1 was 77, whilst those with a ratio (MN)/ratio(NC) <1 was 25. Kawahara et al provided the first evidence that edited miRNAs have a biological significance in vivo (46). We speculated that miRNA base edits may be involved in the pathogenesis of MN and may provide an innovative method for investigating the mechanisms responsible for the development of MN.
This study is one of the few that have profiled the expression patterns of miRNAs on a genome-wide scale in patients with MN. Our results demonstrated the differential expression of miRNAs and a discrepancy in nucleotide bias between normal and diseased individuals. Base edits showed a clear difference between the two groups, which strongly indicated that dysregulated miRNAs may present an area of research into the pathogenesis of MN. Our study also indicated that miRNAs can serve as a potential clinical target to diagnose and treat MN patients in the future. Whilst the results indicated the significant potential and benefit of miRNAs, further studies are required to provide further insight into the molecular functions of miRNAs. Our data may be used as basic research to support novel methods for the investigation, diagnosis and treatment of MN. We anticipate that miRNA-based genetic therapies will be developed to replace traditional therapies for the future benefit of patients.