Markers of disease and steroid responsiveness in paediatric idiopathic nephrotic syndrome: Whole-transcriptome sequencing of peripheral blood mononuclear cells.

Objective To identify markers of disease and steroid responsiveness in paediatric idiopathic nephrotic syndrome. Methods Whole-transcriptome sequencing was performed of peripheral blood mononuclear cells (PBMCs) from patients with NS. Differentially expressed genes (DEGs) were identified in patients with active NS vs those in remission, and those with steroid-sensitive NS (SSNS) vs steroid-resistant NS (SRNS). Results A total of 1065 DEGs were identified in patients with NS ( n = 10) vs those in remission ( n = 9). These DEGs correlated with cytokine and/or immune system signalling and the extracellular matrix. Comparisons between SSNS ( n = 6) and SRNS ( n = 4) identified 1890 DEGs. These markers of steroid responsiveness were enriched with genes related to the cell cycle, targets of microRNAs, and genes related to cytokines. Conclusions Meaningful DEGs were identified. Additional studies with larger numbers of patients will provide more comprehensive data.

Introduction

The first-line therapy for children with idiopathic nephrotic syndrome (NS) is steroid treatment, which induces remission in most patients.[1-3] The main clinical problems associated with steroid-sensitive NS (SSNS) are frequent relapse and subsequent drug toxicity.[4] Patients with steroid-resistant NS (SRNS) who do not respond to steroids and other treatments are at risk of the deterioration of renal function leading to end-stage renal disease.[5,6] Both SSNS and SRNS are associated with effacement of glomerular epithelial cell (podocyte) foot processes, a cardinal morphological feature of NS.[7] The aetiology of podocytopathy resulting in NS, reasons for steroid non-responsiveness, and the mechanisms underlying relapse in SSNS remain to be fully established.[8]

It has been speculated that the pathophysiology of SSNS involves disturbance of the immune system, especially T cells. This speculation is based on findings including the association between NS and lymphoma in some cases, relapse coinciding with infection, response to various immunosuppressive medications, and imbalances of a subpopulation of lymphocytes.[3,9-12] Several cytokines and other soluble plasma components may also be associated with NS.[13,14] A case has been described in which SSNS disappeared after bone-marrow transplantation,[10] suggesting that hematopoietic cells are involved in the pathogenesis of SSNS. SRNS has been shown to recur after kidney transplantation in some patients, suggesting that the pathogenesis of this condition resides outside the kidney; in addition, the efficacy of plasmapheresis in most recurrent cases indicates the presence of circulating factor(s) that cause SRNS.[15,16] However, these contributing factors remain to be identified and validated.[17-20] Although immunosuppressive agents are effective in some patients with SRNS, there are currently no tools to determine the optimal treatment for a patient before a therapeutic trial, or for predicting recurrence after kidney transplantation.[15,21]

Comprehensive information regarding NS would lead to a better understanding of the pathogenesis of the disease, mechanism of relapse, optimal medication choice, and prediction of prognosis. Thus, the present study applied whole-transcriptome sequencing of peripheral blood mononuclear cells (PBMCs) from patients with NS, using a next-generation sequencing (NGS) method of RNA sequencing.[22] PBMCs were used because of the high probability of immune system involvement in the pathogenesis of NS, and their easy accessibility, a prerequisite for a useful biomarker.[23,24] Compared with microarray technologies, RNA sequencing can capture the dynamic range of transcriptomes in terms of both expression profiling and differentially expressed isoforms (DEIs) on a massive scale.[25-27] We report the preliminary results of signature gene sets of NS and steroid responsiveness.

Patients and methods

Study population

The study recruited children aged <18 years who were newly diagnosed with idiopathic NS at Seoul National University Children’s Hospital, Seoul, Republic of Korea, between January 2008 and December 2011. Patients who were on long-term treatment prior to transfer to our hospital were excluded from the study. Pathological diagnosis was obtained only in patients with SRNS.

The study was approved by the Seoul National University Hospital Institutional Review Board (No. 0812-002-264), and the participants’ parents or legal guardians provided written informed consent prior to enrolment.

