Integrative microRNA-gene expression network analysis in genetic hypercalciuric stone-forming rat kidney

Background. MicroRNAs (miRNAs) influence a variety of biological functions by regulating gene expression post-transcriptionally. Aberrant miRNA expression has been associated with many human diseases. Urolithiasis is a common disease, and idiopathic hypercalciuria (IH) is an important risk factor for calcium urolithiasis. However, miRNA expression patterns and their biological functions in urolithiasis remain unknown. Methods and Results. A multi-step approach combining microarray miRNA and mRNA expression profile and bioinformatics analysis was adopted to analyze dysregulated miRNAs and genes in genetic hypercalciuric stone-forming (GHS) rat kidneys, using normal Sprague-Dawley (SD) rats as controls. We identified 2418 mRNAs and 19 miRNAs as significantly differentially expressed, over 700 gene ontology (GO) terms and 83 KEGG pathways that were significantly enriched in GHS rats. In addition, we constructed an miRNA-gene network that suggested that rno-miR-674-5p, rno-miR-672-5p, rno-miR-138-5p and rno-miR-21-3p may play important roles in the regulatory network. Furthermore, signal-net analysis suggested that NF-kappa B likely plays a crucial role in hypercalciuria urolithiasis. Conclusions. This study presents a global view of mRNA and miRNA expression in GHS rat kidneys, and suggests that miRNAs may be important in the regulation of hypercalciuria. The data provide valuable insights for future research, which should aim at validating the role of the genes featured here in the pathophysiology of hypercalciuria.


INTRODUCTION
Kidney stones are commonly found in children and adults (Coe, Evan & Worcester, 2005), and can becaused by multiple factors. Idiopathic hypercalciuria (IH) is an important risk factor for calcium urolithiasis (Worcester et al., 2013;Yoon et al., 2013). Patients with IH have normal serum Ca 2+ , and increased urinary calcium excretion. But the pathophysiological process and molecular mechanism of IH are still unclear. The genetic hypercalciuric stone-forming (GHS) rat, has many pathophysiological characteristics identical to that of IH patients, such as normal serum Ca 2+ , hypercalciuria, elevated intestinal Ca 2+ resorption and a tendency to lose calcium from the bone (Frick et al., 2013;Frick, Krieger & Bushinsky, 2015), which is regarded as an ideal animal model of calcium urolithiasis.
MicroRNAs (miRNAs) are a group of small, non-coding RNAs that regulate proteincoding gene function at the post-transcriptional level by binding to complementary sites on target mRNAs in the 3 UTR (Ambros, 2004). Meanwhile, miRNAs have been shown to regulate a wide range of biological processes including cell growth, metabolism, differentiation, proliferation and apoptosis, which have important implications in diseases processes (Ambros, 2001). Dysregulation of miRNAs has been associated with many human diseases. However, miRNA expression patterns and their biological functions in urolithiasis remain unknown. Understanding the relevance of miRNA and mRNA expression patterns in GHS rat kidneys is important to better elucidate the relationship between pathophysiological process and genes.
In this study we therefore conducted mRNA and miRNA expression profiling of three pairs of GHS and normal Sprague-Dawley (SD) rat kidneys. A subset of differentially expressed genes was validated by qPCR in 12 pairs of kidneys. Bioinformatic analysis was further performed to construct an integrative regulatory network of altered miRNA-mRNA transcriptsin GHS rat kidney.

Animals
The colony of GHS rats were created by selective breeding of male and female Sprague-Dawley (SD) rats with the highest 24-hour urine calcium excretion as previously described (He et al., 2015). By the 28th generation, GHS rats stably excreted significantly higher levels of urine calcium than wild-type normal SD rats. Normal SD rats were purchased from the Experimental Animal Center, Tongji Medical College, Huazhong University of Science and Technology, China. Twelve GHS rats with a body weight of 250-280 g were used for this study. A total of 12 normal SD rats were matched with GHS rats with respect to body weight and age, and served as control rats. All rats were fed 13 g/day of a normal calcium diet (1.2% calcium, 0.65% phosphorus, 0.43% chloride, 0.4% sodium, and 0.24% magnesium per gram of food). All animal procedures were approved by the Ethical Committee of Tongji Hospital, Tongji Medical College, Huazhong University of Science and Technology (No. TJ-A20141211). All surgeries were performed under sodium pentobarbital anesthesia.

