An In Vitro Human Segmentation Clock Model Derived from Embryonic Stem Cells.

Defects in somitogenesis result in vertebral malformations at birth known as spondylocostal dysostosis (SCDO). Somites are formed with a species-specific periodicity controlled by the "segmentation clock," which comprises a group of oscillatory genes in the presomitic mesoderm. Here, we report that a segmentation clock model derived from human embryonic stem cells shows many hallmarks of the mammalian segmentation clock in vivo, including a dependence on the NOTCH and WNT signaling pathways. The gene expression oscillations are highly synchronized, displaying a periodicity specific to the human clock. Introduction of a point of mutation into HES7, a specific mutation previously associated with clinical SCDO, eliminated clock gene oscillations, successfully reproducing the defects in the segmentation clock. Thus, we provide a model for studying the previously inaccessible human segmentation clock to better understand the mechanisms contributing to congenital skeletal defects.

INTRODUCTION

The periodic development of symmetrical somite pairs during vertebrate embryogenesis is controlled by a molecular oscillator known as the segmentation clock (Palmeirim et al., 1997; Dale and Pourquié, 2000; Jiang et al., 2000; Gomez et al., 2008; Tam, 1981) that operates in the presomitic mesoderm (PSM). In humans, the first somite begins to form during Carnegie stage 9 (between 19 to 21 days), eventually resulting in a total of 42–44 somite pairs (Müller and O’Rahilly, 1986; Tam, 1981). The NOTCH signaling pathway plays a crucial role regulating gene oscillation during somitogenesis (Dale and Pourquié, 2000; Jiang et al., 2000; Dequéant et al., 2006; Lewis et al., 2009; Gibb et al., 2010; Hubaud and Pourquié, 2014). In particular, mutations in the NOTCH signaling pathway components and downstream targets, including DLL3, MESP2, LFNG, and HES7, all cause congenital vertebral malformation, likely due to defects occurring during somitogenesis (Bulman et al., 2000; Whittock et al., 2004; Sparrow et al., 2006, 2008, 2010, 2012). In mice, Hes7 null mutation causes neonatal lethality, whereas steady-state expression of Hes7 eliminates the oscillatory expression pattern (Bessho et al., 2001b, 2003). Remarkably, accelerating the tempo of Hes7 oscillation results in both faster somite formation and growth of additional numbers of vertebrate (Takashima et al., 2011; Harima et al., 2013). These data strongly argue that Hes7 is a key driver of the mouse segmentation clock. However, whether human HES7 displays a cyclic expression pattern in PSM cells is currently unknown. Comparisons of the transcriptomic data between model organisms display divergence of the segmentation clock genes as well as critical signaling differences (Xi et al., 2017; Chal et al., 2015; Krol et al., 2011). Given these disparities, we reasoned that a segmentation clock system derived from human embryonic stem cells (ESCs) might serve as a more relevant model to understand the human segmentation clock and to elucidate mechanisms of developmental disorders.

There are several published protocols for differentiation of human pluripotent stem cells (PSCs) into somites and their derivatives (Chal et al., 2015; Loh et al., 2016; Xi et al., 2017; Nakajima et al., 2018; Russell et al., 2018). A recent study using single-cell RNA-sequencing (scRNA-seq) analyses found that the cells pass through a transitory state that displays gene expression signatures similar to somitomeres before specifying into somite cells (Loh et al., 2016); however, no oscillatory gene expression pattern has been reported. Our prior studies found that species-specific developmental timing is conserved even in the in vitro environment (Barry et al., 2017), thus we hypothesized that the segmentation clock would remain operative in vitro. To test the hypothesis, we performed RNA-seq on human ESC-derived PSM cells every 30 min and captured a transient expression signature similar to somitogenesis. Next, we generated a segmentation clock reporter and detected a robust, synchronized HES7 oscillation with a constant human specific periodicity of ~5 h. We demonstrated that chemical inhibition and conditional transgene expression could be directly employed to further dissect the signaling interplay during the initiation and propagation of HES7 oscillation. To demonstrate the utility of our system, we introduced a C to T transition in exon 2 of the endogenous HES7 coding region (Sparrow et al., 2008). This single sub-stitutional mutation (R25W) leads to a congenital vertebrae malformation condition known as spondylocostal dysostosis-4 (SCDO4; OMIM 608059) (Sparrow et al., 2008, 2010, 2012, 2013). In cells homozygous for the mutation, we observed a complete disruption of HES7 oscillation in PSM cells. Altogether, we present a system to further understand the nature of the human segmentation clock as well as demonstrate the system’s potential as a platform to model developmental disorders.

RESULTS AND DISCUSSION

RNA-Seq Analyses Identified a Transient Somitogenesis Program

We set out to derive human PSM cells from ESCs by adapting previously described protocols to induce a somite cell state (Nakajima et al., 2018; Loh et al., 2016; Chal et al., 2015; Xi et al., 2017). Human ESCs were stepwise differentiated in chemically defined medium, first to mesendoderm by culturing cells in “mesendoderm medium” (which activates WNT, transforming growth factor β [TGF-β], and fibroblast growth factor [FGF] signaling pathways), then to PSM by culturing cells for the second day in “PSM medium” (which activates WNT and FGF signaling but inhibits TGF-β and BMP4 signaling), and lastly to somite cells by culturing cells for the third day in “somite medium” (inhibition of WNT, FGF, BMP [bone morphogenetic protein], and TGF-β signaling pathways) (Figures 1A and S1A; see STAR Methods for further details). Under these conditions, the expression of paraxial mesoderm and PSM markers (TBX6, HES7), as well as somite markers (MEOX1, DLL1, FOXC2, TCF15/PARAXIS) increased hundreds- to thousands-fold compared to undifferentiated cells (Figure S1B). Immunofluorescent staining and RNA-seq results confirmed the downregulation of pluripotency markers (POU5F1, NANOG, DNMT3B, ZFP42, LIN28A), followed by sequential and uniformed upregulation of markers representing mesendoderm (T, MIXL1, NODAL, GSC), paraxial mesoderm (TBX6, MSGN1), and somite (MEOX1, FOXC2, PAX3, TCF15/PARAXIS, RIPPLY1) (Figures 1A and 1B). Thus, the progression of human ESC differentiation toward somite cells are highly aligned with the known marker expression during mammalian somitogenesis in vivo (Hubaud and Pourquié, 2014; Oates et al., 2012; Chal et al., 2018; Chal and Pourquié, 2017; Hicks and Pyle, 2015).