Sample collection

Peripheral blood samples were collected from patients and PBMCs were isolated using Ficoll-Hypaque density gradient centrifugation, then stored at −80℃ until RNA extraction. Nephrotic samples were collected at the time of onset or relapse of NS, before commencing any treatment. Remission samples were collected from patients with SSNS during remission, when having not been taking steroids for >2 months.

Whole-transcriptome sequencing

Total RNA was extracted from PBMCs using a QIAamp RNA mini kit (Qiagen, Austin, TX, USA). Libraries were prepared based on the Illumina protocol according to the manufacturer’s instructions, and 54 bp of paired-end RNA sequencing data were generated using the Illumina Genome Analyzer IIx (Illumina, San Diego, CA, USA). The prepared libraries were quantified using quantitative polymerase chain reaction (PCR) according to the quantification protocol guide in the manufacturer’s instructions. The read quality was checked, then the differentially expressed gene (DEG) sequences were identified using R package DEGseq (version 1.10.0),[28] through counting the reads and assessing the distribution of count differences between samples. Raw read quality scores and read counts were summarized.

For annotation, RNA sequence reads were aligned to the human reference genome (University of California, Santa Cruz [UCSC] hg19; 20 October 2011) using TopHat software (version 1.4.0)[30] and Bowtie software (version 1.12.5),[31] with the supplied annotations, a set of gene-model annotations and known transcripts, and the –no-novel-juncs option to disable mapping for novel splice junctions.[29-31] The aligned reads were quantified with Cufflinks (version 1.3.0) to obtain the fragments per kilobase of exons per million fragments mapped (FPKM) values for the genes or gene transcripts, and then merged into an expression table for the next analysis step, outlined in Figure 1 and conducted as described.[31]

Figure 1.

Workflow of the RNA sequencing data analysis in a study investigating disease markers of paediatric idiopathic nephrotic syndrome (NS) and steroid responsiveness. First, a pipeline was built to identify differentially expressed genes (DEGs) based on mRNA expression levels. Functional annotations were applied to the DEGs, including pathway enrichment analysis, functional annotation clustering, and gene set enrichment analysis.

Expression profiling and functional annotation

The average number of reads produced from each sample was 74 million. Only those of protein coding genes listed on the UCSC Genome Browser[32] were analysed. Loci with low variance in FPKM values or zero reads across all samples were removed. Variance-stabilizing normalization and upper-quartile normalization were applied to the boost sensitivity without a loss of specificity.[33]

The DEGs were obtained from one-way analyses of variance (ANOVA) for each group, and false discovery rate (FDR) multiple testing corrections were applied. Post-hoc analyses were performed to detect the relationships between groups via the Tukey’s honest significance test. Analyses of DEIs were performed similarly, but no significant DEIs were obtained.

The DEGs of the groups of interest were obtained by t-tests. For functional annotation and clustering, the Gene Set Enrichment Analysis (GSEA) program (version 2.0.8) with the Molecular Signatures Database (version 3.1)[34] and the Database for Annotation, Visualization and Integrated Discovery (DAVID, version 6.7) were used to enhance understanding of the underlying biological relevance.[35,36] Clustering analysis was performed using the kmeans function in R 3.0.2, which performs k-means clustering (K = 10 clusters specified) on a given expression profile for DEGs. The hypergeometric distribution is used to compute P-values for Gene Ontology (GO) annotation for clusters with the Molecular Signatures Database (version 5.1).[34] For upstream analysis of DEGs, gene-sets of microRNA targets (n = 221) and transcription factor targets (n = 615) from Molecular Signatures Database (version 5.1) were downloaded and compared with DEGs.[34]

Results

