Introduction

Primary cilia are highly conserved signaling organelles that protrude from the surface of most vertebrate cells. Constructed from a microtubule scaffold, the ciliary axoneme is anchored in the plasma membrane via the basal body, a modified mother centriole. A specific structure above the basal body, the ciliary transition zone (TZ), regulates the entry and exit of proteins, thereby acting as a “ciliary gate” that regulates the composition of cilia and creates an independent cell compartment. By modulating several signaling pathways such as SHH, WNT and PDGFR, cilia play essential roles in embryonic development and adult tissue homeostasis. Consequently, defects in primary cilia can lead to a large group of rare genetic diseases called ciliopathies1,2, which can affect the development of multiple organs, particularly the retina, kidneys, lungs, skeleton, heart and the central nervous system (CNS). Ciliopathies are characterized by considerable phenotypic and genotypic overlap and, rather than being considered separate clinical entities, they are now recognized as a spectrum of genetic diseases. Additionally, there is extensive phenotypic and genetic heterogeneity amongst ciliopathies, and mutations in the same gene can cause strikingly different phenotypes, without necessarily any genotype-phenotype correlation.

Although this is not completely deterministic, mutations in genes encoding ciliary TZ proteins often lead to neurodevelopmental ciliopathies such as Joubert Syndrome (JS) and the fetal lethal Meckel Gruber syndrome (MKS)1,3. Most of the current knowledge about neurodevelopmental defects in ciliopathy patients and their underlying mechanisms was gained by in vivo studies on animal models and in vitro analyzes. Using these approaches, it was demonstrated that primary cilia are present throughout the neuronal lineage, from neural stem cells to progenitor cells, mature neurons and glia, where they play an essential role in CNS formation by modulating signaling pathways and cell cycle progression4,5,6,7,8,9,10. We and others have demonstrated that the ciliopathy causal gene RPGRIP1L encodes a large scaffolding protein acting as a keystone of TZ construction11,12,13,14,15,16,17,18,19. In addition, it regulates proteasomal activity at the ciliary base and autophagy in murine cells20,21 and is involved in planar cell polarity in zebrafish and mice22. Furthermore, an essential function for Rpgrip1l in SHH signaling was demonstrated in mouse embryonic tissues and cultured cells11,15,16,23. However, our understanding of the role of RPGRIP1L and other ciliary proteins in human neural fate is still very partial and has been hampered due to the paucity of human models, limiting our ability to understand the role of cilia in human development, and the origin of neural developmental defects in ciliopathies. Studies using increasingly diverse hiPSC-based organoid models in which the role of ciliary proteins can be investigated are therefore of growing interest to the field. These include studies that investigate ciliary proteins in retinal24, cortical25,26 or renal27 organoids.

Moreover, evidence suggests that primary cilia function could be only partially overlapping between human and animal models, partly due to specificities in human CNS development. For example, while Rpgrip1l mutant mouse embryos, like other ciliary mutants, display a severe disorganization of the ventral diencephalic and hypothalamic regions and of the ventral spinal cord, no such phenotype was described in ciliopathy patients6,7,28,29,30.

In this study, we took advantage of pluripotent stem cell-based organoid approaches to study the role of the TZ proteins RPGRIP1L and TMEM67 in human neural progenitor (NP) fate. Both proteins are involved in the gating function of the TZ and have been shown to regulate SHH-dependent spinal patterning in mouse models8,31,32,33. However, they are structurally and functionally different: RPGRIP1L is a scaffolding protein and a cornerstone of TZ assembly18,19, while TMEM67 is a Frizzled-like transmembrane receptor8,31, and the two proteins belong to different TZ protein complexes34. Since the role for ciliary proteins in SHH-dependent patterning of the ventral spinal cord has been largely studied35, we used wild-type (WT) and Rpgrip1l-deficient mouse embryonic stem cell (mESC)- and human induced pluripotent stem cell (hiPSC)-derived spinal organoids to compare the functions of RPGRIP1L in ciliogenesis and spinal progenitors between these two models. We show that mouse Rpgrip1l-/- organoids exhibit altered dorso-ventral patterning in response to SHH, similar to mutant mouse embryos, emphasizing the cell-intrinsic requirement for Rpgrip1l in spinal progenitors for cilia-dependent SHH pathway activation and commitment into ventral neuronal cell types. In contrast, RPGRIP1L-deficient human iPSCs submitted to a similar differentiation protocol are still able to form primary cilia and adopt SHH-dependent ventral spinal fates. Interestingly, this phenotype is not unique to RPGRIP1L deficiency, but is also observed in TMEM67-deficient human ventral organoids. Remarkably, in addition to the species-specific differences of TZ proteins in cilia formation and signaling between mouse and human, we discovered an unexpected role for the human RPGRIP1L protein in the precise assignment of antero-posterior (rostro-caudal) identities of human spinal progenitors. Indeed, mutant motoneurons (MNs) derived from early neuromesodermal progenitors (NMPs) shift their identity from brachial to more rostral/hindbrain identities. This is accompanied by a temporary loss of cilia in early RPGRIP1L-deficient axial progenitor cells. With this study, we bring a proof of concept that inter-species differences can be dissected by dynamically comparing mouse and human ex vivo organoids and we highlight human-specific functions for the ciliopathy gene RPGRIP1L in antero-posterior patterning of axial progenitors.

Results

Defective SHH-induced motoneuron specification in Rpgrip1l KO mouse spinal organoids

In vertebrates, graded SHH signaling along the dorso-ventral (DV) axis of the spinal cord is necessary for the specification of the floorplate (FP) and the four different ventral progenitor domains p3, pMN, p2 and p136,37,38,39,40,41 (Fig. 1a). The study of mouse ciliary mutant embryos demonstrated that primary cilia are required for SHH transduction and the correct establishment of spinal DV patterning15,35,42,43. In Rpgrip1l-/- (KO) mouse embryos, brain and spinal cord progenitors lack a functional cilium, leading to altered SHH pathway activation and repression6,15,16. As a consequence, in the spinal cord, the FP and p3 territories are not established, fewer pMNs are specified and the dorsal domains are shifted ventrally6,15. In order to test whether this phenotype can be reproduced ex vivo, even under conditions where SHH activation is forced by the addition of a constant dose of SHH agonist, we established WT and Rpgrip1l KO mESC lines from mouse blastocysts that we differentiated into ventral spinal organoids44,45,46. Briefly, mESC were differentiated into embryoid bodies (EBs) in N2B27 neuralizing medium for 2 days. EBs were then targeted to a neural spinal fate by combined exposure to the WNT agonist CHIR99021 and to retinoic acid (RA), and further ventralized by addition of SAG (Smoothened AGonist; SHH activation). By day 4, embryoid bodies had organized into spinal organoids that comprised spinal neuronal progenitors from different ventral identities, in particular from the MN lineage. Organoids were let to differentiate until day 7, when most progenitors underwent differentiation (Fig. 1b). WT and Rpgrip1l KO spinal organoids were collected at different time points of differentiation, sectioned and immunostained to assess the identity of the different progenitors obtained. We found that on average, 40% of spinal progenitors from control organoids committed to a pMN fate by day 4 of differentiation, as evidenced by Olig2 staining (Fig. 1c, d). pMNs progressively differentiated into post-mitotic MNs until day 6, a stage when only Islet1-positive MNs could be detected (Fig. 1c). Control organoids were also composed of progenitors from the more ventral domains p3 and FP, since we could detect the presence of Nkx6.1 + /Olig2-, Nkx2.2+ and Foxa2+ cells (Fig. 1d). Moreover, control organoids nicely recapitulated the sequential activation of ventral neural identities observed in vivo in the mouse spinal cord, as Olig2 expression preceded the initiation of Nkx2.2 and Foxa2.

Fig. 1: Rpgrip1l-deficient mouse spinal progenitors fail to adopt SHH-dependent ventral fates.
figure 1

a Diagram depicting the main transcription factors spanning the different dorso-ventral progenitor (light blue) or neuronal (light green) domains of the developing spinal cord. FP: floor plate. p0-p2, pMN and p3: ventral spinal progenitor domains. V0-V2, MN, V3: ventral spinal neuronal domains. b Schematic summary of the spinal 3D differentiation approach. The main neural progenitor and neuronal populations obtained at each time point with their molecular markers are indicated at the bottom. Small molecules used to induce or repress signaling pathways are indicated on top. CHIR: WNT agonist CHIR99021; SAG: SHH agonist; RA: Retinoic Acid. Modified after Duval et al.45. c Immunofluorescence for Olig2 (pMN) and Islet1 (post-mitotic MNs) in sections from WT and Rpgrip1l-/- organoids at indicated time points. d: day of differentiation. d Immunofluorescence for Nkx6.1 (p3-pMN-p2), Olig2 (pMN), Nkx2.2 (p3), Foxa2 (FP) and Pax6 (pMN-dp) on sections from WT and Rpgrip1l-/- organoids at indicated time points. In (c and d) nuclei are stained with DAPI. Scale bars: 50 μm. eh Violin plots quantifying the percentage of nuclei per organoids stained for the selected markers on sections from WT and Rpgrip1l-/- organoids at indicated time points. The median and quartiles are represented by dotted lines. Statistics: two-sided unpaired t tests with Welch’s correction (e p = 0.0003; fh p < 0.0001). N = 4 (e, f) or N = 3 (g, h) independent experiments; n = 1 clone per genotype. i, j Heatmap from bulk RNAseq data based on mean values of log2(tpm+1) for selected genes in WT and Rpgrip1l KO spinal organoids over time. N = 3 independent experiments for each time point; n = 1 clone per genotype. k, l qPCR analysis for Olig2 and Dbx1 in WT and Rpgrip1l-/- organoids at day 4. Data shown as mean ± SEM. Statistics: two-sided Mann-Whitney test (p = 0.0079). N = 5 independent experiments, n = 1 clone per genotype.

