Small RNA sequencing evaluation of renal microRNA biomarkers in dogs with X-linked hereditary nephropathy

Dogs with X-linked hereditary nephropathy (XLHN) are an animal model for Alport syndrome in humans and progressive chronic kidney disease (CKD). Using mRNA sequencing (mRNA-seq), we have characterized the gene expression profile affecting the progression of XLHN; however, the microRNA (miRNA, miR) expression remains unknown. With small RNA-seq and quantitative RT-PCR (qRT-PCR), we used 3 small RNA-seq analysis tools (QIAGEN OmicSoft Studio, miRDeep2, and CPSS 2.0) to profile differentially expressed renal miRNAs, top-ranked miRNA target genes, and enriched biological processes and pathways in CKD progression. Twenty-three kidney biopsies were collected from 5 dogs with XLHN and 4 age-matched, unaffected littermates at 3 clinical time points (T1: onset of proteinuria, T2: onset of azotemia, and T3: advanced azotemia). We identified up to 23 differentially expressed miRNAs at each clinical time point. Five miRNAs (miR-21, miR-146b, miR-802, miR-142, miR-147) were consistently upregulated in affected dogs. We identified miR-186 and miR-26b as effective reference miRNAs for qRT-PCR. This study applied small RNA-seq to identify differentially expressed miRNAs that might regulate critical pathways contributing to CKD progression in dogs with XLHN.

. Canine kidney tissue sequencing data with averaged read length distribution (bar chart) and genome mapping results (pie charts). The pie charts indicate almost all reads obtained from 23 kidney tissues mapped to the dog genome. Among the mapped reads, the vast majority (90.8%) belonged to miRNA, and a small population mapped to non-coding RNA (Rfam database; 7.7%) or mRNA (0.5%). The bar chart shows that most sequencing reads were between 21 and 25 nucleotides in length, consistent with miRNAs. The results of CPSS 2.0 were used. www.nature.com/scientificreports/ time point. The results of the geNorm analysis showed stable expression for 3 of the 4 promising internal controls (miR-186, miR-26b, and miR- 16) and indicated that the 2 most stable miRNAs were miR-186 and miR-26b ( Supplementary Fig. S2a,b). In the "reference target stability" quality control of qbase +, we further verified the stability of miR-186 and miR-26b and confirmed that they are reliable for normalization using qRT-PCR (Supplementary Fig. S2c). The upregulation of miR-21, miR-142, and miR-147 in the kidney tissues of affected dogs at T2 and T3 was detected by small RNA-seq using CPSS 2.0 (Fig. 5) as well as using qRT-PCR ( Supplementary Fig. S3). The expression of miR-21, miR-142, and miR-147 based on qRT-PCR showed upward trending as kidney disease progressed in the affected dogs whereas expression in the control dogs was similar at all time points (Supplementary Fig. S3). The qRT-PCR results, along with the higher genome mapping rate, prompted us to use the data generated by CPSS2.0 for the remaining analysis.
Target prediction, gene ontology (GO), and pathway analysis for differentially expressed miRNAs. The differentially expressed miRNAs identified using CPSS2.0 are listed in Table 2 (adjusted P-value < 0.05, log2 fold change > 1 or < − 1). Among them, the top 10 differentially expressed miRNAs were selected as inputs for miRDB 13,14 . To better characterize the changes in target genes, upregulated and downregulated differentially expressed miRNAs were used as 2 separate inputs. In miRDB, only putative target genes with target scores ≥ 90 were included (Supplementary Table S2) and considered satisfactory for the subsequent overrepresentation test in the PANTHER gene ontology (GO) and pathway analyses 15,16 .   Table S3). "Cellular process" and "RNA metabolic process" and subfamilies of "MAPK cascade" and "regulation of transcription from RNA polymerase II promoter" were the main GO terms enriched in the current study. Specifically, "regulation of transcription from RNA polymerase II promoter" was upregulated at T1 and T3, and "intracellular signal transduction" and "regulation of phosphate metabolic process" were upregulated at T2 and T3.
For Reactome pathway analysis, the "signal transduction" pathway was enriched in affected dogs throughout the progression of kidney disease, and its subfamily pathway the "signaling by TGF-β receptor complex" pathway was also enriched approximately fivefold at T1. The PANTHER pathways, such as the gonadotropin-releasing hormone receptor pathway and Wnt signaling pathway, were enriched in affected dogs at T1 and T2. Ten other  Table 1 for a complete list of the consistently differentially expressed miRNAs in the center of the Venn Diagrams. Only miRNAs with an adjusted P-value < 0.05 and a log twofold change > 1 or < -1 were considered differentially expressed miRNAs. (Pink: T1; Blue: T2; Green: T3).    www.nature.com/scientificreports/ enriched PANTHER pathways were identified at T2, but no PANTHER pathway was identified at T3, presumptively because there were fewer miRNA targets (Supplementary Table S3).

