Identification of novel targets of diabetic nephropathy and PEDF peptide treatment using RNA-seq

P78 drug delivery in diabetic Ins2 ^Akita mice

Kidney tissue samples were obtained from a recently published study that tested physiological effects of P78 on kidney function and pathology [25]. Animals were treated as previously described [25] prior to extracting RNA for RNA-seq analysis. The animal studies were approved by the Penn State University College of Medicine Institutional Animal Care and Use Committee, and performed in strict accordance with the recommendations in the Guide for the Care and Use of Laboratory Animals of the National Institutes of Health. All experiments were conducted using male D2.B6-Ins2 ^Akita /MatbJ diabetic animals and their wild type (WT) littermate mice (DBA/2 J background), recommended by the Animal Models of Diabetes Complications Consortium (AMDCC) as a model of DN [26, 27]. The diabetic D2.B6-Ins2 ^Akita mice develop hyperglycemia at 3 weeks of age and all treatment carried out when the mice were either 6 weeks (3 weeks hyperglycemic exposure; early stage treatment) or 12 weeks (9 weeks hyperglycemic exposure; late stage treatment) of age. Only mice with blood glucose levels > 350 mg/dl (measured using Accu-Chek glucometer, Boehringer Mannheim, Indianapolis, IN) were considered diabetic and used in the study.

The drug tested was P78, a small PEDF active peptide [22, 23], generated by methods previously described [25, 28]. Briefly, P78 peptide at a dose of 0.3 μg/g/day or vehicle (phosphate-buffered saline; PBS) was administered by continuous subcutaneous infusion via the osmotic minipump (no. 2006; Alzet, Durect, Palo Alto, CA), implanted dorsally between the shoulders of the animals as previously described [25, 29–31]. Transcriptome analysis of wild-type and diabetic kidney samples were performed at two stages of diabetes where treatment was initiated at an early stage (6 weeks of age; 3 weeks hyperglycemic) and late stage (12 weeks of age; 9 weeks hyperglycemic). Age, gender, and weight matched diabetic Ins2Akita and wild-type non-diabetic controls were used in the study. All animals including wild-type were implanted with an osmotic minipump infused with either vehicle (wt and diabetic controls) or the P78 peptide (diabetic mice). Duration of treatment was 6 weeks with either peptide or vehicle. One group received treatment at the early stage of diabetes (ET, early treatment) at 6 weeks of age and the experiment terminated at 12 weeks of age. Treatment in the second group was initiated at late stage diabetes (LT, late treatment) at 12 weeks of age and terminated at 18 weeks of age. Mice were provided ad lib access to food and water and were euthanized at the end of the experimental period. Kidney samples for RNA extraction were immediately harvested and frozen in liquid nitrogen at the termination of the experiment.

Tissue samples preparation and RNA isolation