In contrast, Rpgrip1l-deficient progenitors failed to commit to ventral spinal identities. Indeed, the percentage of Olig2+ and Islet1+ cells of the MN lineage as that of Nkx6.1+ cells, was strongly reduced in mutant organoids (Fig. 1c–g), while Nkx2.2+ (p3) and Foxa2+ (FP) progenitor cells were nearly absent from mutant organoids (Fig. 1d). Instead, 80% of the progenitors were positive for Pax6, indicating a shift towards more dorsal identities (Fig. 1d, h). To gain a deeper understanding of the changes in cell fate occurring in Rpgrip1l-deficient organoids, we performed bulk RNAseq on WT and mutant organoids from day 2 to day 7, obtained from 3 independent differentiation experiments (Supplementary Fig. 1a). Our transcriptomic analysis confirmed the shift towards dorsal identities in Rpgrip1l mutant organoids, and showed that SHH-dependent ventral cell types were replaced by progenitors expressing markers characteristics of the p0-p2 domains such as Dbx1 (Fig. 1i and Supplementary Fig. 1b–c). qPCR performed on day 4 organoids confirmed these results and showed a more than 10-fold decrease of Olig2 expression while Dbx1 expression was increased 5 times in mutant compared with control organoids (Fig. 1k–l). Consistently, the molecular signature of post-mitotic neurons changed from V3-MNs in the control, to V1-V0 in Rpgrip1l KO organoids (Fig. 1j). The temporal analysis of the pan-neural markers Sox1 and Tuj1 expression in the time course of the differentiation (Supplementary Fig. 1e–f) showed that the reduction in MNs correlated with a transient decrease in neurogenesis between days 3 and 7 in Rpgrip1l KO compared with control organoids. Since MNs are the first neurons to differentiate in the spinal cord, this result is in line with the dorsalization occurring in the mutant organoids.

In Rpgrip1l mutant mice, the defective patterning of the ventral spinal cord has been attributed to an absence of primary cilia in spinal progenitors, and the consequent defect in SHH pathway activation6,15. Ensemble Gene Set Enrichment Analysis (EGSEA) performed on differentially expressed genes in day 3 and day 4 organoids showed an enrichment in genes related to spinal patterning (Supplementary Fig. 1c) and SHH signaling (Supplementary Fig. 1d). Consistently, qPCR performed on control and mutant organoids showed that while expression of direct targets of the pathway, Gli1 and Patched1 was upregulated in controls from day 2 upon SAG addition into the EB culture medium, their expression barely increased in the mutants, showing impaired activation of SHH transduction in absence of Rpgrip1l (Fig. 2a, b and Supplementary Fig. 1g, h). To test whether this phenotype was correlated with defective ciliogenesis we analyzed the ciliary marker ADP Ribosylation Factor Like GTPase 13b (Arl13b) by immunostaining on organoid sections. We found that, similar to what was observed in vivo, Rpgrip1l deficiency led to a drastic reduction of Arl13b+ cilia (Fig. 2c, d). This phenotype might reflect a ciliary gating defect, or an absence of ciliary axoneme. To discriminate between these scenarios, we labeled the axonemal protein Ift88 and measured the length of Ift88 staining to evaluate the axonemal length. We found that the axoneme of mutant progenitors was abnormally stunted, with a mean length of 0.2 µm compared to 1 µm in the control (Fig. 2c, e), suggesting that only the TZ might be present above the basal body. This ciliary phenotype was correlated with the complete absence of Rpgrip1l at the ciliary TZ in mutant organoids, consistent with the absence of Rpgrip1l protein in the KO mouse ES cells as shown by western-blot analysis (Fig. 2c, f). Overall, our results demonstrate that, as in vivo, Rpgrip1l is required in vitro in mouse spinal progenitors for cilia structure, SHH pathway activation and the subsequent acquisition of spinal ventral identities. Previous studies showed that spinal organoids can be applied to dissect the activity of the Hedgehog pathway in NPs46. Our results provide a proof of concept that they also provide a relevant model for the study of ciliary functions in neural cell fate.

Fig. 2: Rpgrip1l-deficient mouse spinal organoids fail to activate the SHH pathway and lack functional cilia.
figure 2

a, b Relative fold-change for Gli1 and Ptc1 expression in WT and Rpgrip1l-/- organoids at indicated time points. Data shown as mean ± SEM. Statistics: multiple two-sided Mann-Whitney tests (p = 0.028571). N = 4 independent experiments; n = 1 clone per genotype. c Immunofluorescence for the indicated ciliary markers on WT and Rpgrip1l-deficient day 5 organoids. Cilia are present at apical sites pointing outside or into inner cavities of the organoids. Scale bar: 2 μm. d Percentage of Arl13b positive cilia. Data shown are mean ± SEM. Statistics: two-sided unpaired t tests with Welch’s correction (p < 0.0001). N = 2 independent experiments; n = 1 clone per genotype. e Length of Ift88 staining. Data shown are mean ± SEM. Statistics: two-sided unpaired t tests with Welch’s correction (p = 0.0002). N = 2 independent experiments; n = 1 clone per genotype. f Western-Blot for Rpgrip1l and Actin on WT and Rpgrip1l KO mouse ES cells. N = 2 independent experiments.

Human RPGRIP1L-deficient neural progenitors differentiate into MNs in spinal organoids

To tackle the function of primary cilia and the TZ protein RPGRIP1L in human spinal development and SHH transduction, we generated human iPSC clones deficient for RPGRIP1L (RPGRIP1L-/-) via CRISPR/Cas9-mediated genome editing that lead to premature stop codons in two independent hiPSC lines, PCIi033-A and UCSFi001-A (Supplementary Fig. 2a, b). After genome editing, two RPGRIP1L+/+ (WT) and two RPGRIP1L-/- (KO) clones of each hiPSC line (WT: n = 4; KO: n = 4) were validated (Supplementary Fig. 2b–e) and used for differentiation into ventral spinal organoids following a well-established 3D differentiation protocol47,48 (Fig. 3a). Starting the differentiation of ventral spinal organoids from single cell hiPSCs, EBs were formed after one day in a floating culture in neural induction medium (dual SMAD inhibition) containing the GSK-3β antagonist CHIR99021 for WNT activation and axial progenitor specification47,48 (Fig. 3a). From day 4 onwards, axial progenitors were treated with RA and high concentrations of SAG in order to induce a ventral neural fate (Fig. 3a). At day 6 of differentiation, the majority of WT and RPGRIP1L-deficient cells expressed SOX2 and PAX6, both markers of early NPs, indicating that control cells as well as RPGRIP1L-deficient cells adopted a neural fate (Supplementary Fig. 3a). In response to SAG stimulation, the pMN marker OLIG2 started to be expressed between day 4 and day 6 of differentiation with an expression peak reached at day 9 where 75% of the cells within organoids expressed OLIG2 (Fig. 3b and Supplementary Fig. 3b). Combined to that, cells expressed the broad ventral marker NKX6.1 at day 9, further evidence of ventral spinal fate specification (Supplementary Fig. 3b). Surprisingly, no difference in OLIG2 and NKX6.1 pattern was observed between control and RPGRIP1L-deficient organoids (Fig. 3b and Supplementary Fig. 3b). In line with this, control and RPGRIP1L-deficient cells started to express ISLET1 and ISLET2, markers for postmitotic MNs, from day 9 onwards (Supplementary Fig. 3c–e). A final treatment of organoids with the γ-Secretase Inhibitor DAPT was used to inhibit NOTCH signaling and accelerate the neural differentiation into postmitotic MNs from day 9 onwards (Fig. 3a). Consequently, ISLET RNA and protein levels were strongly upregulated in WT and RPGRIP1L-deficient organoids between day 11 and day 14 (Fig. 3b and Supplementary Fig. 3c–f). Taken together, these results show that RPGRIP1L-deficient hiPSCs are able to process differentiation from stem cells to NPs and further from pMNs to postmitotic MNs. Nevertheless, our data reveal 15% less ISLET+ MNs at the end of differentiation in RPGRIP1L-deficient hiPSC-derived spinal organoids compared to controls (Fig. 3b and Supplementary Fig. 3f).

Fig. 3: RPGRIP1L-deficient hiPSC-derived spinal organoids adopt the MN fate.
figure 3

