RNA-seq of serial kidney biopsies obtained during progression of chronic kidney disease from dogs with X-linked hereditary nephropathy

Dogs with X-linked hereditary nephropathy (XLHN) have a glomerular basement membrane defect that leads to progressive juvenile-onset renal failure. Their disease is analogous to Alport syndrome in humans, and they also serve as a good model of progressive chronic kidney disease (CKD). However, the gene expression profile that affects progression in this disease has only been partially characterized. To help fill this gap, we used RNA sequencing to identify differentially expressed genes (DEGs), over-represented pathways, and upstream regulators that contribute to kidney disease progression. Total RNA from kidney biopsies was isolated at 3 clinical time points from 3 males with rapidly-progressing CKD, 3 males with slowly-progressing CKD, and 2 age-matched controls. We identified 70 DEGs by comparing rapid and slow groups at specific time points. Based on time course analysis, 1,947 DEGs were identified over the 3 time points revealing upregulation of inflammatory pathways: integrin signaling, T cell activation, and chemokine and cytokine signaling pathways. T cell infiltration was verified by immunohistochemistry. TGF-β1 was identified as the primary upstream regulator. These results provide new insights into the underlying molecular mechanisms of disease progression in XLHN, and the identified DEGs can be potential biomarkers and therapeutic targets translatable to all CKDs.

characterized. Furthermore, dogs with the same mutation causing XLHN display substantial variation in the rate of disease progression such that some dogs reach end-stage disease by 6 months of age and others at 12 months of age or later. Although varied times of onset and rates of progression are common among different types of mutations in people with AS 6 , disease progression may also vary among members of an AS family with an identical mutation 7 , as seen in dogs 8 .
While several studies have characterized gene expression in humans with CKD and in animal Alport models using microarrays 9,10 or PCR 8,9,11,12 , studies that have incorporated high-throughput RNA sequencing (RNA-seq) with the objective of identifying differentially expressed genes (DEGs) and upstream regulators are lacking. Compared with traditional approaches in gene expression analysis, RNA-seq provides unprecedented flexibility in the discovery of DEGs 13 while preserving accuracy and strong correlation with PCR [14][15][16][17][18] , even considering fold change levels 19 .
The objective of this study was to compare the gene expression between dogs with rapid versus slow disease progression phenotypes at 3 stages of the disease. We conducted Gene Ontology (GO) and pathway analyses to characterize DEGs among sample groups at specific time points. Since all CKDs share common pathways that lead to end-stage kidney disease 20 , the results help elucidate the molecular basis of CKD progression and thus may benefit canine patients and indicate potential therapeutic targets for AS patients.

