Identiﬁcation and Characterization of the Wilms Tumor Cancer Stem Cell

A nephrogenic progenitor cell (NP) with cancer stem cell characteristics driving Wilms tumor (WT) using spatial transcriptomics, bulk and single cell RNA sequencing, and complementary in vitro and transplantation experiments is identiﬁed and characterized. NP from WT samples with NP from the developing human kidney is compared. Cells expressing SIX2 and CITED1 fulﬁll cancer stem cell criteria by reliably recapitulating WT in transplantation studies. It is shown that self-renewal versus diﬀerentiation in SIX2 + CITED1 + cells is regulated by the interplay between integrins ITG 𝜷 1 and ITG 𝜷 4. The spatial transcriptomic analysis deﬁnes gene expression maps of SIX2 + CITED1 + cells in WT samples and identiﬁes the interactive gene networks involved in WT development. These studies deﬁne SIX2 + CITED1 + cells as the nephrogenic-like cancer stem cells of WT and points to the renal developmental transcriptome changes as a possible driver in regulating WT formation and progression.


Introduction
Normal human nephrogenesis involves the reciprocal interaction between two embryonic layers: the branching ureteric bud (UB) and the surrounding cap mesenchyme (CM). [1] Induction by the UB initiates condensation of the CM and the beginning of further renal development characterized by the formation of peritubular aggregates, followed by renal vesicles, Cshaped and S-shaped bodies, leading to full maturation into a functional nephron. [1] Many studies confirm the presence of selfrenewing, uncommitted nephrogenic progenitors (NP) in a specific subdomain of the CM, co-expressing the transcription factors SIX2 and CITED1, the master genes regulating nephrogenesis. [2,3] Importantly, proper nephrogenesis is regulated by these uncommitted progenitors expressing SIX2 and CITED1 and by the committed progenitors that lose CITED1 but still maintain SIX2 expression. [1][2][3] Loss of CITED1 primes the cells to renal differentiation; thus, the tight regulation of SIX2 and CITED1 expression is critical during normal kidney development. In humans, kidney development ceases around 34-36 weeks gestational age (WGA) [4] and these NP are absent in the mature postnatal kidney.
Growing evidence links Wilms tumor (WT), which accounts for 95% of all pediatric kidney cancers, [5] to aberrant nephrogenesis, in which normal prenatal depletion of NP does not occur, and NP persists in the neonatal kidney. Histologically, WT resembles an embryonic kidney, with tumors containing multiple developing renal structures. [6,7] All three lineages of the developing kidney (blastema, epithelium, and stroma) can be identified in classic triphasic WT histology. [8] The combination of surgery, radiation, and chemotherapy has increased the survival rate for many WT patients. [9] However, patients with relapsed disease or initial unfavorable histopathology have a poor prognosis; thus, a better understanding of molecular mechanisms that regulate tumor development and progression is necessary to improve treatments for these WTs. [9] Many studies have focused on characterizing the complicated genetic causes of WT, and ≈40 genes that drive diversity in the genetic landscape of WT have been identified. [10] Even if many of these WT genes, including WT1 and -catenin, are known to play important roles in NP of the developing kidney, [11] and the presence of cells characterized by the expression of SIX2 and CITED1 has been described in WT samples, [12][13][14][15] very few studies have compared NP of WT to NP of the human fetal kidney (hFK) to understand molecular changes responsible for the dysregulation of normal nephrogenesis in WT.
Based on the premise that NP are central to WT development, we isolated and characterized, for the first time, a pure population of live NP co-expressing SIX2 and CITED1 from WT and normal hFK. We demonstrated that SIX2+CITED1+ cells from WT exhibited nephrogenic cancer stem cell (CSC) traits and, upon transplantation into NOD/SCID mice, formed consecutive xenografts with histological features similar to the WT of origin, as well as metastatic capabilities. By comparing SIX2+CITED1+ NP from different WT samples versus hFK, we showed that SIX2+CITED1+ cells from all these sources displayed NP characteristics, but important expression-level and phenotypic differences distinguished WT and hFK-derived cells. Trajectory inference analysis generated from scRNA-seq data integration of SIX2+CITED1+ cells and cells dissociated from the WT tissue (from which SIX2+CITED1+ cells were isolated) and spatiotemporal mapping suggested that SIX2+CITED1+ cells drive WT development. We further showed that the interaction between these cells and the surrounding extracellular matrix (ECM) through integrin signaling influences NP self-renewal and renal specification. Using spatial transcriptomic analysis, we defined specific gene map networks of the WT SIX2+CITED1+ cells and correlated these with the histological setting, thus analyzing the role of the microenvironment in defining their molecular and cellular properties. Our studies identified critical differences within the nephrogenic population between WT and hFK and highlighted the disruption of specific signaling pathways driving WT development.

Isolation and Characterization of SIX2+CITED1+ Cells in WT and hFK
Histologic analysis was used to determine morphological differences and structural organization of developing hFK at different WGA and WT samples ( Table 1). Nephrogenesis through gestation can be followed in distinct renal compartments surrounded by organized stroma (Figure 1A, Figure S1A-E, Supporting Information). By 10 WGA, in the hFK, all nephrogenic structures are identifiable, including the nephrogenic niche, mature glomeruli, and tubules. WT, in contrast, appeared histologically disorganized, without recognizable elements of normal renal architecture, and with structural heterogeneity based on the WT subtype ( Figure 1B, Figures S1F-J and S2A-I, Supporting Information).
In hFK, SIX2+CITED1+ cells are found exclusively in the nephrogenic zone in CM near the branching UB ( Figure 1C). As we have previously reported, [2] CITED1 expression is not detected in renal vesicles, C-shape or S-shape bodies, where some cells still express SIX2. Mature hFK glomeruli and tubules also do not contain SIX2+CITED1+ cells. The distribution of SIX2+CITED1+ cells in WT depends on WT subtype and differs from hFK. These cells are not restricted to a developmental niche in WT, but rather are dispersed in multiple "blastema foci" and around abortive glomeruli and tubules ( Figure 1D; Figure S2A-I, Supporting Information). SIX2+CITED1+ cells were detected in all WT samples included in this study. Using our validated Smart-Flare method for live cell sorting [2,16] (Figure 1E, Figure S3A-E, Supporting Information), we isolated SIX2+CITED1+ cells from hFK and WT samples. We established that the abundance of SIX2+CITED1+ cells was the same in hFK samples of similar WGA but was highly variable among WT samples ( Figure 1F and Figure S3F-H, Supporting Information), likely dependent on the number of blastema foci in each WT.