a Schematic summary of the spinal 3D differentiation approach. Main neural progenitor and neuron populations obtained at each time point with their molecular markers are indicated at the bottom. Drugs used to induce or repress signaling pathways are indicated on top. SB and LDN: SMAD antagonists SB-431542 and LDN-193189; CHIR: WNT agonist CHIR99021; SAG: SHH agonist; DAPT: NOTCH agonist; RA: Retinoic Acid. Modified after Maury et al.47. b Immunofluorescence analysis of pMN (OLIG2) and MN (ISLET1/2) production in WT and RPGRIP1L-/- spinal organoids. Quantifications show relative OLIG2 and ISLET1/2 positive areas per organoid over time. Data are shown as mean ± SEM. Asterisks denote statistical significance according to unpaired t tests with Welch’s correction (p = 0.0267). c Heatmap of selected gene expressions from bulk RNAseq data based on log2(tpm+1) values of WT and RPGRIP1L-deficient spinal organoids over time. d Log2(tpm+1)-graphs illustrate the dynamic expression levels of DBX2, OLIG2, NKX6.1 and NKX2.2 in WT and RPGRIP1L-/- organoids. Data are shown as mean ± SEM. e, f Log2(tpm+1) data from bulk RNASeq analyzes show TUJ1 and MAP2 expression in WT and RPGRIP1L-deficient organoids during differentiation. Data are shown as mean ± SEM. g Immunofluorescence analysis of TUJ1 on day 14. Data are shown as median with quartiles. Unpaired t test with Welch’s correction was performed for statistical analysis. h Immunofluorescence analysis of OLIG2-positive MN progenitors and NKX2.2-positive spinal progenitors at day 9 and day 11. Data are shown as mean ± SEM. Unpaired t tests with Welch’s correction was performed for statistical analysis (p = 0.0141). bh N: number of independent experiments; n: number of different clones analyzed per experiment. For b, g and h, 2 WT clones and 2 KO clones from each line (n = 4 for each genotype). For cf, 2 WT clones and 1 KO clone from each line (n = 4 for WT and n = 2 for KO). N = 3 in b, g; N = 2 in h; N = 1 in (cf). Scale bars: 50 µm in (b and h), 100 µm in j.

By analyzing the expression of various transcription factors that are distinctly expressed in specific DV regions of the developing spinal cord using bulk RNAseq data, we tested whether, despite the normal pMN production, human iPSC-derived organoids deficient for RPGRIP1L exhibited a shift in DV patterning towards the dorsal domains, as it was demonstrated for mESC-derived Rpgrip1l KO spinal organoids (Fig. 1). The unaltered expression of the broad ventral marker NKX6.1, the intermediate markers DBX1 and DBX2, and the dorsal spinal markers PAX3, PAX7, and ZIC1 (Fig. 3c, d) during differentiation, indicated that, contrary to mouse Rpgrip1l-deficient spinal organoids, human RPGRIP1L-deficient spinal organoids did not display a more dorsal progenitor type specification. We noted minor changes among the markers MSX1 and OLIG3 at day 2 and day 4 in RPGRIP1L KO organoids (Fig. 3c), likely reflecting subtle differences in axial progenitors, before DV neural specification.

To investigate whether the mildly reduced number of MNs in RPGRIP1L-deficient organoids could be caused by a delayed neurogenesis, we analyzed the expression of the pan-neuronal markers TUJ1 and MAP2. As expected, expression of both markers increased during differentiation with the highest expression of TUJ1 and MAP2 between day 11 and day 14 (Fig. 3e, f) and high protein levels of TUJ1 at day 14 (Fig. 3g). Comparing WT and RPGRIP1L-deficient organoids, no significant differences between MAP2 and TUJ1 expression on RNA or protein level was observed. In conclusion, neurogenesis appears unaltered and seems not to be responsible for the slightly reduced number of MNs observed in RPGRIP1L-deficient spinal organoids. Surprisingly, NKX2.2 RNA and protein expression was significantly increased at day 11 of MN differentiation in RPGRIP1L-deficient cells (Fig. 3d, h). Since NKX2.2 is expressed in various neural progenitors during differentiation, its increased expression could be caused by generation of different cell types in RPGRIP1L-deficient organoids. For instance, NKX2.2 is expressed in visceral pMNs of the hindbrain that do not express OLIG249,50 as well as in a specific subpopulation of ventral spinal pMNs in humans, where it is co-expressed with OLIG251,52. Furthermore, NKX2.2 is expressed in the most ventral NPs (p3) and neuronal subtype (V3) within the neural tube (Fig. 1a), where its expression plays a primary role in ventral neuronal patterning and requires high SHH signaling53. To test different hypotheses, we first performed co-labeling of OLIG2 and NKX2.2 at day 9 and day 11 of differentiation. We did not detect any increase in OLIG2 + /NKX2.2+ cells in RPGRIP1L-deficient organoids compared to WT controls (Fig. 3h and Supplementary Fig. 3g). Consequently, the increased NKX2.2 amount in RPGRIP1L-deficient spinal organoids is unlikely to be caused by increased numbers of ventral and human-specific MN subtypes. However, this does not rule out the possibility of an increased production of hindbrain pMNs in RPGRIP1L-deficient spinal organoids. Next, to test whether increased levels of NXK2.2 could be related to the formation of more ventral NPs and spinal cord neurons, we determined SHH pathway activity in WT and RPGRIP1L-deficient spinal organoids by analyzing the expression profiles of SHH target genes GLI1 and PTCH1 along spinal differentiation.

RPGRIP1L deficient spinal neural progenitor possess cilia and transduce SHH signaling

In both WT and RPGRIP1L-deficient spinal organoids, an increase in GLI1 and PTCH1 expression in response to SAG stimulation was observed between day 4 and day 9 of differentiation, followed by an expected decrease in their expression in agreement with the removal of SAG from the differentiation medium (Fig. 4a and Supplementary Fig. 3h, i). No significant difference in SHH target gene expression was detected between WT and RPGRIP1L-deficient organoids. This result is in line with the functional differentiation of RPGRIP1L-deficient hiPSCs into OLIG2+ pMNs and ISLET+ postmitotic MNs that highly depends on SHH signaling transduction. In contrast, the earlier described increase in the number of NKX2.2+ cells in RPGRIP1L-deficient spinal organoids (Fig. 3d, h and Supplementary Fig. 3g) cannot be explained by an increased SHH pathway activity and a ventralization effect of NPs. The remaining hypothesis that could explain the increased amount of NKX2.2 in RPGRIP1L-deficient spinal organoids is therefore the production of visceral NKX2.2+ pMNs.

Fig. 4: RPGRIP1L-deficient hiPSC-derived neural progenitors harbor cilia and transduce SHH signaling.
figure 4

a Log2(tpm+1) data from bulk RNASeq analysis show expression profiles of GLI1 and PTCH1 as mean ± SEM. b, c Immunofluorescence analysis of ciliary proteins in WT and RPGRIP1L-deficient spinal organoids at day 6. b Cilia are labeled by ARL13B and basal bodies by γ-TUBULIN. c Cilia are labeled by IFT81 in green. d, e Cilia length and density measurements in WT and RGPRIP1L-/- organoids at day 6. Data are shown as median with quartiles. Mann-Whitney test was performed for statistical analyzes. f Stimulated-Emission-Depletion (STED) images of WT and RPGRIP1L-deficient cilia at day 6. Cilia are labeled by ARL13B and RPGRIP1L. g Quantification of the ciliary RPGRIP1L amount based on confocal images. Data are shown as median with quartiles. Statistics: unpaired t tests with Welch’s correction (p < 0.0001). h, i Immunofluorescence of WT and RPGRIP1L-deficient cilia on spinal organoids at day 6. Cilia are labeled by ARL13B and AC3 (h) or by ARL13B and INPP5E (i). Basal bodies are labeled by γ-TUBULIN. j, k, l Quantifications of ciliary ARL13B (j), AC3 (k) and INPP5E (l) amounts in WT and RPGRIP1L-deficient organoids. Data are shown as median with quartiles. Statistics: unpaired t tests with Welch’s correction (p < 0.0001). a-l N: number of independent experiments; n: number of different WT or KO clones analyzed per experiment. For a, 2 WT clones and 1 KO clone from each iPSC line (n = 4 for WT and n = 2 for KO). For d, e, g and jl, 2 WT clones and 2 KO clones from each line (n = 4 for each genotype). N = 3 for (d, e, g, j-l); N = 1 for a. Scale bars: 300 µm in b; 0.5 µm in c; 1 µm in h, i.

Based on the results of functional SHH transduction and MN differentiation in human RPGRIP1L-deficient spinal organoids, we expected RPGRIP1L-deficient spinal organoids to possess cilia. We analyzed cilia formation at day 6 of MN differentiation, a time point where organoids were stimulated with SAG in order to activate SHH signaling and started to adopt the pMN fate. In WT organoids, short primary cilia were detectable on cells lining the inner lumens of the organoids, thereby pointing into internal cavities (Fig. 4b). These cilia were IFT81-positive and displayed a mean length of 1.3 µm (Fig. 4c, d). RPGRIP1L was present at the ciliary TZ in WT but not RPGRIP1L-deficient spinal organoids (Fig. 4f, g and Supplementary Fig. 2c). Despite the loss of RPGRIP1L, cilia were present lining apical surfaces (Fig. 4b) and displaying an unaltered length and density (Fig. 4d, f) in RPGRIP1L-deficient organoids. However, while WT cilia accumulated the ciliary membrane markers ARL13B, Inositol Polyphosphate-5-Phosphatase E (INPP5E) and Adenylate Cyclase 3 (AC3) (Fig. 4h, i), reduced amounts of ARL13B and AC3 (Fig. 4j, k) but not INPP5E (Fig. 4l) were observed in RPGRIP1L-deficient cilia. This indicates a disrupted ciliary gating in human spinal NPs upon RPGRIP1L deficiency. These results are in line with data from C.elegans, different mouse tissues and human cell lines, in which a disrupted TZ organization and in consequence an altered ciliary gating has been described under Rpgrip1l deficiency17,18,19,54,55,56. Since ciliary gating can be involved in SHH signaling function2,54,57 e.g., by regulating ciliary shuttling of pathway components such as SMO (Smoothened) thereby regulating ciliary GPR161 (G Protein-Coupled Receptor 161) amounts58,59, we tested whether GPR161-trafficking was disrupted in RPGRIP1L mutant spinal organoids. GPR161 was enriched in the ciliary membrane in SHH pathway-inactive WT and RPGRIP1L-deficient cells (Supplementary Fig. 3j). In WT and RPGRIP1L mutant NPs stimulated with SAG, GPR161 did not accumulate within primary cilia and thus did not inhibit SHH pathway activation (Supplementary Fig. 3j). These results are in perfect agreement with the previously described functional SHH pathway activation in WT and RPGRIP1L-deficient spinal organoids and show that ciliary gating defects did not affect SHH signaling transduction in RPGRIP1L-deficient NPs.