For RNAseq we used 13 kidney tissue samples from wild-type mice, 7 from the diabetic Ins2 ^Akita mice, 8 from early P78 treatment of diabetic Ins2 ^Akita mice, and 7 from late P78 treatment the diabetic Ins2 ^Akita mice [25]. Total RNA was extracted using mirVana kit (Life Technologies) with some modifications. Briefly, a bead mill homogenizer (Bullet Blender, Next Advance) was used to homogenize the tissue using a safe-lock microcentrifuge tube (Eppendorf) and a mass of stainless steel beads (Next Advance, cat# SSB14B) equal to the mass of the tissue. Immediately after two volumes of lysis buffer were added to the tube, samples were mixed in the Bullet Blender for 1 min at a speed of six. Samples were visually inspected to confirm desired homogenization and then incubated at 37 °C for 5 min. The lysis buffer was added up to 0.6 ml, and samples were mixed in the Bullet Blender for 1 min. Optical density values of extracted RNA were measured using NanoDrop (Thermo Scientific) to confirm an A₂₆₀:A₂₈₀ ratio above 1.9. RNA integration number (RIN) was measured using BioAnalyzer (Agilent) RNA 6000 Nano Kit to confirm RIN above 7.

Library preparation and sequencing

The cDNA libraries were prepared using SureSelect Strand Specific RNA Library Preparation Kit (Agilent) as per the manufacturer’s instructions. Briefly, polyA RNA was purified from 1000 ng of total RNA using oligo (dT) beads. Extracted RNA was subjected to fragmentation, reverse transcription, end repair, 3’-end adenylation, adaptor ligation and subsequent PCR amplification and SPRI bead purification (Beckman Coulter). The unique barcode sequences were incorporated in the adaptors for multiplexed high-throughput sequencing. The final product was assessed for its size distribution and concentration using BioAnalyzer High Sensitivity DNA Kit (Agilent) and Kapa Library Quantification Kit (Kapa Biosystems). 12 libraries were pooled and diluted to 2 nM in EB buffer (Qiagen) and then denatured using the Illumina protocol. The denatured libraries were diluted to 10 pM by pre-chilled hybridization buffer and loaded onto TruSeq SR v3 flow cells on an Illumina HiSeq 2500 and run for 50 cycles using a single-read recipe (TruSeq SBS Kit v3) according to the manufacturer's instructions (Illumina).

Quality control, mapping and quantification of RNA-Seq reads

Illumina CASAVA pipeline Version 1.8 was used to extract de-multiplexed sequencing reads. FastQC (version 0.11.2) (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) was used to validate the quality of the raw sequence data. Additional quality filtering used FASTX-Toolkit (http://hannonlab.cshl.edu/fastx_toolkit) using a quality score cutoff of 20. Next, alignment of the filtered reads to the mouse reference genome (mm10) was done using Tophat (version 2.0.9) [32] allowing 2 mismatches. Picard (version1.102) (https://github.com/broadinstitute/picard) was used to assess proportion of mapped bases to coding, UTR, intronic, and intergenic regions, respectively. Picard was used to find coverage across gene body to determine 5’- or 3’- bias. FPKM (Fragments Per Kilobase of Exon Per Million Fragments Mapped) values were calculated using Cufflinks Version 2.0.2 [33] as provided with the Ensembl gene annotation (release 78).

Differential Gene Expression (DEG) analysis

Only reliably expressed genes were included in the analysis, defined to be those with at least 2 samples with FPKM > = 1. The ComBat function of the sva v3.10.0 R package was used for batch normalization of the FPKM values. The DEGexp function of the DEGseq v1.18.0 R package [34] was used to identify differentially expressed genes (DEG) between diabetes and control (DvC), treatment and diabetes (TvD), early treatment and diabetes (ETvD) and late treatment and diabetes (LTvD), using the Likelihood Ratio Test method. Significantly DEG were defined to be those with q-value < 0.05 calculated by the Storey et al. 2003 method.

Visualization

The batch adjusted RNA-seq FPKM values of the reliably expressed genes were averaged for samples in each category: C, D, ET and LT. A heatmap of the Z scaled resulting FPKM values was generated using the heatmap.2 function of the gplots R package, with parameters “average” for clustering method and 1-correlation for distance. The Z scaling was performed with the genefilter v1.46.1 R package. A volcano plot was generated of the –log10(p-value) vs. log2(normalized FC) of the DvC and TvD DEG analyses, with points with q < 0.05 and abs(log2(normalized FC)) > = 1 (i.e. a FC threshold of 2) colored red for DvC and green for TvD. Labeled genes have a q-value < 0.05 and abs(log2(normalized FC)) > = 2 (i.e. a FC threshold of 4). Quad Venn diagrams showing the number of significant DEG genes and protein coding subsets for various comparisons were generated using the VennDiagram v1.6.9 R package.

Functional analysis

The Tox Analysis function of Ingenuity Pathway Analysis (IPA) (QIAGEN, California) was used to create gene sets based on biological processes and toxicological responses to xenobiotic insult using the mammal filter (human, rat, mouse) and the following tissues in turn: adipose, brain, heart, kidney, liver, retina and all tissues. IPA Build 2015-03-23.

Statistical analysis

Statistical tools embedded in the previous instruments and statistical packages were applied. Significant DEG were defined to be those with q-value <0.05 calculated by the Storey et al. 2003 method.

RT-PCR

Global transcriptome changes

Results obtained from DEGseq analysis applied on RNA-seq FPKM (Fragments per kilobase of exon per million fragments mapped) showed significant differences in coding and non-coding genes expressed in the kidney among control (C, wild-type), diabetic (D), early (ET), and late (LT) P78 treatment in the Ins2Akita mouse model of diabetic nephropathy (DN). Below we provide results for protein coding and non-coding genes that were regulated in the diabetic kidney relative to P78 treated and wt controls: The comparisons are (1) diabetic relative to wild type (DvC), (2) P78 treated relative to non-treated diabetic (TvD), (3) diabetic returned to normal levels by treatment (early and late stage treatment), and (4) P78 early stage treatment (ET) compared to late stage treatment (LT). Early (ET) and late (LT) P78 treatments are grouped in the TvD analyses to determine all gene targets of the treatment whether given at early or late stages of diabetes. These are later separated in the study to determine those targets that were unique to ET and those unique to LT.

Global expression changes

The heatmap in Fig. 1 and volcano plot in Fig. 2 show global transcriptome changes in control, diabetic, and P78 early and late treatments in kidney tissues of the Ins2 ^Akita mouse model of diabetic nephropathy. 43,168 RNA transcripts were identified in the kidney. 25,254 had an FPKM value > 0 and 14,316 were considered reliably expressed, defined to be those with FPKM > =1. Of the 14,316 reliably expressed sequences, 12,328 were protein coding genes and 1,988 non-coding RNA transcripts.

The distribution of all genes (protein-coding and non-coding) with respect to significance (y axis) versus fold changes (x axis) is shown in the scatter plot in Fig. 2 (–log10(p-value) vs log2(normalized fold change (FC)). Genes of interest are outside the midline where absolute normalized FC > =(+/-)4 are colored red for DvC (diabetic versus wt controls) and green for TvD (early and late stage P78-treated diabetic mice versus diabetic controls). A few genes including Nhej1, Ept1, Cyp4a14, and Ugt1a2, which are significantly regulated by diabetes and treatment with an absolute normalized FC > =(+/-)4, are labeled as examples to show their location and relationship to other points in the plot and to the midline. Many protein-coding and non-coding genes showed small changes in their expression levels in either direction but only those that met our selection criteria of FC > =(+/-)1.5, q < 0.05 were of interest in this study and rank-ordered in all tables according to fold changes in either direction.

All transcripts regulated by diabetes and P78

Gene expression differences in the kidney between the diabetic compared to control (DvC) and P78 treatment compared to diabetic (TvD, early and late) groups are shown in the Venn diagram in Fig. 3 and in Table 1. These represent all protein-coding and non-coding RNA transcripts that showed statistically significant differential expression (DEG) (q < 0.05). Changes in expression levels were seen in 1237 transcripts in the kidney (a total of all numbers in the Venn diagram). 208 (46 + 51 + 87 + 24) of these were regulated by both diabetes and treatment. 1,129 RNA transcripts (coding plus non-coding) were regulated by diabetes with log2 fold changes (FC) ranging between -86.2 and +86.0 (q < 0.05) (Table 1). 164 were upregulated and 965 downregulated by diabetes. P78 treatment (early and late; TvD group) altered expression of 316 transcripts with FC between -93.3 and +50.9 (q < 0.05). Of these, 183 were upregulated and 133 downregulated. A total of 138/1129 (12.2%) transcripts whose levels were either increased (51 transcripts) or decreased (87 transcripts) by diabetes were reversed in the direction of their expression to normal levels by P78 treatment (Table 1). While the regulation of most transcripts by diabetes were not reversed by treatment, the effects of P78 in reducing pathology and progression of DN that we have previously shown [24, 25] may be due, in part, to reversing diabetes-induced expression changes in this subset (138) of coding and non-coding RNA transcripts.

Group	# of Genes
Upregulated in D	164
Downregulated in D	965
Upregulated with T	183
Downregulated with T	133
Upregulated in D, reversed by T	51
Downregulated in D, reversed by T	87

Protein-coding genes (mRNA) regulated by diabetes and P78

Figure 4 and Table 2 show the number of protein coding genes (mRNA) regulated by diabetes and P78. Changes in gene expression occurred in 934 protein-coding genes (all numbers in the Venn diagram) in the kidney of the Ins2 ^Akita mouse model of DN, with 88 of these (33 + 15 + 34 + 6) regulated by both diabetes and treatment. Diabetes altered expression of 901 genes (q < 0.05), of which 71 were upregulated and 830 downregulated by the disease (Table 2). P78 (early and late treatment) modulated expression of 121 coding genes with 65 of these showing increased and 56 showing decreased expression levels. P78 treatment reversed direction of expression of 49/901 (5.44%) of the protein coding genes regulated in diabetic kidneys towards normal levels (Table 2).

Group	# of Genes
Upregulated in D	71
Downregulated in D	830
Upregulated with T	65
Downregulated with T	56
Upregulated in D, reversed by T	15
Downregulated in D, reversed by T	34

Table 3B lists genes (FC > =(+/-)1.5 in diabetic animals that showed changes in expression by P78 treatment (early and late) of which the top 5 include Nhej1 (-37.94), Tceanc2 (+5.76), Ept1 (-4.45), Ugt1a2 (+3.03), and Tmsb15l (-3.0) with FC > (+/-)3.0 (q < 0.05).

P78 Treatment normalized gene expression altered by diabetes

From this analysis (Table 3A and 3B), we identified a small set of genes where the changes in diabetes were substantially (~50%) or completely reversed by P78 treatment. These are listed in Table 3C and rank-ordered by FC > =(+/-)1.5 (q < 0.05). Of these, non-homologous end-joining factor 1 (Nhej1) had the highest absolute fold change in expression with both diabetes and treatment. Diabetes induced expression levels of Nhej1 by +32.04 while P78 complete reversed expression by -37.94 fold changes. Treatment also completely reversed diabetes-induced expression changes of Mamdc4, Kdm4b, Tmem252, Selm, and Hpd while others were returned in the direction of normal levels from their diabetes-induced state by ~50% (Table 3C).

Regulation of genes with FC > =(+/-)2.0 (q < 0.05) by diabetes, or treatment (P78 early (ET), and P78 late (LT)) (Table 3C) was subsequently cross-validated for trend in expression by RT-PCR and QRT-PCR (Fig. 5) with the exception of Mamdc4, whose expression level was below the detection threshold using these PCR conditions. While regulation trends remain the same, some differences in levels of gene regulation between the RNAseq and PCR data were noted, which were possible due to primer specificity, alternatively spliced forms of the genes, RNA integrity, and/or PCR conditions.

This small subset of genes regulated in opposite directions by both diabetes and P78, a treatment previously shown to reduce DN pathology [24, 25] may be critically involved in diabetes-induced mechanisms that control the development of the disease.

Transcripts regulated by early and late stage P78 treatment

Non-coding small RNAs (miRNA)

Biological function - pathways Ingenuity Pathway Analysis (IPA)

IPA analysis shows that protein-coding genes regulated by diabetes and returned to near normal levels by treatment (Table 3C) are associated with key biological functions including lipid metabolism, post translational modification, endocrine and hematological functions, cell death and survival, and protein synthesis, all important processes associated with the disease (Fig. 6a). Many of these are also associated with other pathologies including organismal injury and abnormalities, cancer, developmental and hereditary disorders (Fig. 6b)

IPA gene ontology algorithms and KnowledgeBase mining identified several important canonical pathways enriched for the differentially expressed genes (DEG) in our dataset. DEG fell into one of several functionally relevant canonical pathways associated with diabetes and kidney pathology. Figure 7 shows the top pathways represented with FC > =(+/-)1.5 (q > 0.05) and Fig. 8, the percentage of genes regulated relative to those present in the IPA knowledgebase that are assigned to a given pathway and the direction of their regulation (red: upregulated; Green: downregulated) for diabetic (DvC) and P78 treatment (TvD) groups (relative to wt controls). The total number of genes assigned to a pathway in the IPA KnowledgeBase is given to the right of each bar. The open area of each bar represents genes in the pathway that are not represented in our dataset. Pathways were rank-ordered by IPA base on the significance of the directional change (–log(p-value) and are limited to those enriched in IPA KnowlegeBase.

The graph in Fig. 8a shows a striking trend in diabetes-induced DEGs with mitochondrial dysfunction in the diabetic group relative to wild type controls. The largest functional cluster of genes modulated by diabetes was represented in the mitochodrial dysfunction canonical pathway. Approximately 50% of the 123 genes assigned to this pathway were differentially expressed by diabetes and ~48% of these were downregulated. The overall trend observed from this analysis was that diabetes decreased expression of more kidney genes clustered in these pathways (12/16) when compared to controls. In three of these pathways associated with chronic renal injury, renal failure, and renal inorganic phosphate homeostasis, all genes were downregulated by diabetes that were present in our kidney dataset.

P78 targeted some of the same pathways as diabetes but fewer genes showed altered expression levels after treatment relative to controls (Fig. 8b). For example, ~50% of the genes in the mitochondrial pathway showed changes in expression in the diabetes group relative to controls while only ~5% showed altered expression in the treated diabetes group relative to controls. Genes in the P78 treated diabetic mice represented in the cytochrome P-450 and reversible glomerulonephritis associated pathways were all downregulated compared to those in the same pathway for non-treated mice. While those assigned to the recovery from ischemic acute renal failure pathway showed increased levels of expression after treatment compared to the non-treated diabetic mice. Since we previously showed that P78 significantly reduced physiological and pathological changes in diabetic nephropathy [24, 25] these target pathways and their associated genes may give mechanistic clues to the disease.