Novel insights into the disease transcriptome of human diabetic glomeruli and tubulointerstitium.

BACKGROUND: Diabetic nephropathy (DN) is the most common cause of end-stage renal disease, affecting ∼30% of the rapidly growing diabetic population, and strongly associated with cardiovascular risk. Despite this, the molecular mechanisms of disease remain unknown.
METHODS: RNA sequencing (RNAseq) was performed on paired, micro-dissected glomerular and tubulointerstitial tissue from patients diagnosed with DN [n = 19, 15 males, median (range) age: 61 (30-85) years, chronic kidney disease stages 1-4] and living kidney donors [n = 20, 12 males, median (range) age: 56 (30-70) years].
RESULTS: Principal component analysis showed a clear separation between glomeruli and tubulointerstitium transcriptomes. Differential expression analysis identified 1550 and 4530 differentially expressed genes, respectively (adjusted P < 0.01). Gene ontology and Kyoto Encyclopedia of Genes and Genomes enrichment analyses highlighted activation of inflammation and extracellular matrix (ECM) organization pathways in glomeruli, and immune and apoptosis pathways in tubulointerstitium of DN patients. Specific gene modules were associated with renal function in weighted gene co-expression network analysis. Increased messengerRNA (mRNA) expression of renal damage markers lipocalin 2 (LCN) and hepatitis A virus cellular receptor1 (HAVCR1) in the tubulointerstitial fraction was observed alongside higher urinary concentrations of the corresponding proteins neutrophil gelatinase-associated lipocalin (NGAL) and kidney injury molecule-1 (KIM-1) in DN patients.
CONCLUSIONS: Here we present the first RNAseq experiment performed on paired glomerular and tubulointerstitial samples from DN patients. We show that prominent disease-specific changes occur in both compartments, including relevant cellular processes such as reorganization of ECM and inflammation (glomeruli) as well as apoptosis (tubulointerstitium). The results emphasize the potential of utilizing high-throughput transcriptomics to decipher disease pathways and treatment targets in this high-risk patient population.

KEY LEARNING POINTS

What is already known about this subject?

Diabetic nephropathy (DN) is the most common cause of end-stage renal disease, with a fast disease progression and limited treatment options.

The molecular causes of DN are unknown.

RNA sequencing (RNAseq) profiling provides insight into the active disease transcriptome.

What this study adds?

We provide high-quality data from one of the largest numbers of micro-dissected human kidneys with histologically verified diagnosis of DN ever reported.

The study elucidates the complexity of DN pathogenesis by reporting large differences in the gene expression profiles of both glomerular and tubulointerstitial fractions.

Together, results highlight disease- and compartmental-specific changes in the DN kidney.

What impact this may have on practice or policy?

Our detailed molecular characterization of the human DN kidney provides an important foundation for future translational research.

Since our RNAseq data are drawn also from early DN stages, the study illuminates early molecular changes that are of interest for biomarker discovery and better tools for early diagnosis and follow-up.

This opens up for novel targets and mechanisms-targeted treatment strategies of DN.

INTRODUCTION

Diabetes mellitus (DM) is a growing global health problem [1] and a major cause of chronic kidney disease (CKD) [2], with one-third of diabetic patients developing diabetic nephropathy (DN) [3]. DN is the main risk factor for end-stage renal disease (ESRD) [4] and significantly increases the risk of cardiovascular morbidity and mortality [5-7].

While there is no cure for DN, available treatments aim to alleviate symptoms or disease progression and cardiovascular complications. However, renoprotective therapies, including angiotensin-converting enzyme inhibitors (ACEis) and angiotensin II receptor blockers (ARBs) [8, 9], or sodium–glucose cotransporter-2 inhibition [10], have so far been insufficient to reduce ESRD and/or mortality rates.

To date, the diagnosis of DN is based on clinical features, occasionally in combination with a kidney biopsy, typically demonstrating thickening of glomerular basement membrane, mesangial expansion and hyalinosis of afferent and efferent arterioles [11]. However, gaining better insight into the molecular pathophysiology will be a prerequisite for future design of diagnostic and prognostic tools, and new disease intervention targets, which will allow clinicians to distinguish patients at high versus low risk for rapid progression, enabling individualized treatment and guiding selection of patients for clinical trials.

RNA sequencing (RNAseq) has been shown to be a promising avenue to increase molecular understanding of various diseases [12]. In DN research, however, most transcriptomic studies reported to date have been conducted on animal models or human cell lines and show inconsistent results [13-16], partly explained by limited accessibility to human kidney tissue material [17], difficulties in generating adequate animal models or use of biased microarray techniques. Thus, the primary aim of the present study was to, for the first time, provide high-throughput RNAseq data on compartment-specific disease-associated transcriptomic profiles in the human DN kidney. This was accomplished by sequencing paired glomerular and tubulointerstitial fractions from patients with biopsy-verified DN and living donors (LDs). The findings in the present study not only provide novel molecular insight but also reveal potential targets for future clinical and pharmacological intervention.

MATERIALS AND METHODS

Kidney biopsies