In total, 18 patients with idiopathic NS were enrolled (15 males/3 females; mean age 8.2 ± 4.0 years; age range 2.7–16.7 years). The median age at onset of NS was 5.9 years (range 3.0–14.4 years). Nephrotic samples (n = 10) were obtained from six patients with SSNS and four with SRNS. Pathological diagnosis was obtained only in those with SRNS, and was focal segmental glomerulosclerosis in all cases. Of the four patients with SRNS, two responded to cyclosporine treatment (calcineurin inhibitor [CNI] responders [CRs]), and two responded to neither steroids nor CNI (nonresponders [NRs]). A total of nine remission samples were collected from patients with SSNS.

The gene expression profile was determined by analysing 19 samples from 18 patients (one patient provided both a nephrotic sample and a remission sample) and 18 551 genes. Statistical analyses identified 1065 DEGs in the NS group (n = 10) relative to the remission group (n = 9) (Figure 2). Functional annotations of these genes revealed that these DEGs were related to dorsal/ventral pattern formation (enrichment score [ES] 2.05), extracellular matrix structural constituents (ES 1.75), and actin binding (ES 1.36) according to the DAVID functional annotation module. Based on the GSEA, compared with the remission group, the gene-expression profile of the NS group was enriched with genes pertaining to steroid hormones, matrix metalloproteinase (i.e., enzymes that degrade the extracellular matrix)-inducing cytokines, extracellular matrix-receptor interaction, acyl chain remodelling of phosphatidylglycerol, G β:γ signalling through PI3Kγ, CTLA4 inhibitory signalling, the early response to TGFβ1, IL4 receptor signalling in B lymphocytes, pantothenate and CoA biosyntheses, the syndecan 3 pathway, and the mTOR signalling pathway.

Figure 2.

Principal component analysis of peripheral blood mononuclear cell whole-transcriptome sequencing data from children with nephrotic syndrome (NS, red dots) and those in remission (control group, green dots). Groups are segregated according to expression patterns in RNA sequencing, based on 1065 DEGs (P < 0.05).

More stringent criteria (P < 0.01 and >2-fold changes of expression) were applied to identify the highly significant genes in idiopathic NS. A total of 49 genes were found to be significantly upregulated in NS, and 67 genes were found to be downregulated (Table 1). K-means clustering for 116 DEGs revealed 10 clusters of 3–37 genes, with enriched GO terms listed in Table 2 (hypergeometric test, P < 0.005). Upstream analysis revealed that DEGs of NS were enriched with targets of MIR-370 (P = 0.0163, reported in Wilms tumour) and MIR-519E (P = 0.0428, clinical relevance not yet known), as well as targets of transcription factors ATF2, ATF6, EVI1, HMGA1, IRF8, ITGAL, JUN, MEF2A, NFAT, PGR, POU3F2, and STAT6.

Table 1.

Differentially expressed genes (DEGs) in paediatric idiopathic nephrotic syndrome (NS). (nephrotic status vs remission status; P < 0.01; relative change >2-fold).