Figure 1.

Human ESC Differentiation to PSM and Somite Cell States

(A) Schematic of differentiation strategy of human ESCs differentiation toward mesendoderm, PSM, and somite cell states. Immunofluorescence co-staining for POU5F1, T, TBX6, and MEOX1 for characterization of the differentiation protocol. All scale bars represent 100 μm.

(B) Heatmap of RNA-seq data of the somite differentiation. Triplicate samples are shown for each time point. Selected markers are provided to represent the ESC, mesendoderm, PSM, and somite cell states.

(C) PCA of RNA-seq data collected every 30 min for the first 12 h after switching from PSM medium to somite medium. Each time point is collected in duplicates and are indicated by the color key.

(D) Heatmap of selected marker gene expression from the experiment in (C), representing PSM, somitogenesis (blue font), and somite cell states. All expression values (normalized expect counts [nECs]) are scaled minimum to maximum expression per gene row, indicated as a horizontal bar.

To investigate the potential of dynamic gene expression pattern during differentiation, we performed RNA-seq collecting samples every 30 min during the first 12 h immediately following the switch to the somite medium. Principal component analysis (PCA) was employed to probe the transcriptomic transitions by reducing the dimensionality of RNA-seq data. Using the top 500 most variable genes across the entire time course, PCA revealed a highly coordinated, continuously progressive gene expression pattern over time (Figure 1C). Remarkably, a group of genes involved in somitogenesis and the segmentation clock were captured as a transitory expression signature, rising only after the downregulation of mesodermal markers but declining before the expression of somite markers (Figures 1D and S1C). Using correlation analysis, we further characterized these genes into three expression groups: (1) “early downregulated” (100 genes), (2) “transitory” (415 genes), or (3) “late upregulated” (162 genes) (Figure S1C). Gene Ontology (GO)-term enrichment analysis showed that early downregulated genes are enriched for terms including mesoderm morphogenesis and anterior-posterior axis specification; transitory genes are enriched for the NOTCH signaling pathway and mesenchymal cell development; and late upregulated genes are enriched for skeletal system development and spinal cord development (Figure S1C; Table S1).

The transitory group of genes also displayed a relatively synchronized gene expression pattern. Specifically, at a time window from 300–450 min (5–7 h), after switching to somite medium, a group of transiently expressed somitogenesis genes peaked once before downregulation (Figure 1D). This transitory, co-activation pattern could also be detected at the single-cell level when the previously published scRNA-seq data characterizing DLL1-positive somitomere cells was reanalyzed (Loh et al., 2016). Temporal reordering of single cells using Wave-Crest followed with Trendy peak gene analyses (Chu et al., 2016; Bacher et al., 2017, 2018; Vermillion et al., 2018) (Figures S1D and S1E) highlighted a subset of somitogenesis genes (i.e., LFNG, SNAI2, DLL1, HEYL) displaying a co-expression pattern concentrated around a window of ~108–132 reordered cells (Figure S1E). Interestingly, the HES7 expression displayed a flat pattern in this analysis (Figure S1E). Nonetheless, these results are consistent with our data that a transitory signature of somitogenesis gene expression can be detected during the transition from the PSM to the somite cell states (Figure 1D). To confirm our results, we further examined the expression of the segmentation clock genes (HES7 and LFNG) and somite markers (MEOX1 and FOXC2) by qRT-PCR. We found that HES7 and LFNG expression peaked once when the cells were switched to the somite medium (Figure S1F) but peaked twice when the cells were maintained in the PSM medium (Figure S1G). As expected, MEOX1 and FOXC2 were highly induced in the somite medium (Figure S1F) but remained minimally expressed in PSM medium (Figure S1G). These results suggest that a more sustained gene oscillation profile (see further) is likely to be detected when cells are maintained in the PSM state.

Generation of Segmentation Clock Reporter and Characterization of HES7 Oscillation

To directly observe gene oscillation in real time, we tagged the endogenous HES7 locus with NanoLuc (NLuc) Luciferase and 2A-tdTomato cassettes via CRISPR-Cas9-mediated gene targeting (Hou et al., 2013; Chu et al., 2016; Zhang et al., 2017) (Figure 2A). A normal karyotype was confirmed after expansion of the correctly gene-targeted clones (Figures S2A–S2F). Human ESCs differentiated to the PSM state displayed cobblestone-like morphology and were highly proliferative and mobile (Video S1). At the PSM state, we detected robust expression of HES7 reporter in the majority of cells (>95%~99%, Figure 2A; Video S1). To examine HES7 gene oscillation, the HES7 reporter line was differentiated toward PSM without switching to the somite medium (Figure 2B). At day three (~72 h post differentiation), five to six cycles of HES7 gene oscillation could be detected over ~32 h of real-time live cell luminescence recording. The initial number of cells seeded for differentiation is crucial for detecting HES7 oscillation. In serial dilution experiments, we tested cell seeding density across ~30-fold range (1.88 × 103 to 6.0 × 104 cells per well seeding) and identified that an optimal density of 7.5 × 103 cells per well (~2.4 × 104 cells/cm2 in 96-well plate) provided the most robust and reproducible HES7 oscillation profile (Figure 2B). Strikingly, the sinusoidal oscillation profile is highly consistent across independent differentiation experiments, displaying comparable peak-to-peak (311 min ± 15.1 SD) and valley-to-valley (308 min ± 17.6 SD) cycle time (Figure 2C). Thus, our experimental observation of ~5-h periodicity of HES7 oscillation supports the estimated 4–6 h of human somitogenesis periodicity (Müller and O’Rahilly, 1986; Hubaud and Pourquié, 2014; Sparrow et al., 2007).

Figure 2.

Detection and Quantification of HES7 Oscillation

(A) Schematic of knock-in HES7-NLuc-2A-tdTomato-PGK-NEO cassette to the endogenous HES7 locus. The arrow indicates the position of the SpCas9/guide RNA (gRNA) cut site. Grey boxes indicate the exons. The lower panel shows analysis by flow cytometry of HES7-reporter expression during the differentiation from human ESCs to PSM cells. The y axis indicates the tdTomato detection channel. The x axis indicates the GFP channel.

