Identification of housekeeping genes for microRNA expression analysis in kidney tissues of Pkd1 deficient mouse models

Polycystic kidney disease is a complex clinical entity which comprises a group of genetic diseases that leads to renal cyst development. We evaluated the most suitable housekeeping genes for microRNA expression by RT-qPCR analyses of kidney tissues in Pkd1-deficient mouse models from a panel of five candidates genes (miR-20a, miR-25, miR-26a, miR-191 and U6) and 3 target genes (miR-17, miR-21 and let-7a) using samples from kidneys of cystic mice (Pkd1flox/flox:Nestincre, CY), non-cystic controls (Pkd1flox/flox, NC), Pkd1-haploinsufficient (Pkd1+/−, HT), wild-type controls (Pkd1+/+, WT), severely cystic mice (Pkd1V/V, SC), wild-type controls (CO). The stability of the candidate genes was investigated using NormFinder, GeNorm, BestKeeper, DataAssist, and RefFinder software packages and the comparative ΔCt method. The analyses identified miR-26a as the most stable housekeeping gene for all kidney samples, miR-20a for CY and NC, miR-20a and miR-26a for HT and WT, and miR-25 and miR-26a for SC and CO. Expression of miR-21 was upregulated in SC compared to CO and trends of miR-21 upregulation and let-7a downregulation in CY and HT compared to its control kidneys, when normalized by different combinations of miR-20a, miR-25 and miR-26a. Our findings established miR-20a, miR-25, and miR-26a as the best housekeeping genes for miRNA expression analyses by RT-qPCR in kidney tissues of Pkd1-deficient mouse models.