Gene symbol	Official gene name	Relative change	Statistical significance
RWDD1	RWD domain-containing 1	6.3	P = 0.000073
IVD	Isovaleryl coenzyme A dehydrogenase	2.2	P = 0.00019
ZNF48	Zinc finger protein 48	6.6	P = 0.0002
FAM65B	Family with sequence similarity 65, member B	4.8	P = 0.00038
USP2	Ubiquitin-specific peptidase 2	3.0	P = 0.00038
DNMT3B	DNA (cytosine-5-)-methyltransferase 3 beta	3.4	P = 0.00038
DUSP23	Dual specificity phosphatase 23	4.6	P = 0.00110
ZNF229	Zinc finger protein 229	3.9	P = 0.0012
C11orf74	Chromosome 11 open reading frame 74	2.1	P = 0.0014
CAPN6	Calpain 6	2.0	P = 0.0018
LARP6	La ribonucleoprotein domain family, member 6	3.7	P = 0.0018
BIRC6	Baculoviral IAP repeat containing 6	2.2	P = 0.0019
TMEM134	Transmembrane protein 134	3.9	P = 0.0019
RPS15A	Ribosomal protein S15a pseudogene 17	2.3	P = 0.0021
ETV4	ETS translocation variant 4	5.6	P = 0.0022
SLFN5	Schlafen family member 5	2.2	P = 0.0023
ROR2	Receptor tyrosine kinase-like orphan receptor 2	6.5	P = 0.0024
GADD45G	Growth arrest and DNA-damage-inducible, gamma	7.8	P = 0.0025
HAUS4	HAUS augmin-like complex, subunit 4	2.3	P = 0.0025
SNTA1	Syntrophin, alpha 1 (dystrophin-associated protein A1, 59-kDa, acidic component)	2.2	P = 0.003
SIRT6	Sirtuin (silent mating type information regulation 2 homologue) 6 (Saccharomyces cerevisiae)	2.4	P = 0.0037
CD8A	CD8a molecule	3.6	P = 0.0038
OTOP2	Otopetrin 2	2.6	P = 0.0038
AMIGO1	Adhesion molecule with Ig-like domain 1	2.4	P = 0.004
C15orf48	Chromosome 15 open reading frame 48	2.1	P = 0.0041
MPP6	Membrane protein, palmitoylated 6 (MAGUK p55 subfamily member 6)	5.4	P = 0.0044
EPT1	Selenoprotein I	3.7	P = 0.0046
TTLL12	Tubulin tyrosine ligase-like family, member 12	3.5	P = 0.0048
WDR5	WD-repeat domain 5	2.3	P = 0.0049
IL1RAP	Interleukin 1 receptor accessory protein	2.2	P = 0.0053
WDR27	WD-repeat domain 27	2.4	P = 0.0054
CYB5B	Cytochrome b5 type B (outer mitochondrial membrane)	2.9	P = 0.0054
PTK2	PTK2 protein tyrosine kinase 2	4.3	P = 0.0055
LAMA4	Laminin, alpha 4	3.7	P = 0.0055
SGK223	Homologue of rat pragma of Rnd2	2.3	P = 0.0058
DAZAP2	DAZ-associated protein 2	2.0	P = 0.0058
LAX1	Lymphocyte transmembrane adaptor 1	4.3	P = 0.0059
C17orf100	Chromosome 17 open reading frame 100	2.0	P = 0.006
PEBP1	Phosphatidylethanolamine-binding protein 1	3.8	P = 0.006
BAI1	Brain-specific angiogenesis inhibitor 1	5.6	P = 0.0062
RMND5A	Required for meiotic nuclear division 5 homologue A (S. cerevisiae)	4.6	P = 0.0064
OSCP1	Chromosome 1 open reading frame 102	3.2	P = 0.0065
WNT5A	Wingless-type MMTV integration site family, member 5A	9.9	P = 0.0067
LMX1B	LIM homeobox transcription factor 1, beta	2.3	P = 0.0075
TYMP	Thymidine phosphorylase	4.6	P = 0.0075
HIST1H2BN	Histone cluster 1, H2bn	2.8	P = 0.00780
MYH2	Myosin, heavy chain 2, skeletal muscle, adult	2.8	P = 0.0087
IFNA5	Interferon, alpha 5	4.3	P = 0.0096
MMP24	Matrix metallopeptidase 24 (membrane-inserted)	2.4	P = 0.0099
NDST2	N-deacetylase/N-sulfotransferase (heparanglucosaminyl) 2	−2.1	P = 0.00332
ZNF670	Zinc finger protein 670	−2.1	P = 0.00045
VNN1	Vanin 1	−2.1	P = 0.000231
CTIF	CBP80/20-dependent translation initiation factor	−2.2	P = 0.000469
NEDD1	Neural precursor cell expressed, developmentally down-regulated 1	−2.2	P = 0.00907
ALKBH8	AlkB, alkylation repair homolog 8 (Escherichia coli)	−2.3	P = 0.00594
CDKN3	Cyclin-dependent kinase inhibitor 3	−2.3	P = 0.00890
TNNC2	Troponin C type 2 (fast)	−2.3	P = 0.00225
IL21	Interleukin 21	−2.4	P = 0.00711
LCN2	Lipocalin 2	−2.4	P = 0.0055
PELI3	Pellino homologue 3 (Drosophila)	−2.4	P = 0.00426
DIRAS3	DIRAS family, GTP-binding RAS-like 3	−2.5	P = 0.00238
PRKD2	Protein kinase D2	−2.5	P = 0.00909
SLC26A9	Solute-carrier family 26, member 9	−2.7	P = 0.00517
AGXT2 L1	Alanine-glyoxylate aminotransferase 2-like 1	−2.8	P = 0.00712
TMEM81	Transmembrane protein 81	−2.8	P = 0.0023
REEP2	Receptor accessory protein 2	−3.0	P = 0.00516
WNT8B	Wingless-type MMTV integration site family, member 8B	−3.1	P = 0.00381
SPZ1	Spermatogenic leucine zipper 1	−3.1	P = 0.00892
HLCS	Holocarboxylasesynthetase (biotin-[proprionyl-coenzyme A-carboxylase (ATP-hydrolysing)] ligase)	−3.2	P = 0.0000898
SGK2	Serum/glucocorticoid-regulated kinase 2	−3.3	P = 0.00829
SLC10A6	Solute carrier family 10 (sodium/bile acid cotransporter family), member 6	−3.4	P = 0.00696
ERVV-2	Endogenous retrovirus group V, member 2	−3.4	P = 0.00302
OPLAH	5-Oxoprolinase (ATP-hydrolysing)	−3.4	P = 0.00323
PXMP4	Peroxisomal membrane protein 4, 24 kDa	−3.5	P = 0.00689
ANO10	Anoctamin 10	−3.5	P = 0.00627
ST8SIA1	ST8 alpha-N-acetyl-neuraminide alpha-2,8-sialyltransferase 1	−3.6	P = 0.00548
RAD51AP2	RAD51-associated protein 2	−3.6	P = 0.00591
PRPF38B	PRP38 pre-mRNA processing factor 38 (yeast) domain containing B	−3.8	P = 0.00388
RFPL4B	Ret finger protein-like 4B	−4.0	P = 0.007
PSRC1	Proline/serine-rich coiled-coil 1	−4.1	P = 0.00206
FAM55C	Family with sequence similarity 55, member C	−4.5	P = 0.00554
TRIML2	Tripartite motif family-like 2	−4.5	P = 0.00399
MATN1	Matrilin 1, cartilage matrix protein	−4.5	P = 0.00207
MED29	Mediator complex subunit 29	−4.6	P = 0.00543
STH	Saitohin	−4.6	P = 0.00238
CCDC64	Coiled-coil domain containing 64	−4.6	P = 0.00665
KIAA0922	KIAA0922	−4.9	P = 0.00479
LY75	CD302 molecule; lymphocyte antigen 75	−4.9	P = 0.00985
CYP11B1	Cytochrome P450, family 11, subfamily B, polypeptide 1	−5.3	P = 0.00456
TSC2	Tuberous sclerosis 2	−5.6	P = 0.00627
CASQ2	Calsequestrin 2 (cardiac muscle)	−5.7	P = 0.00237
MAGEA9	Melanoma antigen family A, 9; melanoma antigen family A, 9B	−5.8	P = 0.00104
ZNF358	Zinc finger protein 358	−5.9	P = 0.00897
TBX20	T-box 20	−6.1	P = 0.00765
CYP2C9	Cytochrome P450, family 2, subfamily C, polypeptide 9	−6.3	P = 0.00386
FYN	FYN oncogene related to SRC, FGR, YES	−6.3	P = 0.00701
LRRIQ4	Leucine-rich repeats and IQ motif-containing 4	−6.4	P = 0.00827
CLOCK	Clock homolog (mouse)	−6.6	P = 0.00433
SLC7A10	Solute carrier family 7 (neutral amino acid transporter, y + system), member 10	−6.8	P = 0.00315
HOXA11	Homeobox A11	−7.0	P = 0.00935
SPINT4	Serine peptidase inhibitor, Kunitz type 4	−7.1	P = 0.00225
PAGE4	P antigen family, member 4 (prostate-associated)	−7.2	P = 0.00475
OPTC	Opticin	−7.7	P = 0.00442
CEACAM20	Carcinoembryonic antigen-related cell adhesion molecule 20	−7.8	P = 0.00209
ZFHX2	Zinc finger homeobox 2	−8.1	P = 0.0068
CHRD	Chordin	−8.5	P = 0.00949
TGM7	Transglutaminase 7	−8.5	P = 0.00491
CHIA	Chitinase, acidic	−8.9	P = 0.00159
SVIP	Small VCP/p97-interacting protein	−9.0	P = 0.00863
AOX1	Aldehyde oxidase 1	−9.5	P = 0.00905
FAM75A6	Family with sequence similarity 75, member A6	−9.7	P = 0.00053
LIPI	Lipase, member I	−10.1	P = 0.00344
KCNA10	Potassium voltage-gated channel, shaker-related subfamily, member 10	−10.7	P = 0.00366
SLC34A2	Solute carrier family 34 (sodium phosphate), member 2	−11.4	P = 0.0048
LRP1B	Low-density lipoprotein-related protein 1B (deleted in tumours)	−11.6	P = 0.00936
S100A9	S100 calcium-binding protein A9	−21.0	P = 0.000391