Discussion
Dogs with XLHN have been studied as an example of canine CKD and used as an animal model for human Alport syndrome. The gene expression in XLHN dogs has been partially characterized using qRT-PCR 17 and microarrays 18 . Previously, our group used mRNA-seq to investigate the gene expression linked to rapid CKD progression in dogs with XLHN 3 . In the current study, we performed small RNA-seq on 23 renal biopsies collected serially from 5 affected dogs and 4 healthy littermates to characterize miRNA expression during CKD progression. The majority of differentially expressed miRNAs in the kidney tissue of affected animals were upregulated, as in previous studies using CKD mouse models 19,20 . We identified up to 23 miRNAs that were differentially expressed at specific clinical time points including 5 miRNAs (miR-21, miR-146b, miR-802, miR-142, miR-147) that were consistently upregulated at all 3 time points (Table 2). Meanwhile, we compared 3 miRNA analysis tools and identified promising miRNA internal controls (miR-186 and miR-26b) for qRT-PCR normalization using canine kidney tissue. Several studies have documented increased renal miR-21 in different models of CKD, including, but not limited to, Alport syndrome mouse model 21 , B6.MRLc1 mouse model 19 , type 2 diabetes mouse model 22 , and IgA nephropathy patients 23 . MiR-21 is regulated by TGF-β1/Smad3 24-26 and contributes to renal fibrosis by silencing metabolic pathways 27 . In the current study, we used qRT-PCR to examine the upregulation of renal miR-21 at T2 and T3 in affected dogs. We found that the expression of renal miR-21 did not display an upward trending at the onset of azotemia compared with more advanced azotemia. This is similar to our previous finding in a larger group of XLHN dogs where, based on qRT-PCR, the expression of renal miR-21 did not increase significantly until affected dogs became azotemic 28 . Our results also corroborate our previous finding that TGF-β1 is the top upstream regulator of CKD progression in XLHN dogs 3 .
MiR-146b also has been previously described in the context of CKD. MiR-146b was one of the upregulated miRNAs in the kidney tissue of 12-month-old mice with spontaneous CKD (B6.MRLc1) compared with healthy controls (C57BL/6) 19 . The upregulation of renal miR-146b has also been associated with 4 additional kidney conditions in mouse models (folic acid-induced kidney injury, unilateral ureteral obstruction, bilateral renal ischemia/reperfusion, and cisplatin-induced kidney injury) with peak expression associated with fibrosis 20 . In the current study, miR-146b reached peak expression at T2, when XLHN dogs became azotemic and had visible fibrotic change on histopathological evaluation 3 . In dogs, mice, and humans, miR-146b and miR-146a are highly homologous miRNAs that differ by only 2 nucleotides. Therefore, they share similar mRNA targets and biological functions 19 . Upregulation of miR-146a has been seen in the kidney tissues of IgA nephropathy patients 29 , a CKD mouse model 19 , a diabetic nephropathy rat model 30 , and human glomeruli in lupus nephritis 31 and membranoproliferative glomerulonephritis 32 . In the current study, both miR-146a and miR-146b showed an upward trending at T2 in dogs with XLHN, indicating that the function of miR-146a/b is worth further investigation.
Other miRNAs identified in our study, including miR-802, miR-142, and miR-147, have been only rarely described in the context of CKD. In our study, the expression of renal miRNA-802 was upregulated at all 3 time points in the dogs with XLHN. This increased expression of miR-802 has been observed in the kidneys, specifically in the cortical collecting ducts, of mice exposed to high-potassium diets 33 . In vitro, miR-802 targets caveolin-1 (CAV1), decreasing caveolin-1 expression, which in turn increases the surface expression of the renal outer medullary potassium channel and facilitates potassium excretion 33 . In dogs with XLHN, we found CAV1 gene expression increased at T2 and T3 3 , although we would expect the putative target for miR-802 to be downregulated. More studies on canine XLHN miR-802 expression are needed to resolve this discrepancy.
Renal miR-142 (hsa-miR-142-3p) is upregulated in human patients with acute rejection of a renal allograft 34 . Another study suggests the upregulation of miR-142 in renal allografts is due to an influx of lymphoid cells with acute T-cell mediated rejection 35 . Previously, we have demonstrated that T cells are the predominant lymphoid cells infiltrating the kidneys of dogs with XLHN 3 , which could explain the upward trending of miR-142 in the current study. Upregulation of renal miR-147 was detected by RNA-seq in a 3-chloro-1, 2-propanediol (3-MCPD)-induced acute kidney injury mouse model, but its significance has not been discussed 36 . Overall, research on miR-147 and its role in the development of CKD is lacking. To our best knowledge, this study is the first to document upregulation of miR-147 in CKD, which reemphasizes the importance of an unbiased approach using sequencing technology for miRNA discovery 37 .
Using GO terms and pathway analyses, we identified "MAPK cascade" and "regulation of transcription from RNA polymerase II promoter" as the enriched subfamilies of biological processes. Both of these were enriched in patients with CKD and other diseases in previous studies 38,39 . The Reactome pathway analysis supported involvement of the "signal transduction" pathway, particularly the "signaling by TGF-β receptor complex, " consistent with our previous finding that TGF-β1 is the top upstream regulator of CKD progression in XLHN dogs 3 . Additionally, based on PANTHER pathway analysis, several pathways identified (e.g., "gonadotropin-releasing hormone receptor" and "Wnt signaling") were also related to signal transduction. Indeed, these pathways contained several genes that were downregulated in dogs with XLHN in our previous mRNA-seq study (e.g., ACTG2, ADCYAP1R1, and CACNA1C at T1; MAP3K3, FZD4, ACVR1B, PER1, GATA2, MYH15, and GATA4 at T2) 3 .
Although several tools are available for analysis of RNA sequences, there is a lack of comparisons among these tools. We therefore used 3 analysis tools (CPSS 2.0, OmicSoft Studio, and miRDeep2) and the same version of the canine genome and miRNA annotations with the default settings. Among the 3 analysis tools, miRDeep2 is most widely used for miRNA identification; however, both CPSS 2.0 and OmicSoft Studio surprisingly detected more miRNAs than miRDeep2. With our dataset, the use of CPSS 2.0 resulted in a significantly higher genome mapping rate than that of OmicSoft Studio and miRDeep2. There are several possible reasons for the differences  Fig. S1). And, CPSS 2.0 detected similar numbers of miRNAs as OmicSoft Studio but outperformed miRDeep2 ( Supplementary Fig. S1). Also, at each time point, CPSS 2.0 identified more or equal numbers of differentially expressed miRNAs compared to OmicSoft Studio and miRDeep2 (Fig. 3). In our hands, CPSS 2.0 appears to be the preferred analysis tool for small RNA-seq analysis. However, further evaluation is needed as the conclusion cannot be made using a single dataset. In all 3 analysis tools, miR-21 was repeatedly detected at all time points, whereas both miR-142 and miR-147 were detected by only 2 of the 3 analysis tools. Therefore, we incorporated qRT-PCR to detect the expression of miR-21, miR-142, and miR-147 ( Supplementary Fig. S3). The qPCR results support the miR-142 and miR-147 expression detected by CPSS 2.0 but additional samples are needed for verification. In general, PCR is often used to verify sequencing results. It has been argued that qRT-PCR suffered from GC-bias 42 and inconsistency among different assays 43 . Also, the expression level of a particular miRNA depends on the normalization method used. Although small nuclear RNA (snRNA) (e.g., U6 snRNA) are frequently used for qRT-PCR normalization, snR-NAs differ structurally and functionally from miRNAs. The difference in nucleic acid composition, length, and secondary structure could potentially introduce variation between a snRNA control and miRNAs of interest 44 . In order to identify miRNAs for use in normalization, we first applied NormFinder 45 to the small RNA-seq data and then selected 4 miRNAs that appeared to be stably expressed (miR-186, miR-26b, miR-16, and miR-99a) for further assessment. Based on our qRT-PCR results, miR-186, miR-26b, and miR-16 had similarly low M values in the geNorm analysis ( Supplementary Fig. S2), and all of these have been proposed as internal controls in previous renal studies 10,46 . In particular, miR-186 was identified as one of the most stably and ubiquitously expressed miRNAs across 16 types of canine tissues, including kidney 10 . Also, miR-26b has been described as a suitable internal control for glomerular miRNA quantification in IgA nephropathy patients 46 . Our data support miR-186 and miR-26b are more suitable for reference miRNAs for future studies in canine kidneys.
One limitation of the current study was the inability to perform microscopic evaluation and RNA isolation on the same kidney biopsy. If we were able to evaluate the sample used for small RNA-seq by light microscopy, we might have found that the unexpected similarity in renal miRNA expression patterns between a few of the affected dogs at T1 and T2 and the controls (Fig. 2) could be explained by having sampled an area of the kidney that was less affected by the disease than expected, for example. Another limitation of our study is the implementation of in silico analysis and the lack of verification of the miRNA targets and performance of functional studies. To minimize false positive results, we have chosen the most credible tool, miRDB 13,14 , and adapted the most stringent rule of a target score higher than 90, out of the range of 50-100. If samples were available, the study could potentially benefit from in situ hybridization to confirm the downregulation of predicted miRNA targets. Lastly, the miRNA expression profile used whole kidney cortex and therefore mostly represents the tubulointerstitium. Future studies could use laser capture microdissection or single-cell RNA-seq to help narrow down the cell types of interest to achieve a more specific miRNA expression.
In the current study, we applied 3 analysis tools (QIAGEN OmicSoft Studio, miRDeep2, and CPSS 2.0) to identify differentially expressed miRNAs in dogs with CKD using small RNA-seq and qRT-PCR. Based on our data, CPSS 2.0 seems to have advantages over the other two tools in its high genome mapping rate and the comprehensive identification of differentially expressed miRNAs. We identified 5 miRNAs (miR-21, miR-146b, miR-802, miR-142, miR-147) that were consistently upregulated in dogs with CKD compared to age-match controls. The upregulation of miR-21 and miR-146b correspond to the onset of azotemia and could be associated with fibrosis, whereas miR-802, miR-142, and miR-147 were novel miRNAs that warrant further investigation in the context of CKD. The putative targets of these miRNAs and the enriched pathways aligned well with our mRNA-seq study 3 and further characterize the development and the progression of CKD that involves the signal transduction and TGF-β receptor complex. In both studies, we acquired biopsies at defined clinical stages of disease in the same dogs, which was helpful to accurately compare miRNA expression between controls and dogs with XLHN throughout disease progression. MiRNAs found in the current study may be diagnostic biomarkers for CKD and potential therapeutic targets for CKD in both dogs and humans.

