Proximal tubule cell maturation rate and function are controlled by PPARα signaling in kidney organoids

hiPSC culture and the differentiation of kidney organoids

Kidney organoids were generated from the hiPSC line CRL1502.C32 according to our previously published protocol33 with some modifications (Fig. 1a). Briefly, hiPSCs were treated with 8 μM CHIR99021 (TOCRIS, 4423) in APEL2 (StemCell Technologies, ST-05275) supplemented with 1% Protein Free Hybridoma Medium II (PFHM II, GIBCO, 12040077) (basal medium) for 5 days, followed by treatment with FGF9 (200 ng/mL, R&D, 273-F9) and heparin (1 μg/mL, Sigma, H4784) for another 3 days. The cells were then collected and dissociated into single cells using TryPLE select (GIBCO, 12563011); 1 × 105 collected cells were seeded into PrimeSurface® 96M plates (Sumitomo Bakelite) and were cultivated for 1 day in a CO2 incubator to make spheroids. On the following day, these spheroids were transferred onto transwells with a 0.4 μm pore polyester membrane (Corning, 3450) and cultured at the liquid-air interface. Spheroids were treated with 10 μM CHIR99021 in basal medium for 3 h and then cultured with FGF9 (200 ng/mL) and heparin (1 μg/mL) for another 6 days, followed by another 16 days in basal medium with changing the medium every other day. For the PPARα agonist or antagonist treatment, kidney organoids were treated with 0.3 μM CP775146 (PPARα agonist, TOCRIS, 4190) and 1 μM 9-cis-retinoic acid (RXRα agonist, Abcam, ab141023) as agonist treatment group or 3 or 10 μM GW6471 (PPARα antagonist, TOCRIS, 4618) and 1 or 10 μM HX531 (RXR antagonist, TOCRIS, 3912) as antagonists from day 15 in basal medium with medium change performed every other day.

Single-cell isolation from kidney organoids