SIX2+CITED1+ Cells from WT Meet the criteria for Cancer Stem Cells (CSC)
CSCs are capable of self-renewal, differentiation, and tumorigenicity in animal hosts and are defined by a specific set of criteria: small numbers of transplanted cells should generate a tumor in vivo, cells should be resistant to chemotherapeutics, and can be identified by specific markers reflecting the tumor of origin. [27][28][29][30] WT SIX2+CITED1+ cells tumor formation capabilities were assessed through a series of in vivo experiments ( Figure 2C and Figure S5, Supporting Information). Freshly isolated SIX2+CITED1+ cells from different WT samples (Figure 2D) that were injected subcutaneously and intrarenally into NOD/SCID mice generated tumors ( Figure 2E, upper panel, and Figure S6, Supporting Information). Limiting dilution experiments also demonstrated the capability of xenograft formation even at low concentrations if injected intrarenally ( Figure  S6A, Supporting Information). Importantly, culture-expanded WT SIX2+CITED1+ cells (6-12 passages) injected subcutaneously also generated tumors in vivo ( Figure 2E, lower panel). Serial transplantation experiments (2 nd and 3 rd xenograft generation) also demonstrated the capability of these cells to generate in vivo WT ( Figure 2C, Figure S6B, Supporting Information). Transplanted WT SIX2-CITED1-cells and WT SIX2+CITED1-did not generate tumors within 30 weeks (not shown). These data indicate that the cells responsible for tumor formation in our experimental conditions are the SIX2+CITED1+ cells.
Tumors generated by subcutaneous and intrarenal injection of freshly isolated or cultured WT SIX2+CITED1+ cells had classic WT morphology, including blastema, stroma, and epithelial structures containing SIX2+CITED1-cells and SIX2+CITED1+ cells ( Figure 2E, Figure S6A,B, Supporting Information). Injected WT#8 SIX2+CITED1+ cells, which had deletions in 7p and 11q (common in relapsed patients), also generated a WT that metastasized to the liver detected 6 months after injection ( Figure S6C, Supporting Information). SIX2+CITED1+ cells from hFK did not form tumors, but in one case, the cells proliferated and generated a "mass" (<0.2 cm at 5 months) without histologic features of WT and no expression of SIX2 and/or CITED1 ( Figure S6D, Supporting Information).  favorable stage III) show the nephrogenic zone (white dotted line) and differentiating structures (second panel: ureteric bud, UB; cap mesenchyme CM; tubule, glomerulus, and stroma) of hFK, and unorganized WT histology with triphasic components (second panel, stroma, blastema, and epithelial structures including abortive glomeruli and tubules). 10× images acquired and composed using Photoshop DC (Adobe) for whole images, right panels of 20X images. C,D) SIX2 (red) and CITED1 (green) immunofluorescence staining of C) hFK 10 WGA and D) WT#4. SIX2+CITED1+ co-expression in hFK (C, second panel) in the nephrogenic niche (uninduced cap mesenchyme, UCM) but absent within developing (renal vesicle, C-shape, S-shape) and mature (glomerulus and tubule) structures. SIX2+CITED1+ expression is dispersed throughout the WT (D, second panels) in blastema but not in stroma or abortive structures (glomerulus and tubule). Nuclei stained with DAPI (blue), 10× images acquired and composed using Photoshop DC (Adobe) for whole images, right panels of 20× images. E) SmartFlare technique validation by flow cytometry. SIX2-Cy5 and CITED1-Cy3 probes (top left and right panel respectively) were individually used to isolate cells from hFK (  Histologically, the "mass" appeared to contain primitive normal tubular structures. Some post-transplantation proliferation of hFK NPs is expected since the cells have intrinsic proliferative potential [31,32] before committing to differentiation. To evaluate resistance to chemotherapeutic drugs, we treated NOD/SCID mice injected with WT SIX2+CITED1+ cells with vincristine, commonly used clinically for WT. [33][34][35][36] Standard vincristine dosing protocols for mice [36] did not delay or prevent tumor growth ( Figure 2F), suggesting that WT SIX2+CITED1+ cells also manifest chemotherapy resistance. We recognize that different dosing regimens [37] of vincristine might result in variations in response rates. Here we showed proof of the principle that xenografts generated from SIX1+CITED1+ cells (and not from WT tissue, as reported in previous publications [33][34][35][36][37] ) show resistance to this specific vincristine regimen.