(B) Upper left panel: a schematic of the experimental strategy to measure HES7 oscillation in PSM medium. Three subpanels indicated the cell density for differentiation to measure gene oscillation. The y axis indicates the luminescence signals (scaled from minimum to maximum from representative samples in arbitrary units [a.u.]). The x axis indicates 32 h of live cell luminescence detection.

(C) Upper panel: a schematic of the quantification of peak-to-peak (blue) and valley-to-valley (red) periodicity. The lower panel shows the quantification results. The y axis indicates minutes. The x axis indicates quantifications from the first to the fourth peaks.

(D) Left panel: an example frame of NLuc luminescence imaging experiments. Overlay shows the locations of 25 regions of interest (ROIs) indicated as white boxes. Right panel: the quantification results of 25 ROIs over 20 h of time-lapse imaging. The y axis indicates the luminescence signals (scaled from minimum to maximum from representative samples in a.u.). The time-lapse video is shown as Video S2. Insert shows the corresponding bright field cell image. All scale bars represent 100 μm.

Using luminescence imaging analysis, we observed predominantly synchronized HES7 oscillatory expression across a field of PSM cells (Figure 2D; Video S2). As a control for imaging analysis, gene oscillation was never observed when NLuc expression was driven by a constitutive promoter (EF1a) (Figure S2G). To confirm synchronization at the single cell level, we mixed parental H1 human ESCs with HES7 reporter cells (at a ratio of 100:1) and examined their differentiation to PSM using luminescence imaging. Quantification of diluted reporter cells revealed that single cells distantly located from the larger clusters of cells still oscillated in a synchronized manner (Figure S2H). We therefore concluded that HES7 oscillation is a synchronized behavior among a population of PSM cells. A recent study using mouse ESC-derived PSM cells detected Hes7 reporter oscillation propagating from central-to-peripheral directions, accompanied with segmentation-like features in vitro (Matsumiya et al., 2018). While our imaging analysis detected robust HES7 oscillation (Figures 2A and S2H; Video S2), we did not observe morphological features mimicking segmentation formation in vitro. We reason that our fully chemically defined medium provides an environment to “lock in” the PSM state and should be advantageous for dissecting the signaling requirement for propagating gene oscillation (see further).

Following this rationale, we next determined the roles of key signaling pathways on HES7 oscillation using small molecule inhibitors or conditional transgene expression (Jiang et al., 2000; Lewis et al., 2009; Gibb et al., 2010). As expected, reduced intracellular NOTCH signaling by DAPT (a γ-secretase inhibitor) treatment quickly diminished HES7 oscillation (Figure 3A, left panel). Overexpression of NOTCH signaling by NOTCH Intracellular Domain (NICD) or HES7 transgene induction (+doxycycline [DOX]) also eliminated HES7 oscillation (Figure 3A, middle and right panels). These results support previous studies that constitutive expression of NOTCH signaling pathway components will render a non-oscillatory profile (Harima et al., 2013). WNT signaling pathways are also known to contribute to the NOTCH-driven segmentation clock (Dubrulle et al., 2001; Aulehla et al., 2003; Dequéant et al., 2006; Niwa et al., 2007; Sonnen et al., 2018; Naiche et al., 2011). Inhibition of conical WNT signaling via AXIN2 stabilization (using the small molecules IWR1 or XAV939) all dampened HES7 oscillation (Figure S3A). In addition, disrupting MYC-MAX dimerization or interfering with the ability of MYC to regulate its transcriptional targets (using the small molecules 10058-F4 or (+)-JQ1) also eliminates HES7 oscillation (Figure S3B), corroborating a recent report that Myc activity is required for segmentation clock gene oscillation in mice (Mastromina et al., 2018). Together, these data indicate that HES7 oscillation in human ESC-derived PSM cells recapitulates many features of segmentation clock in vivo. We thus conclude that the maintenance of the human segmentation clock in vitro requires a delicate balance among multiple signaling pathways (Dubrulle et al., 2001; Hubaud and Pourquié, 2014; Mastromina et al., 2018; Matsumiya et al., 2018; Sonnen et al., 2018).

Figure 3.

Characterization of HES7 Oscillation and Modeling SCDO4 with Genome Editing

(A) The impact of modulating NOTCH signaling pathways on HES7 oscillation. Left panel: DAPT treatment (magenta). Middle panel: inducible NICD transgene expression (blue). Right panel: inducible HES7 transgene expression (green). The controls are shown in black. Blue arrowhead indicates when the NLuc luciferin and small molecules or DOX were applied. All luminescence oscillation data presented is scaled from minimum to maximum signal and shown as mean ± SEM (shaded area), obtained from 4 to 5 replicates in each experiment.

(B) A schematic for measuring the propagation of HES7 oscillation. The lower panels show the comparison between PSM medium (black) and with base medium supplemented with A83-01 and CHIR (yellow) only. Right panel: PSM medium (black) compared with only withdrawal of A83-01 (green) or the withdrawal of CHIR (magenta).

(C) Strategy for genome editing to introduce a R25W substitution. The sequences for gRNA protospacer adjacent motif (PAM) domain (magenta) and targeted sequences (blue) are indicated. The codon to be mutated is underlined, CGG (Arg) to TGG (Trp). Right panel: Sanger sequencing results of the parental line (wild type [WT]) and a homozygous R25W clone.

(D) Representative oscillation profile of WT (black) overlay with a R25W clone (red). Right panel: quantification of HES7 expression level between WT and R25W clones. ****p < 0.0001, Student’s t test. Data presented as the mean ± SD.

(E) A proposed working model of human segmentation clock and HES7 oscillation. NOTCH signaling is upstream of HES7 activation and oscillation. During the transition from mesendoderm to PSM state, WNT/β-catenin signaling activation and TGF-β signaling inhibition are crucial to propagate HES7 oscillation. Luminescence is measured in arbitrary units (a.u.).

Activation of WNT and TGF-β Inhibition Are Required to Propagate HES7 Oscillation