Given the mild phenotypes observed so far during spinal differentiation of RPGRIP1L-/- hiPSCs compared to the results obtained in the murine system, we wondered whether mutations could lead to residual RPGRIP1L activity. To validate the phenotype of our RPGRIP1L KO clones, we generated an additional RPGRIP1L-deficient hiPSC clone, this time carrying a full deletion of RPGRIP1L (Δexon3-exon27) (Supplementary Fig. 4a, b). The phenotype of this RPGRIP1L-deletion clone was consistent with our previous results obtained from clones carrying InDel mutations in RPGRIP1L. We observed the loss of RPGRIP1L at the TZ of cilia in NPs (Supplementary Fig. 4c, d) with cilia displaying an unaltered length (Supplementary Fig. 4e) and density (Supplementary Fig. 4f). Furthermore, SHH signaling was unaltered as SHH target gene expression of GLI1 and PTCH1 was induced as efficiently in the full deletion RPGRIP1L clone as in WT controls (Supplementary Fig. 4g, h). ISLET1 expression was highly upregulated in control and RPGRIP1L-deficient organoids at day 14 (Supplementary Fig. 4i) and quantification of ISLET+ cells at day 14 validated the efficient generation of MNs in control and RPGRIP1L-deletion clones (Supplementary Fig. 4j, k). Taken together, the consistency of the results demonstrates the authenticity of the phenotype along our RPGRIP1L mutant spinal organoids.

To test whether the results obtained in RPGRIP1L KO spinal organoids could be extended to other neurodevelopmental ciliopathy genes, we generated TMEM67-deficient hiPSCs via CRISPR/Cas9-mediated genome editing in the PCIi033-A hiPSC line. One control (TMEM67+/-) and one KO (TMEM67-/-) clone were used for further analyzes (Supplementary Fig. 5a-d). Cilia were present on control and TMEM67-deficient organoids at day 6 of differentiation where they are supposed to transduce SHH signaling (Supplementary Fig. 5c, e). No alterations of cilia length or quantity have been observed (Supplementary Fig. 5f, g), nonetheless gating defects might be present indicated by lower ciliary ARL13B intensities (Supplementary Fig. 5c). As observed in RPGRIP1L-deficient NPs, gating defects in TMEM67-deficient organoids did not affect SHH transduction demonstrated by unaltered SHH target gene expression of GLI1 and PTCH1 (Supplementary Fig. 5h, i). Furthermore, control and TMEM67 KO organoids started to express ISLET1 between day 9 and day 14 of differentiation (Supplementary Fig. 5j) without significant differences in the percentage of ISLET+ cells per organoid at day 14 (Supplementary Fig. 5k). Thus, human iPSCs mutant for two distinct ciliopathy genes encoding TZ proteins, RPGRIP1L and TMEM67, are able to form cilia, transduce the SHH pathway and differentiate into MNs, strongly suggesting that human and mouse spinal progenitors react differently to the loss of TZ proteins.

RPGRIP1L-deficient spinal organoids show antero-posterior patterning defects with MNs adopting hindbrain identities

To extend the analyzes of WT and RPGRIP1L KO organoids beyond consideration of MN differentiation, we analyzed the bulk RNAseq data in greater detail. Principal Component Analysis (PCA) demonstrated that WT and KO clones clustered closely together at each analyzed time point, indicating low variation between the genotypes (Supplementary Figs. 6a, b). Further testing for overrepresented gene categories via gene ontology (GOSEQ) did not reveal significant changes in biological processes, molecular functions or cellular components between WT and RPGRIP1L-deficient samples at different time points. Among the few genes that were differentially expressed between WT and RPGRIP1L-deficient organoids, expression of HOX genes such as HOXC6 and HOXA7 were significantly downregulated in RPGRIP1L-deficient organoids at day 11 and day 14 (Supplementary Fig. 6c, d, Supplementary Table 3). In addition, ensemble of gene set enrichment analyses (EGSEA) on GeneSetDB gene ontology (GO) at day 2 revealed significant downregulation of the biological process anterior/posterior pattern specification (GO:0009952) involving downregulation of HOX genes as well as the CDX genes CDX1, CDX2 and CDX4 (Supplementary Fig. 6e). In addition to the pairwise comparison of WT and RPGRIP1L KO gene expressions at each selected time point, we performed differential temporal profiles analyzes on our data. These analyzes revealed additional genes whose expression varied significantly over time (Supplementary Fig. 6f), and many of which were related to antero-posterior pattern specification, such as HOXC8 and HOXC9 (Supplementary Figs. 6f, g).

HOX genes are expressed in four clusters along the rostro-caudal axis of bilateral organisms: HOXA, HOXB, HOXC and HOXD. In humans, HOX proteins are encoded by 39 genes whose regionalized expression is realized by their sequential and collinear activation in axial progenitors, which gradually give rise to mesodermal and neuroectodermal structures with rostral and later caudal identity, each expressing a specific HOX code48,60,61,62,63,64. Importantly, upon differentiation HOXB genes become restricted to dorsal spinal neurons and are not expressed in MNs, therefore they were not considered in our analyzes62,65. By analyzing the dynamic expression pattern of the remaining 3 HOX gene clusters over time, we identified early differences in the expression of posterior HOX genes in RPGRIP1L-deficient organoids compared to controls. In WT, early axial progenitors at day 2 started to express anterior HOX genes such as HOXA1 followed by activation of later and more posterior HOX genes such as HOXA7 at day 4 (Fig. 5a), as it was previously reported48. At later stages, WT organoids expressed anterior as well as more posterior HOX genes (Fig. 5a, d and Supplementary Fig. 7a) as HOXA genes from HOXA1 to HOXA7 (Fig. 5a–c and Supplementary Fig. 7b, c) and HOXC genes from HOXC4 to HOXC8 (Fig. 5d-f and Supplementary Fig. 7d, e). Consistently to the HOX gene expression data, MNs within WT organoids adopted a caudal brachial antero-posterior identity at day 14, that could be visualized by immunofluorescence labeling of HOXA7 and HOXC8, markers for brachial MNs62 (Fig. 5g, h), whereas more anterior HOX proteins such as HOXA5, a cervical spinal marker48,62, were not expressed (Fig. 5i). During early stages of differentiation, RPGRIP1L-deficient progenitors expressed anterior HOX genes comparable to the WT situation while expression of the more posterior HOX genes HOXA7 and HOXC8 was already reduced at day 4 (Fig. 5a-f, Supplementary Table 3). This difference was even more striking at later stages, where the latest most caudal HOX genes expressed in RPGRIP1L-deficient MNs were HOXA5 instead of HOXA7 (Fig. 5a, b, c and Supplementary Fig. 7b, c) and HOXC4/HOXC6 instead of HOXC8 (Fig. 5d, e, f and Supplementary Fig. 7d, e). These results were confirmed by immunofluorescence analyzes, which showed that RPGRIP1L-deficient MNs do not express HOXA7 and HOXC8 (Fig. 5g, h), but significantly increased amounts of HOXA5 (Fig. 5i). Taken together, the differences in HOX gene expression indicate a shift in antero-posterior patterning in RPGRIP1L-deficient organoids in which MNs adopt more cervical identities than the caudal brachial identity observed in WT organoids (Fig. 5j) [57, 65, 67].

Fig. 5: RPGRIP1L-deficient spinal organoids show antero-posterior patterning defects with MNs adopting hindbrain identity.
figure 5

a, d Heatmaps showing HOXA (a) and HOXC (d) gene expression in WT and RPGRIP1L-deficient spinal organoids over time. The graphs were generated from log2(tpm+1) files of bulk RNASeq analysis. b, c, e, f Log2(tpm+1) graphs show the expression profiles of selected HOX genes HOXA5 (b), HOXA7 (c), HOXC4 (e) and HOXC8 (f) over time. Data are shown as mean ± SEM. g, h, i Immunofluorescence of HOX proteins in WT and RPGRIP1L-deficient spinal organoids at day 14. Quantifications show the percentage of HOXA7 (g) HOXC8 (h) and HOXA5 (i) positive nuclei per organoid. Data are shown as median with quartiles. Statistics: unpaired t tests with Welch’s correction (p < 0.0001). j Schematic overview about the three HOX gene clusters expressed in MNs from anterior to posterior positions. The black and gray-dotted squares indicate the genes expressed in WT and RPGRIP1L-deficient spinal organoids, respectively. k, l Log2(tpm+1) graphs generated from bulk RNASeq analyzes show the expression profiles of PHOX2B and TBX20. Data are shown as mean ± SEM. m Immunofluorescence of PHOX2B in WT and RPGRIP1L-deficient spinal organoids at day 14. Quantifications show the percentage of PHOX2B positive nuclei per organoid. Data are shown as median with quartiles. Statistics: unpaired t tests with Welch’s correction (p = 0.0004). am N: number of independent experiments; n: represents the number of different WT or KO clones analyzed per experiment. For (af), k and l, 2 WT clones and 1 KO clone from each line (n = 4 for WT and n = 2 for KO). For gi and m, 2 WT clones and 2 KO clones from each line (n = 4). N = 3 in gi, m; N = 1 in af, k, l. Scale bars: 500 µm in g, h, i and m.