Table 2.

Enriched Gene Ontology (GO) terms from K-means clustering of differentially expressed genes (DEGs) in paediatric idiopathic nephrotic syndrome (NS). Clusters listed based on hypergeometric test of P < 0.005.

Cluster	Genes	Enriched GO term(s)	Statistical significance
Cluster 1, n = 4	AMIGO1, CYB5B, HIST1H2BN, MMP24	Heterophilic cell adhesion	P = 0.0022
		Homophilic cell adhesion	P = 0.0034
		Regulation of action potential	P = 0.0037
		Cell recognition	P = 0.0041
Cluster 2, n = 37	ALKBH8, CCDC64, CDKN3, DAZAP2, DUSP23, EPT1, FAM65B, HAUS4, HLCS, IVD, LCN2, LMX1B, MATN1, MPP6, NDST2, NEDD1, OPLAH, OPTC, OSCP1, OTOP2, PELI3, PSRC1, PXMP4, REEP2, RPS15A, SIRT6, SLFN5, SNTA1, ST8SIA1, TMEM81, TNNC2, TYMP, VNN1, WNT8B, ZFHX2, ZNF48, ZNF670	Transferase activity transferring pentosyl groups	P = 0.0007
		Extracellular matrix structural constituent	P = 0.0013
Cluster 3, n = 4	CASQ2, FAM75A6, MED29, STH	Striated muscle contraction	P = 0.0030
Cluster 4, n = 23	ANO10, BAI1, BIRC6, C17orf100, CAPN6, CTIF, DNMT3B, ETV4, FAM55C, GADD45G, IFNA5, IL1RAP, LARP6, LAX1, MYH2, PEBP1, PTK2, RMND5A, ROR2, SGK223, TMEM134, USP2, WDR27	SH2 domain binding	P = 0.0002
		Cysteine type endopeptidase activity	P = 0.0011