We next attempted to determine the minimal number of signaling pathways required to induce and propagate HES7 oscillation from the mesendoderm state. We examined the impact of combinations or drop-outs of the four major signaling components or growth factors, namely A83-01 (TGF-β inhibition), CHIR99021 (CH, WNT activation), FGF2, and NOGGIN (BMP4 inhibition) in the PSM medium on the second day (~24 h post differentiation) (Figures 3B and S3C). As a control, the first peak of HES7 oscillation was detected as early as 6–8 h after switching to PSM medium (~30–32 h post differentiation) (Figure 3B, left panel). Interestingly, none of the single small molecules or growth factors were sufficient to induce HES7 oscillation even though HES7 was expressed in a significant portion of cells (~40% to >90%) (Figure S3C, conditions 2–6). Dual combinations revealed that supplying only A83-01 and CH in the base medium was sufficient to induce HES7 oscillation (Figure 3B, middle panel; Figure S3C, condition 7). When either A83-01 or CH was removed from the PSM medium, HES7 oscillation was disrupted (Figure 3B, right panel; Figure S3C, conditions 13–14). Importantly, antagonizing TGF-β signaling alone renders a ~5-fold higher activation of HES7 compared to control (Figure S3E). These results indicate that TGF-β signaling is crucial to suppress HES7 activation. Applying small molecules RepSox or SB-431542 (SB), which also target TGFBR1 (also known as ALK5) (Inman et al., 2002; Laping et al., 2002; Ichida et al., 2009; Tojo et al., 2005), had similar effects on HES7 activation (Figure S3E). Combining CH with RepSox or SB was sufficient to induce HES7 oscillation (Figure S3D). We therefore conclude that both WNT activation and TGF-β inhibition are required for oscillations from the mesendoderm state. Interestingly, TGF-β signaling inhibition alone renders a higher expression of total luminescence signals, while the PSM medium induced a more uniformed population activation of HES7 (Figures S3C and S3E). Collectively, these results corroborate with a recent report that observed the downregulation of the TGF-β signaling pathway in developing, nascent human somites (Xi et al., 2017). Our data further suggest that reducing TGF-β signaling is required for HES7 oscillation in the PSM state.

Introducing HES7 R25W Mutation to Model SCDO4

Mutations in HES7 cause congenital abnormalities of vertebrae that manifest as an irregular formation of the vertebrae and ribs, known as spondylocostal dysostosis-4 (SCDO4) (Sparrow et al., 2008, 2010, 2013). This disorder is likely due to complications during somitogenesis (Sparrow et al., 2012), but how the mutations disrupt the process of somitogenesis is unknown. To examine the consequences of HES7 mutation, we employed CRISPR-Cas9-mediated scarless gene editing (Liang et al., 2015, 2017; Richardson et al., 2016; Yang et al., 2013) to introduce a C to T transition (73C-T) in exon 2 of the HES7 locus (Sparrow et al., 2008). This mutation results in a missense mutation in HES7 protein (R25W) (Figure 3C). In this cohort of experiments, clones carrying both heterozygous (73 C/T) and homozygous (73 T/T) mutations were obtained (Figure S3F). In our conditions, HES7 oscillation was completely eliminated in R25W homozygous clones but displayed a normal oscillation profile in heterozygous clones (Figures S3G and S3H) and controls (Figure 3D, left panel). This result is consistent with the observation that siblings of individuals with SCDO4 carrying only heterozygous mutations appears to be unaffected (Sparrow et al., 2008). We also observed a significant >3-fold increase of HES7 expression in the R25W clones as measured by the luminescence signal intensity (Figure 3D, right panel). This is consistent with the notion that R25W mutation lies within the DNA binding domain and results in the inability of the HES7 protein to repress its target genes, including its own promoter (Bessho et al., 2001a, 2003). Together, these results demonstrate that HES7 oscillation can be tractably employed to report on the function of the human segmentation clock (Figure 3E). We therefore anticipate that our in vitro system can serve as a modeling platform to further understand the mechanistic link between human gene mutations and congenital vertebral defects.

In conclusion, we report a human segmentation clock system in vitro. Our HES7 reporter can be modulated by small molecule inhibitors or transgene expression, allowing us to dissect human-specific signaling dynamics (Figures 2 and 3). Specifically, the human clock displayed a ~5-h periodicity and is sensitive to TGF-β signaling. Furthermore, the in vivo segmentation clock displays unidirectional, posterior-to-anterior traveling waves (Masamizu et al., 2006), while our system exhibits a radial or predominantly synchronized oscillation profile within a two-dimensional PSM colony (Figures 2A and S2H; Video S2). This may be due to the lack of an anterior-posterior signaling axis in vitro. Future efforts combining our system with a tunable signaling gradient (Sonnen et al., 2018) or with a cell culture device for patterning cells (Warmflash et al., 2014) may allow for better modeling of the spatial context of the human segmentation clock.

One of the most fascinating features of somitogenesis is the species-specific periodicity of gene oscillation, which sets the tempo for each somite pair formation. Recently, a study reported mouse ESC-derived PSM cells that recapitulate a mouse segmentation clock of ~2.5-h periodicity (Matsumiya et al. 2018). This study further advances the idea that PSC-derived progenitors largely obey a species-specific developmental timing in vitro (Ebisuya and Briscoe, 2018; Barry et al., 2017). Improper human segmentation leads to congenital disorders such as SCDO or congenital scoliosis (Sparrow et al., 2006, 2008, 2010, 2012, 2013; Bulman et al., 2000; Whittock et al., 2004). We suggest that human ESC and/or patient-specific induced PSC (iPSC)-derived PSM cells can serve as a practical platform for modeling these disorders during this critical stage of human development.

STAR★ METHODS

LEAD CONTACT AND MATERIALS AVAILABILITY

Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact, James A. Thomson (JThomson@morgridge.org). Plasmids generated in this study have been deposited to Addgene, 130932-130936.

EXPERIMENT MODEL AND SUBJECT DETAILS