Nineteen patients were enrolled from two sites: Sahlgrenska University, Hospital, Gothenburg (SU) and Karolinska University Hospital, Stockholm (K), Sweden. All patients underwent percutaneous kidney biopsies, which were evaluated by renal pathologists. Inclusion criteria were clinically and histopathologically verified DN and with age >18 years. Additionally, 20 living kidney donors from K were included as controls, and the biopsy tissues were obtained immediately post-donation to minimize ischaemic damage. Surplus tissues from diagnosis were immersed in RNAlater (Qiagen, The Netherlands) and stored at −80°C (Table 1;Supplementary data, Methods).

Table 1

DN patient and LD characteristics

Variable	DN patients (n = 19)	LD (n = 20)	P-value
Age, years	61 (30–85)	56 (30–70)	0.527
Sex (male/female)	15/4	12/8	0.301
BMI, kg/m2	28.9 (22.3–33.2)	25.9 (21.1–28.8)	0.079
Systolic blood pressure, mmHg	140 (110–177)	135 (114–158)	0.177
Diastolic blood pressure, mmHg	80 (60–101)	75 (56–89)	0.226
Diabetes duration, years	15 (0–30)a	NA	NA
Diabetes type (DM1/DM2)	5/14	NA	NA
HbA1c, mmol/mol	64 (40–96)b	NM	NA
S-Creatinine, µmol/L	142 (64–233)	70 (53–102)	<0.001
GFR, mL/minc	40 (21–104)d	90 (81–118)	<0.001
eGFR, mL/min (CKD-EPI)	40 (28–114)	95 (78–115)	<0.001
P-Albumin, g/L	30 (21–37)	39 (34–51)e	<0.001
UACR, mg/mmol	212 (33.4–785)f	0.68 (0.37–14.1)	<0.001
P-Sodium, mmol/L	141 (134–147)	141 (140–145)	0.691
P-Potassium, mmol/L	4.4 (2.8–5.6)	4.0 (3.3–4.8)	0.061
P-Calcium, mmol/L	2.23 (2.00–2.44)	2.3 (2.05–2.46)e	0.231
P-Phosphate (mmol/L)	1.2 (0.8–1.5)f	0.9 (0.4–1.2)e	0.001
B-Haemoglobin, g/L	123 (91–152)	147 (125–165)	<0.001
P-hsCRP, mg/L	2.0 (0.4–59.9)	0.7 (0.2–4.9)	0.011

All numerical values are presented as median range. Statistical comparison between patients and controls was performed with Wilcoxon rank-sum test.

Diabetes duration is shown for 17/19 patients.

HbA1c at time of biopsy is presented for 15/19 patients.

GFR is presented as either Cr-EDTA or iohexol clearance.

GFR is presented for 14/19 patients.

Data presented in 19/20 LDs.

Data presented for 18/19 patients.

BMI: body mass index; B: blood; CKD-EPI: Chronic Kidney Disease Epidemiology Collaboration; DM1/DM2: diabetes mellitus Type 1/Type 2; hsCRP: high-sensitivity C-reactive protein; NA: not applicable; NM: not measured; p: plasma; S: serum.

RNA preparation and full-length complementary DNA generation

The tissues were thawed and micro-dissected into the glomerular fraction (5–10 glomeruli per sample) and the remainder referred to as tubulointerstitial fraction. RNA preparation and full-length complementary DNA (cDNA) generation are described in Supplementary data, Methods.

Library generation and sequencing

Next-generation sequencing libraries were made using the NexteraXT transposon-based library generation kit (Illumina Inc.). The DNA libraries were amplified (12 cycles) and purified using AmpureXP beads. Paired-end HiSeq2500 flow-cells were used for sequencing on a HiSeq2500 instrument. For detailed description, see Supplementary data, Methods.

Sequencing data processing and quality control