The Role of ITG 1 and ITG 4 in Self-Renewal and Renal Specification in SIX2+CITED1+ Cells
LAM511 is important in renal development, but together with LAM332 and LAM111, it can promote carcinogenesis. [38] LAM511-cell interaction is largely mediated by integrins ITG 3 1, ITG 6 1, and ITG 6 4. [39] Though the role of ITGs in cancer progression is complex and still not fully elucidated, [40,41] induction of drug resistance in renal carcinoma www.advancedsciencenews.com www.advancedscience.com cells is associated with changes in ITG 1 expression. [42,43] ITG 4 has been identified in aggressive breast CSC [44] and ITG 4-targeted immunotherapies showed some efficacy against both CSCs and non-CSCs in breast and colon carcinoma mouse models. [41,45] To examine the role of LAM511-ITG binding in WT SIX2+CITED1+ cells, we first determined the localization of ITG 1 and ITG 4 in hFK and WT. ITG 1 and ITG 4 were expressed in various developing structures in the hFK, including UB and CM, in proximity to SIX2+ and CITED1+ cells, while expression of ITG 1 and ITG 4 in WT was more diffuse ( Figure  3A,B; Figure S7A, Supporting Information). Transcriptomic and protein analysis revealed different patterns of expression of these two ITG ( Figure 3C,D), with ITG 1 being more highly expressed in hFK SIX2+CITED1+ cells compared to WT SIX2+CITED1+ cells.
We next tested if LAM511-ITG 1 and LAM511-ITG 4 signaling regulate the NP state in cultured hFK SIX2+CITED1+ cells. Antibody neutralization of ITG 1 increased the percentage of cells co-expressing SIX2 and CITED1 but did not change the percentage of cells expressing SIX2 (at 72 h, 5 d, 28 d, Figure 3E-G). Similar results were achieved using treatment with an antibody specific to the active form of ITG 1 (9EG7), [46,47] while no changes in SIX2+CITED1+ cells were detected by activating ITG 1 with MnCl2, which induces localized conformational changes mimicking ligand-receptor occupancy [48] (Figure S7B, Supporting Information).
Neutralization of ITG 4 did not promote an increase in the percentage of SIX2+CITED1+ cells but did stimulate expression of LEF1, one of the first genes expressed in committed NP that determine cell cycle exit to differentiation. [2,3,49,50] On the contrary, LEF1 expression was significantly reduced by anti-ITG 1 treatment ( Figure 3H). Blockage of ITG 4 also activated components of WNT signaling (CSNK2A1, WNT4, CSNK2B, WNT10A, CSNK1D, and CSNK1G2), which are required for commitment ( Figure S7C, Supporting Information). [51] Finally, in mouse recipients of WT SIX2+CITED1+ cells, the speed of tumor growth tended to be enhanced by anti-ITG 1 ( Figure 3I) versus mice treated with anti-ITG 4.
By integration analysis, we identified 12 distinct cell clusters ( Figure 4C,D). Except for cluster 7 (almost equally shared between the two samples), the clusters were exclusively or preponderantly represented by either hFK or WT cells ( Figure 4C-E). WT SIX2+CITED1+ cell clusters were highly enriched in WT1 ( Figure S12A,B, Supporting Information) but devoid of H19 (tumor suppressor gene) expression ( Figure S12C-E, Supporting Information), confirming the previous findings. [6,59,60] Since SIX2 and CITED1 expression ( Figure S12F,G, Supporting Information) was very low and deeper sequencing was prevented by saturation of the samples, qPCR was performed to confirm CITED1 and SIX2 expression ( Figure S12H, Supporting Information).
Using the cell characterization of the hFK by Lindström et al., [61] we stratified clusters based on the expression of specific nephrogenic genes ( Figure 4D,E). We identified 7 different cell types, and interestingly, the clusters identifying early podocyte maturation (clusters 4 and 12) were present only in hFK, suggesting that WT SIX2+CITED1+ cells might not be able to initiate podocyte commitment correctly. Only two clusters (8 and 9) appeared committed, but they lacked a specific differentiation signature defined in; [61] these clusters were enriched in angiogenesis and mesenchymal stem cell maintenance pathways, possibly reflecting a stromal phenotype. These results were also confirmed by GO enrichment, IPA analysis, and hierarchal clustering ( Figure 4F, Figure S12I-L and Dataset S#3-4, Supporting Information) and are further discussed in Discussion S2 (Supporting Information). We also detected by scRNA-seq differences in ITG 1 and ITG 4 expression between hFK and WT SIX2+CITED1+ cells, thus further supporting their role in regulating these cells ( Figure S12M-O, Supporting Information).
Cell cycle regulation is critical to balancing proliferation versus differentiation, FACS, and pathways involved in cell cycle regulation stratified distinctly between hFK and WT. Cluster 2 (hFK) expressed fewer G1/S checkpoint genes, while cluster 11 (WT) expressed fewer G2/M checkpoint control genes. To further  understand the role of the cell cycle in SIX2+CITED1+ cells, we investigated the expression of phase-specific genes [62] within clusters ( Figure 4F; Figure  To decipher the renal commitment signature of the SIX2+CITED1+ cells, we perform trajectory analysis on the SIX2+CITED1+ cells together with the tumor from which they were originally selected (WT total, WT-TOT) and the xenograft that they generated in vivo, identifying 3 distinct states. The distribution of the integrated samples along the trajectory is shown in Figure 4G. To assess the temporal progression to differentiation, we performed pseudotime analysis and determined the pseudotemporal ordering of the cells along the trajectory ( Figure 4H), which produced a major trajectory (right side) bifurcating into two major branches identified as upper and lower. Interestingly, when samples are visualized along the pseudotime, WT SIX2+CITED1+ cells segregated toward the earlier stages of the developmental trajectory (right branch) with some cells scattered along the left lower bifurcation (Figure 4G, top panel; WT-NP: blue, WT-TOT: red). On the other hand, WT-TOT cells were more scattered along the pseudotime scale, with a small cluster of cells localized together with WT SIX2+CITED1+ along the right-side branch and a larger portion of the cells populating the two developmentally more mature branches (left side, upper and lower).
To further study the tumorigenicity of SIX2+CITED1+ cells from WT, we compared their gene expression pattern with that of WT-TOT. scRNA-seq data (Figure S14A-D, Supporting Information) revealed a strong overlap between the tumor of origin and the SIX2+CITED1+ cells ( Figure S14A,B, Supporting Information). By PCA, WT SIX2+CITED1+ and WT-TOT samples varied minimally ( Figure S14C, Supporting Information). An exception was cluster 4, predominantly composed of WT-TOT cells ( Figure  S14D, Supporting Information) and characterized by enrichment for genes involved in the cellular response to metal ions ( Figure  S14E, Supporting Information), including MT1A metallothionine and SOD1 (Figure S14F,G and Dataset S#5, Supporting In-formation). While metallothionines have known roles in carcinogenesis, including in some renal carcinomas, [63] this is the first association of these genes with WT and may suggest an acquired resistance of the tumor to metal toxicity. We also saw a consistent expression of superoxide dismutase 1 (SOD1, involved in detoxification of reactive oxygen [64,65] ) across WT clusters ( Figure S14G and Dataset S#5, Supporting Information). Additionally, scRNAseq analysis showed enrichment of drug resistance genes such as ABC genes in WT versus hFK clusters ( Figure S14H, Supporting Information), possibly explaining the drug resistance of the xenografts originated from these cells in our in vivo experiments ( Figure 2F).
Integration of scRNA-seq data ( Figure 4I,J, GSE175698) showed distinct differences in gene expression between SIX2+CITED1+ cells from hFK and WT, WT xenografts and the tumor of origin (WT-TOT); similarities between WT SIX2+CITED1+ cells and their tumor of origin (WT-TOT) are evident. The transcriptomic analysis confirmed the high similarity of the xenografts generated from freshly isolated versus expanded SIX2+CITED1+ cells ( Figure S15A,B, Supporting Information). Tumors generated from freshly isolated or cultured SIX2+CITED1+ cells exhibited minimal transcriptional differences from the primary WT, as visualized by PCA ( Figure  S15C,D, Supporting Information).

Spatial Transcriptomic (ST) Analysis Reveals Unique Transcriptomic Maps of WT versus hFK
To correlate histological organization with gene expression, we characterized WT samples using spatial maps of gene expression patterns with the Visium 10X Genomics Platform, using normal hFK as a reference (Figure S16, Supporting Information and Methods for details, GSE178349, GSM5388190, GSM5388191, and GSM5388192). A detailed spatial characterization of the samples individually is reported in Discussion S2, Figure S17, and Dataset S#6 (Supporting Information). Briefly, this analysis showed that, unlike the hFK in which genes reflecting all stages of nephrogenesis can be recognized in distinct clusters that histologically identified the nephrogenic niches and more differentiated renal structures, in WT these nephrogenic genes are spatially scattered without relation to clearly recognizable structures and mixed with areas characterized by gene expression of aberrant (particularly muscle) differentiation.
ST was then used to study differences between clusters in WT versus hFK (Figure 5A-E  identified. Hierarchical clustering showed that clusters representative of each sample grouped together, with WT#3 clustering closer to hFK ( Figure 5C) than WT#12.
Clusters 0, 1, 6, and 8 (mostly representing WT#12, histologically defined as favorable, Figure 5B) corresponded to histological compartments of stroma with marked rhabdomyomatous differentiation, blastema, connective tissue, and mixed stroma/blastema, respectively. By GO analysis ( Figure S18 and Dataset S#7, Supporting Information), these clusters were highly enriched in non-renal development pathways (specifically muscle and bone differentiation, heart, and vasculature development), metabolism, and cell cycle regulation but with a clear paucity of kidney development genes (glomerulus and tubule morphogenesis, UB and metanephric mesenchyme development). These data showed a lack of nephrogenesis specification genes in WT expressing aberrant non-renal differentiation pathways, a feature shared by many tumors other than WT. [66,67] Cluster 4 (spatially a blastema, specific of WT#3, histologically defined as unfavorable) was highly enriched in kidney development genes (including those associated with metanephric mesenchyme, and UB branching), with downregulation of other tissue/organ differentiation, cell cycle, and metabolic pathways. Cluster 7 (WT capsule with fibrous tissue, also specific to unfavorable WT) highly expressed immunological response genes (adaptive, innate, and humoral) with downregulation of metabolic, IL12, and cell cycle pathways.
Cluster 2 (expressed by all samples with many spots from WT#3, with histology mixed between blastema and stroma) was highly represented by GO sets specific for mesonephric and ureteric development, possibly representing a "common" nephrogenic cluster shared by hFK and both types of WT (Figure S18, Supporting Information).
Clusters 3 and 5 (largely specific to hFK, with cluster 5 representing the CM and cluster 3 more differentiated structures) highly expressed GO sets typical of all stages of nephrogenesis ( Figure S18, Supporting Information). As in Figure S19 (Supporting Information), though clusters 3 and 5 expressed multiple GO sets related to kidney development, specific expression patterns and tissue localization suggested that cluster 3 is more induced toward nephrogenic differentiation (expression of NPHS1 and PTPTO), while cluster 5 is in a more uncommitted state (expression of PAX2 and SALL1).
Next, we investigated the differences between the blastema of the WT patient samples (clusters 4 and 6) against the nephrogenic zone of the hFK (cluster 5). Even if these clusters shared some similarities ( Figure 5F), distinct gene expression patterns were noted in WT blastema versus the hFK nephrogenic zone. GO analysis, identified by DE genes between the groups (upregulated in Figure S20A,B, Supporting Information; downregulated in Figure S20C,D, Supporting Information) showed enrichments of renal development pathways with reduced cell death regulation, response to heat and peptide metabolic processes in the hFK. WT#3 blastema showed enrichment of immune related pathways and the reduction of regulation of cell motility, extracellular matrix organization, and blood vessel formation, while WT#12 blastema showed enrichment of cell cycle, proliferation, and hypermethylation pathways and reduction in kidney development pathways. This data highlighted that although histologically similar, the blastema portions of the tumors differ from the hFK nephrogenic zone by a distinct set of genes and molecular pathways (Dataset S#8, Supporting Information).
We then compared the histologically represented blastema in both WT patient samples (clusters 6 and 4, Figure S19, Supporting Information). DE gene analysis between these 2 clusters revealed important differences in gene expression between WT samples, allowing us to identify specific blastema transcriptomic maps. We identified potential targets of interest specific to different WT blastema ( Figure S21, Supporting Information), such as CLEC4M ( Figure S21B, Supporting Information); these data were also confirmed by immunohistochemistry in multiple patient samples, as shown in Figure S21E-L, and Dataset S#7 (Supporting Information) and discussed in Discussion S2 (Supporting Information).

Spatial Analysis of the SIX2+CITED1+ Cells Reveals Intrinsic Transcriptomics Differences in WT versus hFK
We next used ST data to identify the spatial localization of these cells in the WT samples versus hFK. Figure 6A,B shows that spots expressing SIX2 and CITED1 were present mainly in cluster 5 in the hFK (as expected, since this cluster represents the nephrogenic zone), while in WT, expression of these two genes was present in multiple clusters that histologically identified the blastema (4 and 6) but also in clusters that represented connective tissue (cluster 1) and stroma (cluster 0); histological regions in which presence of uncommitted progenitors is not expected during normal nephrogenesis. Analysis of the DE and GO sets of the SIX2+CITED1+ spots in the hFK versus WT revealed that the SIX2+CITED1+ spots of the WT (specifically in the WT blastema clusters 4 and 6) were highly enriched in pathways related mainly to muscle differentiation and regulation of the immune system (Dataset S#9, Supporting Information). They also showed clear downregulation of pathways related to renal development (Figure 6C), thus highlighting that even if these spots express SIX2 and CITED1, they are intrinsically different in the hFK compared to WT.
Since regulation of SIX2 and CITED1 expression is key to guiding proper renal development, [2,3] we investigated the correlation between SIX2 and CITED1 expression in our samples ( Figure 6D). Using the Pearson correlation coefficient, we determined that in the hFK (cluster 5) and WTs (clusters 4 and 6) there is a significant correlation between the expression of SIX2 and CITED1. In the hFK and WT#12, the correlation was positively linear (clusters 5 and 6). However, although in WT#3 (cluster 4) we identified a significant inverse correlation, the number of spots detected for this sample was limited, therefore the interpretation of this correlation requires further investigation.
Next, we also used ST to identify spots representing the committed NP population (SIX2+CITED1-cells ( Figure S22 and Dataset S#9, Supporting Information). In the hFK, the spots with committed progenitors were mainly present in cluster 3, characterized by developing nephrogenic structures, while in WT, these spots were present in different clusters (in addition to the blastema clusters 4 and 6). Analysis of DE and GO sets showed that the committed progenitors in the hFK are indeed committed towards renal differentiation in contrast to WT, where the nephrogenic differentiation pathways were not expressed ( Figure S22, Supporting Information). Of note, even if the majority of the SIX2+CITED1-spots were detected in cluster 3, ST also determined the presence of these committed progenitors in cluster 5 in regions that marked the early nephrogenic structures (renal vesicles) and not within the cap mesenchyme, where only SIX2+CITED1+ cells reside. Interestingly, the GO set analysis showed that even if the SIX2+CITED1-spots of clusters 3 and 5 share a nephrogenic signature, the cells in cluster 5 showed enrichment of earlier nephrogenic differentiation pathways (like metanephric mesenchyme regulation and mesenchymal-epithelial transition) versus cluster 3 representing more differentiated pathways like blood vessel development or anatomical structure morphogenesis ( Figure S22 and Dataset S#9, Supporting Information).
These ST observations highlighted that even if spots were identified by the same pattern of genes (SIX2 and CITED1), their  Table summarizing the distribution of SIX2+CITED1+ spots per cluster and the percentage of total SIX2+CITED1+ spots compared to the total number of spots per sample (right side column) C) List of top gene ontology (GO) biological process of hFK nephrogenic zone (cluster 5, SIX2+CITED1+ spots) compared to WT#3 blastema (cluster 4) and WT#12 blastema (cluster 6) SIX2+CITED1+ spots (left panels) and in WT#3 blastema (cluster 4) or WT#12 blastema (cluster 6) SIX2+CITED1+ spots compared to hFK nephrogenic zone (cluster 5, SIX2+CITED1+ spots, right panels). Fold enrichment of each gene set is shown, P < 0.05; upregulated DE genes were used for each comparison. D) Pearson correlation between the expression of CITED1 and SIX2 was performed on SIX2+CITED1+ spots found in clusters of WT#3 (cluster 4), hFK (cluster 5), and WT#12 (cluster 6), with dots representing spots; dots color intensity indicating the number of spots superimposed with a similar expression of SIX2 (x-axis) and CITED1 (y-axis). www.advancedsciencenews.com www.advancedscience.com transcriptomic profile was different based on their morphological context.
Next, we used ST to identify spots characterized by uncommitted NP (SIX2+CITED1+) and committed (SIX2+CITED1-) in relation to ITG 1 and ITG 4 expression ( Figure S23 and Dataset S#9, Supporting Information). In hFK, the SIX2+CITED1+ ITG 1+ or ITG 4+ spots are localized in the CM (cluster 5) and the SIX2+CITED1-ITG 1+ or ITG 4+ spots are present in both cluster 5 and 3; possibly suggesting that expression of these ITGs could be correlated with SIX2 and CITED1.
DE gene expression in the hFK versus the WT samples of the SIX2+CITED1+ spots or SIX2+CITED1-ITG 1+ spots was represented in the hFK by pathways related to mainly to renal nephrogenesis in addition to the regulation of cell cycle, cell division, mitosis, and ECM reorganization. A very similar signature was also evident in the hFK in the spots (SIX2+CITED+ or SIX2+CITED1-) characterized by the presence ITG 4+. Transcriptomics (a pathways analysis) of these spots in the WT patient sample was different: in WT#12 gene expression identified pathways related to a commitment towards differentiation while in WT#3 gene expression was more representative of pathways of immune-response and self-renewal; but none of these spots (both in WT#3 and WT#12) were committed towards nephrogenesis.
One possible interpretation of this pattern in WT spots could suggest that in WT cells are not committed toward a mature (terminally differentiated) renal fate, despite the expression of SIX2 and CITED1. This interpretation is supported by the absence of expression of the major regulator of renal development, WNT, [3,66] and the presence of FGF14 (member of the FGF family, a major signaling pathway of renal development involved in self-renewal, [2,3] Figure S24A-C, Supporting Information) in WT#3, for example. Furthermore, high expression of SIX2 in WT#12 could also explain the prevalence of muscle-like tissue since SIX2 is a well-recognized factor during muscle differentiation. [66] The inability to lose SIX2, accompanied by a push toward differentiation caused by the loss of CITED1, may lead to aberrant maturation toward skeletal muscle phenotypes. [68] On the contrary, in the WT#3, the unbalanced presence of CITED1 prevents differentiation, although a few cells appear to acquire a cardiac-like phenotype, possibly driven by the high expression of CITED1, which normally promotes cardiac differentiation. [69] In summary, ST data (based on gene expression) suggested (as the in vitro data which are based on ITG activity) that ITG 1 and ITG 4 are crucial to prime SIX2+CITED1+ and SIX2+CITED1cells toward differentiation during normal development. In WT this process is altered: NP is maintained in a self-renewal state or, if primed towards a differentiation process, they are not pushed toward nephrogenesis.

Discussion
This study provides an extensive characterization of the WT SIX2+CITED1+ cells and, for the first time, their spatial gene profiles. The presence of NP and the expression of SIX2 or CITED1 in WT tumor cells has been confirmed [12][13][14][15] but unlike previous studies of human WT NP/CSC in which limited subtypes of WT were examined, [70] or analysis was performed in chemotherapy-treated samples, [70,71] or focused only on SIX2+ committed NP, [13] our study incorporated several unique design features. First, uncommitted NP expressing both SIX2 and CITED1 were derived from several WT samples. Second, naïve SIX2+CITED1+ cells from hFK served as a critical reference (control). Third, the WT samples analyzed were not exposed to any chemotherapy, which is well known to alter CSC's expression pattern and differentiation status.
Our analysis, based on the developmental biology of the human kidney [2] shows that SIX2+CITED1+ cells present with WT CSC characteristics. Limiting dilution experiments showed that when these cells are injected into immunodeficient mice, they form xenografts that recapitulate histologic features of the tumor of origin, demonstrate metastatic potential, and resistance to a chemotherapeutic drug. Further, pseudotime trajectory analysis points to SIX2+CITED1+ cells as the root cells that give rise to WT. We also demonstrated that WT SIX2+CITED1+ cells differ from hFK in the expression of genes regulating self-renewal, commitment/differentiation, and proliferation, likely changing the ultimate fate of the cells and preventing normal depletion of NP from the nephrogenic niche before birth. We also found a signature of muscle/stroma differentiation in two clusters composed mainly of WT SIX2+CITED1+ cells. The same signature was found in more differentiated cells in the tumor of origin and the tumor generated by the cells after transplantation.
Other methods of selecting nephrogenic CSC (for example, NCAM1+ALDH1+ cells [70,72] ) have not clearly defined their uncommitted states, such as coexpression of SIX2 and CITED1. NCAM1+ALDH1+ cells can generate WT xenografts with triphasic morphology (even at lower concentrations) only if isolated from propagating WT fragments but do not reliably generate xenografts if freshly isolated from WT tissue and transplanted (30%-10% success rate). [36][37]70,72] Moreover, these cells lose xenograft generation capacity completely when cultured. [70] We used ST data to identify spots positive for NCAM1 and ALDH1 (transcript variant ALDH1A2 was used for analysis based on high expression of this variant in many WT as indicated in https://target-data.nci.nih.gov/Public/WT/mRNA-seq. ALDH1A2 is also expressed in high-risk WTs [73] ) and compared these spots with the SIX2 and CITED1 positive spots shown in Figure 6A ( Figure S25, Supporting Information). In the hFK, NCAM1+ALDH1+ spots were highly expressed in the developing structures (cluster 3); in WT, these spots were not always present in the blastema portion if compared with the SIX2+CITED1+ cells. Even if the comparison of these putative classes of CSC needs more studies (not the focus of this work), our preliminary ST analysis suggests that the NCAM1+ALDH1+ cells might represent a more committed population of cells compared to the SIX2+CITED1+ cells, in addition to the observation that NCAM1 is expressed in cells mainly expressing SIX2 during the first stages of epithelization and not in uncommitted NP. [74] This may be why the generation of xenografts from NCAM1+ALDH1+ cells obtained from primary digested tumors is difficult.
We defined that LAM511 and ITG 1 and ITG 4, the major LAM511 binding integrins, are crucial to prime uncommitted NP (SIX2+CITED1+ cells) and committed NP (SIX2+CITED1cells). LAM511 was a rational candidate to explore since it is involved in the maintenance of pluripotency in other systems. [55] When SIX2+CITED1+cells are cultured on LAM511, expression of SIX2 and CITED1 is maintained for prolonged periods, both facilitating in vitro study of normal nephrogenesis and offering, for the first time, a stable in vitro system for studying the WT CSC.
We also confirmed that hFK regions with uncommitted SIX2+CITED1+ NP cells also express ITG 1 (or ITG 4) and are mainly found in the nephrogenic zone, while histologic areas with committed (SIX2+CITED1-) cells expressing ITG 1 are found in developing structures. In WT, areas with uncommitted and committed cells and ITG 1 expression are found across blastema and stroma with no specific histological characteristics suggestive of a nephrogenic niche. Unlike hFK, none of the WT clusters express WNT9b or WNT7b, the most important genes for initiating renal differentiation. [3,75] WNT9b secreted from cells at the tip of the UB induces the formation of the renal vesicles, the first step of nephrogenesis. [3] These results support the idea that proper nephrogenic differentiation requires a specific pattern of ITG 1 and ITG 4 expression in both uncommitted and committed NP. This nuanced regulation of differentiation is lost in WT CSC, leading to the loss of the signaling to receive proper instruction to complete differentiation, as evidenced by the dispersed irregular renal structures in WT histology.
This study also correlated gene expression patterns with histological features using spatial transcriptomics. Analysis of the ST data, supported by bulk and sc-RNAseq data, clearly revealed that the transcriptional profiles of the WT samples have a paucity of genes involved in later steps of renal differentiation, correlating with the histologic absence of mature kidney structures like tubules or glomeruli. WT areas outside blastemic foci also contained SIX2+CITED1+ cells and included clusters primed for non-renal fate, suggesting that some WT CSC can deviate from normal renal differentiation. Based on our analysis, we think that in some instances in WT, NP seems capable of commitment to differentiation but incapable of the typical mesenchymalto-epithelial transition of renal differentiation. Alternatively, induction of NP to undergo nephrogenesis could be deficient, and in the absence of WNT signaling and the presence of high FGF expression, cells seem fixed in an uncommitted pre-renal state. [1,3,76] We believe these different unbalanced NP gene expression patterns have important implications in guiding our understanding of tumor initiation and progression since different mutations can lead to the development of different specificsubtype of WT. [10] Most importantly, our studies show that histologically similar blastema regions characterized by the expression of SIX2 and CITED from WTs had very different gene expression profiles. We identified genes predominantly expressed in unfavorable WT blastema (like CLEC4M, histologically validated in multiple WT of the same subtype classification) that may ultimately be useful in deeper molecular characterization of WT for purposes of staging and prognosis.
In conclusion, highlights of this study are that i) SIX2+CITED1+ cells are the nephrogenic CSC at the origin of WT, ii) an interplay between ITG 1 and ITG 4 regulates NP state, and iii) important differences in gene expression distinguish WT samples even in histologically similar regions. Given the heterogeneity of WT from sample to sample, we recognized that more studies are needed to characterize WT fully. Nevertheless, we believe that our extensive transcriptomic studies, together with in vivo and in vitro experiments, might help the development of novel strategies, such as manipulating CSC-ECM-niche interaction signaling or targeting the nephrogenic signature in WT CSC that could limit cancer development or spread.

Experimental Section
Sample Collection and Cell Suspension: Wilms tumor (WT) surgical explants were collected after informed consent under a protocol approved by the IRB of CHLA. Human fetal kidney (hFK) samples from elective terminations at 10-20 weeks gestational age (WGA) were obtained from the CHLA Tissue Bank and approved by the CHLA and USC IRB. WT samples were provided to the lab after examination by the CHLA Pathology Department (Table 1). WT staging was determined by the CHLA Pathology department per Children's Oncology Group (COG) staging guidelines [9] and standard CHLA clinical practice.
Cell suspensions were prepared as previously published. [2] After mechanical dissociation, cells were digested with 125 U mL −1 collagenase I (Worthington, # LS004197) in RPMI-1640 (Gibco, #11875093) at 37°C for 35 min, then passed through a 100 μm cell strainer and a 40 μm cell strainer (Corning, #352360, 352340) with washes of 1x PBS (Gibco, # 14190144). Undigested glomeruli and tubules from hFK samples did not pass through the filters. The filtrate suspension was then centrifuged at 1500 rpm for 5 min and erythrocytes were eliminated using a red blood cell lysis kit (Miltenyi Biotec, #130-094-183). Aliquots of isolated cells were stored in CryoStor cell cryopreservation medium (Sigma-Aldrich, #C2874) for later use and the rest were used for flow cytometry analyses, cell isolation using Smartflare RNA probes, scRNA-seq, bulk RNA-seq, or transplantation into NOD/SCID mice (xenograft studies) and in vitro experiments.
Western Blotting: Western blotting for protein quantification was performed as previously published. [16] Membranes were blocked for 1 h at RT in either 5% skim milk (w/v) or 5% BSA (w/v) in tween-20 trisglycine buffer solution (TBS-T), depending on the primary antibody used, and probed with primary antibodies at 4°C O/N. Horseradish peroxidase (HRP)-conjugated secondary anti-rabbit and anti-mouse antibodies (Sigma-Aldrich) were used at 1:20000 and 1:30000 dilutions respectively in 2.5% skim milk in TBS-T. Membranes were developed using the Super-Signal West Femto Western Blotting detection reagents (Thermo Scientific, #34096) or SuperSignal West Pico (Thermo Scientific, #34577) and impressed on Biomax Light Films (GE Healthcare, #66302). Nuclear fraction protein isolation was performed using a Nuclear Extraction Kit (Abcam, #113474) following manufacturer's protocols. Data from 3 independent experiments were quantified by densitometry (all normalized against a housekeeping gene, -actin or H100).
qRT-PCR: qRT-PCR was performed as published. Propidium Iodide Staining: hFK and WT samples were used to assess cell cycle phase distribution of SIX2+CITED1+ and SIX2+CITED1-cells. Briefly, total digestates of hFK and WT samples were suspended in 1 mL PBS (1 million cells) and incubated for 15 min in -20°C ethanol while vortexing. Next, cells were centrifuged and incubated for 15 min in 5 mL PBS. Cells (10 7 cells/100 μL) were blocked in a human IgG 1× solution for 10 min and incubated with antibodies against SIX2 and CITED1 conjugated with Zenon labels (Life Technology, #647R-Z25308 and 488M-Z25002,) for 1 h. Cells were then washed and counterstained with 1.5 × 10 −6 m propidium iodide -FluoroPure Grade solution (Thermo-Fisher, #P21493). Flow cytometry analysis of propidium iodide staining and costaining for SIX2 (Alexa-Fluor APC) and CITED1 (Alexa-Fluor 488) was performed on a FacsCanto BD Biosciences instrument.
Bulk RNA-Seq: RNA extraction was performed immediately after sorting of SIX2+CITED1+ cells (passage 0) using the RNeasy Micro Kit (Qiagen, #74004) following manufacturer's recommendations. After cDNA production and construction of DNA libraries, the samples were run on an Illumina NextSep500 sequencer. Differential gene expression was analyzed using ERCC ExFold probes with the Remove Unwanted Variation R/Bioconductor software package combined with edgeR 2 . WT datawere deposited in Gene Expression Omnibus (GEO) under accession number (GEO for WT: #GSE176342 and GEO for hFK: #GSE74450) and processed as previously described. [2] Differentially expressed (DE) genes between samples were determined by calculating fold change using the Reads Per Kilobase of transcript, per Million mapped reads (RPKM) values. DE genes with fold change >1.5 or ←1.5, for upregulated and downregulated genes respectively, were considered for downstream analysis using the Gene Ontology (GO) resource and Ingenuity pathway analysis (IPA).
Final sequencing libraries were run and analyzed on an Illumina HiSeq 4000 sequencing system, PE150, and 625 M total reads/lane (QuickBiology). Approximately, 300 million reads/samples were sequenced. Using Partek Flow software, cells that contained fewer than 1600 expressed genes or >10% mitochondrial transcripts were removed from analyses. For each cell, the expression of each gene was normalized to the sequencing depth of the cell, scaled to a constant depth (10 000), and log-transformed. To filter out immune cell sequences, following k-means (k = 5) clustering, the cluster containing immune cells (based on CD45 expression) was filtered out. However, for analysis of WT-Xe and WT-Xe from cultured NP, mouse cells were removed prior to processing, so the elimination of CD45 cells was not necessary for those samples ( Figure S11, Supporting Information). After appropriate filtering, unsupervised graph-based clustering was performed and dimensionality reduction and visualization using the UMAP algorithm. [77] DE genes for each cluster versus all other clusters were analyzed by ANOVA using the Partek software. Genes with a fold change >2 or ←2, for upregulated and downregulated genes respectively, were considered for further analysis on the GO platform or IPA. Trajectory inference analysis was performed using Monocle2 algorithm within the Partek Flow scR-NAseq toolbox. The root of the trajectory for the pseudotime analysis was based on expression of the nephrogenic signature (as described by Lindström et al. [61] ).
Transcriptomic Analysis: Gene Ontology Enrichment Analysis, Ingenuity Pathway Analysis and Heatmap Generation: The GO resource platform (http://geneontology.org) was used to determine the enrichment of biological processes in samples of interest. [78][79][80] DE genes (refer to bulk RNAseq, scRNA-seq, and spatial gene expression analysis for gene expression cutoff) were uploaded onto the GO platform and GO sets were generated, powered by the PANTHER (Protein analysis through evolutionary relationships) classification system. GO sets with an enrichment score >2 and p-value < 0.05 were considered significant and biologically relevant.
IPA software (Qiagen) was used to investigate biologically significant processes in samples of interest and determine the activation of canonical pathways based on bulk and scRNA-seq data (http://www.ingenuity.com). Fold-change values were uploaded onto the IPA software and overlaid with the global gene network in the Ingenuity Knowledge Base. Specifically, the pluripotency pathway, nephrogenic development, and WT-associated genes were explored between samples from the bulk RNA-seq experiments (refer to the bulk-RNA seq section for data processing). Significant differentially regulated functional clusters or single pathways were further grouped by the indicated functional classes and compared by the enrichment score. Pathways of interest were selected as the representative for each sample and cluster if significantly enriched (p-value < 0.05) or folddifference >1.5 patterns emerged from analysis. Z-score >2 was considered to predict potential activation/inactivation of pathways.
Morpheus software (https://software.broadinstitute.org/morpheus) was used to generate heatmaps for visualization of either total gene expression or specific gene sets. Heatmaps were generated by uploading raw values [reads per kilobase of transcript/million reads (RPKM) for bulk RNAseq data, least squares mean (LSMean) for scRNA-seq data]. Hierarchical clustering was applied to rows (gene lists) and columns (samples) based on a "one minus Pearson correlation" metric. Gene expression values are mapped in colors using the minimum (blue) and maximum (red) of each row (gene) independently. No cutoffs were applied when generating heatmaps.
Spatial Gene Expression: OCT embedded 16 WGA hFK, WT#12, and WT#3 samples were sectioned to 10 μm thickness (Leica, Rotary Cryostat CM1510). RNA was isolated using the Qiagen RNeasy Mini kit (Qiagen, Hilden, Germany) and RNA quality was assessed using the 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA). All RNA tissues had RNA integrity (RIN) above 9. Tissues were then cryosectioned onto preequilibrated Visium tissue optimization or gene expression slides (10× Genomics, Pleasanton, CA). These sections were then fixed in chilled methanol, stained with hematoxylin and eosin, then imaged with an Aperio AT Turbo (Leica Biosystems, Buffalo Grove, IL) at 20× magnification. Automated real-time stitching of tiled images yielded a final image of the whole slide which was imported to image analysis software (ImageScope, Leica Biosystems).
Using visium spatial gene expression reagent kits -tissue optimization user guide (10× Genomics), tissues on optimization slides were permeabilized in a time course experiment and reverse transcription was performed using fluorescently-labeled nucleotides, resulting in fluorescent cDNA bound to the capture areas. Tissues were enzymatically removed, and fluorescence imaging was performed using a Zeiss Axio Scan.Z1 (Carl Zeiss Microscopy, White Plains, NY) digital slide scanner with a Texas Red filter set. Whole slide scanning was carried out at 20× with 150 ms exposure per image frame. Sequentially imaged frames were automatically stitched by the data acquisition software (Zen 3.0).
After H&E staining and imaging of tissue sections on gene expression slides, the sections were permeabilized for 18 min to release polyadenylated mRNA from overlying cells onto the capture areas of the slide. (A permeabilization time of 18 min resulted in the maximum fluorescence signal with the lowest signal diffusion.) Following the manufacturer's user guide (10× Genomics), the bound mRNA was then reverse transcribed, resulting in spatially barcoded, full-length cDNA. Second-strand synthesis was performed, followed by denaturation and transfer of the cDNA from the slide to a PCR tube. qPCR (KAPA SYBR FAST qPCR Master Mix, Roche Sequencing and Life Science, KAPA Biosystems, Wilmington, MA) was used to determine the number of cDNA amplification cycles required. After cDNA amplification, fragment analysis of the cDNA was performed on a 4200 TapeStation (Agilent Technologies). Library construction, performed on a portion on the cDNA, consisted of enzymatic fragmentation of the cDNA, end-repair, and A-tailing, followed by double-sided size selection. Libraries were sequenced on a NovaSeq 6000 using a custom paired-end sequencing protocol, consisting of: read 1, 28 cycles; index 1, 10 cycles; index 2, 10 cycles; read 2, 90 cycles.
Using the Visium 10× Genomics Platform we generated spatial maps of gene expression in favorable and unfavorable WT and compared them with the hFK map. Three samples (WT3, WT12 and hFK) were sequenced at the sequencing depths of 288 M, 696 M and 352 M, and detected 1766, 3735, and 1774 spots in the samples, respectively. This corresponded to a median of unique molecular indices (UMIs) per spot of 6805, 15 369, and 12 442, and a median of genes per spot of 3123, 4826, and 5326.
ST Data Analysis: Raw sequencing data were demultiplexed and converted to fastq format by using bcl2fastq v2.20. Space Ranger software v1.0.0 (10× Genomics) was used for read alignment, tissue detection, fiducial detection, and barcode/UMI counting with pre-defined default parameters. Briefly, raw reads were aligned to the human reference genome and transcriptome, then reads with the same barcode, UMI and gene annotation were grouped for UMI counting. For image processing, the slide barcoded spot pattern was manually aligned to the input slide image. Then, tissue and background in the slide image were discriminated. Finally, a UMI count matrix was generated consisting of the gene identities as rows and spatial barcodes as columns.
Seurat v3.2 in R v3.3.1 was used for advanced downstream data analysis. The matrices were loaded into Seurat to create a Seurat object, then normalized using sctransform [81] in order to account for variance in sequencing depth across data points. To obtain 2D projections of the population dynamics, principal component analysis (PCA) was first run on the normalized gene-barcode matrix of the top 5000 most variable genes, to reduce the number of features. These genes were identified based on their mean and dispersion as described by Macoscko et al. [82] After PCA, dimensionality was reduced by uniform manifold approximation and projection (UMAP) technique which enabled further visualization of cells in a 2D space. Based on their relative positions in the UMAP plot, an unsupervised graph-based clustering was performed to group cells for the clustering analysis. [83] To identify molecular features that correlate with spatial location within a tissue, differential expression was performed based on pre-annotated anatomical regions. Alternatively, features were searched exhibiting spatial patterning which models spatial transcriptomic data as a mark point process and computed a "variogram" to identify genes expression levels dependent on their spatial location. Data from each sample were visualized for exploration via Loupe Cell Browser, including generation of heatmaps. Significantly DE genes (p-adj value < 0.05 and avg_log fold-change >0.25 or ←0.25) were considered biologically relevant and used for downstream GO analysis. For Venn analysis of the ST data a p-adj value < 0.05 and avg_log fold-change >0.5 or ←0.5 was used.
Statistical Analysis for In Vitro and In Vivo Experiments: All comparisons between groups were analyzed using either a two-tailed unpaired Student's t-test or Mann-Whitney U-test depending on the normality of distribution. Comparisons between more than two groups were performed by one-way ANOVA with multiple comparisons or two-way analysis of variance with Bonferroni correction for multiple comparisons, as applicable. All results are represented as the mean ± standard error of the mean. Data from in vivo experiments comparison of the probability of survival were analyzed using the Gehan-Breslow-Wilcoxon test. Data analysis was performed with the Prism Software (version 8.0; GraphPad). Significance was expressed as p-value (*p < 0.05; **p < 0.01; *** p < 0.001; **** p < 0.0001).

Supporting Information
Supporting Information is available from the Wiley Online Library or from the author.