Before testing the expression stability of the five candidate housekeeping genes, all cDNA samples were normalized at the RNA level. After RT-qPCR cycling, the median Ct values of triplicate reactions were acquired, representing raw expression data. Expression of the miR-20a, miR-25, miR-26a, miR-191 and U6 candidate housekeeping genes is shown in Fig. 2 [24.28 (3.95)], placing it as the highest expressed gene (Fig. 2).
Here, we included RT-qPCR analyses of kidney tissues from three mouse models with distinct patterns of Pkd1 deficiency and correlated them with their respective control tissues. The CY model presented some smaller variations in Ct values compared to its control NC (Fig. 2). miR-26a varied by more than 0.5 cycles in CY group but did not exceed 1.0 and Ct varied less than that for the other genes, except for U6. The HT model showed higher Ct variability compared to its control WT (Fig. 2). U6 Ct values varied by more than 1.0 cycle between HT and WT. Expression differences were far lower between SC and its control (CO kidneys); only the miR-20a difference exceeded 0.3 cycle.
Expression stability analysis of the housekeeping miRNAs. The stability values of the candidate housekeeping genes were obtained applying six software packages (Supplementary Table 2). The top-ranked genes (associated with the smallest stability value) are the most stably expressed ones in the set of kidney tissues samples from the three Pkd1-deficient mouse models. Apart from U6 in the HT group, all candidate housekeeping miR-NAs presented M values below 1.5, the GeNorm set threshold, findings that are consistent with stability 40 . The Bestkeeper software points out inconsistency when SD is higher than 1.0. The SD value was higher than 1.0 only for U6 expression in the All, CY, NC, HT and WT sample groups (Supplementary Table 2).
Based on the different utilized algorithms and a visual inspection of all ranks generated by these analyses, miR-26a seems to be the best housekeeping gene for the All, HT and CO groups; miR-20a for the CY, NC and WT groups; and miR-25 for the SC group (Table 1).
NormFinder recommends a SD value lower than 0.5 for genes to be considered relatively stable. Only miR-20a, miR-25 and miR-26a had an SD value below 0.5, while miR-191 and U6 showed it above 0.5 in the All, CY and NC groups. However, miR-25 presented SD higher than 0.5 in HT and WT and miR-191-related SD lower than 0.5 in SC and CO. While these data represent only a selection of possible tissue pairs (Pkd1 deficiency kidney tissues versus controls), they illustrate that optimal housekeeping genes can significantly vary between kidneys of mouse models with distinct patterns of Pkd1 deficiency (Tables 1, 2 and Supplementary Table 2).
The evaluation of relative expression stability by the BestKeeper software defines the genes that display SD higher than 1.0 as unstable. Again, miR-20a, miR-25 and miR-26a presented SD lower than 1.0 and CV below 3.0 in the All, CY, NC, HT and WT groups (Supplementary Table 2). However, U6 showed SD below 1.0 only for the SC and CO groups, thus it is not a suitable housekeeping gene for the All, CY, NC, HT and WT sample groups (Supplementary Table 2). Based on these results, miR-20a, miR-25 and miR-26a were ranked as the most stable candidate housekeeping genes, whereas miR-191 and U6 were deemed least stable. www.nature.com/scientificreports www.nature.com/scientificreports/ Analysis of the best combination of housekeeping miRNAs. The GeNorm software package recommends at least two genes as combination of housekeeping genes. Table 2 shows the best combination of housekeeping miRNAs for each model/control group pair, based on the different software packages and a visual inspection of all ranks generated by such analyses. The comparisons of Pkd1-deficiency kidney tissues versus their respective controls identified miR-20a for the CY and NC group pair, miR-20a and miR-26a for HT and WT, and miR-25 and miR-26a for SC and CO as the most stable housekeeping genes. The comparison including all groups, in turn, revealed the miR-20a + miR-26a pair as the most stable housekeeping gene selection ( Table 2).
To evaluate the effects of the best candidate housekeeping genes determined by the different algorithms, the expression levels of the top three candidate miRNAs (miR-20a, miR-25 and miR-26a) were normalized by each other (Fig. 3). miR-20a, miR-25 and miR-26a did not differ between the groups ranking for the best ones when normalized by each other (Fig. 3): All comparisons showed no statistically differential expression, but when miR-20 and miR-26 were used as housekeeping genes we observed a more cohesive distribution and equivalent expression in the CY vs NC, HT vs WT and SC vs CO group pairs These results indicate that miR-20 and miR-26a are the most suitable genes to be used together as housekeeping genes in the assessed experimental scenario.
Determination of the suitable number of housekeeping genes. After ranking the candidate housekeeping genes according to their stability, a second approach was applied to determine the optimal number of reference genes to be used in each dataset. This analysis was performed using the GenEx software package. The optimal number of reference genes was calculated using the Acc.SD for the five candidate housekeeping genes (Fig. 4). Estimating the use of the ideal number of genes for normalization, one gene seems to be the best number in HT (mir-26a), WT (mir-20a), SC (mir-25), and CO (mir-26a) groups. The analysis showed that two (mir-20a + mir-26a) is the optimal number of references to be considered for normalization of miRNAs gene expression when using all the samples and the NC group. Three (mir-20a + mir-25 + mir-26a), in turn, was the optimal number of housekeeping genes to be applied to the CY group. We did not observe difference among the numbers of housekeeping genes to be used for the CO group.  www.nature.com/scientificreports www.nature.com/scientificreports/ Correlation between the top three candidate housekeeping miRNAs expression. Correlation analyses were performed using the miRNAs expression data from all evaluated kidney samples. The expression level of the three best candidate housekeeping genes demonstrated a strong correlation between miR-20a and miR-25 (ρ = 0.80, p < 0.001, Fig. 5). Additionally, a moderate correlation was observed between miR-20a and miR-26a (ρ = 0.62, p < 0,001, Fig. 5), as well as between miR-25 and miR-26a (ρ = 0.61, p < 0,001, Fig. 5). These results suggest that, in addition to miR-25 and miR-26a having shown a moderate correlation, they are still correlated in all the samples herein evaluated and can be used together as suitable housekeeping genes. Our results  www.nature.com/scientificreports www.nature.com/scientificreports/ showed that CT data dispersion of miR-20a and miR-26a increases substantially after normalization by U6. These results, were represented in Supplementary Fig. 3.
Validation of the best candidate housekeeping genes for normalizing miR-17, miR-21 and let-7a target genes. In order to validate the three best candidate housekeeping genes stability, the relative expression of miR-17, miR-21 and let-7a target genes was assessed using different combinations of miR-20a, miR-25 and miR-26a (Fig. 6). The expression levels of miR-21 were consistent with upregulation in SC compared to its control CO when normalized by different combinations of the best candidate housekeeping genes (Fig. 6C), trends of upregulation was observed in the CY and HT groups when compared with their respective controls (Fig. 6A,B). Trends of let-7a downregulation in CY, HT and SC compared to its control kidneys, was observed when their expression was normalized by different combinations of miR-20a, miR-25 and miR-26a (Fig. 6D,F). www.nature.com/scientificreports www.nature.com/scientificreports/ Moreover, miR-17 did not present expression statistically different in the CY vs NC, HT vs WT and SC vs CO groups (Fig. 6G,I).
When U6 was used as housekeeping gene, the target genes did not differ between CY vs NC and HT vs WT, except for miR-21 that showed downregulation in CY vs CO, HT vs CO and SC vs CO (Supplementary Fig. 2). These results also suggest that the SC group and its respective control CO have an expression profile slightly different from the other studied ADPKD-related mouse models. Taken together, the use of miR-25 and miR-26a showed to be the most suitable pair of housekeeping genes among all considered sample groups.
To evaluate whether the severity of renal cystic phenotype was a significant modifier of expression of potential housekeeping genes, we investigated potential correlations between the expression levels of miR-20a, miR-25 and miR-26a and kidney weight/body weight (KW/BW) in the CY group. The expression levels of miR-20a, miR-25 and miR-26a were not correlated with the KW/BW (ρ = −0.52, p < 0,18, Supplementary Fig. 3A; ρ = −0.24, p < 0,57, Supplementary Fig. 3B; and ρ = −0.14, p = 0.74, Supplementary Fig. 3C, respectively). The lack of correlation observed between the expression of the evaluated candidate housekeeping genes and kidney weight suggests that the cystic burden does not significantly influence the expression level of the analyzed housekeeping genes. These results provide support to use these miRNAs as controls in studies involving animals with different severities of renal cystic phenotype.