Preprocessing of the sequencing data was done using BCBIO v1.0.1 (https://github.com/chapmanb/bcbio-nextgen). Hisat2 v2.0.5 [18] was used to align the data against the human hg38 assembly. MultiQC v0.9 [19] was used to summarize pre- and post-alignment quality control statistics. Sailfish 0.10.1 [20] was used to quantify gene-level read counts. For subsequent analyses, counts of protein-coding genes with more than two counts in ≥80% of the samples were used.

Multivariate analysis

The data were variance-stabilizing transformation (VST) normalized using DESeq2 R package [21]. Principal components analysis (PCA) on the top 1000 most variable genes selected by median absolute deviation was performed using MixOmics R package [22]. As for the subsequent analyses, we handled the research centre adjustment, either by applying limma batch correction [23] or by using the design formula ‘.∼Center+Group’.

Differential expression analysis

Differential expression analysis (DEA) was performed with DESeq2 [21] accounting for the research centre factor. Shrunken fold changes were obtained to minimize the impact of the genes with low counts. P-values were adjusted for multiple testing using Benjamini–Hochberg method, and genes with adjusted P <0.01 (transcriptome-wide analyses) and <0.05 (targeted analyses) were considered differentially expressed. Hexbin [24] was used for hexagon binning of overlapping points in the volcano plots. Eulerr [25] was used to visualize the gene lists overlap.

Enrichment analyses

DAVID v6.8 [26, 27], Enrichr [28, 29] and GOplot [30] were used for functional annotation of the differentially expressed genes and visualization of the enriched terms.

Weighted gene co-expression network analysis

Glomerular VST-transformed and limma batch-adjusted [23] counts were used, excluding 10% of low-variance genes. Weighted gene co-expression network analysis (WGCNA) [31] was used to construct an unsigned bicor gene co-expression network and detect gene modules and their association to clinical variables. The soft thresholding power was determined based on the criterion of approximate scale-free topology and set to 10. Minimal module size was 30, and module merging cut-off 0.2. Cytoscape [32] was used for visualization.

Validation with quantitative PCR

Unconsumed RNA from glomerular and tubulointerstitial samples from five patients and five controls were used for validation analyses with quantitative PCR (qPCR). In brief, the Smart-seq2 and KAPA HiFi HotStart (Roche) protocols were used for cDNA synthesis and amplification. COL1A1, CCL13, SFRP2, C3, Nephrin (NPHS1) and Podocin (NPHS2) were analysed in glomerular samples, whereas COL1A1, CCL13, SRFP2, C3, TCL1A, CD3D, FOS and EGR1 were analysed in the tubulointerstitial samples using SYBR Green (BioRad) protocol and primers (Supplementary data, Table S1). The data were analysed by the 2−ΔΔCT method.

Urinary analyses

For urinary analyses, see Supplementary data, Methods.

Statistical analyses of clinical data

Clinical variables were expressed as median (range), and comparisons between two groups were assessed with Wilcoxon rank-sum test or Chi-square test using JMP v11.2 and SPSS Statistics 24. Statistical significance was set at the level of P < 0.05.

Ethics statement

All patients gave informed consent before participating and the procedures were in accordance with the ethical standards of the regional ethical board (Stockholm and Gothenburg) and with the Declaration of Helsinki 1975, as revised in 2013.

Data availability

Included RNAseq data are available at http://karokidney.org/rna-seq-dn/.

RESULTS

Patient disposition

Patients and controls differed with regard to glomerular filtration rate (GFR) (P < 0.001), serum creatinine (P < 0.001) and urine albumin-to-creatinine ratio (UACR) (P < 0.001), but not age and sex (Table 1). All DN cases and 2 LD were on ACEi/ARB; 16 DN patients were treated with other antihypertensive drugs; and 13 DN cases and 2 LD were on statins. Only DN patients used anti-diabetic therapeutics: 16 patients were on insulin; 1 on metformin; 1 used both insulin and metformin; and 1 DN patient used repaglinide. The predominant underlying diabetes diagnosis was DM Type 2 (14 out of 19 patients).

Sequencing quality control

Sequencing data were of overall good quality (Supplementary data, Table S2) with 87% aligned reads on average. Three samples from two individuals were flagged due to >10% rRNA content but were kept for the subsequent analysis since the depth of the remaining sequences was deemed sufficiently high. Normalized expression levels of 359 housekeeping genes [33] were comparable across samples (Supplementary data, Figure S2A). Densities of podocyte- and tubular-specific marker expression (Supplementary data, Figure S2B) visualize separation of the micro-dissected compartments.

Unbiased exploration of the global transcriptome

Unsupervised multivariate analysis

We performed PCA to explore the global trends and reveal the inherent patterns in the transcriptomics data. The PCA scores plot (Figure 1A) shows that the samples are clearly separated by tissue compartment (glomeruli versus tubulointerstitium) along the first principal component 1 (PC1), which explains 64% of variance (Figure 1C). The glomerular samples also clustered according to the research centre, indicating presence of a batch effect. Furthermore, PC2 separated the samples by disease status (DN and LD) (Figure 1A and B).

FIGURE 1

Unbiased analyses of renal transcriptome. (A and B) PCA scores plots: PC1 and PC2 (A), PC2 and PC3 (B). The samples are separated by the tissue compartment (glomeruli and tubulointerstitium) and the disease status (LD and DN). The shape of the plotting characters designates the research centre (GU/SU, University of Gothenburg/Sahlgrenska University Hospital and KI/K, Karolinska Institutet/Karolinska University Hospital). (C) PCA scree plot showing explained variance per PC, clearly illustrating that the first three PC are the most important ones in our dataset. (D and E) The volcano plots visualize log-fold changes in gene expression (x-axis) in relation to the negative logarithm of the adjusted P-value (y-axis) of the DN versus LD kidney transcriptomes comparison. To avoid overplotting, the individual values were grouped by similarity in bins, and the colour gradient reflects the gene count per bin. The blue horizontal dotted line shows the significance threshold at adjusted P = 0.01. The black vertical dotted line indicates zero-fold change. The glomerular compartment was characterized by more genes being down-regulated, while the opposite was true for the tubulointerstitial compartment. (F) Counts of differentially expressed genes (DEGs) at different fold change levels. (G) Venn diagram visualizing the overlap of differentially expressed genes between the renal compartments. Out of all the genes, 727 were differentially expressed in both compartments, indicating that also similar processes occur in the two compartments.

Differential expression analysis

DESeq analysis [21] revealed major differences in DN versus LD kidney transcriptomes as visualized by volcano plots (Figure 1D and E). As a general trend, most differentially expressed DN genes in the tubulointerstitium were up-regulated, while glomerular transcriptome tended to be down-regulated (Figure 1F). In total, we identified 1550 significantly differentially expressed genes in the glomerular fraction (576 up-regulated and 974 down-regulated) and 4530 significantly differentially expressed genes (2308 up-regulated and 2222 down-regulated) in the tubulointerstitium (Supplementary data, Table S3). We observed a substantial overlap in differential expression from both compartments as well as presence of compartment-specific modulated genes (Figure 1G). The majority of the shared genes (672/726) were modulated in a directionality-consistent fashion, while 54 (7.4%) showed opposite patterns (Supplementary data, Table S4) with XKR4, FZD2, DCN, KCNB1, PRRT2, KIRREL1, DNAH9 and XYLT1 up-regulated more than 2-fold in the tubulointerstitium while being more than 2-fold down-regulated in the glomeruli.

Data validation by qPCR

Eight differentially expressed genes and two podocyte-specific markers were cross-validated using qPCR and showed a good agreement between the technologies with respect to case–control relative gene expression levels (Supplementary data, Figure S3).

Enrichment analyses

Gene ontology (GO) enrichment analysis showed that in DN glomeruli (Figure 2A), the top enriched biological pathway (BP) categories were cell adhesion, inflammatory response, complement activation and extracellular matrix (ECM) organization, while the top enriched cellular component (CC) was ECM. In the DN tubulointerstitium (Figure 2B), angiogenesis, apoptosis and oxidation–reduction processes were enriched, and mitochondria was highlighted as a top subcellular component. Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis confirmed and further refined these findings, featuring ECM-receptor interaction, infectious and inflammatory pathways, complement and coagulation cascades, haematopoiesis, Phosphoinositide 3-kinase-Akt (PI3K-Akt) and Hippo signalling pathways (in the DN glomeruli), and cytokine–cytokine receptor interaction, immunity, phagosome, mineral absorption and mitogen-activated protein kinase (MAPK) pathway (in the DN tubulointerstitium) (Table 2).

FIGURE 2

GO enrichment results. (A and B) GO enrichment results in glomeruli (A) and tubulointerstitium (B) are presented as bubble plots visualizing the z-score (x-axis) in relation to the negative logarithm of the adjusted P-value (y-axis). The area of the displayed circles is proportional to the number of genes assigned to the term and the colour corresponds to the GO category. The terms legend on the right provides description of the enriched GO terms above the significance threshold (orange horizontal line). BP: biological process; CC: cellular component; MF: molecular function.

Table 2

Significantly changed pathways and genes

Group of differential expression	KEGG term (ID)	Overlap	P-value, adjusted	Combined score	Genes
Up-regulated in glomeruli	ECM-receptor interaction (hsa04512)	15/82	2.94E-09	42.7510506	TNXB; ITGA4; LAMC3; ITGB4; LAMA4; FN1; HMMR; THBS2; COL1A1; COMP; COL1A2; COL6A3; CD36; ITGB6; CD44
	Complement and coagulation cascades (hsa04610)	14/79	1.11E-08	41.1019959	C1QB; C1QA; ITGB2; F13A1; CLU; C2; C4B; C3; C4A; C8G; C7; VSIG4; CFB; C1QC
	Staphylococcus aureus infection (hsa05150)	11/56	2.41E-07	34.0328219	C4B; C1QB; C3; C4A; C1QA; SELPLG; ITGB; ITGAL; CFB; C1QC; C2
	Amoebiasis (hsa05146)	13/100	1.3E-06	31.2981641	LAMC3; ACTN1; LAMA4; ITGB2; FN1; CXC1; PIK3CG; COL1A1; COL3A1; COL1A2; C8G; PLCB4; CTSG
	PI3K-Akt signalling pathway (hsa04151)	22/341	5.17E-06	28.4396763	CSF1R; TNXB; ITGA4; LAMC3; ITGB4; HGF; LAMA4; LPAR1; FN1; IL2RG; THBS2; PIK3CG; COL1A1; COMP; KITLG; COL1A2; CCND2; GNG4; IL2RB; COL6A3; ITGB6; PIK3AP1
	Haematopoietic cell lineage (hsa04640)	12/88	2.09E-06	28.0313106	CD2; CSF1R; KITLG; ITGA4; CD7; CD3G; CD36; CD3E; MS4A1; CD3D; CD33; CD44
	Focal adhesion (hsa04510)	17/202	3.47E-06	27.6027344	TNXB; ITGA4; LAMC3; ITGB4; ACTN1; HGF; LAMA4; FN1; THBS2; PIK3CG; COL1A1; COMP; COL1A2; CCND2; RAC2; COL6A3; ITGB6
	Pertussis (hsa05133)	11/75	2.93E-06	25.6397118	C4B; PYCARD; C1QB; C3; C4A; CXCL6; C1QA; ITGB2; IRF8; C1QC; C2
	Systemic lupus erythematosus (hsa05322)	14/135	3.47E-06	25.2587947	CD86; C1QB; C1QA; ACTN1; C2; C4B; C3; C4A; GRIN2A; C8G; C7; CD28; CTSG; C1QC
	Cytokine–cytokine receptor interaction (hsa04060)	18/265	2.85E-05	22.3852033	CX3CR1; CCL13; CSF1R; CXCL6; CCL21; CD70; IL10RA; HGF; TNFRSF11B; CXCL1; IL2RG; KITLG; IL2RB; CD27; FAS; TNFRSF17; LTB; CCR2
Down-regulated in glomeruli	Hippo signalling pathway (hsa04390)	19/153	5.72E-05	26.0480563	TCF7L1; FZD2; SERPINE1; PARD6G; FZD8; FGF1; MPP5; WTIP; BMP7; AREG; NKD1; FRMD1; LATS2; DLG2; DVL1; CTNNA3; AMH; WNT3; TEAD3
Up-regulated in tubulointerstitium	Staphylococcus aureus infection (hsa05150)	31/56	1.89E-14	67.963694	CFD; C1QB; C1QA; SELPLG; ITGAM; CFH; C1S; C1R; ITGB2; PTAFR; FPR3; ITGAL; C2; C3; FCGR3A; HLA-DMA; HLA-DMB; C3AR1; HLA-DOA; FCGR1A; HLA-DOB; HLA-DQA1; HLA-DPA1; IL10; SELP; FCGR2A; HLA-DPB1; FCGR2B; FCGR2C; CFB; C1QC;
	Cytokine–cytokine receptor interaction (hsa04060)	69/265	2.19E-11	54.4638677	CNTFR; CXCL6; CXCL9; CSF3R; TNFRSF13B; IL23R; IL24; IL5RA; FASLG; CXCL1; TNF; EDA2R; TNFSF13B; TNFSF1; CCR7; CCR6; CCR5; TNFRSF4; PF4V1; CCR2; IL10; PDGFRB; IFNAR2; PDGFRA; HGF; TNFRSF18; IL19; OSMR; TNFRSF1B; IFNG; XCL2; XCL1; LTB; CX3CR1; CCL13; CSF1R; FIGF; TNFRSF11B; CSF2RB; CXCR6; IL2RG; CSF2RA; CXCR3; CCL5; IL21R; TNFRSF17; CCL19; IL12RB1; IL12RB2; NGFR; TGFB2; XCR1; CCL22; TGFB1; TNFSF14; CCL21; CD70; IL10RA; TNFRSF9; TNFRSF10C; GDF5; CD40LG; IL2RA; IL2RB; CD27; TNFSF8; IL7R; CCL28; CCL26
	Primary immunodeficiency (hsa05340)	20/37	3.61E-09	39.3871413	CIITA; TNFRSF13B; TAP1; IL2RG; CD3E; CD3D; ZAP70; CD79A; CD4; PTPRC; CD40LG; CD8B; LCK; IGLL1; CD8A; CD19; BTK; IL7R; ICOS; JAK3
	Haematopoietic cell lineage (hsa04640)	31/88	1.76E-08	37.3935788	CSF1R; CSF3R; ITGAM; ITGB3; IL5RA; CD3G; CD1E; CD1D; CD3E; CD1C; TNF; CD3D; CSF2RA; CD19; CD38; CD37; CD36; FCGR1A; CD33; ITGA4; CD2; CD4; CD8B; CD5; CD8A; IL2RA; CD7; IL7R; MS4A1; CD22; CD44
	Leishmaniasis (hsa05140)	27/73	4.72E-08	35.2730309	ITGAM; NCF1; NCF2; ITGB2; NCF4; TNF; C3; FCGR3A; HLA-DMA; HLA-DMB; HLA-DOA; FCGR1A; HLA-DOB; HLA-DQA1; HLA-DPA1; IL10; TGFB2; MARCKSL1; TGFB1; ITGA4; STAT1; PRKCB; FCGR2A; IFNG; HLA-DPB1; PTPN6; FCGR2C
	Cell adhesion molecules (hsa04514)	41/142	2.45E-08	34.954997	CD86; SELPLG; ITGAM; SDC3; ITGB2; ICAM3; ITGAL; CLDN1; SPN; HLA-DMA; CDH3; HLA-DMB; CTLA4; NCAM1; HLA-DOA; ICOS; TIGIT; HLA-DOB; HLA-DQA1; HLA-DPA1; CADM3; VCAM1; NLGN4X; ITGA4; PDCD1LG2; CD2; SELP; CD4; VCAN; PTPRC; CD40LG; CD6; CD8B; SELL; CD8A; CD28; HLA-DPB1; CD226; PDCD1; SIGLEC1; CD22
	Osteoclast differentiation (hsa04380)	38/132	8.2E-08	34.0189588	CSF1R; SPI1; NCF1; NCF2; ITGB3; SIRPG; NCF4; TNFRSF11B; TREM2; LILRA1; LILRA2; TNF; SIRPB1; PIK3CG; PIK3R5; FCGR3A; CTSK; TNFSF11; FCGR1A; IFNAR2; TGFB2; TGFB1; STAT1; CYBB; LILRB1; LILRB2; LILRB4; LILRB5; OSCAR; TYROBP; FCGR2A; IFNG; LCK; CAMK4; BTK; LCP2; FCGR2B; FCGR2C
	Chemokine signalling pathway (hsa04062)	44/187	2.42E-06	26.8161614	CX3CR1; ITK; CXCL6; CCL13; CXCL9; NCF1; WAS; ADCY3; ADCY2; CXCL1; CXCR6; RASGRP2; ADCY7; PIK3CG; PIK3R5; PREX1; GNG2; CXCR3; CCL5; RAC2; CCR7; CCR6; CCL19; CCR5; JAK3; PF4V1; CCR2; XCR1; CCL22; CCL21; PRKCB; STAT1; VAV1; FGR; HCK; PLCB4; ELMO1; GNB4; XCL2; XCL1; DOCK2; CCL28; PLCB2; CCL26
	Phagosome (hsa04145)	39/154	2.19E-06	25.3724165	COLEC12; ITGAM; NCF1; NCF2; C1R; ITGB3; NCF4; ITGB2; TCIRG1; THBS2; CORO1A; CTSS; COMP; C3; MRC2; FCGR3A; HLA-DMA; TUBA1A; HLA-DMB; TUBB3; CLEC7A; MRC1; CD36; HLA-DOA; FCGR1A; HLA-DOB; HLA-DQA1; HLA-DPA1; MSR1; TAP1; CYBB; TUBB4A; MARCO; FCGR2A; CD209; HLA-DPB1; TLR6; FCGR2B; FCGR2C
	ECM-receptor interaction (hsa04512)	26/82	2.31E-06	23.8605696	TNXB; LAMC3; ITGB4; LAMA4; ITGB3; TNC; LAMC2; THBS2; COMP; RELN; SV2B; SV2A; CD36; ITGB6; ITGA4; VWF; FN1; COL1A1; COL1A2; COL4A2; COL4A1; COL6A2; ITGA11; COL6A3; COL6A6; CD44
	Intestinal immune network for IgA production (hsa04672)	19/48	2.31E-06	23.5968068	IL10; CD86; PIGR; TGFB1; ITGA4; TNFRSF13B; TNFSF13B; HLA-DMA; CD40LG; HLA-DMB; CD28; HLA-DPB1; TNFRSF17; ICOS; HLA-DOA; CCL28; HLA-DOB; HLA-DQA1; HLA-DPA1
	Inflammatory bowel disease (hsa05321)	22/65	4.65E-06	23.0034505	IL10; TGFB2; TGFB1; STAT1; IL23R; NOD2; IL2RG; FOXP3; TNF; HLA-DMA; HLA-DMB; IFNG; TBX21; IL21R; STAT4; HLA-DPB1; IL12RB1; HLA-DOA; HLA-DOB; HLA-DQA1; IL12RB2; HLA-DPA1
	Malaria (hsa05144)	19/49	2.84E-06	22.9538282	IL10; TGFB2; TGFB1; VCAM1; GYPC; LRP1; KLRB1; HGF; ITGB2; ITGAL; THBS2; TNF; COMP; SELP; CD40LG; KLRK1; IFNG; CD36; ACKR1
	T cell receptor signalling pathway (hsa04660)	28/104	2.13E-05	20.2904308	ITK; CD3G; CD3E; TNF; CD3D; PIK3CG; PIK3R5; GRAP2; CTLA4; PAK3; ICOS; IL10; VAV1; ZAP70; CD4; PTPRC; CD40LG; CD8B; IFNG; CD8A; LCK; CD28; PTPN6; LCP2; PDCD1; CD247; CARD11; LAT

Hypothesis-based exploration of targeted genes

Kidney-specific genes

We tested a list of 334 genes with preferential renal expression [Human Protein Atlas (HPA); www.proteinatlas.org] [34] for differential expression in DN. HPA defines three categories of genes with elevated expression in kidney compared with other organs: kidney-enriched (at least 5-fold higher mRNA levels as compared with all other tissues, n = 54), group-enriched (at least 5-fold higher mRNA levels in a group of two to seven tissues, n = 119) and kidney-enhanced (at least 5-fold higher mRNA levels as compared with average levels in all tissues, n = 161). We observed dramatic down-regulation of all categories of kidney-specific transcriptome in the DN versus LD tubulointerstitium (Figure 3A).

FIGURE 3

Hypothesis-driven analyses of renal transcriptome. (A–D) Volcano plots for biological groups of genes with renal relevance. The plots visualize log2-fold changes in gene expression (x-axis) in relation to the negative logarithm of the adjusted P-value (y-axis) for the DN versus LD transcriptomes comparison. The red horizontal dotted line shows the significance threshold at adjusted P = 0.05. The black vertical dotted line indicates zero-fold change. Overall, the gene expression of marker genes and Mendelian glomerulopathy genes found in DN patients resembles other known kidney diseases. (E) Biomarkers: case–control tubulointerstitial gene expression (iand iii) and corresponding creatinine-adjusted urinary concentration (iiand iv) of proximal tubular injury biomarkers: KIM-1 and NGAL showing a strong correlation between gene and protein enrichment in diseased kidneys.

Renal cell type-specific genes

A collection of bona fide cell-specific markers of resident renal cells [35] was tested for differential expression in our data, shown in Figure 3B and C. Specifically, we evaluated markers of podocytes, mesangial and endothelial cells, proximal and distal tubular epithelial cells, fibroblasts, vascular endothelial and smooth muscle cells. We saw a marked down-regulation of many differentiated podocyte (e.g. MAGI2, SYNPO, NPHS1, WT1) and tubular genes (e.g. SLC12A3, PVALB, BHMT, UMOD) in DN patients (Figure 3B and C). On the other hand, mesangial (PTGR, COL4A1, ACTA2) and fibroblast (NGFR, PDGFRB, S100A4, VIM) markers were up-regulated.

Mendelian kidney disease genes

From Parsa et al. [36], Mendelian genes responsible for various glomerulopathies (n = 172) and tubulopathies (n = 55) were extracted. Thirty-four glomerulopathy genes were differentially regulated in DN versus LD (P < 0.05), which was a significant enrichment as compared with the global transcriptome (P = 0.04). The majority of the genes (26/34) were down-regulated (Figure 3D). Correspondingly, 14 tubulopathy genes were differentially expressed (12 down- and 2 up-regulated) in our tubulointerstitial samples (data not shown), however they did not differ from the overall transcriptome-wide trends (P = 0.675).

Urinary biomarker genes

We also tested a set of genes encoding known protein urinary biomarkers [37, 38] for differential expression in our data. Filtration markers originating from extra-renal tissues were excluded from this analysis as we focused only on the intra-renal gene products along various segments of the nephron. Tubular damage markers genes LCN2, HAVCR1, HGF, CLU and TIMP2 [encoding proteins NGAL, KIM-1, hepatocyte growth factor (HGF), Clusterin and TIMP-2, respectively] were significantly up-regulated in the DN tubulointerstitium as compared with LD (Supplementary data, Table S5). Accordingly, creatinine-normalized urinary concentrations of both NGAL and KIM-1, assayed immediately before the biopsy, were also elevated in DN patients (Figure 3E) and significantly correlated to LCN2 and HAVCR1 mRNA levels (data not shown). In contrast to the damage markers, expression of EGF, a marker of tubular integrity, was substantially decreased in DN.

Clinical variables and WGCNA

On univariate analyses, no correlations between expression of individual genes and clinical variables were significant after multiple testing adjustment. Instead, we employed a modular, as opposed to univariate, approach using WGCNA of the glomerular transcriptome and identified 55 gene modules (Figure 4A). The obtained modules were summarized into eigengenes representing the weighted average of the expression profile for each module and were tested for correlations with relevant clinical parameters (Figure 4B). Two distinct association patterns emerged: Darkolivegreen, Lightcyan1, Steelblue, Yellow, Turquoise, Cyan and Darkgreen gene modules were associated with better renal function (positive correlation with GFR, serum albumin, haemoglobin; negative correlation with serum creatinine, UACR, high-sensitivity C-reactive protein, blood pressure), while many others showed the inverse tendency. The grey module, which contained all the genes that were not assigned to any module and was kept in the analysis as a negative control, showed no significant correlations. The Darkred module (225 genes), exemplifying a module with strongly negative association with renal function, was chosen for subsequent functional analysis. The subnetwork of the Darkred module, overlaid with the DEA results (Figure 4D), shows the predominance of disease-upregulated genes in its core. A list of genes for each WGCNA module is shown in Supplementary data, Table S6. KEGG pathway enrichment analysis revealed Wnt signalling, complement and coagulation cascades, as well as chemokine signalling (Figure 4E).

FIGURE 4

Glomerular WGCNA. (A) Gene dendrogram obtained by hierarchical clustering of adjacency-based dissimilarity. The colour row below the dendrogram indicates module membership identified by dynamic tree cut followed by merging of highly similar modules. (B) Module–trait associations. Rows correspond to eigengenes from respective modules, columns to clinical variables. Each cell contains the corresponding correlation coefficient and P-value. Colour reflects direction of the correlation (red: positive, blue: inverse), with colour intensity being proportional to the correlation strength (correlation coefficient). Red box shows the Darkred module as an example of a strongly renal function-associated module which positively correlates to creatinine and UACR and negatively correlates to GFR and albumin. (C) Whole network view with the Darkred module genes highlighted. (D) Subnetwork of the Darkred module. The colour reflects log2-fold changes from the DEA (red: up-regulated genes, blue: down-regulated genes). (E) Top enriched KEGG pathways.

DISCUSSION

In this comprehensive RNAseq study, we present novel information on DN-associated transcripts and pathways in human glomerular and tubulointerstitial tissue fractions. This extensive data resource should be utilized to illuminate DN biology and to identify potential biomarkers relevant to the diagnosis and prognosis of DN patients.

Despite DN being the most common cause of ESRD in the Western populations, with a considerable elevated risk for cardiovascular morbidity and mortality, there is insufficient knowledge about the molecular mechanisms and pathways responsible for onset and progression of disease. This has several explanations, including infrequent biopsy collection, poor translatability from animal models to the human setting and a likely heterogeneity of underlying molecular pathogenesis, as the clinical course of DN varies tremendously. Limitations of the analytical methods have also placed restrictions on the ability to generate high-quality data with low-quantity input material. The current study addresses these obstacles by accessing 39 paired glomerular and tubulointerstitial tissue samples from patients with biopsy-verified DN and LD with matched urine collected around the time of biopsy. The use of LD tissue as control is superior to other more commonly used control tissues from nephrectomies [39-41], most importantly because they provide less ischaemic tissue damage. A further strength of our study is the carefully chosen methodologies, including RNA amplification and cDNA library synthesis protocols that were refined and optimized specifically for low amounts of starting material, as well as the RNAseq technology that allows unbiased analyses of both the known and the unknown transcriptome. Together, this provides higher data specificity and sensitivity with a wider dynamic range in comparison with transcriptional profiles generated by microarray, which are restricted to the set number of gene probes available on the array.

Our study population ranged from CKD Stages 1–4, with most patients having renal function corresponding to CKD Stages 3 and 4. At this stage, early morphological changes, such as deposition of ECM causing both mesangial matrix expansion and glomerular baseline membrane thickening, have already occurred [42]. Consistent with this, we found that the top enriched CC in the glomerular compartment was ECM, along with upregulation of ECM pathways and organization, as well as cell adhesion. The KEGG enrichment analysis of the glomerular compartment further highlighted ECM–receptor interaction, but also infectious and inflammatory conditions, as well as complement and coagulation cascades. Both inflammatory response and complement activation were also displayed as top BP categories altered in DN glomeruli. As our RNAseq data were generated on histologically verified DN kidneys in earlier disease stages (likely with more intact glomeruli with less fibrotic tissue), using larger sample size and LDs as controls, than previous transcriptomic studies [39, 43–45], these results corroborate and add robustness to the available literature. Together, this confirms the importance of inflammation in DN, a feature that is known to increase the risk for cardiovascular events [46]. Genes in the tubulointerstitial compartment were mainly up-regulated; activated GO pathways included apoptosis and oxidation–reduction. Mitochondria was highlighted as a top subcellular component, which is consistent with the oxidative stress known to occur in DN.

Whereas several transcriptomic studies on human DN kidneys focused on a single compartment, i.e. either the glomeruli [40] or the tubuli [41, 45, 47], or even used undissected tissue [44], our study clearly demonstrates disease-specific gene expression profiles in both tissue compartments. This observation, which is also in support of data from Berthier et al. [48] and Woroniecka et al. [39], underscores the complexity of the disease and the importance to study not only the diabetic glomeruli but also the tubules. Indeed, the role of the tubular system in disease progression has been emphasized previously [49, 50]. Hypertrophy of proximal tubules is observed early in DN and is thought to promote increased glucose and sodium re-uptake, contributing to glomerular hyperfiltration [51]. As DN progresses, tubules become atrophic and interstitial fibrosis and glomerulosclerosis increases. Changes in the DN transcriptome reported here are consistent with this. We believe this is an example of how transcriptomic profiles reflect tissue morphology, and this highlights the usefulness of integrating transcriptomics with the histopathological assessment in the future.

The detailed insight into the compartmental changes taking place also revealed an overlap of differentially expressed genes between the glomerular and tubulointerstitial fraction. Since many of the processes in these cell types are likely to be shared, this is not unexpected. While it is necessary to bear in mind that this may indicate technical difficulties in the micro-dissection procedure, our testing of tubular and glomerular-specific transcripts indicates that this was unlikely. Also, kidney-enriched as well as cell-specific genes were down-regulated in DN, potentially reflecting a loss of cellular differentiation or changes in cell numbers and composition—or both. This is in line with previous human studies [39, 40] and is a known feature of progressing glomerular disease [52]. Single-cell RNAseq studies, with side-by-side comparison of specialized cell types, will likely enhance the resolution of DN pathology even more and help to disentangle the issue on cell loss versus cell type-specific changes.

Our analyses of a curated gene set with Mendelian genes responsible for various glomerulopathies and tubulopathies [36] showed that these genes were down-regulated, which is directionally consistent with the hypothesis of a phenocopy and a loss-of-function mutation effect.

Fourteen tubulointerstitial genes were significantly differentially regulated in DN cases versus LD controls, suggesting that DN may share several common pathogenic mechanisms with inherited nephropathies. Among non-Mendelian genes, the genes corresponding to tubular biomarkers u-NGAL and u-KIM-1 (gene symbols LCN2 and HAVCR, respectively) were significantly increased in DN patients. Since both u-NGAL and u-KIM-1 are increased in early stages of DN, indicating disease state even before elevation of UACR [53], we also analysed the urine protein levels, which fully reflected the gene expression changes. This is the first time a correlation between protein urine levels and tubular gene expression of these biomarkers is reported, corroborating how the use of non-invasive biomarkers mirrors the kidney molecular profile. This result serves as a clear example of the potential of our dataset to find new biomarker candidates.

Some limitations of the present study should be acknowledged. Although this is the largest DN cohort studied with RNAseq to date, it may still be underpowered to detect significant changes. This limitation is mainly a result of the fact that biopsies are collected for diagnostic purposes only from a minority of all diabetic kidney disease patients, also introducing a risk for selection bias. Additionally, some site-specific differences in transcriptomic profile were seen. These site/batch effects were corrected for in all downstream analyses to not influence the final results. Finally, as this study has a cross-sectional design, it does not provide information on causality or prognostic factors. Our cohort includes both DM Types 1 and 2 patient samples, and since the DN disease has different primary causes this may affect the results, although our dataset is too small for stratification.

To conclude, the present study adds important detailed information at an unprecedented resolution on the DN-specific transcriptomic profile of both glomeruli and tubulointerstitial fractions. The demonstrated sensitivity of the RNAseq technique and the potential usefulness for transcriptomic signatures to be included alongside traditional histopathological methods for diagnostic purposes pave the way for future studies on biomarkers and treatment targets.

Click here for additional data file.