Urine and serum calcium and phosphorus determination
Two successive 24-hour urine samples were collected before rats were killed, then blood samples were taken after rats were killed. Urine calcium, serum calcium, and phosphorus were measured using an Abbott Aeroset AutoAnalyzer (Abbott Diagnostics, Chicago, IL, USA).

RNA extraction
Total RNA was extracted from kidney tissue using RNeasy Fibrous Tissue kit (Qiagen, Dusseldorf, Germany) according to the manufacturer's protocol. RNA purity and concentration were assessed by NanoDrop ND-2000 spectrophotometer (Thermo Fisher Scientific Inc., Waltham, MA, USA) and electrophoresis of RNA on agarose gel containing formaldehyde was used to evaluate the integrity of RNA.

Microarray
Three GHS rats and three SD rats were randomly selected for microarray analysis. The Affymetrix GeneChip miRNA 4.0 Array and Affymetrix Gene 1.0 Array for rats were used to compare miRNA and mRNA expression profiles, respectively, in GHS and control rat kidneys. Microarray analysis was performed by GMINIX Informatics Ltd. Co, Shanghai, China. The data have been deposited in the NCBI Gene Expression Omnibus and are accessible through GEO Series accession number GSE75543.

Strategy
As shown in Fig. 1, we used a multi-step strategy to identify genes dysregulated in GHS rats relative to the control group. First, significantly differentially expressed mRNA and miRNA were identified using a random variance model (RVM) corrective analysis of variance (ANOVA) (Wright & Simon, 2003). Second, the miRanda database was used to predict putative miRNA targets, and the overlap between these target genes and differentially expressed mRNAs was established. Third, these intersecting genes were classified according to their biological functions using the Gene Ontology System. Similarly, pathway analysis was used to identify affected KEGG pathways. Fourth, genes that were present in both the enriched GO terms and significant KEGG pathways were used to construct a miRNA-gene network (Joung et al., 2007;Shalgi et al., 2007) and signal-net (Spirin & Mirny, 2003;Zhang & Wiemann, 2009). The center of the network is represented by a degree, which indicates the predicted interaction of a given miRNA with its target genes. All of the data mentioned above were analyzed by GCBI working platform (GMINIX Informatics Ltd. Co, Shanghai, P. R. China), and the principles and methodologies of data analysis were described in File S1. Finally, a subset of the predicted miRNA-mRNA pairs was selected for validation by qRT-PCR in an extended cohort of kidney tissues.

Real-time RT-PCR
For mRNA, total RNA was extracted and 500 ng of RNA was used for cDNA synthesis using the Takara reverse transcription kit (Takara, Dalian, China) according to the manufacturer's instructions. PCR was conducted using SYBR Premix Ex Taq (Takara, Dalian, China) according to the manufacturer's instructions on an Mx3000P system (Agilent Stratagene, Santa Clara, CA, USA). The primers were chemically synthesized by Tsingke, Wuhan, China and are listed in Table S1. The All-in-One TM miRNA qRT-PCR Detection Kit (GeneCopoeia, Guangzhou, China) was used for both cDNA synthesis and quantitative detection using miRNA specific primers (GeneCopoeia, Guangzhou, China). GAPDH and U6 were used as internal controls to determine the relative expression of target mRNA and miRNA. All reactions were performed in triplicate.

Statistical analysis
Continuous variables were expressed as means ± standard deviation. For each triplicate of microarray data, the geometric mean was used. Student's t -test was applied for comparisons between two groups, and ANOVA for comparisons between more than two groups. A P-value of <0.05 was considered statistically significant.

Serum calcium and phosphorus levels, and urine calcium excretion
Serum calcium and phosphorus levels were not significantly altered in GHS compared with SD control rats (Figs. 2A and 2B). GHS rats did, however, excrete significantly more urine calcium (mg/24 h) than SD control rats (Fig. 2C).