We further confirmed the anterior shift in MN identities upon loss of RPGRIP1L by analyzing hindbrain specific pMN and MN markers in WT and RPGRIP1L KO organoids. Expression of PHOX2B, a marker of visceral pMNs and MNs of the hindbrain66, and TBX20, a late hindbrain-specific MN marker in mice67, was increased in RPGRIP1L KO compared to WT spinal organoids (Fig. 5k, l and Supplementary Fig. 7f). In line with this, the percentage of PHOX2B+ nuclei per organoid was significantly elevated in RPGRIP1L-deficient organoids (Fig. 5m). The results were confirmed in the full deletion RPGRIP1L KO clone (Supplementary Fig. 4l-p), where an increased expression of HOXA5 (Supplementary Fig. 4l) combined with decreased expressions of HOXA7 and HOXC8 was observed (Supplementary Fig. 4l-o). In line with this, the percentage of PHOX2B+ cells per organoid was significantly increased in the full deletion clone of RPGRIP1L (Supplementary Fig. 4p).

Interestingly, the results emphasize the previously discussed hypothesis that the increase of NKX2.2+ cells in RPGRIP1L-deficient organoids (Fig. 3d, h) could be explained by the formation of visceral PHOX2B+ MNs.

Analyzes of early hindbrain and spine-specific markers further confirmed the shift in antero-posterior patterning upon RPGRIP1L loss. The expressions of the hindbrain patterning markers EGR2/KROX2068,69,70,71, FGF371,72,73, PHOX2A66,68,74 and the visceral pMN marker NKX2.8/NKX2.950,75 were strikingly increased in RPGRIP1L-deficient organoids at day 2, day 4 and day 6 (Supplementary Fig. 7g–j). In contrast, expression of WNT5A, involved in axis elongation in mammals76, was decreased in early RPGRIP1L KO progenitors (Supplementary Fig. 7k). In line with this, expression of the spine-specific markers CDX1 and CDX4 was strongly reduced in RPGRIP1L-deficient organoids (Supplementary Fig. 7l, m), further emphasizing the caudal to cranial shift in antero-posterior identities of RPGRIP1L KO progenitors77,78,79,80,81,82.

RPGRIP1L controls axial progenitor fate specification and ciliogenesis in early human spinal organoids

The sequential and collinear expression of HOX genes that leads to specification of rostral and later caudal identities is described to be set-up within neuromesodermal progenitors (NMPs)60,61,62,63,83,84,85,86,87,88. It has previously been shown that cells in the human spinal organoid model used in this study transit through an early NMP-like identity, characterized by the co-expression of CDX, BRACHYURY and SOX247,48.

To validate that WT and RPGRIP1L-deficient cells adopted the axial NMP-like fate during differentiation, we analyzed the co-expression of BRACHYURY and SOX2 in organoids at day 2 and day 4. SOX2 was equally expressed in WT and RPGRIP1L KO organoids (Supplementary Fig. 8a) with nearly 100% of SOX2+ cells per organoid at day 2 and day 4 (Supplementary Fig. 8c, d, g, h). In contrast, while BRACHYURY expression was equally induced in WT and RPGRIP1L-deficient organoids at day 2 (Supplementary Fig. 8b, c, e, f), its expression dropped as soon as day 4 in RPGRIP1L-deficient cells (Supplementary Fig. 8g, i, j) indicating an altered induction of the NMP-like fate upon RPGRIP1L loss. In addition to BRACHYURY and SOX2, axial NMP-like progenitors in WT organoids expressed the caudal marker CDX2 at day 2 and day 4 of differentiation (Fig. 6a-d). Analyzes of CDX2 gene expression profiles showed its strong upregulation between day 0 and day 2 of differentiation, followed by its downregulation from day 4 onwards (Fig. 6e). In RPGRIP1L-deficient organoids, CDX2 expression was reduced (Fig. 6e, Supplementary Table 3) and significantly less CDX2+ cells per organoid were detectable at day 2 and day 4 of differentiation (Fig. 6a–d), with clusters of CDX2-negative cells in RPGRIP1L-deficient but not WT organoids (Fig. 6c). Together, the results of reduced BRACHYURY and CDX2 expression in KO organoids point to a defective induction of the axial NMP-like cell fate upon RPGRIP1L loss.

Fig. 6: Altered axial patterning and reduced ciliogenesis in early RPGRIP1L-deficient progenitors.
figure 6

a, c Immunofluorescence of CDX2 in WT and RPGRIP1L-deficient spinal organoids at day 2 (a) and day 4 (c). b, d Quantifications show the percentage of CDX2 positive nuclei per organoid at day 2 (b) and day 4 (d). Data are shown as median with quartiles. Statistics: unpaired t tests with Welch’s correction with (b) p = 0.0007 and (d) p = 0.0003. e Log2(tpm+1) graph generated from bulk RNASeq analyzes show the dynamic expression profile of CDX2. Data are shown as mean ± SEM. f, i, j Immunofluorescence of cilia in WT and RPGRIP1L-deficient spinal organoids at (f) day 2 and (i, j) day 4. Cilia are labeled by (f, i) ARL13B and INPP5E or by (j) IFT81 and ARL13B. Magnified areas are indicated by yellow rectangles and magnified images are displayed on the right. Scale bars: 5 µm. g, h, k, l Cilia length and cilia density measurements in WT and RPGRIP1L-/- organoids at day 2 and day 4. Data are shown as median with quartiles. Statistics: g, k Mann-Whitney (p = 0.0007) and h, l unpaired t tests with Welch’s correction (p < 0.0001). a-l N: number of independent experiments; n: number of different WT or KO clones analyzed per experiment. For b, d, g, h, k and l, 2 WT clones and 2 KO clones from each line (n = 4). For e, 2 WT clones and 1 KO clone from each line (n = 4 for WT and n = 2 for KO). N = 2 (except for 1 KO N = 1) for b; N = 3 for d; N = 1 for e, g, h, k, l. Scale bars: 150 µm in a, c, f, i and j.

We next analyzed the ciliary status of axial progenitors at day 2 and day 4 of differentiation. At day 2, cilia were present in WT organoids in an unpolarized fashion (Fig. 6f–h) whereas cilia density was significantly reduced in RPGRIP1L KO organoids (Fig. 6f, h). The same difference was observed at day 4 of differentiation. Here, cilia were present at the apical sites of WT organoids (Fig. 6i-l), but absent from RPGRIP1L KO cells (Fig. 6i-l).

Since the loss of cilia in early RPGRIP1L-deficient organoids correlated with the stage of NMP-like fate induction and antero-posterior specification of progenitors,we hypothesize that cilia loss on early RPGRIP1L-deficient progenitors is responsible for the reduced specification of the NMP-like axial fate that lead to the antero-posterior patterning defects observed in RPGRIP1L KO organoids. To further investigate this correlation, we analyzed TMEM67-deficient organoids. Cilia were present on control and KO progenitors at day 4 of differentiation (Supplementary Fig. 9a, b), with TMEM67 KO cilia showing the same length and density as controls (Supplementary Fig. 9c, d). Furthermore, the NMP-like fate was correctly acquired in control and TMEM67 KO cells, demonstrated by BRACHYURY and SOX2 expression at day 4 (Supplementary Fig. 9e). Consequently, no antero-posterior patterning defect was observed in TMEM67-deficient organoids at day 14. Expression of the HOX genes HOXA5, HOXA7, HOXC6 and HOXC8 was induced between day 4 and day 14 in control and TMEM67 KO organoids (Supplementary Fig. 9f) and the percentage of HOXA7+ and HOXC8+ cells per organoid was unaltered in TMEM67-deficient organoids compared to controls (Supplementary Fig. 9g). In both control and KO organoids a few HOXA5+ nuclei were present (Supplementary Fig. 9g). Furthermore, low expression of PHOX2B was detected in control and TMEM67 KO organoids at day 14 (Supplementary Fig. 9h), consistent with a few PHOX2B+ cells in control and KO organoids (Supplementary Fig. 9i). These results strengthen the correlation between cilia loss, altered NMP-like identity and perturbed antero-posterior patterning of spinal progenitors.

Interestingly, no spinal-to-cranial shift of progenitor identity was observed in Rpgrip1l-/- mouse organoids, as evidenced by the unaltered expression of Cdx1-2-4, and Cyp26a1 in mutant progenitors compared to controls (Supplementary Fig. 10a-d). Furthermore, no difference in the efficiency of axial progenitors induction was observed in mutant organoids, where at day 3 similarly to control organoids, the majority of cells was positive for Cdx2 (Supplementary Fig. 10e-f). To gain more insight into the antero-posterior identity of control and Rpgrip1l-/- mouse spinal progenitors present in our organoids, we also analyzed Hox gene expression by bulk RNAseq and qPCR analyzes. We found that in both genotypes, organoids expressed predominantly brachial and anterior thoracic Hox genes (Supplementary Fig. 10g-k). Moreover, although we noticed a slight decrease in Hox10 gene expression, no major alterations were observed in Hox gene expression that might indicate severe antero-posterior specification defects in the mutants (Supplementary Fig. 10g-k). Since the antero-posterior specification defect in RPGRIP1L KO human organoids was correlated with a reduction in NMP-like cells co-expressing BRACHYURY and SOX2, we analyzed this population in mouse organoids. At day 3, although the percentage of progenitors co-expressing Brachyury and Sox2 was only 20% in average, it did not change between control and mutant organoids (Supplementary Fig. 10l-m), nor did the levels of Brachyury expression (Supplementary Fig. 10n), suggesting that at least part of the progenitors migrate through a bi-potent NMP-like state in both genotypes. Finally, analysis of cilia in organoids at day 3 revealed that mouse axial progenitors deficient for Rpgrip1l presented abnormally shortened cilia (Supplementary Fig. 10o-p), similar to later stages.