Results
Histopathological evaluation of kidney biopsies. Figure 1 presents representative cortical fields of the kidney biopsies and the mean interstitial fibrosis scores comparing the rapid versus slow groups. Clinical time points in affected dogs were defined as: T1 -onset of proteinuria (the earliest time point that clinical disease can be detected); T2 -onset of azotemia (serum creatinine ≥ 1.2 mg/dL); and T3 -end-stage disease (serum creatinine ≥ 5 mg/dL). At both T2 and T3, mean fibrosis and chronic inflammation scores were significantly higher in diseased dogs than in controls (Supplementary Tables S1 and S2). The degree of fibrosis in the rapid group was more severe at T2 and T3 than that in the slow group, despite the two groups being clinically indistinguishable; however, statistical significance between the scores was reached only at T3 (Supplementary Table S1 and Fig. S1). No statistically significant difference was observed for the chronic inflammation score between the rapid and slow groups (Supplementary Table S2 Table S3). Because of the variable quality of RNA, the proper library preparation kit was used to compensate for the low-input samples according to the best practice for RNA with variable qualities 21 (see "RNA isolation and sequencing" section in Methods). After performing quality control, we obtained an average of over 30 million paired-end reads from each sample (n = 24). Overall, 91-96% of reads were mapped to the canine genome (CanFam 3.1) by HISAT2 22 . Among them, 70-78% of reads were uniquely mapped (Supplementary Table S3 and Fig. S3). Based on the union setting of HTSeq 23 , ambiguous reads that mapped to multiple genes were not included in our analysis.
Principal component analysis (PCA) and hierarchical clustering analysis. We performed PCA at each time point to determine whether samples in each group clustered with each other or other groups. First, we used HTSeq 23 to count reads that uniquely aligned to one gene, and these data were then imported into DESeq2 24 to generate PCA plots. At T1, the PCA results demonstrated that most samples clustered together, regardless of the grouping (Fig. 2). Except for one dog, rapid and slow groups became separated from controls at T2. Although the slow group tended to be closer to the controls than the rapid group, there was no clear distinction between rapid and slow groups at T2 or T3. Furthermore, PCA scree plots confirmed that principal components 1 (PC1) and 2 (PC2) accounted for 76-90% of the total variation in gene expression at each time point ( Supplementary  Fig. S4). To further investigate the time-dependent nature of the DEGs, we performed hierarchical clustering of the top 100 DEGs (i.e., those with the smallest q-values identified in the time course analysis in DESeq2). In agreement with the PCA plots, this analysis demonstrated clustering of almost all sample groups at T1 (Fig. 3a). At time points T2 and T3, the rapid and slow groups clustered together (Cluster 1 in Fig. 3a) and were distinctly separated from the control group for all but one T2 sample (Cluster 2 in Fig. 3a).
Differentially expressed genes (DEGs). In average, 20,090 genes were mapped by at least one read in each of the kidney biopsy samples (Supplementary Table S3). Overall, 1,947 DEGs with a q-value < 0.05 were detected over the 3 time points in the time course analysis of DESeq2 (Supplementary Table S4). We applied the plot counts function in DESeq2 to visualize the top 10 genes with the smallest q-values (Fig. 3b). While these genes were not differentially expressed at T1, group-specific changes were observed over time, and expression in the slow group was consistently closer to that in the control group for each gene at T2 (Fig. 3b).
We also compared the rapid and slow groups with the control group, both individually and combined as a single "affected" group (Fig. 4a,c, and d). In these comparisons, the number of DEGs increased with advancing disease, with the largest number of DEGs identified at T3. This phenomenon indicates that the DEGs are disease-dependent as they are more differentially expressed in the later time points (T2 and T3) than T1 (Fig. 4).
A substantial overlap of DEGs was present when comparing rapid and slow groups with the control group at T2 and T3 (Fig. 5). The overlapping DEGs between the rapid and slow groups were more numerous at T3 (2,952 DEGs) than at T2 (190 DEGs), supporting that the two groups behave similarly at the end-stage disease, as expected based on Figs 2-4. Furthermore, the number of DEGs (1,189 DEGs) identified in both T2 and T3 in the rapid group was higher than the number of DEGs (171 DEGs) identified in both T2 and T3 in the slow group. This supports the theory that rapidly-progressing dogs express end-stage DEGs at a young age. The complete lists of DEGs from the time course analysis and all pairs of comparisons appear in Supplementary Table S4.
Gene Ontology (GO) and pathway analysis of DEGs. To characterize the GO terms, including molecular functions, biological processes, cellular components, and functional pathways of DEGs, we conducted over-representation tests for all pairs of comparisons in PANTHER version 11.1 (released on October 24, 2016) ( Fig. 6 and Supplementary Table S5). We used the GO-Slim PANTHER annotation data set, which represents phylogenetically inferred annotations 25 .
Overlaps of GO terms among comparisons were commonly seen in the current study. We focused on "biological process" since it is the most characterized GO term. Within this category, the "immune system process" family was upregulated in all 10 comparisons presented in Fig. 6, and the "immune response" family was upregulated