Cluster 7, n = 5	AOX1, CHIA, CHRD, CLOCK, LY75	Inflammatory response	P = 0.0005
		N acetylglucosaminemetabolic process	P = 0.0032
		Cellular polysaccharide metabolic process	P = 0.0043
Cluster 8, n = 9	CEACAM20, FYN, LAMA4, S100A9, SLC7A10, SPZ1, SVIP, TGM7, TSC2	Neutral amino acid transport	P = 0.0048

Gene expression patterns differed significantly between SSNS and SRNS (Figure 3), with 1890 DEGs identified (P < 0.1). These DEGs were enriched with genes related to the microtubule organizing centre and regulation of the response to biotic stimuli based on the GO terms. Based on the GSEA, compared with the SRNS group, the gene expression profile of the SSNS group was enriched with genes pertaining to TGFβ1 signalling, the cell cycle and p53 signalling, Y branching of actin filaments, FoxP3 targets in T lymphocytes, cytokines IL6 and IL4, and targets of MIR106B (related to renal cell carcinoma[37]) and MIR16 (expressed in the kidneys[38]).

Figure 3.

Principal component analysis and heat map of peripheral blood mononuclear cell whole-transcriptome sequencing data from children with steroid sensitive nephrotic syndrome (SNNS, green dots) and steroid resistant NS (SRNS, red dots). Groups are segregated according to expression patterns in RNA sequencing, based on 1890 DEGs (P < 0.1).