Discussion

Previous analysis of Rpgrip1l and Tmem67 mutant mice demonstrated that deficiency for these TZ proteins led to DV mis-patterning of the neural tube marked by the strong decrease of ventral neuronal subtypes like V3 and MNs6,15,32. These phenotypes were correlated to the absence of primary cilia and functional SHH signaling, and were also described for Ift mutant mice15,43,89. Since these mice are lacking or displaying a reduced FP, a question remained as to whether altered patterning of the neural tube resulted only from a reduced SHH secretion ventrally, or from impairment of SHH signal transduction within ventral spinal progenitors themselves. Here, by exposing control and Rpgrip1l mutant spinal organoids to high dose of SHH agonist, we show that primary cilia are intrinsically required within mouse spinal progenitors to transduce SHH. Upon inactive SHH signaling, mutant progenitors switch on intermediate-dorsal identities instead of ventral ones.

Strikingly, we find that contrary to the mouse, human RPGRIP1L and TMEM67 are not required for the establishment of ventral spinal neuronal cell types in organoids. Despite a clear ciliary gating defect, with decreased ARL13B or AC3 ciliary localization, human mutant spinal NPs were still able to grow primary cilia and transduce SHH signaling similarly to controls, adopting MN identities. Ciliogenesis in the absence of Rpgrip1l and Tmem67 has also been observed in other cellular contexts, like in mouse limb bud progenitors, MEFs, HEK cells18,19,20,32 or human fibroblasts derived from RPGRIP1L and TMEM67 mutated JS patients17,90. Whether the difference in ciliary phenotype observed between mouse and human spinal progenitors stems from differences in ciliary composition or stability, and accounts for the difference in SHH activation in these two systems remains to be addressed. In addition to the species-specific differences in the ciliary phenotype, our results highlight a context-specific role for ciliary TZ proteins in human neural progenitors. In RPGRIP1L-deficient human organoids, NPs at day 6 display cilia of normal length and a moderately altered ciliary content, while at day 2 and day 4, progenitors display a strong reduction of the cilia density. Differences in ciliary structure and composition between tissues, already suggested by the high heterogeneity of phenotypes in human ciliopathies, have been described in other models18,91,92,93. However, here we uncover a change in cilia sensitivity to the loss of a TZ protein along the process of neural progenitor specification in humans, at a key step of cell fate transition.

Whether this change parallels differences in ciliary structure, and which signaling and molecular cascades lead to these changes in cilia stability during the differentiation process are important questions raised by our study. In this regard, recent data show that cell fate can regulate ciliogenesis by different mechanisms such as transcriptional regulation94, niche condition-specific regulation95,96, and several non-canonical signaling pathway activities such as SHH97 and WNT98. It will be a key future task to examine cell type specific cilia features and explore whether different modes of ciliogenesis regulation play a role in diverse neural progenitors to explain at least in part the regional specificity of defects found in neurodevelopmental ciliopathies.

Quite unexpectedly, we found that although they do not display SHH-dependent dorsoventral patterning defects, human RPGRIP1L mutant spinal organoids shift to more anterior (hindbrain/cervical spinal) identities as assessed by altered HOX expression patterns and activation of genes specific for hindbrain neurons such as PHOX2B. This switch occurs at early stages of HOX activation, at the axial progenitor stage (day 2 in our spinal organoids) and is paralleled by the above discussed transient loss of cilia and a reduction in the expression of CDX genes. Later (at day 4), reduced numbers of BRACHYURY+/SOX2+ NMP-like cells are found in RPGRIP1L mutants. In TMEM67 mutants, which do not lose cilia at early stages, no perturbation of AP patterning nor NMP-like fate establishment is observed; Thus, although this is not a proof of causal relationship, we find a correlation between cilia loss at axial progenitor and NMP-like stage and the activation of a more anterior progenitor program in spinal organoids. Why TMEM67 depletion does not lead to cilia loss at day 4, while RPGRIP1L depletion does, remains to be understood. As mentioned earlier, these two TZ proteins have different structures and functions and participate in distinct complexes at the TZ8,18,31,34. Thus, it is possible that they contribute differently to TZ stability and ciliary gating. The analysis of other ciliary gene mutants will be essential to understand the contribution of different TZ proteins in cilia stability.

Another remaining question is the role of Rpgrip1l and more generally cilia, in early antero-posterior specification in the mouse embryo. Our study suggests that, unlike in humans, Rpgrip1l is not required for the specification of the antero-posterior identity of murine axial progenitors. This could reveal an intrinsic difference in the role of Rpgrip1l in ciliary construction and function between human and mouse at the level of axial progenitors. In support of this hypothesis, the ciliary phenotypes observed in axial progenitors appear to differ between the two species as in human organoids, a decrease in cilia density is observed, while in murine organoids cilia are present but truncated. Furthermore, we cannot exclude that variations in the protocols used to generate human and murine organoids may also be at the root of the antero-posterior phenotype differences. In mouse organoids, NMP-like progenitors were observed at a much lower frequency compared to humans, probably because of early addition of retinoic acid at day 2. Such differences in NMP-like cell kinetics might mask the phenotype in mouse organoids. Re-analysis of the phenotype of Rpgrip1l KO mice and other ciliary mutants in the light of these results will be essential to tackle the role of cilia in antero-posterior regionalization of the embryonic axis in an in vivo context. The altered antero-posterior patterning that was observed in human RPGRIP1L KO organoids could arise from several mechanisms. WNT, RA, Fibroblast Growth Factor (FGF), and more recently NOTCH signaling, have been involved in the determination of axial progenitors and thereby regulation of gradual HOX gene expression47,64,77,99,100,101,102,103. From our transcriptome data, we could not observe any indication of perturbations in FGF or NOTCH pathways upon RPGRIP1L deficiency. Analysis of Wnt3a- and Cdx-deficient mouse embryos78,104,105 or experimental manipulation of CDX or BRACHYURY activity106,107,108,109 lead to a working model in which early BRACHYURY and CDX expression activates and fine-tunes HOX gene expression profiles along the antero-posterior axis downstream of WNT/β-Catenin signaling in axial progenitors77,79,83,85,99,100,104,106,109,110,111,112. Interestingly, it was demonstrated recently that WNT signaling during early differentiation of hiPSCs into spinal MNs allows progressive, collinear HOX activation and that a switch to RA treatment at any time point during this process stops collinear HOX activation and favors transition into neuroectodermal progenitors with a discrete HOX gene profile48,100. Thus, we hypothesize that the shift in antero-posterior patterning that we describe in RPGRIP1L-deficient human spinal organoids results from reduced WNT transduction during early stages of axial progenitor differentiation, leading to a reduction in BRACHYURY, CDX2, CDX1 and CDX4 expression and the consequent decreased expression of posterior HOX genes, so that MNs adopt hindbrain or cervical rather than caudal brachial identities. Whether RPGRIP1L or cilia act directly on maintaining WNT signaling in axial progenitors remains to be demonstrated. While cilia and TZ proteins, and specifically RPGRIP1L and TMEM67, have been involved in modulating WNT pathways, their effect is complex and depends on the tissue and cell type22,113,114,115,116,117. Further analyzes will show if and how cilia regulate WNT signaling and consequently BRACHYURY, CDX and HOX gene expression during axial progenitor fate specification, or whether other pathways such as FGF and RA are implicated.

Another important question raised by our study is whether the defect in antero-posterior patterning upon loss of RPGRIP1L could be related to phenotypic features of neurodevelopmental ciliopathies. Since the antero-posterior patterning defect of MNs described here most likely lead to a mild shift of boundaries between the posterior hindbrain region and the anterior spinal region with an increase in anterior identities, a corresponding phenotype could be present in ciliopathy patients. Unfortunately, analyzes of such patterning defects in patients are not available. A first important step would be therefore to test the differentiation capacity of patient-derived hiPSCs and analyze their phenotype related to antero-posterior patterning. Apart from that, our study provides interesting entry points for further investigations where it will be important to analyze antero-posterior patterning more carefully and in additional CNS regions such as the hindbrain, as it is highly sensitive to antero-posterior patterning defects and one of the major targets of neurodevelopmental ciliopathies. Here, the increased expression of PHOX2B appears to be particularly interesting, since its deregulation was shown to cause breathing pattern regulation defects118,119,120,121,122, a feature that can also be observed in JS patients123.

Methods

Ethics

The work performed in this manuscript complies with all relevant ethical regulations. Mouse breeding, mating and embryo collection were performed under ethical agreement number APAFIS-1382 from the French Ministry of Higher Education and Research. The obtaining and manipulation of human iPSC lines were performed under ethical agreements from the French Ministry of Higher Education and Research DC-2021-4734 (for the UCSFi001-A line), DC-2015-2595 and DC-2016-A00773-480 (for the PCIi033-A line) and GMO authorization number DUO-9509.