For preparing single cells from kidney organoids at days 15, 19, 23, 27 and 31 after the induction of differentiation, experiments were performed using 30–36 organoids from the same experimental run (36 organoids on day 15, 30 organoids on all other days). Kidney organoids were incubated for 30 min at 6 °C in cold active protease (DPBS(+), 5 mg/mL Bacillus Licheniformis Protease (Sigma, P5380), 125 U/mL DNase (Roche, 11284932001)34 with pipetting 30 times every 5 min. Each cell suspension was diluted with fluorescence-activated cell sorting (FACS) buffer (DPBS(−), 0.5% FBS, 2 mM EDTA) + 125 U/mL DNase, centrifuged at 300×g for 10 min at 4 °C after which the supernatant was removed. This step was repeated twice, after which the suspension was passed through a 35 μm filtered tube (Falcon) to prepare a cell suspension containing single cells. Dead cells and Doublets were removed from these cells by BD FACSAria™ II & III using an Annexin V-FITC Apoptosis Detection Kit (Nakalai, 15342-54).

Library preparation, sequencing and data pre-processing

Cell suspensions were adjusted to 1000 cells/μL, and libraries were prepared using Chromium Single Cell 3’ reagent Kits (V3) according to the manufacturer’s instructions. Libraries were sequenced using an Illumina NovaSeq 6000 Sequencing system with an average depth of 60,653 reads per cell, then mapped to the human genome (build GRCh38) and demultiplexed using Cell Ranger pipelines (10X Genomics, v.3.0.2). After performing a cellranger count for each sample, all samples were aggregated with cellranger aggr to create a gene expression matrix.

Data processing and cell clustering

Pre-processed data from samples were further processed and analyzed individually using the R package Seurat (v.3.1.1) on RStudio (v.3.6.1)35. Clustering was performed using Seurat. For quality control of the data, UMI counts of cells, mitochondrial gene detection rate, and gene detection rates were prepared in violin plots, and then nCount_RNA > 5000, percent.mt < 30, nFeature_RNA < 7500 were filtered. The number of reads per droplet was then aligned to 10000 and normalized on a log scale. For clustering analysis, 2000 highly variable genes were extracted using the variance in expression levels and were Z-scored per gene using data from the entire sample. Notably, gene dependence on UMI counts, mitochondrial gene detection rate, S-phase score and G2M-phase score were corrected. PCA was performed using the Z-scored expression levels, and from JackStrawPlots and ElbowPlots, the number of principal components for clustering was set to the top 20 and the resolution to 0.4. Cell types were identified by examining the expression of marker genes and cluster-specific genes of putative cell types.

Sub-clustering

To identify detailed cell types that could not be classified by the parameters used for the whole sample, sub-clustering was performed for clusters (hereafter CL) 2, CL12, CL13 (NPC/PTA), CL4, CL5, CL9 (Tubules), and CL0, 1, 7 (Glomeruli). Clustering analysis was performed in the same way as the method for clustering all cells described above after extracting the clusters of interest, normalization, extraction of high variation genes, calculation of cell cycle score, Z-scoring, PCA, selection of a number of principal components, clustering, extraction of marker genes and cluster-specific genes, and annotation was performed for the sub-clusters.

Pseudotime cell trajectory in kidney lineage cells

Samples of kidney cell lineages CL2, CL12, CL13 (NPC/PTA), CL4, CL5, CL9 (Tubules), and CL0, 1, 7 (Glomeruli) were extracted and differentiation trajectories were calculated using Monocle2 (ver.2.14.0)36. Using the extracted cells, PCA was performed on 1624 genes with large variance, and the top 10 principal components were used to compress the data into 2 dimensions using tSNE. Using tSNE density and peak data, clustering was performed with rho > 100 and delta > 15, which resulted in 8 clusters. Specific genes were detected among the clusters and a Pseudotime cell trajectory plot was created using 2021 genes with q-values of 0. To check whether the genes were aligned on the differentiation axis, a branched heatmap was created using known marker genes.

Pseudotime cell trajectory in PT cells

Cells in the PT state of the pseudotime cell trajectory (PT cells in Fig. 2b) and cells annotated in the developmental lineage of the proximal tubule by UMAP sub-clustering (Supplementary Fig. 1f: cells in RV, Medial SSB, PT1 and PT2 clusters) were extracted. These cells were analyzed similarly to the analysis of pseudotime cell trajectory in kidney lineage cells to detect specific genes, and the top 320 genes with the lowest q-values were extracted. Pseudotime cell trajectory plots were created excluding mitochondrial genes and ribosomal protein-coding genes from the 320 genes. Maturation markers were obtained by extracting DEGs from the latter part of the pseudotime cell trajectory plot (Pseudotime 14-19) and other cells and selecting those that were more highly expressed in proximal tubule cells.

DEG analysis and IPA analysis

DEGs were extracted using Seurat’s function FindMarkers, using Poisson’s method. Upstream regulators were analyzed by Ingenuity Pathway Analysis (IPA, Qiagen) using DEGs extracted by the above method.

Isolation of LTL-positive cells from kidney organoids

LTL-biotin conjugated (Vector Laboratories, B-1325-2, 1:50) was added to the prepared single cell suspension and incubated at 4 °C for 15 min. The samples were then centrifuged at 300×g for 5 min at 4 °C, and the supernatant was removed. The pellet was washed twice with FACS buffer. Anti-biotin micro Beads Ultra Pure (Miltenyi) were then added and incubated at 4 °C for 15 min, after which the samples were centrifuged at 300×g for 5 min at 4 °C and the supernatant was removed. The pellet was washed twice with FACS buffer. The cell suspension was purified by the MACS method using an MS column (Miltenyi) fixed in a magnetic field to obtain LTL-positive cells.

qPCR and heatmap

Total RNA was extracted from the LTL-positive proximal tubule cells obtained using a Purelink RNA micro kit (ThermoFisher Scientific). The total RNA obtained was used as a template for cDNA synthesis using Prime Script RT master mix (Perfect Real Time) (TaKaRa). The synthesized cDNA was quantified by real-time PCR using TB Green Premix Ex Taq II (Tli RNaseH Plus) (TaKaRa) for expression analysis. The sequences of primers used for qRT-PCR are listed in Supplementary Table 1. Heatmaps were created using Heatmap2 in R packages, and the tree was analyzed using Ward’s method.

IF staining

Organoids were fixed with 4% paraformaldehyde in DPBS (−) for 30 min at 4 °C followed by 3 washes with DPBS (−). The organoids were then blocked with 10% donkey serum, 0.3% Triton X/DPBS (−) for 3 h at room temperature and incubated with primary antibodies overnight at 4 °C. After 5 washes with 0.3% Triton X/DPBS (−), secondary antibodies were incubated overnight at 4 °C. The following antibodies and dilutions were used: mouse anti-ECAD (1:300, 610181, BD Biosciences), sheep anti-NPHS1 (1:300, AF4269, R&D Systems), LTL-biotin-conjugated (1:300, B-1325, Vector Laboratories), rabbit anti-Megalin (1:300, ab76969, Abcam), mouse anti-γH2AX (1:100, DNA Damage Detection Kit γH2AX Red, G266, Dojindo Laboratories). Images were taken using a Zeiss LSM800 confocal microscope or CQ1 (Yokogawa). All immunofluorescence analyses were successfully repeated more than three times and representative images are shown.

Quantification of megalin fluorescence intensity

Fluorescence imaging of stained kidney organoids was performed using a Zeiss LSM800 confocal microscope. Z-stack images were acquired with a 10x objective at 3.4 µm intervals. Maximum intensity projection (MIP) was applied to both megalin and DAPI fluorescence channels to generate 2D images from the 3D datasets. Megalin-positive areas and fluorescence intensities were quantified from the MIP images using Fiji software. Megalin fluorescence intensity per unit area was calculated by dividing the total megalin fluorescence intensity by the corresponding megalin-positive area. These values were then normalized to the average intensity of the day 20 control group, which was set as a baseline value of 1.

scRNA-seq analysis of PPARα agonist-treated kidney organoid samples

Kidney organoids were differentiated according to the previously described protocol until day 15. From day 15 to day 23, the control group was cultured in base medium supplemented with 0.2% DMSO, while the PPARα agonists-treated group received base medium supplemented with 0.3 μM CP775146 and 1 μM 9-cis-retinoic acid. ScRNA-seq data were generated as described for the time-course experiment. Libraries were sequenced on an Illumina NovaSeq 6000 system, yielding an average depth of 165,692 reads per cell. Quality control filtering was applied to retain cells with mitochondrial gene content <25% (percent.mt, <25) and gene counts between 2500 and 9500 (2500 <; nFeature_RNA < 9500). Control and PPARα agonists-treated group datasets were integrated using the Seurat package’s anchor-based CCA integration workflow (FindIntegrationAnchors and IntegrateData functions). Integrated expression data were Z-score normalized, and PCA was performed. Clustering was subsequently conducted using the top 20 principal components and a resolution of 0.2.

Integration of kidney organoids and Kidney Cell Atlas (KCA) fetal kidneys

Cells annotated as “Cap mesenchyme,” “Proliferating cap mesenchyme,” “Proximal renal vesicle,” “Distal renal vesicle,” “Proliferating distal renal vesicle,” “Proximal S-shaped body,” “Medial S-shaped body,” “Distal S-shaped body,” “Proximal UB,” “Podocyte,” “Proximal tubule,” “Loop of Henle,” “Endothelium,” “Stroma progenitor,” and “Proliferating stroma progenitor” were extracted from the KCA fetal metadata, assuming their presence in the kidney organoids.

The KCA fetal data from six donors and two kidney organoid samples were integrated using Seurat’s anchor-based CCA integration workflow (FindIntegrationAnchors and IntegrateData functions). Integrated expression data were Z-score normalized, and PCA was performed using the top 20 principal components for clustering with a resolution of 0.2.

Pseudotime generation

Pseudotime was generated using the KCA metadata. Differentially expressed genes between proximal tubule and mSSB cell groups in the KCA fetal samples were identified (FindMarkers: logfc.threshold = 0.5, min.pct = 0.5, test.use = “poisson”). PT cells from KCA fetal cluster 4 (CL4) and three KCA child samples (Wilms1, Wilms2, Wilms3) were merged, and PCA was performed using the identified gene set.

Pseudotime extrapolation

The generated pseudotime was used to extrapolate and compare the maturity of proximal tubules in organoids and adult kidneys. Pseudotime of organoid proximal tubule cells was calculated from the product of normalized expression in CL4 and the loading of the first principal component in the PCA.

Fetal kidney data for extrapolation was obtained from GSE114530. Proximal tubule cells (classified as CL4 using TransferData) were extracted, and pseudotime was calculated as described for organoids.

Adult kidney count data for extrapolation was obtained from GSE202109, with metadata from the UCSC Cell Browser. Proximal tubule cells were extracted, and pseudotime was calculated similarly.

Transmission electron microscopy (TEM)

The tissue samples for TEM were fixed in phosphate-buffered 2% glutaraldehyde and subsequently post-fixed in 2% osmium tetra-oxide for 3 h in an ice bath. The specimens were then dehydrated in a graded ethanol series and embedded in epoxy resin. Ultrathin sections were obtained using an ultramicrotome. Ultrathin sections were stained with uranyl acetate for 10 min and a lead staining solution for 5 min and then observed by TEM (HITACHI H-7600 at 100 kV). The acquired images were used for quantification by Fiji to measure the area of the cytoplasm of the proximal tubule cells and the number and size of endosomes. As per the analysis, 24 images were used for the control group and 23 images for the PPARα agonists-treated group, with 26 cells extracted from each.

Dextran uptake assay

Kidney organoids on day 26 of differentiation cultured at the air-liquid interface in transwells were added to medium with pHrodo™ Red Dextran, 10,000 MW, for Endocytosis (Invitrogen, P10361) at 10 mg/mL and incubated at 37 °C in a 5% CO2 incubator for 24 h. Organoids were then washed 3 times with ice-cold DPBS (−) and loaded into DPBS (−) containing LTL-FITC ((1:100) Vector Laboratories, FL-1321) for 2 h at 4 °C. Stained kidney organoids were observed for fluorescence using a CV8000 (Yokogawa). Images were taken with an X20 lens, and after acquiring 2 × 2 tiled images, images were taken every 5 µm of 100 µm thickness as a Z-stack. Average Intensity Projection (AveIP) was used for LTL fluorescence and Sum Intensity Projection (SumIP) was used for Dextran fluorescence to create 2D images from 3D images. Data analysis was performed using CellPathfinder (Yokogawa). The fluorescence intensity of dextran and LTL-positive region was measured. Then, the fluorescence intensity of dextran divided LTL-positive region. These values were then normalized to the average value of the control group, which was set as 1.

Flow cytometry analysis of kidney organoids treatment with cisplatin

Kidney organoids on day 23 of differentiation were cultured for 6 days in medium containing 5 μM cisplatin (cisplatin intravenous infusion 10 mg “Marco” (Nichi-Iko Pharma)). The medium was changed every other day. The kidney organoids were dissociated into single cells by the same method used for single-cell isolation from kidney organoids, and then the cells were suspended in 4% PFA/DPBS (−) and incubated at 37 °C for 10 min. The cells were then centrifuged at 300×g at 4 °C for 5 min, and the supernatant was removed. The cells were resuspended in 90% of cold methanol and incubated on ice for 30 min. The cells were resuspended in FACS buffer to make 106 cells/100 μL. LTL-FITC (FL-1321, Vector Laboratories) was added at a concentration of 1:100, and the mixture was incubated for 30 min at 4 °C with a rotator. The cell suspension was centrifuged at 300×g at 4 °C for 5 min, the supernatant was removed and resuspended in FACS buffer and passed through a 35 µm filtered tube (352235, Falcon). The cell suspensions were measured using a flow cytometer (NovoCyte, Agilent). The data were analyzed in Flowjo (Ver 10.6.2). The gating strategy for LTL-FITC was performed as follows: the population with the lowest fluorescence was gated to be negative, and all experiments were performed using the same gating. The gating strategy for FSC/SSC is shown in Supplementary Fig. 7.

Quantification of γH2AX-positive cells in proximal tubules of kidney organoids treated with cisplatin

Stained kidney organoids were observed for fluorescence using a CQ1 (Yokogawa). Imaging was performed under the following conditions to capture the entire organoid. Images were taken with an X10 lens, and after acquiring 3 × 3 tiled images, images were taken every 10 µm as a Z-stack. Maximum intensity projection (MIP) was applied to LTL, γH2AX and DAPI fluorescence channels to generate 2D images from the 3D datasets. Data analysis was performed using CellPathfinder (Yokogawa), and the number of γH2AX-positive cells was calculated by dividing the total number of γH2AX-positive cells present within the LTL-positive region of the organoid by the total area (µm2) of the LTL-positive region. These values were then normalized to the average value of the cisplatin-untreated group, which was set as 1.

Statistics and reproducibility

Statistical significance was determined by the unpaired two-tailed Student’s t-test, Welch’s t-test or paired t-test using SAS software for Windows, release 9.4 (SAS Institute Japan Ltd.). The mean, standard error, and statistical analysis results are shown in Supplementary Data 2. The number of replicates is defined in each figure legend.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.