Down-regulated DEGs in the rapid group
Continued in all except for the rapid versus slow comparison at T2 and the time course analysis. Both immune-related GO terms appeared to be more upregulated at T2 than T3. At T2, the "biological adhesion" family, especially the cell-cell adhesion subfamily, was expressed in the rapid group more than in the slow group (Supplementary Table S5). The most common pathway represented by the DEGs within the various comparisons was the "integrin signaling pathway" (Fig. 6). This was the only pathway identified in the comparison between rapid and slow groups. It was also the top upregulated pathway at T2 comparing rapid and control groups. Interestingly, analyzing the overlapping DEGs between T2 and T3, the integrin signaling pathway was identified as the top upregulated pathway within the rapid group, but not within the slow group, as compared with control.
The "T cell activation pathway" was another frequently detected pathway in our study (Fig. 6). It was upregulated in all comparisons with control at T2 and in overlapping DEGs between T2 and T3 within both the rapid and the slow groups. Another upregulated pathway included the "inflammation mediated by chemokine and cytokine signaling pathway" that was identified only when comparing all affected dogs with controls at T2. Because only limited numbers of DEGs were discovered, no enriched pathways were identified at T1 with any comparison.
To further explore the possible biological interaction between orthologous genes in the human, mouse, and rat, we performed Ingenuity Pathway Analysis (IPA) to discover the most prevalent pathways and upstream regulators within each comparison. The "hepatic fibrosis/hepatic stellate cell activation pathway" was identified as the top pathway in multiple comparisons (Supplementary Table S6). Transforming growth factor beta 1 (TGF-β1) was the most activated upstream regulator when the rapid and slow groups were compared at T2 and in the time course analysis (Supplementary Table S6).

Immunohistochemistry (IHC) validation of inflammatory pathways.
Lastly, we aimed to validate the overexpression of inflammatory pathways in kidney biopsies. We chose to identify the predominant lymphocyte subtype present (T versus B cell) since the identified inflammatory pathways are all closely related to the presence of T lymphocytes 26 and antibodies for CD3 (T cell) and CD20 (B cell) are validated for use in dogs 27,28 . As shown in Fig. 7, we confirmed the lymphocyte infiltration in the affected dogs to be composed mostly of T cells rather than B cells. This result correlated with the previously assigned chronic inflammation scores based on the histopathological evaluation (Supplementary Table S2 and Fig. S2) and validated our RNA-seq data.