More stringent criteria (P < 0.01 and >2-fold changes of expression) were applied to identify the markers of steroid responsiveness. Consequently, 23 genes were selected (Table 3; enriched GO terms per k-means clustering Table 4). Upstream analysis did not reveal any significant findings.

Table 3.

Differentially expressed genes (DEGs) in paediatric patients with steroid sensitive idiopathic nephrotic syndrome (SSNS) or steroid resistant nephrotic syndrome (SRNS) (P < 0.01; relative change >2-fold).

Gene symbol	Official gene name	Relative change	Statistical significance
TCF4	Transcription factor 4	11.3	P = 0.0067
BMPR1B	Bone morphogenetic protein receptor, type IB	8.5	P = 0.008
LOC255411	Hypothetical LOC255411	6.7	P = 0.0075
ITGA4	Integrin, alpha 4 (antigen CD49D, alpha 4 subunit of VLA-4 receptor)	6.4	P = 0.00017
XKR6	XK, Kell blood group complex subunit-related family, member 6	6.3	P = 0.0041
PLA2G5	Phospholipase A2, group V	6.1	P = 0.00017
DSG4	Desmoglein 4	5.5	P = 0.004
ODZ4	Teneurin transmembrane protein 4	5.4	P = 0.00011
GAGE2D	G antigen 2A	5.2	P = 0.004
KRTAP5-10	Keratin-associated protein 5-10	5.1	P = 0.0038
OC90	Otoconin 90	4.8	P = 0.0074
GPC1	Glypican 1	4.7	P = 0.0042
FAM48B2	Family with sequence similarity 48, member B2	4.5	P = 0.0041
MSRB3	Methionine sulphoxide reductase B3	4.3	P = 0.0065
ZXDC	ZXD family zinc finger C	2.9	P = 0.000052
ZNF566	Zinc finger protein 566	2.3	P = 0.0051
ZNF251	Zinc finger protein 251	−2.0	P = 0.00085
SAP30BP	SAP30-binding protein	−2.3	P = 0.003
LGALS4	Lectin, galactoside-binding, soluble, 4	−2.8	P = 0.0042
RAD9B	RAD9 homologue B (Schizosaccharomyces pombe)	−3.6	P = 0.0012
ATCAY	Ataxia, cerebellar, Cayman type	−3.7	P = 0.0049
PDF	Peptide deformylase	−4.0	P = 0.0048

Table 4.

Enriched Gene Ontology (GO) terms from K-means clustering of differentially expressed genes (DEGs) in steroid sensitive paediatric idiopathic nephrotic syndrome (SSNS). Clusters listed based on hypergeometric test of P < 0.005.

Cluster	Genes	Enriched GO term(s)	Statistical significance
Cluster 1, n = 6	BMPR1B, ITGA4, LOC255411, ODZ4, PLA2G5, TCF4	SMAD binding	P = 0.0039
		Phospholipase A2 activity	P = 0.0042
Cluster 3, n = 3	MSRB3, ZNF566, ZXDC	Oxidoreductase activity acting on sulphur group of donors	P = 0.0007
Cluster 5, n = 6	ATCAY, LGALS4, PDF, RAD9B, SAP30BP, ZNF251	N-terminal protein amino acid modification	P = 0.0036

Discussion