miRNA and mRNA expression profiles in GHS ratkidney
The miRNA expression profiles of three GHS rats' kidney tissues and the control rats were determined using miRNA microarray analysis. A total of 19 miRNAs were significantly differentially expressed between the two groups (p < 0.05), of which 10 and nine miRNAs were up-or downregulated, respectively, in GHS vs. control rats (Table 1 and Table S2); of these rno-miR-184, rno-miR-21-3p and rno-miR-672-5p had the largest positive fold changes, while rno-miR-484, rno-miR-138-1-3p and rno-miR-201-3p had the largest negative fold changes. The expression levels of these 19 miRNAs are illustrated by the heatmap in Fig. S1. Microarray-based mRNA expression analysis was also conducted and 2,418 genes were identified as significantly differentially expressed in GHS rats (P < 0.05) including 1,057 upregulated genes and 1,361 downregulated genes (Table S3).

miRNA target gene prediction
Target mRNAs for differentially expressed miRNAs were predicted using the miRanda database. Since miRNAs negatively regulate gene expression, upregulated miRNAs result in downregulated target mRNAs, and vice versa. A total of 29,164 miRNA-mRNA pairs (based on 19 dysregulated miRNAs and their 10,521 target mRNAs) were predicted (Table S4).

Integrated analysis of dysregulated miRNAs and mRNAs
The set of intersecting mRNAs between the predicted target mRNAs and differentially expressed mRNAs were selected (Table S5) and those that were negatively correlated with their predicted miRNA matches were used for downstream GOanalysis and KEGG pathway analysis.
There were 417 upregulated and 286 downregulated GO terms (P < 0.05), with the most significant GO terms including negative regulation of apoptotic process (GO:0043066), response to lipopolysaccharide (GO:0032496) and inflammatory response (GO:0006954). Table 2 shows the top 15 up-and downregulated GO terms, with further detail in Tables  S6 and S7. According to KEGG pathway analysis, 93 KEGG pathways were significantly enriched, of which 75 were upregulated and 18 downregulated at P < 0.05 (Tables S8 and S9). The most highly enriched pathways included the Pantothenate and CoA biosynthesis and synthesis and degradation of ketone bodies pathways. The top up-and downregulated pathways are shown in Fig. 3. The darker bars indicate pathways that reported to be related to urolithasis in previous studies.
Next, to illustrate the role of key miRNAs in the regulation of kidney mRNAs in GHS rats, a miRNA-mRNA-Network was built based on the subset of significantly differentially expressed mRNAs that were also members of significantly enriched GO terms and KEGG pathways (Table S10, Fig. 4). In total, 223 mRNAs and 19 miRNAs were included in the network, where box nodes represent miRNAs, and circular nodes represent mRNAs. The degree of connectivity, which represents the number of genes regulated by a given miRNA, is indicated by the size of the node with a higher degree of connectivity represented by larger nodes (Table S11). Rno-miR-674-5p, rno-miR-672-5p, rno-miR-138-5p, and rno-miR-21-3p had high degrees of connectivity and may play crucial roles in this regulatory network. Meanwhile, Sema5a, Lpin2, Gcnt1, Masp1, Olr1 and Traf3 were the most common miRNA targets. Table 3 presented a subset of significantly regulated miRNA-mRNA hybrids. Moreover, to investigate the regulatory relationships between these genes and their potential role in hypercalciuria, we performed a signal-net analysis based on significantly regulated KEGG pathways (Fig. 5, Table S12). Signal-net analysis has shown that NF-kappa B signal pathway, including the members of NF-kappa B1, RelA and NF-kappa B2 et al., might play a core role in the gene regulatory network. And NF-kappa B1 and RelA had the highest degrees of connectivity of 19 and 18, respectively. Box nodes represent miRNAs, circular nodes represent mRNAs. Blue represents downregulation, while red represents upregulation. The higher the degree of connectivity of a gene, the larger the node within the network. In total, 223 mRNAs and 19 miRNAs were included in the network. Rno-miR-674-5p, rno-miR-672-5p, rno-miR-138-5p, and rno-miR-21-3p were found to have the highest degrees of connectivity. Figure 5 Signal-net. Blue, downregulation; red, upregulation; a, activation; a(ind), activation (indirect effect); a(ind)(p), activation (indirect effect)(phosphorylation); a(p), activation (phosphorylation); b, binding/association; disso, dissociation; disso(ex), dissociation (expression); disso(inh), dissociation (inhibition); ex, expression; ind(inh), indirect effect (inhibition); inh, inhibition; inh(s), inhibition (state change); inh(u), inhibition (ubiquitination); p, phosphorylation; u, ubiquitination.