Discussion
Dogs with XLHN have been studied as both an example of progressive canine glomerular disease and an animal model of human AS, which has CKD as a major syndrome component 2 . The genetic cause is well characterized; however, the gene expression and molecular pathways influencing disease progression are incompletely known. In particular, the variable rate of disease progression in dogs with the same mutation and within families affected by AS is intriguing 7 . We, therefore, aimed to evaluate differential gene expression, overrepresented pathways, and upstream regulators by comparing RNA-seq data in dogs that displayed a rapid clinical progression of the disease to those with relatively slow disease progression.
In this study, we examined serial biopsies from rapid and slow groups as well as healthy age-matched littermates. To understand the biological changes during the pathogenesis of CKD, we included the earliest time point at which clinical disease could be detected in these dogs (onset of proteinuria, T1), with the onset of azotemia (T2) and the advent of end-stage disease (T3). By performing renal biopsies when animals reached specific clinical markers of disease progression, we could compare the same clinical stage in the rapid and slow groups. We believe this type of approach provides more confidence in identifying DEGs that are involved in the rate of disease progression than the traditional method of using age-driven time points, as differences detected are likely to be the driving force rather than the consequence of disease progression. Our data demonstrate the dynamic changes  in gene expression at different stages of the disease. It also supports that the biological processes and pathways of fibrosis/adhesion and inflammation are the likely driving forces producing different rates of disease progression in the rapid and slow groups.
While it is known that serum creatinine correlates strongly with tubulointerstitial fibrosis 29 , one of the most intriguing findings when comparing the rapid and slow groups was the increase in fibrosis observed in the rapid group at the same clinical stage of disease (T2 and T3; Fig. 1). This corresponds with many of the 70 DEGs identified using RNA-seq between the rapid and slow groups that are implicated in fibrosis, almost all of which were identified at T2. Several of these genes have been previously described as upregulated in XLHN dogs, Alport  mice, and other kidney diseases 9,11,30,31 . Among these, one of the upregulated genes, CD248 (endosialin or tumor endothelial marker 1, TEM1), has been found to mediate the adhesion and migration of cells through ligand interaction with the upregulated collagen-related genes: collagen type I, collagen type IV, and fibronectin-1 (FN1) 32 . CD248+ stromal cells bind extracellular matrix and have been implicated in kidney 33 , liver 34 , and lung fibrosis 35 . In non-inflamed kidneys, CD248 is expressed by mesangial cells located in glomeruli. In fibrotic kidneys, CD248 is additionally expressed by myofibroblasts and stromal fibroblasts, and the increased expression is closely related to prognostic indicators, such as albuminuria and renal scarring 33 . Of note, one of the downregulated genes in the rapid group, prolactin receptor (PRLR), decreased in association with the extent of interstitial collagen I deposition in kidney transplant rejection 36 , suggesting that PRLR might be a protectant against renal fibrosis. In our study, the decrease in PRLR could be responsible for the more rapid development of fibrosis and consequential faster progression of disease in the rapid group.
Involvement of inflammatory components in the progression of CKD is another major finding of this study. Several inflammatory genes involved in fibrotic changes, such as biglycan (BGN), kielin/chordin-like protein (KCP), and matrix metallopeptidase-2 (MMP2) were upregulated in the rapid versus slow groups. BGN plays a role in bone growth, muscle development and regeneration, and collagen fibril assembly in multiple tissues. BGN is upregulated in renal fibrosis 37 , and BGN protein expression strongly correlated with chronic kidney progression in one study 30 , which may suggest its role in regulating inflammation and innate immunity. KCP expression is stimulated by renal stress, and it enhances the antifibrotic function of BMP7 to attenuate the profibrotic stimulus of TGF-β and to suppress proinflammatory cytokines 38 . In Alport mice, the administration of recombinant BMP7 reduces glomerular and interstitial fibrosis but also upregulates MMP2 39 . This upregulation of MMP2 seems contradictory, as it is associated with renal fibrosis in several animal models, including XLHN dogs 11 . However, the function of MMP2 is specific to the temporal context of fibrosis. At the prefibrotic phase, increased MMP2 induces epithelial to mesenchymal transition, tubular atrophy, and fibrosis. For established fibrosis, inducing MMP2 synthesis by BMP7 promotes proteolytic removal of accumulated extracellular matrix, which is thought to be a potential therapeutic strategy 40 . Since both KCP and MMP2 are upregulated at T2, when fibrosis is already relatively well established in the rapid group, their downstream actions are likely skewed toward anti-fibrotic effects.
Mechanisms other than fibrosis and inflammation can also play roles in the rapid progression of CKD. NAT8, which is almost exclusively expressed by tubular cells in the renal cortex 41 , is a cysteine S-conjugate N-acetyltransferase that is responsible for glutathione-mediated detoxification of nephrotoxic substances. The downregulation of NAT8 in the rapid group suggests a more severe loss of normal renal function in the presence of similar serum creatinine concentrations compared with the slow group. Another downregulated gene, organic anion transporter 3 (OAT3, also known as SLC22A8), is decreased in kidney biopsies from human CKD patients and in a nephrectomized rat model of CKD 42,43 . Reduced protein expression of Oat3 is associated with decreased excretion of an endogenous uremic toxin; meanwhile, the accumulation of this uremic toxin further inhibits Oat3-mediated transportation, accelerating toxin accumulation in serum 42,43 . Therefore, downregulation of OAT3 in the rapid group may result in impaired urinary excretion that is not adequately represented by serum creatinine concentration.
Only 2 DEGs, SCD5 and TK1, were differentially expressed in the rapid group as compared with the slow group at T1. SCD5 was downregulated at both T1 and T2 in the rapid compared with the slow group. And, it was downregulated in the rapid versus control comparison throughout all 3 time points. SCD5 is an isoform of stearoyl-CoA desaturase that is responsible for the formation of monounsaturated fatty acids. Although SCD5 has not previously been described in the CKD literature, it has been proposed as a novel regulator of neural cell proliferation and differentiation, likely through β-catenin-independent (non-canonical) Wnt pathways 44 . In fibrotic kidneys, the canonical Wnt pathway induces myofibroblast differentiation, and the non-canonical Wnt pathway leads to cytoskeleton rearrangement, cell adhesion, and cell movement 45 . Thymidine kinase 1 (TK1) is the only DEG that was exclusively upregulated at T1. Thymidine kinase is responsible for producing dTMP that is later incorporated into DNA. The cytoplasmic isoform of TK1 is cell cycle-dependent, as it substantially increases in the S phase of the cell cycle. Given that unregulated proliferation is the hallmark of neoplasia, TK1 is a valuable serum marker for breast cancer, non-Hodgkin's lymphoma, plasmacytoma, and lung cancer 46 . Its upregulation in our study could indicate increased cell proliferation at T1 in the rapid group. However, further investigation of SCD5 and TK1 is needed to determine their roles in CKD progression.
In addition to abovementioned genes, many genes that are rarely described in CKD progression were found differentially expressed in the rapid versus slow groups at T2 (e. g., LOX, PAMR1, CDCA8, C1QTNF6, FNDC1,  CCDC80, MFSD7, FMO2, SI, SLC26A4). The protein product of the upregulated gene lysyl oxidase (LOX) is an extracellular enzyme that is essential for covalent cross-linking of collagen in irreversible extracellular matrix deposition 47 . Despite reports of its upregulation in liver fibrosis 48,49 and cardiomyopathy 50 , the upregulation of LOX in CKD has only been described in one study using a glomerulonephritis mouse model 51 . Simtuzumab, a monoclonal antibody that inhibits one of the LOX family members, lysyl oxidase homologue 2 (LOXL2), has recently been a focus of research as a possible new treatment for lung, liver, and kidney fibrosis due to a similar pathogenesis 52 . Another gene minimally described in nephrology is the downregulated sucrase-isomaltase (SI). SI is an α-glucosidase that commonly appears on the brush border of small intestinal enterocytes and is involved in glucose digestion. SI is also present in small amounts in non-intestinal cells such as blood leukocytes and kidney cells 53 . Little is known about the function and significance of SI in the renal tubule. However, the decreased expression of SI could indicate renal damage.
To characterize DEGs in functional groups, we performed GO terms and pathway analyses by comparing affected dog groups, both individually and collectively, with control dogs at each time point. GO terms showed that the functions of identified DEGs were associated with "biological adhesion, " "immune system processes, " and "immune response, " representing a common mechanism of disease progression in the early stages of CKD 54 .
SCIeNTIfIC REPORTs | 7: 16776 | DOI:10.1038/s41598-017-16603-y Consistent with the GO terms, the "integrin signaling pathway" was the most upregulated pathway, and the "T cell activation pathway" and the "chemokine and cytokine signaling pathway" were also upregulated in multiple comparisons. The exclusive early expression of the "integrin signaling pathway" in the rapid group and universal expression at a later time point could indicate that it is an essential pathway driving rapid progression of disease in these dogs. The "integrin signaling pathway" consists mainly of collagen and integrin genes. The integrin subunit alpha 2 gene (ITGA2) was increased in the rapid versus control group at both T2 and T3. The COL4A3 −/− /ITGA2 −/− double knockout Alport mouse model has delayed renal fibrosis compared with COL4A3 −/− /ITGA2 +/+ Alport mice, which express significantly higher levels of MMP2, MMP9, MMP12, and TIMP1 55 . Upregulation of MMP2 is consistent with our Alport dog model, suggesting that the "integrin signaling pathway" is involved in matrix accumulation.
To verify the pathway analysis results, we used IHC to identify the infiltrating lymphocyte population. The "T cell activation, " "integrin signaling, " and the "inflammation mediated by chemokine and cytokine signaling" pathways are all closely related to the presence of T lymphocytes 26 . Consistent with this, T lymphocytes were identified as the predominant inflammatory cell population present in the affected dogs, as well as in a study of canine end-stage renal disease 56 . Moreover, T cell infiltration inversely correlates with renal function at the time of renal biopsies in AS patients 57 . Overall, IHC validated the results of GO terms and pathway analyses, which showed that inflammatory pathways and corresponding biological processes are altered during the progression of CKD.
Ingenuity Pathway Analysis (IPA) allows for characterizing orthologous genes, and results identify possible mechanisms that have been validated in humans, mice, and rats. The top upregulated pathway identified in multiple comparisons was the "hepatic fibrosis/hepatic stellate cell activation pathway. " The IPA identified enriched pathways based on the over-represented DEGs, and the principle is the same as that used in the GO terms analysis via Panther. There was an extensive overlap between the DEGs we identified and the genes involved in the "hepatic fibrosis/hepatic stellate cell activation pathway" in the IPA, including collagen genes, cytokine-related genes, matrix metalloproteinases, tissue inhibitor of metalloproteinase, and the TNF receptor superfamily. Therefore, the upregulation of this pathway in our study supports common mechanisms involved in hepatic and renal fibrosis. It could also indicate contributing genes beyond those identified using known canine gene pathways.
The upstream regulator analysis of IPA identified the TGF-β group, especially TGF-β1, as the top upstream regulator, with both the highest activated z-score and the lowest p-value across multiple comparisons. Previous studies in canine 8 and murine 58 models of AS demonstrated the expression of TGF-β mRNA in kidney tissue. However, TGF-β was not differentially expressed in our comparisons, and IHC staining for TGF-β did not demonstrate appreciable differences in XLHN dogs compared to controls in a previous study 8 . The IPA Upstream Regulator Analysis predicts the upstream regulator of gene expression changes based on the knowledge of this regulatory cascade in the literature compiled in the Ingenuity Knowledge Base. Thus, the identified upstream regulator may not be identified as a DEG despite its importance in gene regulation.
A limiting factor of this study is that the kidney biopsies were immediately placed into RNALater to preserve RNA integrity, so microscopic evaluation could not be performed to determine whether the biopsy used for RNA isolation was representative of the cortex as a whole. Thus, one of the samples in the control group at T1 could have represented an area affected by a clinically insignificant insult. Another limiting factor is that expression data represents the mean expression by many cell types. Because the kidney has many cell types, all with different roles in CKD progression, it would be ideal to study gene expression changes in individual cells using laser-capture microdissection to further elucidate the progression of CKD. Last, extensive validation of the RNA-seq results was not performed; however, previous studies have shown RNA-seq to be a robust tool that highly correlates with qPCR results [14][15][16][17][18][19] . RNA-seq may even be more reliable than qPCR due to its higher sensitivity and lower probe bias 16 . We did perform IHC to further characterize the inflammatory population, which supported the pathways identified through the RNA-seq results.
In what appears to be the first RNA-seq study of a canine CKD model, we identified several previously described and novel genes and enriched pathways involved in the pathogenesis and development of CKD. The approach of acquiring biopsies at time points determined by the clinical stage of disease was an attempt to target causative gene expression starting at the clinical onset of proteinuria rather than secondary changes. Regardless of initial insult, CKD has common pathways that lead to end-stage kidney disease 20 . Therefore, many genes found in the current study may serve as predictive or diagnostic biomarkers for early detection of CKD in dogs and people. They may also be potential targets for drug development for this condition.