The experimental procedures for H1 human embryonic stem cells (ESCs) described in this study were approved by the ethics committee with IRB Approval Number: SC-2015-0010. The H1 human ESCs are registered in the NIH Human Embryonic Stem Cell Registry with the Approval Number: NIH hESC-10-0043 (https://grants.nih.gov/stem_cells/registry/current.htm?id=29).

METHOD DETAILS

Human ESCs Maintenance and Differentiation

H1 human embryonic stem cells (ESCs) were routinely maintained in E8 medium prepared in-house (Chen et al., 2011); commercially available as Essential 8 Medium (Thermo Fisher Scientific Cat. No. A1517001) on Matrigel (0.5 μg/6-well plate, BD Bioscience) coated tissue culture plates with daily medium feeding. Cells were passaged every 2 to 3 days with 0.5 mM EDTA in phosphate buffered saline (PBS) at 1:6 to 1:12 ratio for maintenance. Differentiation was performed similarly to previously established protocols (Chu et al., 2016; Zhang et al., 2017; Hou et al., 2015; Loh et al., 2016). In brief, H1 cells were individualized with Accutase (Thermo Fisher Scientific, Cat. No. A1110501), washed once, and then seeded onto Matrigel-coated plates at a density of 1.0–2.0 × 104 cells/cm2 and cultured in various types of differentiation medium supplemented with ROCK inhibitor Y-27632 (10 μM) overnight. For mesendoderm and PSM differentiation, individualized H1 cells were first seeded in E8 supplemented with Activin A (R&D Systems, 100 ng/mL) and CHIR99021 (3 μM) for the first 24 h, then switched to the PSM medium: E6 base medium (E8 medium minus FGF2, minus TGFβ1) supplemented with A83-01 (2 μM), CHIR99021 (3 μM), FGF2 (100 ng/mL), and NOGGIN (100 ng/mL) for another 24 h. The cells were then switched to somite medium: E6 base medium supplemented with PD0325901 (1 μM), IWR1 (2.5 μM), A83-01 (2 μM), and NOGGIN (100 ng/mL). Additional small molecules were used to test their impact on HES7 oscillation in the PSM state: XAV939 (10 μM), (+)-JQ-1 (10 μM), 10058-F4 (10 μM), RepSox (25 μM), and SB (10 μM). All small molecule inhibitors and reagents are listed in the Key Resources Table. All the cell cultures used in our laboratory are routinely tested as negative for mycoplasma contamination.

KEY RESOURCES TABLE

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Antibodies
Anti-POU5F1	Santa Cruz	Sc-5279; RRID: AB_628051
Anti-T (Brachyury) (APC-conjugated)	R&D Systems	IC2085A; RRID: AB_2271455
Anti-TBX6	R&D Systems	AF4744; RRID: AB_2200834
Anti-MEOX1	Atlas Antibodies	HPA045214; RRID: AB_10959299
Chicken anti-Mouse IgG (H+L) Cross-Adsorbed Secondary Antibody, Alexa Fluor 594	Thermo Fisher Scientific	A-21201; RRID: AB_2535787
Donkey anti-Goat IgG (H+L) Cross-Adsorbed Secondary Antibody, Alexa Fluor 350	Thermo Fisher Scientific	A-21081; RRID: AB_2535738
Chicken anti-Rabbit IgG (H+L) Cross-Adsorbed Secondary Antibody, Alexa Fluor 488	Thermo Fisher Scientific	A-21441; RRID: AB_2535859
Chemicals, Peptides, and Recombinant Proteins
CHIR 99021	Tocris	4423
A 83-01	Tocris	2939
DAPT	Tocris	2634
endo-IWR1	Tocris	3532
XAV 939	Tocris	3748
(+)-JQ1	Tocris	4499
10058-F4	Tocris	4406
Y-27632 dihydrochloride	Tocris	1254
RepSox	Tocris	3742
SB 431542	Tocris	1614
Recombinant Human Noggin Protein	R&D Systems	6057-NG-01M/CF
Recombinant Human TGF-beta 1 Protein, CF	R&D Systems	240-B/CF
Human Holo-Transferrin Protein, CF	R&D Systems	2914-HT
Porcn Inhibitor II, C59	Sigma-Aldrich	5004960001
Insulin solution human	Sigma-Aldrich	I9278
Doxycycline hyclate	Sigma-Aldrich	D9891
BD Matrigel Matrix Growth Factor Reduced	BD Biosciences	354230
UltraPure 0.5M EDTA, pH 8.0	Thermo Fisher Scientific	15575020
Dulbecco’s phosphate-buffered saline (DPBS)	Thermo Fisher Scientific	14190250
DMEM/F-12	Thermo Fisher Scientific	11320082
StemPro Accutase Cell Dissociation Reagent	Thermo Fisher Scientific	A1110501
TrueCut Cas9 Protein v2	Thermo Fisher Scientific	A36497
Critical Commercial Assays
GeneArt Precision gRNA Synthesis Kit	Thermo Fisher Scientific	A29377
GeneArt Genomic Cleavage Detection Kit	Thermo Fisher Scientific	A24372
POU5F1 TaqMan Assay (4331182)	Thermo Fisher Scientific	Hs00999634_gH
TBX6 TaqMan Assay (4331182)	Thermo Fisher Scientific	Hs00365539_m1
HES7 TaqMan Assay (4331182)	Thermo Fisher Scientific	Hs03988571_m1
LFNG TaqMan Assay (4331182)	Thermo Fisher Scientific	Hs00385436_g1
DLL1 TaqMan Assay (4331182)	Thermo Fisher Scientific	Hs00194509_m1
MEOX1 TaqMan Assay (4331182)	Thermo Fisher Scientific	Hs00244943_m1
FOXC2 TaqMan Assay (4331182)	Thermo Fisher Scientific	Hs00270951_s1
MSGN1 TaqMan Assay (4331182)	Thermo Fisher Scientific	Hs03405514_s1
GAPDH TaqMan Assay (VIC_PL, 4448486)	Thermo Fisher Scientific	Hs02786624_g1
Nano-Glo(R) Live Cell EX-6829	Promega	CS2025A03
Deposited Data
Raw RNA-seq data	This paper	GEO: GSE128469
Experimental Models: Cell Lines
Human: WA01 (H1) hESC line (NIH approval number NIHhESC-10-0043)	WiCell	WA01: H1 Human ES
Oligonucleotides
HES7 knockin junction PCR, in Table S2	This paper	Synthesized from IDT
HES7 for Southern blotting probe template, in Table S2	This paper	Synthesized from IDT
Primers for make HES7-C73T-gRNA template, in Table S2	This paper	Synthesized from IDT
Oligo for HES7-C73T-ssODN donor, in Table S2	This paper	Synthesized from IDT
Primers for sequencing HES7 Exon2, in Table S2	This paper	Synthesized from IDT
Recombinant DNA
HES7-NLuc-2A-tdTomato knockin vector	This paper	Addgene, 130932
pCRII-Blunt-U6-HES7-STOP-gRNA	This paper	Addgene, 130933
pB-TAG-NICD	This paper	Addgene, 130934
pB-TAG-HES7	This paper	Addgene, 130935
pB-EF1a-NLuc	This paper	Addgene, 130936
NanoLuc gene cassette	Promega	GeneBank: JQ437370
Software and Algorithms
Bandpass and detrending Perl code	This paper	https://github.com/stewart-lab/Detrending
Wave-Crest	Chu et al., 2016	https://github.com/lengning/WaveCrest
Trendy	Bacher et al., 2018	https://rbacher.shinyapps.io/trendy/
Enrichr	Chen et al., 2013	http://amp.pharm.mssm.edu/Enrichr/
Other
scRNA-seq data for somitomere cells	Loh et al., 2016	http://imlspenticton.uzh.ch:3838/conquer/

Gene Targeting and Construction of the HES7-NanoLuc-2A-tdTomato Reporter Cell Line

HES7-NanoLuc-2A-tdTomato donor vector with ~880 bp homology arms on each side of the stop codon was synthesized and cloned into a gene targeting vector backbone. CRISPR-Cas9 mediated gene targeting experiments were performed as previously described (Hou et al., 2015; Chu et al., 2016; Zhang et al., 2017). Briefly, H1 cells were individualized with Accutase, washed once with E8 medium, and resuspended in E8 medium supplemented with 10 mM HEPES buffer (Thermo Fisher Scientific, Cat. No. 15630080) and Y-27632 (10 μM). For electroporation, 2.5 × 106 cells were mixed with 7.5 μg of pCMV-hCas9 plasmid, 7.5 μg of gRNA construct, and 15 μg of donor DNA template in a total of 500 μL of cell/DNA suspension, transferred to a 4-mm cuvette (Bio-Rad, Cat. No. 1652008) and electroporated with a Gene Pulser Xcell Electroporation Systems (Bio-Rad, Cat. No. 1652660). Electroporation parameters used were 250 V, 500 μF, and infinite resistance. Cells were then plated into Matrigel-coated culture dishes in E8 medium supplied with Y-27632 (10 μM) overnight. Fresh medium was changed daily. Geneticin selection (Thermo Fisher Scientific, Cat. No. 10131027, 100 μg/mL) was started three days after electroporation. Geneticin-resistant colonies were picked 10 to 12 days later. Correct targeted clones were first identified by junction PCR and confirmed by Southern blot analysis. In the HES7-NLuc-2A-tdtomato targeting experiments, clone 10 was verified as a correctly targeted clone with normal karyotype and is used for all the subsequent experiments in this study. All DNA oligos were synthesized from Integrated DNA Technologies (IDT). Oligo sequences are listed in Table S2. The pCMV-hCas9 plasmid was a gift from George Church (Addgene plasmid # 41815) (Mali et al., 2013). All plasmids reported in this study and their sequences are available upon request.

PiggyBac Transposon to Generate Inducible Transgenic Reporter Cell Lines

Human gene open reading frame (ORF) cDNA clones were ordered from GeneCopoeia or GenScript. Entry vectors containing specific cDNA clones were subcloned to an all-in-one piggyBac destination vector (PB-TAG-ERP2) by using Gateway LR Clonase II Enzyme mix (Thermo Fisher Scientific, Cat. No. 11791020). PiggyBac constructs were electroporated into the H1-HES7-NLuc-2A-tdTomato reporter cell line. For electroporation, 2.5 × 106 cells were mixed with 1.5 μg of pCMV-hyPBase plasmid, 20.0 μg of PB-TAG-ERP2 piggyBac construct in a total of 500 μL of cell/DNA suspension, with the identical electroporation parameters listed previously. Cells were then plated at clonal density in E8 medium supplied with Y-27632 (10 μM) overnight. Fresh medium was changed daily. Puromycin (Thermo Fisher Scientific, Cat. No. A1113803) selection (1.0 μg/mL) was started two days after electroporation. Single Puromycin-resistant colonies were picked, expanded, and then tested for the efficiency of DOX induction (2000 ng/mL) based on percentages of GFP-positive expression and confirmed by qPCR analysis of the transgene expression. Only clones that reached >95% purity after DOX induction were used for further experiments. The PB-TAG-ERP2 was a gift from Knut Woltjen (Addgene plasmid # 804790) (Kim et al., 2016). All the piggyBac vectors and their sequences are available upon request.

Genome Editing for Point Mutation by HDR in HES7 Reporter Cell Line

A point mutation was introduced into human ESCs via homology-directed repair (HDR) triggered by CRISPR-Cas9 complexed with single-guide RNA (gRNA) as ribonucleoproteins (RNPs) along with single-stranded oligonucleotide donors (ssODNs) harboring the desired mutation as described previously (Liang et al., 2015, 2017). Specifically, assembled RNP and ssODNs (1 μM) were mixed with 1.0 × 105 parental HES7 reporter cells per electroporation using the Neon Transfection System 10 μL Kit, program #16 (pulse voltage: 1400, pulse width: 20, pulse no.: 2) (Thermo Fisher Scientific). Electroporated cells were then seeded at clonal density (500–1500 cells per 10 cm tissue culture dish) with E8 plus Y-27632 (10 μM) for the first three days. Single colonies were then picked 7 to 10 days post electroporation into 96-well plates for further expansion. Picked colonies were lysed, and the genomic region of HES7 exon 2 were PCR amplified and analyzed by Sanger sequencing performed at the University of Wisconsin Biotechnology Center DNA Sequencing Facility. Identified heterozygous and homozygous mutated 73C-T clones were then selected for expansion and oscillation measurements. All DNA oligos and ssODN donors were ordered from IDT, and full sequences are available in Table S2.

Live Cell Detection of NanoLuc Luminescence Activity and Quantification

Cells were seeded on Matrigel (BD Bioscience) coated 96-well optical plates (Thermo Fisher Scientific, Cat. No. 165306) for initiation of differentiation (2.4 × 104 cells/cm2). Standard HES7 oscillation detection experiments were performed on day 3 differentiated PSM cells (~72 hours post differentiation). For experiments described in Figures 3B and S3D, a density of 1.6 × 105 cells/cm2 were used. Prior to luminescence detection, fresh PSM medium supplemented with NanoLuc substrate Nano-Glo® Live Cell EX-6829 (1000X, Promega, Cat. No. CS2025A03) and small molecules and/or DOX were applied. Cells were allowed to equilibrate for at least 30 min in a regular incubator. The optical plate was then transferred to a GloMax Explorer Luminometer/Fluorometer Instrument (Promega, Cat. No. GM3510) set at 37°C and placed inside an environmental controlled incubator (Thermo Fisher Scientific, Cat. No. 3310) maintained at 5% CO2 during luminescence detection. All experiments reported in this study used the default parameter for 96-well Lid-ON setting, integration time 5 s, and each well is recorded with an interval of 5 min with continuous looping for at least 24 and up to 38 h.

Live Cell Imaging of NanoLuc Luminescence Activity and Quantification

Cells were seeded on a Matrigel-coated 35 mm glass-bottomed dish (MatTek, Cat. No. P35G-1.5-4-C) at 2.5 × 105 to 5.0 × 105 cells per dish to induce differentiation as previously described. All the HES7 luminescence imaging detection experiments were performed with PSM medium supplemented with NanoLuc substrate Nano-Glo® Live Cell EX-6829 (1000X, Promega, CS2025A03). After applying the substrate containing medium, cells were equilibrated for at least 30 min in a regular incubator before being transferred to an environment-controlled long-term bioluminescence imager (Olympus LV200). All the images were captured at 5-min intervals for 16 to 36 h with bright field and luminescence channels. Video S2 was processed by MetaMorph software 7.4 (Olympus). Quantification was performed using Fiji (Schindelin et al., 2012). Raw luminescence intensity data were obtained by manually selecting regions of interests (ROIs) across a population of cells or single cells. For the quantification of a population of cells, 10 to 20 identical sized, square-shaped ROIs were randomly placed to capture luminescence intensity in the field of view. For the quantification of single cells, the shapes of ROIs were adjusted across all the frames to include the traveling path of a single cell while avoiding overlap with neighboring cells. Background ROIs were obtained for each of the time series by selecting regions without observable luminescence activity, usually in the corners of the images. Raw luminesce data were then detrended and processed the same manner as described for data obtained with the plate reader.

RNA-Seq Library Preparation, Reads Processing, and Normalization

Cells gown in 6-well or 12-well cell culture plates were directly lysed with RLT Plus lysis buffer (QIAGEN) and stored in −80°C for later processing. Total RNA was purified from the RLT Plus lysates using the RNeasy Plus Micro Kit (QIAGEN) according to the manufacturer’s instructions. To construct the RNA-seq library, total RNA from each batch of samples (using an input of ~100 ng total RNA) were used following the LM-seq protocol (Hou et al., 2015). The reads generated from the Illumina HiSeq 3000 (69 cycles of insert read and 10 cycles of index read) were processed with CASAVA-1.8.2 basecalling software (Illumina). Reads were mapped via Bowtie (Langmead et al., 2009) against the hg19 RefSeq reference, and gene expression estimates were obtained as expected counts (ECs) and transcripts per million (TPM) via RSEM (Li and Dewey, 2011). A total of 58 bulk samples from two experiments were sequenced: (1) a daily time course of H1 human ESCs differentiating to the somite cell state, 12 samples collected as triplicates per time point and (2) differentiation from H1-derived PSM to the somite cell state, samples collected every 30 min in duplicates, 46 samples. All the gene expression expected counts (ECs) were first normalized by median-by-ratio normalization as indicated in the EBSeq package before downstream analysis (Leng et al., 2013). All the RNA-seq data reported in this study is available through GEO: GSE128469.

Immunofluorescence Staining and Confocal Image Analysis

Cells were seeded and cultured on Matrigel-coated glass-bottom culture dishes (MatTek, 12- or 24-well dishes) for differentiation. Cultured cells were then washed with PBS once then fixed with BD Cytofix at 4°C for 15 min. Cells were then permeabilized with 0.2% Triton X-100 (in PBS) at room temperature for 30 min. Cells were then blocked with blocking buffer (2%BSA and 1%FBS in PBS) for 1 h at room temperature followed by staining with primary antibody diluted in blocking buffer at 4°C overnight. The next day, the cells were washed three times with the blocking buffer before being incubated with AlexaFluor secondary antibodies (Thermo Fisher Scientific, 1:1000 dilutions in blocking buffer) and DAPI or Hoechst 33342 (Sigma, Cat. No. B2261) for 1 h at room temperature. Cells were then washed three times with the blocking buffer and ready for imaging. All the primary antibodies used in this study can be found in the Key Resources Table. Immunofluorescence images were collected using a Nikon A1R laser scanning confocal microscope with Plan Apo 10x, Plan Fluor 20x Ph1 DLL, or Plan Apo 20x DIC M objectives. Images were processed using NIS Elements or ImageJ. Z stacks are presented as maximum intensity projection images.

Reverse Transcription (RT) and qPCR Analysis

RNA isolation, RT reactions, and real-time quantitative PCR analysis were performed essentially as described previously (Chu et al., 2011, 2016). Briefly, cells were lysed in 350 uL of RLT plus buffer, and the total RNA were purified using RNeasy kits with either on-column DNase treatment or genomic DNA removal columns under the guidance of the user manual (QIAGEN). Between 100 ng to 500 ng of purified RNA were reverse transcribed with SuperScript IV VILO Master mix (Thermo Fisher Scientific, Cat. No. 11756050). To perform TaqMan qPCR reactions (10 μL total volumes), 1 μL of the diluted cDNA was subsequently used in each of the triplicate qPCR reactions with individual 1 × TaqMan Gene Expression assays and 1 × TaqMan Universal PCR Master Mix II (Thermo Fisher Scientific). The ViiA 7 System was used to perform qPDR, and ExpressionSuite was used to perform data analysis (all from Thermo Fisher Scientific). All TaqMan Gene Expression assays are from Thermo Fisher Scientific and are listed in the Key Resources Table.

QUANTIFICATION AND STATISTICAL ANALYSIS

Oscillation Data Processing – Detrending of Raw Luminescence Data

Time course luminescence signals (obtained from GloMax plate reader or quantified images) were detrended using a bandpass filter implemented in Perl. For each sample, the raw data were first high-pass filtered by subtracting, at each time point, the mean value over a broadly localized window around that time point. These results were then low-pass filtered via Gaussian smoothing over a much smaller localized window. All samples were run using the same window sizes: ±45 time points for the high-pass filter and ± 3 time points for the low-pass filter. The detrended data were then imported to Prism for visualization. The code for detrending raw data as well as example inputs and outputs are available at Github (https://github.com/stewart-lab).

Oscillation Data Analysis – Polynomial Fitting and Quantification of Oscillation Periodicity

The detrended luminescence data were fitted using polynomial regression to eliminate artifact noise and emphasize main peaks and valleys. We used the base function lm in R version 3.5.0 with polynomial degree equals to 20 by default. Based on raw data, 20 is large enough to illustrate the main patterns of all luminescence datasets while maintaining a reasonable computation complexity. The resolution of fitted data becomes per min instead of per 5 min, allowing for more accurate estimations. By determining the derivative and calculating roots, we calculated a lagged difference of fitted data to allow us to identify peaks and valleys. We identified peaks as “+ to −” sign changes (increasing to decreasing) and valleys as “− to +” sign changes (decreasing to increasing). The peak-to-peak time, valley-to-valley time, peak heights, and valley depths were then calculated. Since the luminescence data are usually noisier at the beginning and the end, we truncated the first 100 min and last 10% of data points during analysis. All luminescence data presented is scaled from minimum to maximum signals and shown as mean ± SEM obtained from 4 to 5 replicates of samples in each experiment.

Bulk High Temporal RNA-Seq Data Analysis – Principal Component Analysis, Pearson Correlation Analysis, and GO-Term Enrichment

Principal component analysis was performed using default principal component analysis in R, based on the top 500 variable genes across all the samples. Before Pearson correlation analysis, 5860 low-expressed genes whose maximum expressions among all time points were less than 10 nECs were filtered out. Downregulated genes (100 genes) were identified by the intersection of correlation ≥0.8 with the expression trend of T (326 genes), MIXL1 (120 genes), and MSGN1 (351 genes). Upregulated genes (162 genes) were identified by the intersection of correlation ≥0.8 with the expression trend of MEOX1 (360 genes), RIPPLY1 (267 genes), and TCF15 (262 genes). To identify transitory genes, we focused on seven marker genes (HES7, DLL1, NOTCH1, LFNG, AXIN2, SNAI2, HEYL) based on their strong literature support for being involved in somitogenesis (Bessho et al., 2001b; Bone et al., 2014; Forsberg et al., 1998; Jouve et al., 2000; Dale et al., 2006; Aulehla et al., 2003; Leimeister et al., 2000) and displayed a strong transitory signature in our data (Figure 1C). In Figure S1C, gene expression values (in nECs) were first smoothed by the bandpass filter code mentioned above (±1000 time points for the high-pass filter and ±4 time points for the low-pass filter). Smoothed data were processed with the following criteria: (1) the position of the peak was within the range of peaks of the seven marker genes; (2) the shape around the peak was within the seven marker genes, wherein shape is measured by the second derivative and estimated by the lagged difference of detrended curve; and (3) the variance of the lagged difference of the gene itself was smaller than those of seven marker genes (less noisy). A total of 415 transitory genes was identified based on the previous criteria. GO-term enrichment analysis was performed using Enrichr (http://amp.pharm.mssm.edu/Enrichr/) (Chen et al., 2013; Kuleshov et al., 2016). All of the downregulated, upregulated, and transitory gene lists, as well as Enrichr results, are in Table S1.

Reanalysis of scRNA-Seq Data from Loh et al., 2016

scRNA-seq data were downloaded from the conquer repository (Soneson and Robinson, 2018) (http://imlspenticton.uzh.ch:3838/conquer/). Dataset#40, SRP073808 contains three types of cells collected during human ES differentiation (Koh et al., 2016; Loh et al., 2016). This dataset contains a total of 188 cells: DLL1-positive paraxial mesoderm (PXM, 56 cells); D2.25 somitomere (87 cells), and early somite (45 cells). Expression counts were normalized using SCnorm (Bacher et al., 2017) (v1.5.7) using the default parameters and specifying each cell’s condition. To mitigate the effects of cell cycle heterogeneity in scRNA-seq data, we employed Seurat’s (Butler et al., 2018) (v2.3.4) CellCycleScoring function to score canonical G2M and S phase cell cycle markers and regress out the differences between phases as an additional preprocessing step (Nestorowa et al., 2016). Normalized, preprocessed scRNA-seq data were then reordered by Wave-Crest (Chu et al., 2016) (v0.0.1) with the default setting of degree = 3, and using 12 key markers according to their expression progression over time: FOXC2, LFNG, HEYL, HES7, HOPX, MEOX1, MIXL1, MSGN1, T, TBX6, TCF15, and RIPPLY1 (Figure S1D). To identify peak gene expression patterns in reordered single cells, we applied Trendy (v1.4.7), which characterizes piece-wise expression trends as up, down, or the same (Bacher et al., 2018). We used default parameters except we set minNumInSeg =10 requiring a minimum of 10 cells for a valid trend and the significance threshold for a trend to be considered either up or down to pvalCut=.2. Peaked expression was defined as genes with an up trend followed by a down trend. Among the list of genes with a peak, we manually checked for somitogenesis genes and selected representative genes are shown in Figure S1E.

Other Software

To generate figures and text, the following software packages were used: Microsoft Word, Excel, and PowerPoint for Mac v16.16.6; Adobe Illustrator CSS 2018 v22.1 and Photoshop CSS 2018 v19.1.2; Endnote X8.2; FlowJo 10.5.0; and Prism Graphpad 8.0.2.

DATA AND CODE AVAILABILITY

The RNA-seq data reported in this paper have been deposited in the NCBI Gene Expression Omnibus and are accessible through accession number GEO: GSE128469.