Real-time quantitative PCR validation
To validate the reliability of our microarray-based results, 11 miRNAs and eight target mRNAs were selected for validation by qRT-PCR in 12 pairs of matched GHS/normal SD rat kidneys (Table S13). As demonstrated in Fig. 6A, rno-miR-184, rno-miR-484 and rno-miR-138-1-3p were the most significantly dysregulated miRNAs in GHS rats, and the results were comparable with our microarray data. Moreover, the expression levels of the seven mRNAs measured by qPCR (Gcnt1, Lpin2, Olr1, Sema5a, Nfκ B1, Rela and VDR) coincided exactly with our microarray data except for CaSR (Fig. 6B).
These results demonstrate a strong consistency between the microarray-and qRT-PCRbased results.

DISCUSSION
IH typically manifests as hypercalciuria with normal serum Ca 2+ , increased intestinal Ca 2+ absorption, bone resorption and decreased bone mineral density in addition to decreased renal tubule Ca 2+ reabsorption (Yoon et al., 2013). Thus, IH is one of the most important risk factors for calcium urolithiasis. However, the pathogenesis of IH is not yet fully understood. GHS rats exhibit many features of human IH including increased intestinal Ca 2+ absorption, increased bone resorption, decreased renal tubule Ca 2+ reabsorption and high levels of vitamin D receptor (VDR) protein in Ca 2+ -transporting organs (Frick et al., 2013;Frick, Krieger & Bushinsky, 2015). The research of Tsuruoka, Bushinsky & Schwartz (1997) strongly suggests that decreased tubular Ca 2+ reabsorption plays a key role in hypercalciuria. In general, GHS rats represent an ideal animal model for idiopathic hypercalciuria urolithiasis. MicroRNAs regulate various disease processes by inhibiting the expression of their target mRNAs. Understanding the relevance of miRNA and mRNA expression patterns in GHS rat kidneys is important to better elucidate the relationship between pathophysiological process and genes. In the present study, we used bioinformatics methods to determine the potential role of differentially expressed miRNAs and mRNAs in GHS rat kidneys, and identified specific miRNAs and possible negative regulative mRNAs that may affect the development of hypercalciuria.
We used a multi-step approach to identify mRNA targets of dysregulated miRNAs in GHS rat kidneys whereby we first identified mRNAs (N = 2,418) and miRNAs (N = 19) that were significant differentially expressed between GHS and control SD rats; next, potential miRNA-mRNA pairs were predicted for the 19 miRNAs using miRanda, resulting in 29,164 miRNA-mRNA pairs; the intersecting set of predicted target mRNAs and differentially expressed mRNAs were then selected for further analysis.
We found a significant enrichment in over 700 GO terms, including inflammatory response (GO:0006954) and response to hypoxia (GO:0001666). Interestingly, local hypoxic conditions and inflammatory injuries have been linked to the initiation of urolithiasis (Cao et al., 2006), whereas antioxidants protect renal tubular cells from cellular injury and decrease the formation of calcium oxalate stones (Itoh et al., 2005). Moreover, the GO term response to calcium ion (GO:0051592) was enriched (based on differentially expressed genes such as IL-6 and VDR). Increased levels of VDR protein in the kidney are regarded as a common feature in both GHS rats and IH individuals (Frick, Krieger & Bushinsky, 2015), and VDR intimately affects urolithiasis formation (Arcidiacono et al., 2014;Zhang, Nie & Jiang, 2013). Glybochko et al. (2010) found a significant increase in serum IL-6 in nephrolithiasis patients compared with healthy individuals, and Hasna et al. (2015) showed that IL-6 was significantly increased in patients with diabetes mellitus and urolithiasis compared with patients with diabetes mellitus alone. In addition, skeletal system development (GO:0001501) and skeletal system morphogenesis (GO:0048705) were significantly increased in GHS rats. Our previous studies demonstrated the occurrence of osteochondral differentiation progression in primary renal tubular epithelial cells in GHS rats (He et al., 2015), which may be linked to the pathogenesis of calcium stone development.
By identifying pathway membership of dysregulated mRNAs, we can improve our understanding of underlying disease-related processes. We detected 93 enriched pathways including the MAPK signaling pathway, mineral absorption as well as the TGF-beta signaling pathway. Khan (2013) reported that the p38 MAPK/JNK pathway regulates crystallization modulator production and influences plaque formation as well as calcium oxalate nephrolithiasis. Our previous research also showed that calcium and TGF-β may participate in the pathogenesis of epithelial-mesenchymal transition and lead to stone formation (He et al., 2015).
To better understand the biological processes linked to differentially expressed miRNAs and their predicted target genes, we constructed an interaction network where rno-miR-674-5p, rno-miR-672-5p, rno-miR-138-5p and rno-miR-21-3p were found to be highly connected, which means that they may play crucial roles in the regulatory network. Interestingly, Wang et al. (2011) reported that miR-674-5p was able to stimulate the expression of osteogenic marker genes; miR-138-5p is a risk factor for pancreatic cancer (Yu et al., 2015) and Alzheimer's disease (Lugli et al., 2015); and miR-21-3p increases resistance to cisplatin in a range of ovarian cell lines (Pink et al., 2015). In addition, rno-miR-484, rno-miR-138-1-3p, rno-miR-201-3p and rno-miR-203a-3p are downregulated in the network, nevertheless, previous studies reported these miRNAs to be tumor-associated (Pizzini et al., 2013;Liu et al., 2015;Merhautova et al., 2015;Murray et al., 2014;Yang et al., 2016;Ye et al., 2015). To our knowledge, there has been very little research on the relationship between these miRNAs and urolithiasis or calcium metabolism. Hu et al. (2014) reported serum and urinary levels of miR-155 were significantly elevated in patients with nephrolithiasis, and miR-155 might influence pathophysiology of nephrolithiasis via regulating inflammatory cytokines expression. However, the expressions of miR-155 were not significantly different between the two groups of our present study. Regarding mRNAs in the constructed network, Sema5a, which belongs to the semaphorin gene family, had the highest degree of connectivity of 13. It has previously been reported that Sema5a increases cell proliferation and metastasis and suppresses tumor formation (Lu et al., 2010;Sadanandam et al., 2012). In addition, Ding et al. (2008) suggested that Sema5a is a risk factor for Parkinson's disease. Lpin2, which is associated with fatty acid, triacylglycerol, and ketone body metabolism (Chen et al., 2015), also had a high degree of connectivity of 9. Olr1, a crucial regulator of lipid metabolism (Tejedor et al., 2015), had a degree of connectivity of 8. Increasing numbers of studies are reporting a strong correlation between obesity/dyslipidemia and kidney stones (De, Liu & Monga, 2014;Fujimura et al., 2014). However, the mechanism whereby obesity and kidney stone disease are linked is still unknown. Recent research has indicated that obesity may increase reactive oxygen species and oxidative stress, which would influence the interaction between calcium oxalate/calcium phosphate crystals and renal epithelial cells (Khan, 2012).
The NF-kappa B signaling pathway was also significantly enriched in GHS rats; in addition, signal-net analysis has shown that NF-kappa B family might be a core role of our gene regulatory network including the members of NF-kappa B1, NF-kappa B2 and RelA. NF-kappa B is a transcription regulator that is activated by various intraand extra-cellular stimuli such as cytokines, oxidant-free radicals, and bacterial or viral products (Hoesel & Schmid, 2013). Activated NF-kappa B translocates into the nucleus and stimulates the expression of genes involved in a wide variety of biological functions. Menon et al. (2014) found that the receptor activator of NF-kappa B ligand mediates bone resorption in IH, while Tozawa et al. (2008) reported that oxalate induced OPN expression by activating NF-kappa B in renal tubular cells. But there has been little research focusing on the function of NF-kappa B in the formation of idiopathic hypercalciuria urolithiasis. Furthermore, abnormal fat metabolism would also activate the NF-kappa B signaling pathway (Yang et al., 2015), which might contribute to the increased risk of urolithiasis in obesity populations. However, the molecular mechanism of NF-kappa B in IH disease is in need of further exploration, which may assist in improving the clinical diagnosis and treatment of patients with IH.

CONCLUSIONS
Here, we successfully identify rat hypercalciuria-related miRNAs and their target mRNAs. This comprehensive analysis will provide valuable insights for future research, which should aim at confirming the role of these genes in the pathophysiology of hypercalciuria stone desease. More further work needs to be done to get the whole picture right in the disease context. To the best of our knowledge, this is the first study to focus on the role of miRNA-mRNA interactions in hypercalciuria urolithiasis.