Methods
Animals. The dogs in this study were part of a colony with XLHN maintained at Texas A&M University 5 .
XLHN is caused by a 10-base deletion in the gene encoding the α5 chain of type IV collagen. Affected males develop juvenile-onset CKD that progresses to end-stage renal disease as previously described 5 . Overall, 6 affected dogs and 2 unaffected littermates were studied. All dogs were raised according to standardized protocols, and no treatments were given to these dogs. All protocols were approved by the Texas A&M University Institutional Animal Care and Use Committee.
Clinical phenotypes. For this study, dogs were selected to represent both extremes in the speed of disease progression in this family of dogs (rapid versus slow progression). Clinical progression was determined by serial monitoring of serum and urine biomarkers of kidney disease, which allowed us to establish specific progression time points 29 : T1 (onset of proteinuria: defined as the presence of microalbuminuria for 2 consecutive weeks (E.R.D. HealthScreen Canine Urine Test Strips, Loveland, CO, USA)); T2 (onset of azotemia: serum creatinine ≥ 1.2 mg/dL); and T3 (end-stage disease: serum creatinine ≥ 5 mg/dL). Rapidly-progressing (rapid) dogs (n = 3) reached each time point at an earlier age than slowly-progressing (slow) dogs (n = 3). On average, the rapid group reached T3 at 26.3 weeks of age (range: 26-27 weeks), while the slow group reached the last clinical time point (T3) at 49 weeks of age (range: 46-52 weeks) (Supplementary Table S7).
Tissue collection. Kidney cortex was serially collected from each dog at the aforementioned 3 clinical time points (independent of age). Control dogs (n = 2) were biopsied to correspond with an affected littermate. All samples were collected by ultrasound-guided needle biopsy. This technique was appropriate for the current study as it is unlikely to induce changes that might be confused with those of CKD progression 59 . Samples for pathology evaluation and immunohistochemistry were placed in formalin and embedded in paraffin. Samples for RNA sequencing were immediately placed in RNAlater Stabilization Solution (Life Technologies, Foster City, CA, USA) and stored at −80 °C until RNA isolation.
Histopathological evaluation. Paraffin-embedded samples were processed and stained as previously described 29 . To determine the severity of interstitial fibrosis and chronic inflammation, a board-certified veterinary anatomic pathologist (REC) evaluated 5 or 20 randomly chosen 20x fields of renal cortex based on core size for each biopsy. For interstitial fibrosis, a score of 0 to 3 was assigned for each field based on the degree of tubulointerstitial architecture distortion caused by fibrosis: 0 -no fibrosis, 1 -fibrosis present but no distortion, 2 -moderate distortion, and 3 -severe distortion. For chronic inflammation, 0 -no inflammatory cells, 1 -scattered inflammatory cells, 2 -aggregates of inflammatory cells that separate or replace tubules, and 3 -diffusely distributed inflammatory cells. Statistical analysis comparing the average scores between groups was performed using bootstrap in R (version 3.2.4) to construct simultaneous 95% confidence intervals for all 3 pairwise comparisons of mean fibrosis and chronic inflammation scores at each of the latter 2 disease stages.
RNA isolation and sequencing. The MirVana miRNA Isolation Kit (Ambion, Austin, TX, USA) was used to isolate total RNA from homogenized kidney tissue according to the manufacturer's instructions. The library preparation, sequencing, and initial quality check were performed by the Texas A&M AgriLife Genomics and Bioinformatics Service (http://www.txgen.tamu.edu/). RNA integrity was assessed by the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). The average RIN was 3.4, and the average RNA yield was 86.2 ng/μL (Supplementary Table S3). To compensate for the low-input samples, we use the TruSeq Stranded Total RNA Library Prep Kit with Ribo-zero Gold (Illumina, San Diego, CA, USA), based on the best practice for RNA with variable qualities 21 (to remove both cytoplasmic and mitochondrial rRNA) and its compatibility with canine samples. Samples were then sequenced using the Illumina Genome Analyzer (HiSeq. 2500v4 High Output). Raw sequencing data were submitted to the NCBI sequence read archive (SRA) (Accession: SRP101707; Samples: SAMN06560417, SAMN06560429-51; BioProject: PRJNA378728). Data analysis. FastQC (version 0.11.2) was used for quality control to ensure that the quality value was above Q30. The canine (Canis lupus familiaris) genome FASTA file (ftp://ftp.ensembl.org/pub/current_fasta/ canis_familiaris/dna/Canis_familiaris.CanFam3.1.dna.toplevel.fa.gz) and gene annotation GTF file (CanFam 3.1 assembly; ftp://ftp.ensembl.org/pub/current_gtf/canis_familiaris/) were obtained from Ensembl. Although RNA-seq is a popular research tool, there is no gold standard for analyzing RNA-seq data. Among the available tools, we chose up-to-date open source tools for mapping, retrieving read counts, and differential analysis. We used HISAT2 22 (version 2.0.3-beta) to generate indexes and to map reads to the canine genome. For assembly, we chose SAMtools (version 1.2) and the "union" mode of HTSeq 23 (version 0.6.1), as the gene-level read counts could provide more flexibility in the differential expression analysis. Both HISAT2 and HTSeq analyses were conducted using the high performance research computing resources provided by Texas A&M University (http:// hprc.tamu.edu) in the Linux operating system (version 2.6.32). Differential expression and statistical analysis were performed using DESeq2 (release 3.3) in R (version 3.2.4). DESeq2 24 was chosen as it is a popular parametric tool that provides a descriptive and continually updated user manual (https://bioconductor.org/packages/release/ bioc/vignettes/DESeq2/inst/doc/DESeq2.html). DESeq2 internally corrects for library size, so it is important to provide un-normalized raw read counts as input. We used variance stabilizing transformation to account for differences in sequencing depth. P-values were adjusted for multiple testing using the Benjamini-Hochberg procedure 60 . A false discovery rate adjusted p-value (i.e., q-value) < 0.05 was set for the selection of DE genes.
Gene ontology (GO), pathway, and upstream regulator analysis. Gene ontology and PANTHER pathway analyses were performed with the PANTHER Overrepresentation Test (released on July 15, 2016) in PANTHER 25 version 11.1 (http://www.pantherdb.org/, released on October 24, 2016). This program supports the canine genome. Also, QIAGEN's Ingenuity Pathway Analysis (IPA, QIAGEN Redwood City, www.qiagen. com/ingenuity) was used to provide overrepresented orthologous genes in human, mouse, and rat databases and to identify orthologous pathways and upstream regulators in our data. PANTHER used the binomial test and Bonferroni correction for multiple testing, while IPA used the right-tailed Fisher Exact test and displayed z-scores to indicate whether a potential regulator was activated or inhibited. We used the default settings for statistical analysis in both the PANTHER pathway and IPA. In PANTHER, only pathways and GO terms with fold enrichment > 0.2 were listed. In IPA, p-value < 0.05 and fold change > 2 were set as cutoff values. Immunohistochemistry (IHC). Three-micrometer, formalin-fixed, paraffin-embedded renal cortex sections of both affected and control dogs were stained for CD3 (n = 24) and CD20 (n = 9). After deparaffinization, the sections were placed in citrate buffer (pH 6.0) for antigen retrieval, using a pressure cooker (Decloaking Chamber, Biocare Medical). Endogenous peroxidase activity and non-specific protein binding were blocked with 3% hydrogen peroxide and Sniper protein block (Biocare Medical), respectively. After blocking, the sections were incubated with primary antibodies CD3 (1:300 dilution; Dakocytomation, Carpinteria, CA) and CD20 (1:500