Discussion
In recent years, miRNAs have emerged as key players in tightly-controlled biological processes such as proliferation 37,41,42 , apoptosis [43][44][45][46] and metabolism 47 . MiRNAs are known to be deregulated in numerous kidney diseases including ADPKD [41][42][43] and there is evidence to suggest that miRNA expression profiles may be more accurate in classifying kidney disease progression than mRNA expression profiles 11 . The relative miRNA expression www.nature.com/scientificreports www.nature.com/scientificreports/ quantification is usually compared with a stably expressed housekeeping miRNA from the same sample 33,48 . The high sensitivity of RT-qPCR demands appropriate normalization to correct for non-biological variation and the use of housekeeping genes remains the most commonly used method 49 . However, there are no universally accepted housekeeping transcripts for miRNA data normalization by RT-qPCR analyses. In this scenario, the present study describes the first assessment of candidate housekeeping genes for the normalization of miRNA expression by RT-qPCR in kidney tissue of mouse models with distinct patterns of Pkd1 deficiency. This is a significant problem to be solved, since different animal models orthologous to ADPKD have been and will keep being used to address distinct questions related to ADPKD pathogenesis and interventional studies.
The expression of small nuclear (snRNAs) and nucleolar (snoRNAs) may exhibit tissue-specific and developmental regulation 53 , emphasizing the need for validation of commercially-available control assays. Although U6 is the most commonly used gene to normalize miRNA RT-qPCR data [37][38][39] , we showed that it was the least stable candidate housekeeping gene in the evaluated samples. This result was also observed in other tissues and diseases 35,36,54,55 .
Since miR-20a, miR-25 and miR-26a were the most suitable candidate genes, they were selected for normalization of miR-17, miR-21 and let-7a target genes. The differences in miR-21 expression detected between the tissue groups markedly varied depending on which single housekeeping gene was used for normalization. These results corroborate previous studies that associated miR-21 with increase in cystogenesis and kidney size 20,38 . On the other hand, Gee et al. observed that miR-21 presented variable expression in breast cancer and head and neck squamous cell carcinoma samples when normalized by U6, U44, U48 and U43 54 . In renal cell carcinoma, miR-28, miR-103 and miR-106a were proven to be more stably expressed than U6 56 . Therefore, these findings indicate that the use of unvalidated housekeeping genes may lead to inaccurate and unreliable results, and the use of snRNAs for normalization of expression of miRNAs might introduce bias in the associations between miRNA and the pathology or outcome 50,54 .
The current observed differences in miR-17 expression detected between the tissue groups were at variance with many studies 21,24,37,47 , depending on which single housekeeping gene and animal models with Pkd1 deficiency had been used for normalization. These findings draw attention to the potential effects of the housekeeping gene choice on the outcome of a study and demonstrates the need for validation of candidate housekeeping genes to generate reliable expression data. Even though our analysis suggested that the SC group and its respective control CO have an expression profile slightly different from the other Pkd1-deficient models, the use of miR-25 and miR-26a showed to be the most suitable pair of housekeeping genes among all.
Finally, we assessed the correlation between the best candidate housekeeping genes and KW/BW ratio in the CY group of samples. The detected absence of correlation between the expression levels of the best ones and the KW/BW ratio in the CY group suggests that the level of cystic involvement does not lead to significant changes in housekeeping gene expression. Therefore, miR-20a, miR-25 and miR-26a can be used as housekeeping genes in groups that include different stages of cystic burden in murine models orthologous to ADPKD.
One essential observation is the careful selection and validation of appropriate reference genes as a basis for normalizing the variability between samples in the corresponding study designs to each animal model used in the literature for experimental ADPKD research, amongst them Ksp/cre mice, Pkhd1/cre mice, Kif3a flox/flox mice, Pkd2 flox/flox mice, Pkhd1 −/− mice, Hnf-1β flox/flox mice, miR-17∼92 flox/flox . Based on current findings, we propose an appropriate selection of the best housekeeping genes for each comparison involving one or more mouse models, following the approach adopted in the present study. This procedure consisted in evaluating reference genes in three specific mouse models, namely Pkd1 flox/flox :Nestin cre (CY) Pkd1 +/− (HT) and Pkd1 V/V (SC). The best housekeeping gene selection for each model/control group pair should be used when comparing a specific model with its respective control, while in comparisons including more than one model the best choice should be based on the whole analysis. Such specific approach enables a more proper and reliable normalization for future studies of dysregulated miRNAs within the ADPKD progression cascade. This task obviously implies reliable comparisons between expression data in different stages of the disease, including mild to severe PKD.
One limitation of the present study relies on the selection of only 5 candidate housekeeping genes for validation based on previous studies, but we recognize that more and other genes could also be suitable for accurate miRNA expression normalization by RT-qPCR in ADPKD in in kidney tissue samples of orthologous mouse models and that the analysis of candidate housekeeping genes should be further tested using other models as well.
Previous studies suggested that, unlike mRNAs, the miRNA fraction present in FFPE tissues is relatively unaffected by the fixation process and that miRNAs extracted from these tissues may be accurately profiled using RT-qPCR [57][58][59] . In this context, the housekeeping genes identified in this study may also prove useful for miRNA RT-qPCR analysis of FFPE kidney tissues. Our findings will allow further analysis in kidney tissue of Pkd1-deficient mouse models gene expression to elucidate the role of different regulatory miRNAs in different scenarios of or related to ADPKD.
Normalizing to a suitable housekeeping gene, therefore, can eliminate differences due to sampling and quality of RNA and identify real changes in miRNA expression levels. In the present study, we analyzed a series of candidates and identified the most suitable housekeeping genes to be used for miRNA expression normalization in RT-qPCR studies in kidney tissues of mouse models orthologous to ADPKD and their respective controls. www.nature.com/scientificreports www.nature.com/scientificreports/ The housekeeping gene spectrum and data generated by our work should therefore be employed in miRNA-related studies involving Pkd1-deficient mouse models. Among the genes currently used in this field, appropriate combinations of the miR-20a, miR-25, and miR-26a housekeeping genes offer increased accuracy and resolution in the quantitation of gene expression data, favoring the detection of smaller changes in miRNA expression than otherwise possible. The proper selection of the best housekeeping genes to be used in each of these scenarios should follow the guidelines specific to each comparison. It must be noted that other Pkd1-deficient models distinct from the ones analyzed in the current study may display slight expression differences of reference genes. Even in these cases, however, our guidelines would still be the best available housekeeping controls for such studies.
Methods ADPKD mouse models. We used three mouse models with distinct Pkd1-deficiency profiles in the current study, all generated and maintained on the C57BL/6 background. Only kidneys from male animals were analyzed in order to avoid potential gender-related experimental variability. Two models were evaluated at 10-12 weeks of age, including a renal cystic mouse (Pkd1 flox/flox :Nestin cre , CY, n = 10) and its corresponding non-cystic control (Pkd1 flox/flox , NC, n = 10), and a Pkd1-haploinsufficient mouse (Pkd1 +/− , HT, n = 6) and its respective wild-type control (Pkd1 +/+ , WT, n = 6). The third model was assessed at 15 days of life due to its severely renal cystic phenotype (Pkd1 V/V , SC, n = 7) along with its wild-type control (CO, n = 5). These animals were genotyped using specific PCRs. The CY mouse is homozygous for a Pkd1-floxed allele (Pkd1 flox ) and displays a mosaic pattern of gene inactivation, induced by a Nestin-Cre transgene through excision of exons 2-4 (Pkd1 flox/flox :Nestin cre ) 25,[27][28][29] . The HT model is heterozygous for a Pkd1 null allele, characterized by early transcriptional interruption 25,26 , and develops no renal cysts by 12 weeks of age. The SC model is homozygous for the Pkd1 knockin T3041V allele, which prevents the autoproteolytic cleavage of PC1 at the GPS site 25,30 . Pkd1 V/V animals have no gross phenotype by postnatal day (P) 6 but develop rapid and progressive distal nephron dilatation thereafter. This severe renal phenotype, that eventually leads to renal failure, is apparently responsible for the early mortality that occurs between the 2nd and 6th week 30  SC (Pkd1 V/V ) and its wild-type control (CO) animals were euthanized by adopting cervical dislocation and the other animal groups were euthanized with intraperitoneal thiopental (0.4 mg/g of body weight); their kidneys were appropriately harvested for RT-qPCR analyses. All experiments were conducted in accordance with international standards of animal care and experimentation. Both kidneys were collected and stored at −80 °C for further use.
Housekeeping genes. We selected these genes based on the observation that they were among the miRNAs with the lowest variances in a previous study of ours (TLDA Taqman Array, unpublished data) and that miRNA expression studies showed minor evidence of differential expression 10,20,37-39 . These genes not only behaved stably but also presented expression stability in transcriptomic analysis. Of note, U6 and miRNA191 presented stable expression, also documented in a previous report 33 . All these five miRNAs are constitutively expressed in kidney tissue of mouse models orthologous to ADPKD, have independent cellular functions, and are assumed not to be co-regulated. RNA extraction. Renal tissue lysis was performed using zirconia beads (Interprise, USA) and the Precellys (BioAmerica, USA) homogenizer. TRIzol (Life Technologies, USA) was employed for total RNA extraction according to the manufacturer's protocol. The RNA quantity and purity were determined using the NanoVue spectrophotometer (GE Healthcare Life Sciences, USA) and analyzed with the Agilent 2100 Bioanalyzer 6000 Nanochip (Agilent Technologies Inc., Waldbronn, BW, Germany). Total RNA was stored at −80 °C until further use. cDNA preparation and RT-qPCR. Complementary DNA (cDNA) was synthesized from 1 µg of total RNA with an oligonucleotide pool for each evaluated gene (Supplementary Table 1), using the TaqMan ® microRNA reverse transcription kit (Life Technologies) according to the manufacturer's instructions. Detection of the expression range of the evaluated genes was achieved by using TaqMan ® assays and the QuantStudio ® 7 Flex real-time PCR system (Applied Biosystems, USA). The used primers are shown in Supplementary Table 1. miR-17, miR-21 and let-7a were employed as target genes [20][21][22][23][24] . All reactions were run in triplicate. The expression of the candidate housekeeping genes is represented by the original cycle threshold (Ct) value.
The optimal number of reference genes was selected using the GenEx software package. GenEx is a software for the processing and analysis of RT-qPCR data and provides methods to select and validate housekeeping genes, classify samples, group genes and monitor time-dependent processes. GenEx calculates Accumulated standard deviation (Acc.SD); lower Acc.SD values, in turn, indicate the optimal number of reference genes (https://www. biomcc.com/genex-software.html).