Generation of mESCs lines

Rpgrip1l-/- and Rpgrip1l+/+ mouse ES cell lines were derived from blastocysts according to standard procedures using 2-4 months-aged male and female Ftm-deficient mice produced and genotyped as described previously15 and maintained in the C57BL/6 J background. One ES clone from each genotype was used in this study. Both Rpgrip1l-/- and Rpgrip1l+/+ clones are males, which avoids sex-dependent bias in RNAseq data.

mESCs culture and differentiation into ventral spinal organoids

mESCs were cultured in DMEM-GlutaMAX (Gibco) supplemented with 15% FCS, 100 µM β-mercaptoethanol, 100 U/ml penicillin, 0.1 mg/ml streptomycin, and 103 U/ml LIF (Millipore). To generate mouse spinal organoids, low-passage (inferior to 10) mESCs were dissociated and differentiated into embryoid bodies by seeding them in ultra-low attachment petri dishes (Corning) at a density of 9 × 104 cells ml−1 in N2B27 differentiation medium containing Advanced DMEMF12 and Neurobasal media (1:1, Life Technologies), B27 w/o vitamin A (Life Technologies), N2 (Life Technologies), 2 mM L-glutamine (Life Technologies), 0.1 mM β-mercaptoethanol, penicillin and streptomycin (Life Technologies). On day 2, the medium was changed and supplemented with 3 µM CHIR, 10 nM RA and 500 nM SAG to induce ventral spinal fate. The medium was changed daily, with addition of 10 nM RA only.

hiPSCs culture