Methods
Animals. The dogs evaluated in this study were part of a colony with XLHN maintained at Texas A&M University 2 . XLHN is caused by a 10-base deletion in the gene encoding the ⍺5 chain of type IV collagen. The affected males develop juvenile-onset CKD that progresses to end-stage renal disease, as previously described 2 . Overall, 9 dogs were included in the study. At each clinical time point, 3-5 affected dogs and 4 age-matched unaffected littermates were included (Supplementary Table S1). All dogs were raised according to standardized protocols, and no treatments were given to these dogs. All animal experiments were performed in accordance with the relevant guidelines and regulations set forth by the Texas A&M University Institutional Animal Care and Use Committee represented by the approved study protocol #2010-132 and #2005-129. The study was carried out in compliance with the ARRIVE guidelines (http:// www. nc3rs. org. uk/ page. asp? id= 1357, last accessed on May 24, 2021).

Clinical phenotypes of XLHN dogs.
Clinical progression in the 5 affected dogs was determined by serial monitoring of serum and urine biomarkers of kidney disease, which allowed comparison of dogs at standardized clinical time points 47  Small RNA-seq data analysis. Preprocessing of small RNA-seq reads and quality control were performed by the sequencing facility (BGI). Low quality reads were defined as reads with more than 4 bases with a quality score less than 10 or reads with more than 6 bases with a quality score less than 13. Briefly, low quality reads, adapter contaminants (no insert reads and 5' primer contaminants) and reads less than 18nt were filtered out. Length distribution was inspected to confirm the composition of the small RNAs as shown in Fig. 1 . For all 3 tools, un-normalized raw read counts were used to perform differential expression and statistical analysis with the identical script using DESeq2 48 (release 3.3) in R (version 3.6.2) as previously described 3 . The 3D PCA plots were made using Plotly (https:// plotly. com/r/). We used the Benjamini-Hochberg procedure 49 to adjust P-values for multiple testing. An adjusted P-value < 0.05 was set for the selection of differentially expressed miRNAs.
Target prediction, GO terms, and pathway analysis. For differentially expressed miRNAs at each time point, we selected the top 10 upregulated (fold change > 2) and down-regulated (fold change < 2) miRNAs as input for miRDB 13,14 . miRDB is an online database for miRNA target prediction that hosts predictive miRNA targets in 5 species: humans, mice, rats, dogs, and chickens. There are 453 miRNAs targeting 170,435 genes in the database (http:// mirdb. org/, last accessed on Dec 30, 2020). All miRNA targets are assigned a target score from 50 to 100 determined by the MirTarget algorithm 50 . We used the target score of 90 as a cutoff value since higher scores represent more statistical confidence. The list of miRNA targets was passed on to gene ontology (GO) and PANTHER pathway analysis using the Canis familiaris reference list (20, qRT-PCR. Small RNA-seq data generated by the 3 analysis tools were analyzed using NormFinder 45 to identify suitable miRNAs as candidates for internal controls in qRT-PCR. Samples remaining after small RNA-seq (n = 10) were used for qRT-PCR. Additionally, control samples used in a previous study 3 (n = 2) were recruited. Overall, there were 2 samples from the affected group and 2 from the control group at each time point, for a total of 12 samples. The miRCURY LNA miRNA PCR Assays (Qiagen, Germany) were used for the miRNA targets: miR-21 (Cat. no. YP00204230), miR-142 (Cat. no. YP02102101), and miR-147 (Cat. no. YP00204368). The miRCURY LNA miRNA PCR Assays (Qiagen, Germany) were used for the miRNA internal controls: miR-186 (Cat. no. YP00206053), miR-26b (Cat. no. YP00205953), miR-16 (Cat. no. YP00205702), and miR-99a (Cat. no. YP00204521). First, total RNA isolated from kidney biopsies was diluted to the concentration of 10 ng/μL. Next, the miRCURY LNA RT Kit (Qiagen, Germany) was used in a 10 μL reaction for RT consisting of 2 μl of 5X Reaction Buffer, 5 μL RNase-free water, 1 μL enzyme mix (omitted for no reverse transcriptase control wells), and 2 μL diluted RNA (10 ng/μl). The RT reaction was performed using a T100 Thermal Cycler (Bio-Rad, UK) with a protocol of 60 min at 42 °C (reverse transcription), 5 min at 95 °C (inactivation), followed by storage at -20 °C. Although the cDNA is stable for up to 5 weeks, all PCRs were performed within 18 days from the RT reaction. For PCR, the miRCURY LNA SYBR Green PCR Kit (Qiagen, Germany) was used to make a 10 μL reaction consisting of 5 μL SYBR Green www.nature.com/scientificreports/ System (Eppendorf, Germany). PCR was performed using a CFX384 Touch Real-Time PCR Detection System (Bio-Rad, UK) with a protocol of 2 min at 95 °C (initial heat activation), 40 cycles of 10 s at 95 °C (denaturation) and 60 s at 56 °C (annealing), followed by a melting curve analysis at 60-95 °C. Negative controls included RNase-free water only, no reverse transcriptase (NRT), and no template controls (NTC) to ensure no genomic DNA contamination was present. For the internal controls, the Cq values of the candidate miRNAs were further analyzed by geNorm 51 in the qbase + software 52 . Candidate miRNAs for internal controls were ranked by geNorm 51 based on stability (M value) and coefficient of variation. The algorithm of geNorm 51 then calculates the normalization factor (V value) and determines the optimal number of internal controls. Lastly, the candidate miRNAs were analyzed for "reference target stability" quality control in qbase + 52 . Target miRNAs were normalized based on the 2 −ΔΔCq method 53 , using miR-186 and miR-26b as reference miRNAs.

Data availability
Raw sequencing reads generated from this study are deposited at the NCBI sequence read archive (SRA) under BioProject ID PRJNA664365.