This study used whole-transcriptome sequencing to identify genes that differed in expression in children with idiopathic NS in remission or with nephrotic status. Analysis using t-testing with P < 0.05 revealed 1065 DEGs for NS independent of steroid responsiveness. These DEGs were enriched with extracellular matrix structural constituent/actin binding/cytoskeletal protein binding according to the GO term of molecular function, as well as cytokine and/or immune system signalling related to steroids; CTLA4, TGFβ1, IL4, and mTOR according to GSEA. IL4 is a representative cytokine of Th2 immune reactions, and Th2 immune reactions have been reported to be predominantly associated with childhood NS.[11] Additionally, CTLA4 and TGFβ1 are related to immune regulation, and impaired regulatory T cell function has been reported in idiopathic NS.[39] Upstream analysis showed that DEGs of NS were enriched with targets of MIR-370, which is related to Wilms tumour of the kidneys, suggesting relevance of DEGs affecting the kidneys.[40] Furthermore, among 12 upstream genes, ITGAL, MEF2A, STAT6 are members of steroid responsiveness panel genes in U.S. patents.[41] Therefore, the findings of the present study generally agree with knowledge regarding NS. Further refinement of these results in larger studies will improve our understanding of NS.

Although steroid treatment is the first-line treatment for children with NS, it is associated with significant toxicity.[4] For patients who do not respond to steroid treatment, initial treatment with steroids could be harmful as well as ineffective. Moreover, more aggressive treatments, such as CNI, rituximab and plasmapheresis, could induce remission in many patients if instituted without delay, as seen in recurrent SRNS after kidney transplantation.[15,21] Therefore, the identification of reliable markers for steroid responsiveness would allow more directed treatment of paediatric NS. Patients who are nonresponsive to steroids could be other treatment options without delay. In search of markers for steroid responsiveness in paediatric NS, we identified a total of 1890 DEGs, and selected 23 genes based on more stringent criteria. Interestingly, the DEGs of patients with SSNS (vs SRNS) were enriched in genes pertaining to the cell cycle and the targets of microRNAs MIR106B and MIR16, in addition to those related to cytokines. The emergence of cell cycle-related genes may imply differences in the proliferative properties of SSNS and SRNS, which could be utilized for the development of novel therapeutic options. The 23 genes that were selected as markers of steroid responsiveness seem heterogeneous, but following refining with different sets of samples for validation, list of genes or part of this list can be used as markers of steroid responsiveness. Interestingly, comparison of the signature genes of SSNS with those listed as SSNS in the patent for “Kit and method for identifying individual responsiveness to steroid therapy of nephrotic syndrome”[41] did not reveal any common genes, despite the similarity of the methods, indicating that clinical utilization of this approach requires further study. Notably, previously proposed circulating factors indicative of SRNS (cardiotrophin-like cytokine factor 118 and urokinase-type plasminogen activator receptor[17]) were not found among the DEGs in the present study, possibly due to the heterogeneous nature of our study population. These proposed circulating factors were discovered in patients with recurrent NS after kidney transplantation in which steroid treatment can achieve remission in the majority of patients.

The present study has several shortcomings. First, the sample size was small, limiting the statistical power. Additionally, some relevant DEGs may not have been identified due to this small sample size. The DEGs identified in this study were able to clearly classify the groups, so our approach seems valid and justifies further studies to identify disease/therapeutic response markers for clinical applications. Secondly, although RNA sequencing was used rather than mRNA microarrays, DEIs and alternative splicing pattern differences between groups were not identified. To discover novel splice sites and rare transcripts, deep sequencing of at least 100 million reads of 76 bp in length is required (according to the guidelines of the Encyclopaedia of DNA Elements Project[42]). The insufficient number of reads of this study (mean 77 million reads with up to 75% of the reads properly aligned against a reference genome) could be the reason for the failure in the DEI search, in addition to the small number of samples per group. Finally, the validation of candidate markers of NS or steroid responsiveness was not performed in this study. Clearly, many of the DEGs are not linked to pathogenesis but rather are the results of or surrogate changes due to disease. A validation study may be helpful in discriminating these differences.

In conclusion, whole-transcriptome sequencing of PBMCs found that DEGs of NS were enriched in immune system signalling, and potential therapeutic targets were suggested. Further studies with larger numbers of patients will provide more comprehensive information to enable the application of precision medicine to paediatric NS.