hiPSCs were thawed in presence of Rock-inhibitor Y-27632 (5 µM, Stemgent-Ozyme #04-0012) and cultured under standard conditions at 37 °C in mTeSR+ medium (Stem Cell Technologies #100-0276) on Matrigel (Corning, VWR #354277) coated plates upon confluency of 70-80%. Passaging was performed using ReLeSR (Stem Cell Technologies #05872) and testing for potential mycoplasma contamination was performed regularly by Eurofins genomic (Mycoplasmacheck). Accutase (Stem Cell Technologies #07920) was used for the dissociation of hiPSC colonies into single cells.

Generation of mutant hiPSC lines

RPGRIP1L mutant and wildtype hiPSC clones were produced via CRISPR/Cas9-mediated genome editing in two different hiPSC lines: PCIi033-A (PHENOCELL (PLI), product reference PCi_CAU2) and UCSFi001-A (deposited at the Coriell Institute for Medical Research under the identifier GM25256). Both cell lines are reported to be males. The targeted exon 3 of RPGRIP1L and flanking sequences were sequenced prior to gene editing in both cell lines (see on- and off-target analyzes section). The gRNA targeting exon 3 of RPGRIP1L was designed using the online design tool CRISPOR124. Annealing of crRNA (Alt-R® CRISPR-Cas9 crRNA, 400 nM, IDT) and tracrRNA (Alt-R® CRISPR-Cas9 tracrRNA, 400 nM, IDT #1072533) was performed and 400.000 hiPSCs of each line were transfected with 1 µl Cas9 protein (30 µM; kindly provided by TACGENE, Paris, France) and 1 µl gRNA (200 µM) using nucleofection. The transfected cells were cultured in mTeSR+ medium (Stem Cell Technologies #100-0276) containing Rock-inhibitor Y-27632 (10 µM, Stemgent-Ozyme #04-0012) for one day, replaced by mTeSR+ medium without Rock-inhibitor. Medium was changed every second day and serial dilutions were performed in mTeSR+ medium containing Rock-inhibitor Y-27632 and CloneR (1:1000, Stem Cell Technologies #05889). Rock-inhibitor was removed after two days and the medium containing CloneR was changed every second day until colonies were ready to be picked. Colonies were isolated in 24-well plates in mTesR+ medium containing CloneR. DNA was isolated for genotype analyzes (see on- and off-target analyzes section). Four WT (two clones of each hiPSC line) and four homozygous mutant hiPSC clones (two of each hiPSC line) carrying InDel mutations that lead to premature STOP were used for further analyzes (Supplementary Fig. 2).

A large deletion of RPGRIP1L was induced in the PCli033-A line by combining the gRNA targeting exon 3 of RPGRIP1L with a gRNA targeting exon 27 of RPGRIP1L (Δexon3-exon27). This approach was performed at the ICV-iPS platform of the Brain Institute in Paris (ICM, Paris). One RPGRIP1L KO clone and one heterozygous control clone were used for further analyzes (Supplementary Fig. 3 and Supplementary Fig. 11). No WT clone was recovered in this experiment. For the generation of TMEM67-deficient hiPSC clones from the PCli033-A cell line, the following strategy was used: The gRNA targeting the 5’UTR of TMEM67 was combined with a gRNA targeting intron 27 of TMEM67 to generate a large deletion of TMEM67 (Δexon1-exon27). One TMEM67 KO clone and one heterozygous control clone were used for further analyzes (Supplementary Fig. 5 and Supplementary Fig. 12). No WT clone was recovered in this experiment. For all clones, genomic stability was assessed by detection of recurrent genetic abnormalities using the iCS-digitalTM PSC test, provided as a service by Stem Genomics (https://www.stemgenomics.com), as described previously125.

On- and off-target analyzes

Genomic DNA was extracted from both hiPSC lines (PCli033-A and UCSFi001-A) prior to CRISPR/Cas9 gene editing. Genotyping of on-targets was performed by PCR, and PCR products were sent for sequencing. After CRISPR/Cas9 gene editing, genomic DNA was isolated from individual hiPSC clones and the on-target sites were amplified by PCR. The PCR products were sent for sequencing to determine the specific mutations that were induced. Off-targets were identified using the CRISPOR online design tool126 and ranked according to their CFD off-target score127. Potential off-targets within exons and with CFD scores > 0.02 were selected for PCR and subsequent sequencing analyzes128. The primers used for on- and off-target analyzes are listed in Supplementary Table 1.

Differentiation of hiPSCs into spinal organoids

Differentiation of wild-type and mutant hiPSCs was performed as previously described. After amplification, hiPSC lines were dissociated into single cells using Accutase (Stem Cell Technologies #07920) and resuspended in differentiation medium N2B27 [vol:vol; Advanced DMEM/F-12 (Gibco) and Neurobasal Medium (Gibco)] supplemented with N2 (Thermo Fisher #17502048), B27 without Vitamin A (Thermo Fisher #12587010), penicillin/streptomycin 1% (Thermo Fisher #15140122), β-mercaptoethanol 0.1% (Thermo Fisher #31350010). Cells were seeded in ultra-low attachment dishes (Corning #3261) to allow EB formation. The N2B27 differentiation medium was used throughout the whole differentiation process (Fig. 3a). Rock-inhibitor Y-27632 (5 μM; Stemgent-Ozyme #04-0012) was added from day 0 to day 2, CHIR-99021 (3 µM; Stemgent-Ozyme #04-0004) from day 0 to day 4, SB431542 (20 μM; Stemgent-Ozyme #04-0010) from day 0 to day 3 and LDN 193189 (0.1 μM; Stemgent-Ozyme #04-0074) from day 0 to day 4. The differentiation was preceded by adding Smoothened AGonist (SAG 500 nM; Merck #566660) and Retinoic Acid (100 nM RA; Sigma #R2625) to the N2B27 differentiation medium from day 4 to day 9. γ-Secretase inhibitor DAPT (10 µM; Tocris Bioscience #2634) was added from day 9 to the end of differentiation. Medium was changed every other day and additionally at day 3 of differentiation.

EB embedding and Cryosectioning

Murine and human spinal organoids were collected at different time points during differentiation, rinsed with PBS and fixed in cold PFA (4%) for 7–12 min at 4 °C. EBs were rinsed in PBS and incubated in 30% sucrose in PBS until completely saturated. Cryoprotected EBs were embedded in OCT embedding matrix (Cell Path #KMA-0100-00A) and stored at -80 °C. 12 µm cryostat sections were prepared and stored at -80 °C.

Western Blot

Mouse embryonic stem cells were collected during amplification and lysis was performed using RIPA buffer (150 mM NaCl, 25 mM Tris pH 7.5, 1% NP40, 0.1% SDS, Sodium Deoxycholate 1%) containing proteases inhibitors. Protein samples in 4x sample buffer were loaded on precast gels (Mini-PROTEAN TGX Gels, BioRad #4561093) and separated by SDS-PAGE in the Mini-PROTEAN Tetra Cell system (BioRad). The Spectra Multicolor High Range Protein Ladder (Thermo Fisher #26625) was used as a molecular weight marker. Blotting was performed on PVDF membranes (BioRad #162-0176) using the Mini Trans-Blot Cell tank blotting protocol (BioRad). Primary antibodies against Rpgrip1l and Actin, and secondary antibodies conjugated to horseradish peroxidase were used (Supplementary Table 2). Visualization of protein bands was realized using the ECL detection kit (SuperSignal WestPico PLUS Chemiluminescent substrate, ThermoScientific #34577) and the ChemiDoc MP Imaging System (BioRad).

Immunofluorescence

Cryostat sections of murine and human organoids were washed in PBS/0.1% Triton X-100. Optional steps for permeabilization with PBS/0.5 %Triton X-100 for 10 min, post-fixation with MeOH (100%) for 10 min or antigen-retrieval with citrate-buffer were performed depending on the antibody and specimen. After washing, cryosections were incubated with a blocking solution containing 10% NGS in PBS/0.1% Triton-X100 for at least 2 hours. The sections were incubated with primary antibodies (Supplementary Table 2) diluted in PBS/0.1% Triton-X100 and 1% NGS overnight at 4 °C. Next day, the slides were washed 3 × 5 min with PBS/0.1% Triton-X100 and incubated with the secondary antibodies (Supplementary Table 2) diluted in PBS/0.1% Triton X-100 and 1% NGS for 2 hours. After final washing steps with PBS/0.1%Triton X-100 for at least 1 hour, sections were mounted with Vecta-Shield (Vector #H-1000) or Mowiol (Roth #0713.2) optionally containing DAPI (Merck #1.24653).

Image acquisition

Fluorescence image acquisition of patterning and differentiation markers was performed at room temperature using a Zeiss Observer Z1 equipped with an Apotome module, an Axiocam 506 monochrome cooled CCD camera and a 40x oil objective with a NA 1.3. All images were processed with Zen software (blue edition). Image analyzes and figure preparation were performed using Fiji (ImageJ; National Institutes of Health). For cilia analyzes, fluorescence images were acquired at room temperature using an inverted confocal microscope (LSM 710, Carl Zeiss AG or TCS SP5, Leica), a 63x oil objective with a NA 1.4, and a monochrome charge-coupled device camera. Z-stacks with 0.2 µm steps were acquired. Images were processed with Zen software (black edition). Image analyzes and figure preparation were performed using Fiji (ImageJ; National Institutes of Health).

STED images were acquired on the STEDYCON confocal and STED module (Abberior Instruments GmbH, Germany) connected to the lateral port of an IX83 inverted microscope (Evident Scientific). The STEDYCON SmartControl software was used for data acquisition. Imaging was performed using a UPLXAPO 100x oil immersion objective (Evident Scientific). For STED, a pulsed 775 nm laser was used. STAR RED was excited with a pulsed 640 nm laser, STAR ORANGE was excited with a pulsed 561 nm laser, and DAPI was excited using a continuous wave 405 nm source. Pixel size was 20 nm in xy direction, with a 250 nm step size in the z direction.

Quantification of fluorescent staining and statistical analyzes

The quantitative analysis of nuclear stainings was done either by scoring the percentage of positive nuclei for the selected markers, or by measuring the stained area. The percentage of positive nuclei for selected stainings was determined by first applying automated nuclear segmentation using the cellpose python environment (https://www.cellpose.org/), followed by manual analysis and scoring of the stained nuclei. The positive stained areas were defined and measured via setting common thresholds in Fiji. At least 10 organoids per clone (n = number of clones) and experiment (N = number of experiments) were analyzed for quantifications. Intensity measurements of ciliary proteins were performed on unprocessed images (raw data) as described before18,19,20,21,54. Analyzes of ciliary density were performed by manually counting cilia axonemes over the number of basal bodies (γ-Tubulin). Whenever possible, at least 100 cilia per clone were used for intensity, length and density measurements. After quantifications were performed, representative images were processed by means of background subtraction and contrast settings via Adobe Photoshop CS2.

Data are presented as mean ± SEM or as median with quartiles. Two-tailed t test with Welch’s correction was performed for all data in which two normally distributed datasets were compared. Mann-Whitney tests were performed for data sets that were not normally distributed. Two-tailed t test Analysis of variance (ANOVA) and Tukey honest significance difference (HSD) tests were used for all data in which more than two normally distributed datasets were compared. Kruskal-Wallis with Dunn’s multiple comparison test was used for comparison of more than two datasets with small or not normally distributed data. Two-sided tests with confidence intervals of 5% were used for all statistical analyzes. Following statistical significances were considered: *p < 0.05, **p < 0.01, ***p < 0.001 and ****p < 0.0001. All statistical data analysis and graph illustrations were performed using GraphPad Prism (GraphPad Software).

RNA isolation and qPCR analyzes

mESC- and hiPSC-derived organoids were collected at different stages during differentiation and washed with PBS. RNA was isolated using the RNeasy Kit (Qiagen #74104) and RNAse-free DNase Set (Qiagen #79254). Isolated RNA was transcribed into cDNA using the GoScript Reverse Transcriptase Kit (Promega #A5001). For quantitative real-time PCR, 50 ng of cDNA of each sample was used in a Maxima SYBR Green/ROX qPCR Master Mix 2x (Thermo Scientific #K0222). Reactions were run in a Step One Real-Time PCR System Thermal Cycling Block (Applied Biosystems #4376357). Primer pairs are listed in Supplementary Table 1. The analysis of real-time data was performed using the included StepOne Software version 2.0 and graphs were generated using GraphPad Prism (GraphPad Software).

Bulk RNAseq analysis on mouse and human spinal organoids

mESC- and hiPSC-derived organoids were collected at different stages during differentiation and washed with PBS. For mouse organoids, pooled organoids from Wt and Rpgrip1l KO were collected from 3 independent experiments (N = 3). For human organoids, material from two control clones of the PCIi033-A hiPSC line and the UCSFi001-A hiPSC line each (WT: n = 4) and from one RPGRIP1L-deficient clone of the PCIi033-A hiPSC line and the UCSFi001-A hiPSC line each (KO: n = 2) was collected. RNA was isolated using the RNeasy Kit (Qiagen #74104) and RNAse-free DNase Set (Qiagen #79254). RNA quality was determined by measuring the RNA integrity number (RIN) via the High Sensitive RNA Screen Tape Analysis kit (Agilent Technologies #5067) on the TapeStation system (Agilent Technologies). A RIN above 9.5 was considered as good quality RNA and 250 ng RNA in a total volume of 25 µl was prepared per sample for further procedure. Bulk RNAseq was realized by the Genotyping and Sequencing Core Facility of the Paris Brain Institute (iGenSeq, ICM Paris). RNAseq data processing was performed in Galaxy under supervision of the ARTBio platform (IBPS, Paris).

Paired-end RNA reads were aligned against the M. musculus mm39 and the homo sapiens hg38 genome, respectively, by using the STAR aligner (v2.7.10, Galaxy wrapper v4)129 and the gene annotations gtf files GRCm39.109 and GRCh38.109, respectively. Quality of sequencing was controlled with fastQC (Galaxy wrapper v0.74+galaxy0) and MultiQC (Galaxy wrapper v1.11+galaxy1)130. Gene expression was assessed by counting sequencing reads aligned to gene exons with featureCounts (Galaxy wrapper v2.0.3+galaxy2)131. Raw counts were further converted to normalized gene expression values using the log2(tpm+1) transformation where tpm is the count of transcript aligned reads per length of transcript (kb) per million of total mapped reads (Galaxy tool “cpm_tpm_rpk” v0.5.2). Principal Component Analyzes (PCA) were performed using r-factominer (v.2.9, Galaxy tool id “high_dimensions_visualisation” v4.3+galaxy0)132 and heatmap visualizations were produced with the Galaxy tool “high_dim_heatmap” (v3.1.3+galaxy0) using the normalized gene expression values. Differentially expressed genes (DEGs) were selected from the gene raw read counts using DESeq2 (Galaxy wrapper v2.11.40.8+galaxy0)133 and the Benjamini-Hochberg p-adjusted cutoff 0.01. DESeq2 statistical tables were used for generation of Volcano Plots (Galaxy tool id “volcanoplot” v0.0.5)134, for goseq tests of overrepresented gene ontology categories (Galaxy wrapper v1.50.0+galaxy0)124, as well as for ensemble of gene set enrichment analyzes (EGSEA) (Galaxy wrapper v1.20.0)135.

Analysis of differential temporal profiles on human spinal organoids

We obtained 62710 gene temporal profiles of expression levels [expressed in log2(tpm+1)], each composed of 7 time points in the WT and RPGRIP1L KO mutant contexts. We selected genes whose variance of expression across the 14 WT and KO time points is higher than 0.1, and we scaled each gene time point (WT and KO) using the minmax function on the TSPred R package, following which the scaled expression in each series is E’ = (E - Emin) / (Emax - Emin), the highest expression of series is 1 and the lowest expression of series is 0. Note that since we scaled wild type and mutant expressions together, their relationship is maintained during the minmax scaling.

We next split the data into gene time series of 7 points in WT or KO context, and for each gene, we computed the euclidean distance between the wild type and mutant series. The density function of the euclidean distance was empirically derived from the observed distances. Using this function, we selected genes with the 1% highest distances between time profiles in wild type or mutant context. Time profiles of selected genes were finally plotted using their unscaled expression values. Since the euclidean distance is not able to capture significant differences between all time series, we repeated the above analysis with 11 other distance metrics, including short time series (sts), dynamic time warping (dtw), time alignment measurement (tam), autocorrelation-based dissimilarity (acf), fourier discrete transform (fourier), compression-based dissimilarity (cdm), complexity-invariant (cid), Pearson’s correlation dissimilarities (cor), integrated periodogram-based dissimilarity (int.per), periodogram-based dissimilarity (per) and Frechet (frechet) distances.

Gene time profiles selected using the various distances metrics are compiled in table [distance mode] which indicates for each gene the metrics that returned a significant distance. In total, we selected 542 genes with significant distance between mutant and wild-type time profiles. These profiles are sorted according to the number of metrics returning significant distance (p < 0.01). R code used for this analysis is available here for download.

Reporting summary

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